text:
# Temperature Variations from HST Spectroscopy of the Orion Nebula

## 1 INTRODUCTION

Most observational tests of the chemical evolution of the universe rest on emission line objects; these define the endpoints of stellar evolution and probe the current state of the interstellar medium (ISM). Gaseous nebulae (H II regions & Planetary Nebulae (PNs)) are laboratories for understanding physical processes in all emission-line sources, and probes for stellar, galactic, and primordial nucleosynthesis.

There is a fundamental issue that continues to be problematic – the discrepancy between heavy element abundances inferred from emission lines that are collisionally excited (CELs) compared with those due to recombination/cascading, the so-called “recombination lines”. Studies of PNs contrasting recombination and collisional abundances (Liu et al. 1995, Kwitter & Henry 1998) often find differences exceeding a factor of two. For NGC 7009, Liu et al. (1995) found that the recombination C, N, and O abundances are a factor of $``$5 larger than the corresponding collisional abundances, which was found to be the case also for neon (Luo, Liu & Barlow 2001). For NGC 6153, Liu et al. (2000) found that C++/H+, N++/H+, O++/H+, and Ne++/H+ ratios derived from optical recombination lines are all a factor of $``$10 higher than the corresponding values deduced from CELs. The discrepancy is even larger, a factor of $``$20, for the Galactic bulge PN M 1-42 (Liu et al. 2001) and Hf 2-2, which has the most extreme abundance difference to date, a factor of 84 (Liu 2002).

Studies of H II regions find similar behavior although the differences are generally smaller than is the case for PNs. For the Orion Nebula, Esteban et al. (1998) derived a factor of 1.5 larger for O++/H+ from recombination lines than from CELs. Using several ions for a study of M8, Esteban et al. (1999) found again that the ionic abundances are always somewhat lower as derived from collisional forbidden lines than those from corresponding recombination lines. The mean abundance difference is a factor of $``$2. Recently Tsamis et al. (2002) presented observations of 5 more H II regions – M 17 and NGC 3576 as well as the Magellanic Cloud H II regions 30 Doradus, LMC N11B and SMC N66. They found that the disparity was always in the same direction. For four of their objects, the O++/H+ abundance from O II recombination lines exceeded the corresponding value inferred from the nebular \[O III\] CELs by factors ranging from 1.8 – 2.7, while the factor was $``$5 for LMC N11B. According to the limited statistics, apparently H II regions exhibit less of an abundance dichotomy than some of the PNs.

Most of the efforts to explain the abundance puzzle between collisional and recombination values have attempted to do so by examining electron temperature ($`T_e`$) variations in the plasma. Because $`T_e`$ is not expected to vary dramatically within the (hydrogen) ionized region of a given nebula (e.g., Harrington et al. 1982; Baldwin et al. 1991; Rubin et al. 1991), a method developed by Peimbert (1967) has been extensively used. He developed a useful formalism by expressing the volume emissivity for a given spectral line in a Taylor series expansion around an average temperature defined such that the first order term vanishes. Since the \[fractional\] mean-square $`T_e`$ variation (called $`t^2`$) is expected to be small, only terms to second order need be retained. The resulting influence on elemental abundance determinations began with the work of Peimbert (1967), Rubin (1969), and Peimbert & Costero (1969). There is a vast literature now that measures the discrepancy between collisional and recombination abundances in terms of $`t^2`$.

The studies of H II regions mentioned above also derived the value of $`t^2`$ that forces the two abundance techniques to yield the same result. For the Orion Nebula, Esteban et al. (1998) required $`t^2`$ $``$ 0.024 to reconcile the difference, while for M8, Esteban et al. (1999) found an even larger $`t^2`$ = 0.032 was needed on average for the various ions in common. Tsamis et al. (2002) found for 30 Doradus a similar value, $`t^2`$ $``$ 0.03. For the drastically different O++/H+ abundances found for NGC 7009, Liu et al. (1995) discussed that it would be necessary to invoke $`t^2`$ $``$ 0.1, which would then force agreement close to the higher recombination value – a value more than 2.5 times larger than the solar O/H of 7.41$`\times `$10-4 (Grevesse & Sauval 1998). Such a large $`t^2`$ is not at all predicted by current theory/models (e.g., Kingdon & Ferland 1998).

The current unsettled situation has led to efforts to broaden the study to include other variables besides $`T_e`$ to analyse the effects upon abundance determinations. One promising avenue is to examine abundance derivation considering density variations, abundance variations, and $`T_e`$ variations in combination. Liu et al. (2000) took this approach in their investigation of the PN NGC 6153 with a two-phase empirical model. Péquignot et al. (2002) continued the study using photoionization models including two components with different heavy element abundances.

This is an extension of an earlier paper (Rubin et al. 2002 – Paper I) that used HST STIS and WFPC2 observations to study the variation of $`T_e`$ in NGC 7009 as determined from the \[O III\] (4364/5008) flux ratio. Very low values for $`t^2`$ in the plane of the sky ($`t_A^2`$) were found (always $``$0.01). In this paper, we focus solely on $`T_e`$ variations with the purpose to determine from the observational data the magnitude of $`t_A^2`$ for the Orion Nebula. In section 2, we present the HST observations and data reduction procedures. Section 3 includes a discussion and analysis of extinction. We determine the electron temperature distributions in section 4. In section 5, we analyse the $`T_e`$ distributions in terms of average temperatures and fractional mean-square temperature variations in the plane of the sky. Section 6 provides a discussion and conclusions.

## 2 HST OBSERVATIONS AND DATA REDUCTION

The observations of NGC 1976 described here were taken as part of our HST Cycle 7 program GO-7514. We observed with 4 different STIS long-slits: slits 1, 2, 4, and 5. These are shown in Figure 1. Our slit positions are chosen to cross several features, including an Herbig-Haro (HH) object, a proplyd, and the Orion Bar. Slit 1 passes through the position 1SW which we observed with FOS and GHRS in Cycle 5 (see Figure 1 in Rubin et al. 1997). The centre of 1SW is at $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}14\stackrel{\mathrm{s}}{\mathrm{.}}71`$, $`5^\mathrm{o}`$23′41$`\stackrel{}{\mathrm{.}}`$5 (all positions are equinox J2000), 18.5′′ S and 26.2′′ W of $`\theta ^1`$ Ori C. This slit was also chosen to pass through proplyd P159-350 (O’Dell & Wen 1994). Slit 2 passes through position x2 observed with FOS in Cycle 5 (see Figure 1 in Rubin et al. 1997). x2 is located in one of the most prominent arcs in the \[O III\] WFPC2 images. Slit 2 is parallel to slit 1, which allowed a substantial saving in overhead (OH). The position angle (PA) = 114.555o; the distance between slit centres is 18.1′′. Note that the bottom of slit 1 and the top of slit 2 have a separation of only 0.211′′. The displacement between these two slits in the direction along the slit is 18.056′′.

Slit 4 crosses the Orion Bar and is oriented to point toward $`\theta ^1`$ Ori C (the dominant exciting star), which also places it essentially orthogonal to the Bar. The southern tip passes through the Herbig-Haro object HH 203. Slit 5 passes through a very bright, sharply-defined “rim” of the Bar where a positional bifurcation begins. Slit 5 is parallel to slit 4 in order to take advantage of the same OH savings as is the case for slits 1 and 2. The PA = 139.068o; distance between centres of slits 4 and 5 is $``$32′′.

Observations were made as follows: data for slits 1 and 2 are from Visit 2 (1998 December 7 UT), Visit 52 (2000 December 7) and Visit 72 (2001 December 16); data for slits 4 and 5 are from Visit 5 (1998 December 22). All our data described here were acquired using the STIS/CCD with a 52′′ slit length and 0.5′′ slit width. Each visit comprises 2 orbits. Spectra were taken with gratings G430M (with wavelength settings: 3680, 4451, 4961) and G750M (settings: 5734, 6581). Each exposure was done in accumulation mode and at least two spectra were taken at each setting in order to cosmic-ray (CR) clean. The data sets we processed were those obtained after sufficient time elapsed from the observation dates in order that “best reference files” would be stable/finalized. We requested On-the-Fly Calibration for science files and Best Reference Files. After retrieving the data sets, we then co-added and cosmic ray cleaned images using standard packages in IRAF.111IRAF is distributed by NOAO, which is operated by AURA, under cooperative agreement with NSF. Calibrations to produce 2-dimensional (2D) rectified images were then carried out. From these, we singled out specific emission lines for further investigation. Data for the \[O III\] 4364 Å line as well as H$`\gamma `$ 4342 Å were contained in the G430M/4451 grating setting; for the \[O III\] 5008, 4960 Å lines and H$`\beta `$ 4863 Å lines, G430M/4961 was used; for \[N II\] 5756 Å, G750M/5734 was used; for H$`\alpha `$ 6565 Å and the \[N II\] 6550, 6585 Å lines, G750M/6581 was used. All wavelengths in this paper are vacuum rest wavelengths.

For G430M, the dispersion is 0.28 Å per pixel for a point source and the plate scale 0.05 arcsec/pixel (Leitherer et al. 2001, Chapter 13). For a uniformly filled slit with width 0.5′′, a degradation in resolving power by a factor of 10 is expected to a spectral resolution of 2.8 Å. For G750M, the dispersion is 0.56 Å per pixel for a point source and the plate scale 0.05 arcsec/pixel. For a uniformly filled slit with width 0.5′′, the spectral resolution would be 5.6 Å.

According to Paul Goudfrooij (private communication), the STIS CCD absolute flux calibration that is performed in the pipeline follows methods described in Instrument Science Report (ISR) 97-14 (Bohlin, Collins & Gonnella 1998) available through the STIS web site at www.stsci.edu/hst/stis/documents/isrs. That report considers only the L (low-resolution) modes specifically, but the flux calibration for the M (medium-resolution) modes is done the same way, i.e., by comparing the observed spectrum of a primary standard with a pure-hydrogen white dwarf model. The calibration observations were made with the 52′′$`\times `$2′′ slit; thus relative transmission corrections are necessary to derive the absolute fluxes for the other slits (such as the 52′′$`\times `$0.5′′ slit we used). This procedure is documented in ISR 98-20 (Bohlin & Hartig 1998).

For the analysis presented here, we were interested mainly in the distribution of line flux along the slit spatial direction. This was accomplished with the IRAF routine blkavg in conjunction with specialized software tools that we developed ourselves. Even after applying the standard CR rejection there still remain many bad pixels due to CRs and/or hot pixels. There is considerable danger that including these can corrupt the flux values we seek. The program developed to eliminate these remaining bad pixels is called PIXHUNTER, which has been described briefly in an earlier paper (see Appendix A of Paper I). Once the columns containing the line have been cleaned for deviant pixels, we are ready to subtract an equivalent spectral range of continuum. We do this by using IRAF functions, including blkavg, to operate upon the appropriate sections of cleaned continuum. The 1-dimensional (1D) distribution of line flux versus spatial coordinate for the various emission lines of interest is what we need for our subsequent physical analysis.

We note that there is excellent agreement with a cross check of the 1D results of flux versus spatial direction by comparing with 1D results of flux versus wavelength for a corresponding spatial sample. The latter were measured with the splot package. With this, the underlying continuum is fitted and the integrated line flux determined with the $`e`$-option (area under the line profile), which was preferable to fitting with a Gaussian profile. Because of the spectral impurity introduced by the relatively wide slit used, the line profiles have flatter tops and less extended bases (i.e., they are more “trapezoidal”) than the Gaussian fits. It is also apparent that the Gaussian fit is overestimating the line flux.

Both the \[O III\] 5008 and 4960 Å lines were observed simultaneously with theG430M/4961 grating setting. Because both transitions arise from the same upper level, the intrinsic flux ratio depends only on the transition probabilities (A-values) and wavelengths. As reported previously for similar STIS observations of the planetary nebula NGC 7009 obtained under program GO-8114 (PI RR) (Paper I), what we found was a surprising variation in the F(5008)/F(4960) ratio with position along the slit. This amounts to a variation in the ratio of roughly 3.0$`\pm `$0.1. Furthermore, it appears more-or-less periodic with an $``$3.5′′ cycle. As described in Paper I, according to Ted Gull (private communication), this is an instrumental effect and is a ratio of two fringe patterns. The source of the problem is a thin blocker filter that had to be matched with each grating and the best (and only) location that it could be placed was above the grating in a stable mounting. To attempt to do anything about fringing would probably require a dedicated HST/STIS calibration program. If there were fringing in the F(5008)/F(4364) ratio at the same level as for the F(5008)/F(4960) ratio, the $`\pm `$ 3.3% error would result in only a minor $`T_e`$ error, e.g., $`\pm `$ 100 K at $`T_e`$ = 104 K (see §4.1).

## 3 EXTINCTION and REDDENING CORRECTION

Before deriving the $`T_e`$ distribution from the STIS data, we first correct for extinction. This is calculated by comparing the observed F(H$`\alpha `$)/F(H$`\beta `$) ratio with the theoretical ratio I(H$`\alpha `$)/I(H$`\beta `$). We use a value of 2.88 assuming $`T_e`$ = 8500 K and $`N_e`$ = 5000 cm-3, Case B (Storey & Hummer 1995).222For conditions applicable to NGC 1976, we find using an online program (see Storey & Hummer 1995) that the I(H$`\alpha `$)/I(H$`\beta `$) ratio used here will depart by less than 3% over the range 7500 $``$ $`T_e`$ $``$ 12500 K and 103 $``$ $`N_e`$ $``$ 105 cm-3. The extinction correction is done in terms of $`c`$(H$`\beta `$), given by the relationship

$$log[F(\lambda )/F(H\beta )]=log[I(\lambda )/I(H\beta )]f(\lambda )c(H\beta ),$$
(1)
where $`f(\lambda `$) is the extinction curve. For the H$`\gamma `$ (4342), 4364, 4863, 4960, 5008, 5756, 6550, 6565, and 6585 Å lines, the respective values for $`f(\lambda `$) are 0.0856, 0.082, 0, $``$0.0153, $``$0.0225, $``$0.1237, $``$0.2185, $``$0.2202, and $``$0.2226, (Martin et al. 1996). This leads to,

$$c(H\beta )=4.541log[F(H\alpha )/F(H\beta )]2.086.$$
(2)
The correction for extinction/reddening from observed to intrinsic flux for the 4364 and 5008 lines is then given by,

$$I(4364)=F(4364)10^{1.082c(H\beta )};I(5008)=F(5008)10^{0.9775c(H\beta )},$$
(3)
and for the 5756 and 6585 lines by,

$$I(5756)=F(5756)10^{0.8763c(H\beta )};I(6585)=F(6585)10^{0.7774c(H\beta )}.$$
(4)
For the STIS data, we binned the pixels along the slit into tiles that are 0.5′′ square (matching the slit width). The STIS fiducial bars are excluded. This produced $`c`$(H$`\beta `$) results for from 79 to 96 tiles, depending on the slit/visit observed, which are the same set used later for the $`T_e`$ analysis. The distributions of $`c`$(H$`\beta `$) along the two slits (slit 1 and slit 2) observed in 3 separate visits (V2, V52, and V72) are remarkably similar. Furthermore, the distributions match well for the regions of slit 1 and slit 2 that are adjacent to each other. Recall that the separation between the bottom of slit 1 and the top of slit 2 is a mere 0.211′′ and that slit 2 is shifted 18.056′′ relative to slit 1 along the slit spatial direction. The average $`c`$(H$`\beta `$) values without regard to any weighting for brightness are: for slit 1: 0.680, 0.699, and 0.701 respectively for the three visits; for slit 2: 0.762, 0.798, and 0.766. For slit 4, it is 0.566 and for slit 5, it is 0.611. We comment further on slit 1, which covers the proplyd P159-350 and 1SW. Starting at the SE end of slit 1, the values are roughly flat at $``$0.85. There is a spike to 1.43 in V2 at the tile containing most of the proplyd. The behavior is similar for V72 with a spike to 1.33 at the tile containing most of the proplyd. The adjacent nebula has for both visits a $`c`$(H$`\beta `$) value of $``$0.9. Unfortunately, because P159-350 was centred in the East fiducial bar on V52, we cannot reliably measure $`c`$(H$`\beta `$) there. Where slit 1 crosses 1SW, $`c`$(H$`\beta `$) is $``$0.6, which agrees well with our previous spectroscopic value 0.605 at our position observed with the Faint Object Spectrograph (FOS-1SW, Rubin et al. 1998). At the end of slit 1 toward the NW, $`c`$(H$`\beta `$) has decreased to $``$0.4.

We examined the vicinity of P159-350 in more detail using the original resolution of 1 pixel (0.05′′) for slit 1 in the spatial direction. The increase in flux and thus signal-to-noise (S/N) in the vicinity of the proplyd permits meaningful analysis here. For V2, we find that the peak observed flux occurs in pixel 311 for the H$`\alpha `$ (5.64$`\times `$10-11 erg cm-2 s-1 arcsec-2) H$`\beta `$ (6.83$`\times `$10-12), and 4364 (3.18$`\times `$10-13) lines, while the 5008 line reaches the highest peak (6.71$`\times `$10-12) at pixel 317 and a second relative maximum (5.14$`\times `$10-12) at pixel 305 with a clear trough between these, having a relative minimum (3.58$`\times `$10-12) occurring at pixel 310. The F(H$`\alpha `$)/F(H$`\beta `$) ratio peaks at 8.269 also at pixel 311 and hence the derived $`c`$(H$`\beta `$) reaches a peak value there of 2.08.

Repeating the analysis for V72, we find that the peak observed flux occurs in pixel 343 for the H$`\alpha `$ (2.19$`\times `$10-11) and 4364 (8.11$`\times `$10-14) lines. The H$`\beta `$ peaks in pixel 344 (3.19$`\times `$10-12) and is lower at pixel 343 (2.74$`\times `$10-12). The 5008 line reaches the highest peak (5.77$`\times `$10-12) at pixel 350 and a second relative maximum (4.58$`\times `$10-12) at pixel 336, again with a distinct trough between these, having a relative minimum (2.85$`\times `$10-12) occurring at pixel 344. The F(H$`\alpha `$)/F(H$`\beta `$) ratio peaks at 7.994 at pixel 343 where the derived $`c`$(H$`\beta `$) reaches a peak value of 2.01. The fact that the peak surface brightnesses are much higher for P159-350 in V2 compared with V72 is most likely because the brightest part of the proplyd was better sampled and/or aligned with the 1 pixel $`\times `$ 10 pixel (0.05′′ by 0.5′′) “smallest rectangular aperture” in the former visit than the latter.

## 4 ELECTRON TEMPERATURE DETERMINATION

### 4.1 \[O III\] Electron Temperature

The electron temperature $`T_e`$ is derived from the intrinsic ratio I(5008)/I(4364) using the following relation,

$$T_e=32966/[ln(I(5008)/I(4364))1.701].$$
(5)
Effective collision strengths are from Burke, Lennon & Seaton (1989) for $`T_e`$ = 104 K. Transition probabilities (A-values) are from Froese Fischer & Saha (1985). Note that this holds in the low-$`N_e`$ limit, which should be valid for Orion where $`N_e`$ values are less than the critical densities ($`N_{crit}`$) for these lines. The lowest $`N_{crit}`$ $``$6.4$`\times `$105 cm-3 for the 5008 line, which is well above $`N_e`$ values determined (e.g., Peimbert & Torres-Peimbert 1977; Osterbrock, Tran & Veilleux 1992; Esteban et al. 1998).

We continue with the analysis using the tiles along the slit described above; the number of usable tiles varies with the slit/visit and ranges from 79 to 96. Equations (2), (3), and (5) are applied to the four emission lines to derive $`T_e`$. Figure 2(a) shows the distribution of $`T_e`$ versus position along slit 1. For slit 1, position along the slit is relative to the location of peak surface brightness (S.B.) in the H$`\alpha `$ line, which occurs at P159-350. We set our zero corresponding to coordinates, $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}15\stackrel{\mathrm{s}}{\mathrm{.}}94`$, $`5^\mathrm{o}`$23′50$`\stackrel{}{\mathrm{.}}`$04 (Bob O’Dell, private communication). As mentioned in the last section, there is actually a relative minimum in the \[O III\] 5008 Å S.B. close to that position. The lower curve in Figure 2(a) shows the observed 5008 S.B. (erg cm-2 s-1 arcsec-2), which is displayed (unsmoothed) at the pixel level providing 0.05′′ spatial resolution along the slit. The feature at the very left (SE) end of this slit is associated with the prominent arc, part of which we used to define our position x2 and which was better observed in slit 2. The open circles on the top ($`T_e`$) curve represent the individual tiles plotted at their midpoint. The dashed lines are a linear interpolation across the tiles that were deemed to have unreliable measurements because of proximity to the fiducial bars. There is a remarkably flat $`T_e`$ distribution with the notable exception of an upward spike to 13600 K at the tile containing the bulk of the proplyd H$`\alpha `$ emission. For slit 1 and V2, the behavior is similar with a sharp upward spike to 15050 K at the tile containing the bulk of the proplyd H$`\alpha `$ emission. It is most unlikely that these high temperatures are realistic (these tiles are omitted in the statistical analysis to follow) and are the result of equation (5) not accounting for the higher electron densities, probably in excess of 106 cm-3, in P159-350. Because the 5008 line will then suffer considerable collisional deexcitation, $`T_e`$ derived using equation (5) will be overestimated (e.g., Viegas & Clegg 1994). In the next section, we will evaluate the distribution of $`T_e`$ along slit 1, as well as along the other slits, in terms of $`T_e`$ variations.

Figure 2(b) shows the distribution of $`T_e`$ versus position along slit 2, where position along the slit is relative to the location of the local peak S.B. in the 5008 line that occurs at x2. Using the position of P159-350 and the specified offset from slit 1 to slit 2 of $`+`$1.08s in RA and $``$8.15′′ in Dec., we find x2 coordinates: $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}16\stackrel{\mathrm{s}}{\mathrm{.}}96`$, $`5^\mathrm{o}`$23′57$`\stackrel{}{\mathrm{.}}`$73. The lower relative maximum near offset +14′′ is where slit 2 passes through the “downstream” tail of P159-350. Again there is a remarkably flat $`T_e`$ distribution.

Figure 3(a) plots the distribution of $`T_e`$ versus position along slit 4, where position along the slit is relative to the location of a feature in the bar seen in the 6585 S.B. distribution (Fig. 3(c)) discussed in the next section. Recall that this slit is aligned to point toward $`\theta ^1`$ Ori C. While there is an underlying decrease in the 5008 S.B. with increasing distance from $`\theta ^1`$ Ori C, there is a leveling off (plateau) at the SE end of slit 4 of 5008 emission that remains substantial well beyond the bar. This is consistent with the WFPC2 images of the region; the composite “drizzled” image by Walsh (1998) that includes filter F502N (which covers the 5008 line) depicted in green, shows a green hue in a parabolic shape that appears to be a wake emanating from HH 203 and HH 204. Walsh’s image is available at stecf.org/newsletter/webnews1/orion/m42col\_drizzle.jpg. The point was also made in an earlier discussion of HH 203 and HH 204 (O’Dell et al. 1997). Along this slit, the $`T_e`$ distribution is also fairly flat, especially in the region of higher S/N (see §4.2). There is a clear increase in the amplitude of the $`T_e`$ fluctuations for the tiles SE of the bar due to poorer S/N in the 4364 Å line.

In Figure 3(b), we plot the distribution of $`T_e`$ versus position along slit 5. Here, the position along the slit is relative to the location of a feature seen in the 6585 S.B. distribution (Fig. 3(d)) discussed below. There is roughly a linear decrease in 5008 emission from the NW to the SE end of slit 5 with a “jump” occurring between 18–22′′ offset position. As was the case for slit 4, here too there remains substantial 5008 emission beyond the bar that is evident in the green hue at this location in the Walsh image. The $`T_e`$ distribution here is flat similar to that for slit 4; the increase in the $`T_e`$ amplitude range at the SE end of the slit is again due to poorer S/N in the 4364 Å line.

### 4.2 \[N II\] Electron Temperature

The $`T_e`$ in the N+ zone is derived from the intrinsic ratio I(6585)/I(5756) (see equations (2) and (4)). Because the critical density for the 6585 line, $``$7.7$`\times `$104 cm-3 (at 104 K), is substantially less than $`N_{crit}`$ for the \[O III\] 5008 line ($``$6.4$`\times `$105 cm-3), we do consider various $`N_e`$ values when deriving $`T_e`$. We also derive temperatures utilizing two different sets of N II effective collision strengths – those calculated by Lennon & Burke (1994) and by Stafford et al. (1994). We use the effective collision strengths for 10000 K; these do not vary much with the $`T_e`$ range of interest in Orion. The A-values used are discussed in Rubin et al. (1998, Appendix A with original references therein).

For non-zero densities, our derivation of $`T_e`$ is done in an iterative fashion, starting with an initial estimate in the low-$`N_e`$ limit. Then the volume emissivities ($`j`$ values) for both the 5756 and 6585 lines are calculated solving the statistical equilibrium equations for the six lowest energy levels. The $`T_e`$ value is then recomputed using the intrinsic ratio I(6585)/I(5756), $`j`$(6585), and $`j`$(5756). The rapidly converging iteration is halted when $`T_e`$ changes by less than 1 K.

In this paper, we are interested in assessing the amount of $`t^2`$ that can occur. Thus, we limit our analysis here by not including a detailed study of density variations, which is the subject of future papers by the present authors. Nevertheless, it is necessary to include some discussion of $`N_e`$ with regard to deriving the \[N II\] $`T_e`$. As is well known, there is an inverse scaling of \[N II\] $`T_e`$ with $`N_e`$. First we perform the calculations for $`T_e`$ in the low-$`N_e`$ limit. This provides a firm upper limit to the $`T_e`$ distributions evaluated for the various slit/visit data. We also repeat the full computations for 4 other $`N_e`$’s: 1000, 2000, 5000, and 10000 cm-3. The last value is more than high enough to bracket expectations for Orion, with the exception of objects like the proplyds as will be discussed below for P159-350. The value of 5000 is our best single $`N_e`$ to cover the region of slits 1 and 2 (e.g., Pogge, Owen & Atwood 1992). Their map of the electron density in the S+ zone indicates that $`N_e`$ is lower in the regions of our slit 4 and slit 5; for these, we adopt a best single $`N_e`$ = 2000. It is these respective densities that are used to depict the distribution along slits 1, 2, 4, and 5 of \[N II\] $`T_e`$ in Figure 2(c), 2(d), 3(c), and 3(d) respectively. For the Figures, we have used the N+ collision strengths of Lennon & Burke (1994) while results from this set, as well as Stafford et al. (1994) will be tabulated later. The lower curve in each of these panels shows the 6585 surface brightness distribution along the slit.

Figure 2(c) shows the tiled \[N II\] $`T_e`$ distribution versus position along slit 1 and is a positional match to Figure 2(a). The lower curve here shows the observed 6585 Å S.B. at the unsmoothed pixel level. There again appears to be a relative minimum in the \[N II\] 6585 S.B. at the centre pixel for P159-350 although it is not as well defined as the dip in the 5008 S.B. there. As with the 5008 curve, the peak S.B. on the NW side of the proplyd exceeds that on the SE side. However, the overall 6585 emission is much narrower in P159-350 than is the 5008 emission. In totality, the structural behavior in the 6585 emission compared with that in 5008 is consistent with the N+ region being more tightly confined than is the O++ region to the low-mass star of P159-350, which is expected if the proplyd is externally ionized by $`\theta ^1`$ Ori C. The prominent absolute maximum of the 6585 S.B. is very close to the 1SW position, 20.26′′ from P159-350, which is somewhat blocked by the NW-side fiducial bar. The narrow secondary peak at $``$13′′ offset appears to be associated with the much broader secondary maximum near offset $``$12′′ in the 5008 S.B. of Figure 2(a). There is perhaps an ionization boundary being observed here, although the slit orientation is far from ideal to test this. A slit pointing toward $`\theta ^1`$ Ori C would be better suited. Once again, there is a notably flat $`T_e`$ distribution with the exception of the proplyd. In Figure 2(c), we have omitted the tile with the highest derived $`T_e`$ of 20815 K. The highest $`T_e`$ in the plot is at 13000 K. Both these tiles are omitted from the subsequent statistical analysis. These high derived $`T_e`$’s are undoubtedly pointing more to high $`N_e`$ values associated with the proplyd than high temperatures. The $`N_e`$ assumed for the calculation of $`T_e`$ in Figure 3(c) is 5000 cm-3, while $`N_e`$ for P159-350 may exceed 106 cm-3. Because the 6585 line emission will suffer enormously from collisional deexcitation at more realistic densities, the derived $`T_e`$ would be much lower. For slit 1, V2, there are two tiles at the proplyd that have an enormous upward spike in $`T_e`$ to 29800 and 35850 K; both tiles are omitted from the statistical analysis in following sections.

Figure 2(d) has the distribution of \[N II\] $`T_e`$ versus position along slit 2. There is a fairly flat $`T_e`$ distribution with some indication of a drop in $`T_e`$ SE of the lower fiducial bar. The maximum 6585 S.B. in the lower curve occurs near offset 28′′. This location is immediately adjacent to the narrow peak at $``$13′′ offset in Figure 2(c), just mentioned. The relatively small peak at x2 (offset 0) further supports the finding in the WFPC2 composite image that the emission arc which includes the position x2 is mainly an \[O III\] feature (see Figure 1). We note that x2 is part of the HH 529 complex and more specifically, the easternmost feature called 170-358 (see fig. 20 in Bally, O’Dell, & McCaughrean 2000). According to their study, including the kinematics, 170-358 is an expanding bow shock seen from the side.

Figure 3(c) plots the distribution of \[N II\] $`T_e`$ versus position along slit 4. Position along this slit is relative to the location of the SE-most peak in S.B. in the 6585 line. We set our zero corresponding to coordinates: $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}21\stackrel{\mathrm{s}}{\mathrm{.}}34`$, $`5^\mathrm{o}`$24′48$`\stackrel{}{\mathrm{.}}`$21 determined from the position of our offset star $`\theta ^2`$ Ori A (HD37041) with $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}22\stackrel{\mathrm{s}}{\mathrm{.}}9`$, $`5^\mathrm{o}`$24′57$`\stackrel{}{\mathrm{.}}`$9. There is a brighter feature that occurs where the slit crosses HH 203 at offset near $`20^{\prime \prime }`$. Comparison with Figure 3(a) shows no enhancement in 5008. This is consistent with the WFPC2 images of the region; the composite image in Figure 1 that includes filter F658N (isolating the 6585 line) shows HH 203 as a reddish-orange object. There is a very flat $`T_e`$ distribution even over the wide range of S.B. and ionization structure along slit 4.

In Figure 3(d), we plot the distribution of $`T_e`$ versus position along slit 5. Here, the zero position along the slit is the SE-most 6585 S.B. peak of where there is a bifurcation in the bar (see Figure 1). This feature appears to mark the boundary between the ionized region and the photodissociation region (PDR). Using the position of $`\theta ^2`$ Ori A and the specified offset from slit 4 to slit 5 of $`1.87`$s in RA and $`15.32^{\prime \prime }`$ in Dec., we find: $`\alpha `$, $`\delta `$ = $`5^\mathrm{h}35^\mathrm{m}19\stackrel{\mathrm{s}}{\mathrm{.}}74`$, $`5^\mathrm{o}`$25′08$`\stackrel{}{\mathrm{.}}`$14. Beyond this position, the $`T_e`$ distribution becomes less reliable because of poorer S/N in the 5756 Å line. Overall, the $`T_e`$ is flat as has been the case throughout.

## 5 FRACTIONAL MEAN-SQUARE TEMPERATURE VARIATIONS

Our STIS analysis above presents results in the plane of the sky. The observations here do not address temperature fluctuation along the line of sight, which may be characterized in terms of the average temperature $`T_0`$ and fractional mean-square $`T_e`$ variation ($`t^2`$) as defined by Peimbert (1967).

$$T_0=\frac{T_eN_eN_i𝑑V}{N_eN_i𝑑V},$$
(6)
$$t^2=\frac{(T_eT_0)^2N_eN_i𝑑V}{T_0^2N_eN_i𝑑V},$$
(7)
where $`N_i`$ is the ion density $`N(N^+`$) or $`N(O^{++}`$). The integration in equations (6) and (7) is over the column defined by each tile, and along the line-of-sight ($`los`$). We are unable to measure the $`t^2`$ along the $`los`$ for any column (cross section 1 tile). If there are $`t^2`$ along the $`los`$, we can say that $`T(4364/5008)`$ $`>`$ $`T_0`$ or $`T(5756/6585)`$ $`>`$ $`T_0`$ (e.g., Peimbert 1967; Rubin et al. 1998).

In the case of \[O III\], for each tile, we have calculated $`T(4364/5008)`$. Then the intrinsic flux I(5008), fully correcting F(5008) for extinction (see equ. 3), in each tile is used in conjunction with $`T_e`$ = $`T(4364/5008)`$ for that tile, and assumed constant along the $`los`$, to derive the following:

$$I(5008)=K(5008)N_eN_iT_e^{0.5}exp(\chi /kT_e)𝑑l=K(5008,T_e)N_eN_i𝑑l.$$
(8)
Here $`\chi `$ is the excitation energy above the ground state for the upper level of the 5008 transition, k is Boltzmann’s constant, $`K(5008)`$ is known from atomic data, and $`K(5008,T_e)`$ has finally incorporated the known $`T_e`$ factor with the atomic constants. Here we again make the safe assumption of the low-Ne limit (negligible collisional deexcitation) discussed earlier.

For the case of \[N II\], for each tile, we have calculated $`T(5756/6585)`$. Then the intrinsic flux I(6585), fully correcting F(6585) for extinction (see equ. 4), in each tile is used in conjunction with $`T_e`$ = $`T(5756/6585)`$ for that tile and the chosen value for $`N_e`$, and assumed constant along the $`los`$, in the following relation:

$$N_eN_i𝑑l\frac{I(6585)}{ϵ(6585)}.$$
(9)
Here $`ϵ(6585)`$ is the normalized volume emissivity which is related to the usual volume emissivity $`ȷ(6585)`$ by $`ϵ(6585)`$ $``$ $`ȷ(6585)`$/($`N_e`$$`N_i`$). $`ϵ(6585)`$ depends on the (fractional) population in a given level obtained by solving the 6-level atom for the specific $`T_e`$ and $`N_e`$. It is not necessary to deal with the constant of proportionality for our purposes.

Following Paper I, we define the average $`T_e`$ ($`T_{0,A}`$) and fractional mean-square $`T_e`$ variation ($`t_A^2`$) in the plane of the sky.

$$T_{0,A}=\frac{T_eN_eN_i𝑑l𝑑A}{N_eN_i𝑑l𝑑A},$$
(10)
$$t_A^2=\frac{(T_eT_{0,A})^2N_eN_i𝑑l𝑑A}{T_{0,A}^2N_eN_i𝑑l𝑑A}=\frac{T_e^2N_eN_i𝑑l𝑑A}{T_{0,A}^2N_eN_i𝑑l𝑑A}1,$$
(11)
where $`dA`$ represents an element of surface area in the plane of the sky and the integration over $`dl`$ is for each tile along the $`los`$. The proper weighting for each tile is provided by either equation (8) or equation (9) depending on whether we are performing an \[O III\] or \[N II\] analysis.

The results using the STIS data for the various slits and visits to determine $`T_{0,A}`$ and $`t_A^2`$ are summarized in Table 1. These quantities for \[O III\] are entered in the first column. For slit 1, Visit 72 (see Figure 2(a)), $`T_{0,A}`$ = 8258 K and $`t_A^2`$ = 0.00925. We have excluded the spike in $`T_e`$ at the proplyd. The value for $`t_A^2`$ is remarkably low, although not as low as the value of 0.0035 that we found for NGC 7009 (Paper I). The results from the other visits (V2 and V52) for slit 1 are close to those for V72. Our results for slit 2 and all three of the visits show $`T_{0,A}`$ values in close accord and similar to those for slit 1. All these averages vary by less than 285 K. By comparing $`t_A^2`$ for slits 1 and 2 in all 3 visits, we see that all six numbers are small, ranging from 0.00682 to 0.0129. For slits 4 and 5 the statistics are roughly comparable and differ from those for slits 1 and 2 in the following ways: $`T_{0,A}`$ is several hundred K smaller and $`t_A^2`$ is somewhat higher, although the largest value, 0.0176, is still notably small. A simple visual comparison of the $`T_e`$ plots in Figs. 2(a), 2(b), 3(a), and 3(b) leads one to conclude that $`t_A^2`$ looks larger for slits 4 and 5, particularly with the larger $`T_e`$ excursions in the regions SE of the bar. We note however that these data with poorer S/N have a minor effect on the statistics. There is natural biasing against these tiles due to weaker emission in both 4364 and 5008 as can be seen from equation (8), which in turn provides a smaller ($`N_eN_i𝑑l`$)–weighting of these tiles in equations (10) and (11).

All the remaining columns in Table 1 pertain to the derivation of $`T_{0,A}`$ and $`t_A^2`$ for \[N II\]. These are displayed for 5 different electron densities and two different sets of effective collision strengths. As described, for slits 1 and 2, our preferred $`N_e`$ is 5000 cm-3, while for slits 4 and 5, it is 2000 cm-3. We will now discuss the resulting numbers in these particular density columns and also using the left-side column with the Lennon & Burke (1994) collision strengths, unless stated otherwise. For slit 1, Visit 72 (see Figure 2(c)), $`T_{0,A}`$ = 10226 K and $`t_A^2`$ = 0.00695. We exclude the high $`T_e`$ near P159-350 (see §4.2). The value for $`t_A^2`$ is even lower than found for this slit/visit in the \[O III\] analysis; $`T_{0,A}`$ is nearly 2000 K higher. There is close agreement among the N+ $`T_{0,A}`$ and $`t_A^2`$ values for all the visits for slit 1; the above intercomparison between N+ and O++ with the corresponding numbers in the respective visits is also true.

For slit 2, the $`T_{0,A}`$ numbers vary by more than 400 K between visits; the respective $`T_{0,A}`$’s are smaller than for slit 1 for V52 and V72 but slightly higher for V2. All of the $`t_A^2`$ values are larger for slit 2 than for slit 1, reaching 0.0146.

For slits 4 and 5, the N+ $`T_{0,A}`$ is several hundred to more than 1000 K smaller than for slits 1 and 2. Compared with the O++ $`T_{0,A}`$, slits 4 and 5 are higher by 1343 and 1749 K, respectively. $`t_A^2`$ at 0.0175 is largest for slit 5 and nearly identical to the corresponding O++ value, which is also the highest in the first column. This larger $`t_A^2`$ may be attributed to the larger excursions in $`T_e`$ in those tiles SE of the bar, where the S/N in the 5756 Å line is poorer. A comparison of Figure 3(c) and 3(d) clearly indicates a much lower 6585 S.B. SE of the bar for slit 5 than for slit 4. Again we emphasize that these data with poorer S/N play only a minor role in contributing to the integrals in equations (10) and (11). When there is weaker emission in both 5756 and 6585, as can be seen from equation (9), there is a smaller ($`N_eN_i𝑑l`$)–weighting of these tiles in equations (10) and (11).

The effect of substituting the Stafford et al. (1994) collision strengths is small for both slit 1 and slit 2 for all 3 visits: $`T_{0,A}`$ decreases by less than 200 K; $`t_A^2`$ decreases by roughly 5–6 percent. The $`T_{0,A}`$ drop for slits 4 and 5 is slightly more, up to a 276 K difference; $`t_A^2`$ decreases by $``$6–7 percent.

On the whole, the results presented in Table 1 allow us to reach two robust conclusions. First, $`t_A^2`$ is always very small ($``$ 0.0182) for either the O++ or N+ analyses. Second, the $`T_{0,A}`$ derived from the \[O III\] lines is always significantly less than $`T_{0,A}`$ derived from the \[N II\] lines. These hold even when we consider for the N+ analysis a wide range in $`N_e`$ that should amply bracket the great bulk of the plasma in Orion and also examine the influence of using alternative collision strengths.

In Appendix A, we present an analysis of the uncertainties in the derived $`T_e`$ values as well as the $`t_A^2`$ values.

## 6 DISCUSSION AND CONCLUSIONS

We extend our earlier study of electron temperature variations from the \[O III\] (4364/5008) flux ratio in the PN NGC 7009 (Paper I) to the Orion Nebula. In addition, we use our STIS long-slit observations to expand the analysis to measure $`T_e`$ variations from the \[N II\] (5756/6585) flux ratio. The observations here do not address $`T_e`$ fluctuation along the line of sight through the specific O++ region or, likewise, through the N+ region. We analyse both the \[O III\] and \[N II\] data sets to derive the average $`T_e`$ and fractional mean-square $`T_e`$ variations in the plane-of-the-sky, which we call $`T_{0,A}`$ and $`t_A^2`$. We assume for each square column (projection of 1 STIS tile 0.5′′ square on the plane-of-the-sky) that the plasma along the line of sight is isothermal at the $`T(4364/5008)`$ in the case of O++ or $`T(5756/6585)`$ in the case of N+. For the latter case, we consider a large range of $`N_e`$ values, which should be sufficient for Orion, in order to produce Table 1. The analysis for each $`N_e`$ assumes that it is constant for the entire STIS slit length and throughout the sheet projected through the nebula along the $`los`$.

Fluctuations in $`T_e`$ (and $`N_e`$) along the $`los`$ are inevitable. We can make some comments about how our results for $`t_A^2`$ might be adjusted by $`T_e`$ variations along the $`los`$. The relationship between $`T(4364/5008)`$ and $`T_0`$ for the O++ region is

$$T(4364/5008)=T_0[1+0.5(91200/T_03)t^2],$$
(12)
and between $`T(5756/6585)`$ and $`T_0`$ for the N+ region is

$$T(5756/6585)=T_0[1+0.5(69000/T_03)t^2],$$
(13)
(e.g., Peimbert 1967; Rubin 1969). With the $`T(4364/5008)`$– or $`T(5756/6585)`$–values for Orion (see Figures 2 and 3) or indeed for H II regions in general, $`T_0`$ will be smaller than these temperatures inferred from the forbidden line flux ratios here. Rather than repeat some examples of how the differences depend on $`t^2`$, the reader is referred to our comments in the earlier paper (§7 in Paper I).

We do not have the data here to characterize $`T_e`$ variations in 3-dimensions (3-D). It is useful to define an overall 3-D average $`T_e`$ ($`T_{0,V}`$) and fractional mean-square $`T_e`$ variation ($`t_V^2`$). These single values apply for the entire source. Equations (6) and (7) define these specific values when the integration is over the entire volume. We note that for a spatially unresolved object (total integrated fluxes observed in the aperture), not the case for Orion, $`t_V^2`$ = $`t^2`$ and a calculation of $`t_A^2`$ is meaningless.

Our measurements of $`T_e`$ reported here are an average along each line of sight. Because each element of area treated in the plane of the sky represents a column which has already created a spatially averaged temperature along the $`los`$ (e.g., see fig. 1 in Rubin 1969), it is likely that the value for $`t_V^2`$ is substantially higher than $`t_A^2`$. Measurements of $`t^2`$ along various sight lines appear to be the most direct way to reliably gauge $`t_V^2`$. Therefore, despite finding remarkably low $`t_A^2`$, we cannot completely rule out much larger temperature fluctuations along the $`los`$. Further work, beyond the scope of this paper, is underway that will use modeling as well as additional observational data in an effort to better determine the relationship between $`t^2`$, $`t_A^2`$, and $`t_V^2`$.

Finding \[N II\] temperatures that are higher than \[O III\] temperatures is not a new result. For Orion, this has been known for years (e.g., Baldwin et al. 1991). What is new/significant here is that we have results from many more sight lines, more than 700 represented by the individual tiles, with improved spatial resolution (0.5 arcsec squares) and excellent spatial registration between the \[N II\] and \[O III\] data sets that HST affords. Of the roughly 700 tiles, about half represent independent sight lines, with the other half having spatial overlap due to the repeated visits for slits 1 and 2.

The fact that the \[N II\] temperatures are higher is also expected on theoretical grounds and again, not something novel that we have uncovered. We enumerate three factors that contribute to the predicted inequality. First, the cooling of the \[O III\] 5008, 4960 Å lines in the O++ region is uniquely efficient for nebular conditions and elemental abundances that prevail in Galactic H II regions, including Orion. On the other hand, in the “singly ionized” region (where the dominant O and N ions are O+, N+) there is no coolant that is nearly as efficient as \[O III\] 5008, 4960 is in the O++ zone. Additionally for Orion, the blister geometry likely results in higher average electron densities in the singly ionized region compared with the O++ region because the former is closer to the PDR. The higher $`N_e`$ in the N+ region would also contribute to less efficient cooling.

Second, there is the predicted hardening of the stellar ionizing photons at progressively larger distances from the exciting star. This is due to the functional form of the H photoionization (predominantly) cross section, which diminishes steeply with higher frequency from its value near threshold. At larger distances from the exciting source, the average energy per photoionization will increase and thus the heating rate will also increase. This causes a rise in $`T_e`$ with an increase in distance (other factors being equal). This has been known for many years (e.g., Rubin 1968).

Third, there is a possibility that there is a contribution to the production of \[N II\] 5756 Å emission (and to a lesser extent 6585) by recombination and cascading. Under some conditions, this may provide significant routes into the upper energy level of the 5756 transition that are not negligible compared with collisional processes. This was examined by Rubin (1986). If there were a significant “recombination” contribution to the observed emission, we would have overestimated the $`T_e`$ derived from the \[N II\] 5756/6585 flux ratio. With our two independent, detailed models for the Orion Nebula (Baldwin et al. 1991 and Rubin et al. 1991), the recombination contribution appears to be small. We discussed this in a subsequent paper where these models were “retrofitted” with more in-common input parameters, including the Stafford et al. (1994) N+ collision strengths, to allow closer intercomparisons (see end of §4 in Rubin et al. 1998).

All three of the above effects are included in our photoionization, plasma simulation models just mentioned. The predicted \[N II\] and \[O III\] temperatures for the entire volume are 8721 K and 7704 K in the Rubin et al. (1991) model and 8649 K and 7692 K in the retrofit (Rubin et al. 1998) using Stafford et al. cross sections. When the Lennon & Burke (1994) N+ collision strengths are used instead of the Stafford et al. set, the respective temperatures become 8706 K and 7701 K.

Finally, we encapsulate the findings here for the proplyd P159-350. We find large local extinction as evidenced by the dramatic increase in the observed F(H$`\alpha `$)/F(H$`\beta `$) ratio along slit 1. For V2, this ratio peaks at 8.27 which implies a $`c`$(H$`\beta `$) of 2.08 from the Balmer decrement. Similar values are found for our V72 observations: F(H$`\alpha `$)/F(H$`\beta `$) peaks sharply at the proplyd position reaching a ratio of 7.99, where the derived $`c`$(H$`\beta `$) = 2.01. Because the adjacent nebular values for $`c`$(H$`\beta `$) are much lower, the extinction must be associated with the P159-350 environs.

A comparison of the \[O III\] 5008 and \[N II\] 6585 surface brightnesses in the vicinity of P159-350 (see Figures 2(a) and 2(c)) shows for both lines that the S.B. on the NW side of the proplyd exceeds that on the SE side. However, the 6585 emission is much narrower than is the 5008 emission which provides evidence of ionization stratification consistent with the N+ region being more tightly confined than is the O++ region to the low-mass star of P159-350. This “inverse H II region” behavior is precisely what is expected due to the external ionization source $`\theta ^1`$ Ori C.

The derived $`T_e`$ distributions were all notably flat except at the position of the proplyd observed in slit 1. The very high derived $`T_e`$’s at the location of P159-350 are no doubt due more to high $`N_e`$ values associated with the proplyd than high temperatures. As stated, the calculations used here to assess $`T_e`$ have not accounted for $`N_e`$ values likely to exceed 106 cm-3 in P159-350. Because both the 5008 line and even more so the 6585 line emission will suffer substantial collisional deexcitation at such high densities, the derived $`T_e`$ for both \[O III\] and \[N II\] would be much lower.

Each of the above mentioned facets of P159-350 deserves further study. The smallest spatial resolution element was 1 pixel (0.05′′) in the spatial direction by 10 pixels (0.5′′) in the dispersion direction. It would be highly desirable to obtain higher spatial resolution of P159-350, which can be readily achieved by observing with a narrower STIS slit. The narrower slit would enable better spectral resolution as well and permit a detailed mapping of P159-350 in the C III\] 1907, 1909 Å lines. The ratio of their fluxes is an excellent diagnostic of $`N_e`$ even at the high $`N_e`$ values expected in proplyds. Improved spectral resolution in the optical lines will also help sort out kinematics from spatial structure, which was a difficulty with the 0.5′′ slit width we used.

This paper is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (ST ScI). Support for proposal #GO-7514 was provided by NASA through a grant from ST ScI. ST ScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Valuable contributions were made by Naman Bhatt, Ephrat Bitton, Brent Buckalew, Bob O’Dell, Chris Ortiz, Aaron Svoboda, and Abby Wong. We thank Bob O’Dell and Janet Simpson for reading the paper and providing useful comments. We much appreciate the exemplary service of our ST ScI contact scientist Paul Goudfrooij and program coordinator Ray Lucas. RHR acknowledges support from the Long-Term Space Astrophysics (LTSA) program, NASA/Ames Research Centre contract NCC2-9018 with Orion Enterprises, and thanks Scott McNealy for providing a Sun workstation.

## Appendix A Error Analysis

The values found for $`t_A^2`$ are very small. It is important to determine whether even these small fluctuations in $`T_e`$ from tile to tile are a real signal or the result of measurement uncertainty.

From equation (5) for O++, it can be shown that the mean square uncertainty in temperature, expressed as $`t_e^2`$, has the following dependence on the uncertainties $`\delta I_n`$ and $`\delta I_a`$ of the nebular ($`I_n`$) and auroral ($`I_a`$) line fluxes, respectively:

$$t_e^2=\left(\delta T_e/T_e\right)^2=\frac{T_e^2}{c_{1}^{}{}_{}{}^{2}}\left[\left(\frac{\delta I_a}{I_a}\right)^2+\left(\frac{\delta I_n}{I_n}\right)^2\right],$$
(A1)
where $`c_1=32966`$. We use this as an approximation for N+ too, with $`c_1=24933`$.

The uncertainty $`\delta I`$ is principally the measurement error in the line flux for a tile. The S/N for the flux of the auroral lines (\[O III\] 4364 and \[N II\] 5756 here) is much lower than that for the nebular lines (\[O III\] 5008 and \[N II\] 6585 here), and so the uncertainty of the weaker auroral line will dominate the uncertainty $`\delta T_e`$. The uncertainty $`\delta I`$ also depends on the (differential) extinction correction, but in the present case this is a small effect compared to the uncertainty in the observed auroral line fluxes (hence, $`I/\delta IF/\delta F=`$ S/N ).

The three independent methods used to calculate the S/N of the auroral lines are discussed in the following sections. We use $`T_eT_{0,A}`$ to calculate the corresponding $`t_e^2`$.

### A.1 Method 1

The flux integrated over the line profile as found using the “blkavg sum” technique in IRAF is

$$F^{line}=\underset{i=1}{\overset{m}{}}DF_i^{line+cont}\frac{m}{n}\underset{j=1}{\overset{n}{}}DF_j^{cont},$$
(A2)
where $`F_i`$ is the monochromatic flux for the tile (erg cm-2 s-1 Å-1), $`D`$ is the dispersion conversion (Å/pixel), $`n`$ is the number of pixels in the spectrum used to define the $`continuum`$ and $`m`$ is the number of pixels used for the $`line+cont`$. Since the auroral line is weak, the rms fluctuation in the $`line+cont`$ can be assumed to be the same as in the $`continuum`$, $`\delta F_j^{cont}`$. Thus, the rms of the line flux is

$$\delta F^{line}=D\sqrt{m}(\delta F_j^{cont})\sqrt{1+m/n}.$$
(A3)
In our measurements we used $`m=n`$. The above three equations are used to create the entries for method 1 (M1) in Table 2. Values are given for both the \[O III\] 4364 and \[N II\] 5756 lines for all slits (for slits with multiple visits, a representative case is given).

### A.2 Method 2

In this method, the line profiles from the STIS spectra are fit with a template (a Gaussian convolved with a slit of width 0.5″) using a combined Gauss-Newton and modified Newton algorithm (Numerical Algorithms Group (NAG) - E04FDF) to find the minimum least-squares solution and the corresponding variances of the variables used in the fit (NAG-E04YCF). The slit width is held constant, but due to anamorphic magnification in the dispersion direction, the plate scale differs from line to line which results in different slit widths expressed as pixels (see Bowers & Baum 1998).

The variables used in the template fitting are the flux $`F^{line}`$, the FWHM and central wavelength of the Gaussian, and the slope and mean brightness of the continuum. The variances reflect how well the model fits the data. Keeping all of these as variables gives one a conservative estimate of the variance of each. Note that a poor template would lead to an overestimate of the error. The templates were tested on the strong nebular lines and generally fit very well. However, there are instances where the surface brightness changes across the slit (in the dispersion direction), so that the implicit assumption of uniform illumination is not perfect. Therefore, the estimated uncertainty from method 2 is an upper limit. As is seen in Table 2 it agrees well with the error estimated using method 1.

### A.3 Method 3

This method utilizes the multiple sets of observations for slit 1 and slit 2 in the three visits V2, V52 and V72 and examines the reproducibility of the line fluxes and $`T_e`$. There are adequate fiducial sources visible in all three visits so that the spectra may be shifted in the spatial direction to achieve alignment to $``$1 pixel. The proplyd P159-350 in slit 1 is particularly useful for this and serves for slit 2 as well because the relative positional HST offsets between slit 1 and slit 2 are not subject to the acquisition uncertainty. There is also positional alignment uncertainty in the dispersion direction which is more difficult to estimate. For all three visits, P159-350 was seen in the 0.5′′ wide slit; although it was blocked for the most part by the East fiducial bar in V52, some of the emission still peeks out both ends. We noted at the end of section 3 that P159-350 has a higher peak surface brightness in V2 than in V72. We estimate that the positional alignment of the tiles perpendicular to the spatial direction between the three visits is within 0.2′′ so that there is substantial area overlap (perhaps 60%) of our tiles between the three visits. Comparison of the nebular line strengths from visit to visit shows excellent reproducibility ($`<3\%`$ per tile for the O++ line). This is much smaller than the auroral line differences found among successive visits. The measurement error of the weaker auroral lines is of course higher. Furthermore lack of registration is a more serious issue for the auroral lines (from a much higher energy level) if there are real $`T_e`$ fluctuations on this spatial scale; in the limit of no overlap, one would expect $`t_e^2t_A^2`$ even if there were no measurement errors. Therefore, having $`t_e^2`$ from method 3 larger than from M1 and M2 is not unexpected.

For each tile in common between V2, V52 and V72, both an average value ($`F_a`$) and the standard deviation ($`\delta F_a`$) of any single determination of $`F_a`$ were calculated. The representative value entered as M3a in Table 2 is the median $`\delta F_a/F_a`$ for each of the slits; as in methods 1 and 2, the $`\delta T_e`$ listed is what is implied by this S/N, through equation A1.

For each in-common tile, we also calculated $`T_e`$ for each visit, then the average $`T_e`$ and the actual standard deviation ($`\delta T_e`$) of any single determination of $`T_e`$. The median $`\delta T_e`$ for each slit is included in Table 2 as M3b. From these median values we calculated the tabulated $`t_e^2`$. Methods 3a and 3b agree closely, as they should given the good reproducibility of the nebular lines (i.e., the effect of $`\delta F_n`$ in equation A1 is negligible).

The analysis summarized in Table 2 shows that the measured $`t_A^2`$ is a real signal not dominated by the noise and/or measurement errors. The two exceptions are for the O++ measurements for slits 4 and 5 which include tiles beyond the Orion bar (SE part of the slit) where the surface brightness for lines of that ion, and hence the S/N, is low. Accordingly, we repeated the entire analysis for only the NW half of these slits. The results recorded in Table 2 show a lower $`t_A^2`$, as anticipated, and a lower $`t_e^2`$ as well which is a smaller fraction of $`t_A^2`$.

Where there is a contribution of measurement error to the apparent $`t_A^2`$, the actual $`t_A^2`$ is smaller, by about $`t_e^2`$. This reinforces our conclusion that $`t_A^2`$ is very small.


format: text

text:
# Mechanical Stability of a Strongly-Interacting Fermi Gas of Atoms

## Abstract

A strongly-attractive, two-component Fermi gas of atoms exhibits universal behavior and should be mechanically stable as a consequence of the quantum mechanical requirement of unitarity. This requirement limits the maximum attractive force to a value smaller than that of the outward Fermi pressure. To experimentally demonstrate this stability, we use all-optical methods to produce a highly degenerate, two-component gas of 6Li atoms in an applied magnetic field near a Feshbach resonance, where strong interactions are observed. We find that the gas is stable at densities far exceeding that predicted previously for the onset of mechanical instability. Further, we provide a temperature-corrected measurement of an important, universal, many-body parameter which determines the stability—the mean field contribution to the chemical potential in units of the local Fermi energy.

Strongly interacting Fermi systems are expected to exhibit universal behavior Heiselberg . In atomic gases, such strong forces can be produced in the vicinity of a Feshbach resonance, where a bound molecular state in a closed exit channel is magnetically tuned into coincidence with the total energy of a pair of colliding particles Verhaar . In this case, the zero energy scattering length $`a_S`$, which characterizes the interactions at low temperature, can be tuned through $`\pm \mathrm{}`$. For very large values of $`|a_S|`$, the important properties of the system (e.g. the effective mean field potential, the collision rate, the superfluid transition temperature, etc.) are predicted to lose their dependence on the magnitude and sign of $`a_S`$, and instead become proportional to the Fermi energy with different universal proportionality constants. For this reason, tabletop experiments with strongly interacting atomic Fermi gases can provide measurements that are relevant to all strongly interacting Fermi systems Heiselberg ; Steele , thus impacting theories in intellectual disciplines outside atomic physics, including materials science and condensed matter physics (superconductivity), nuclear physics (nuclear matter), high-energy physics (effective theories of the strong interactions), and astrophysics (compact stellar objects).

In a gas of degenerate fermions, compression of the gas is resisted by the Pauli exclusion principle, leading to an effective pressure known as the Fermi pressure or degeneracy pressure. The Fermi pressure plays an important role in nature, providing, for example, the outward force which stabilizes neutron stars against gravitational collapse. Because of the Fermi pressure, a confined cloud of degenerate fermions is always larger than a cloud of bosons with an equivalent temperature and particle number. This effect has been directly observed in two elegant experiments Truscott ; Salomon .

An important question is: At what point, if at all, do strong attractive interactions overcome the Fermi pressure? Previously, it was predicted that an atomic Fermi gas could become mechanically unstable for sufficiently large attractive interactions Houbiers . If true, this prediction presents a possible roadblock to attempts to use a Feshbach resonance to produce a high-temperature superfluid Holland ; Timmermans ; Kokkelmans , as the gas should be unstable in the required regime. Recently, however, it has become apparent that the gas may indeed be stable. Heiselberg has estimated the maximum ratio of the attractive potential to the local Fermi energy for a strongly interacting Fermi gas. Using a self-consistent many-body approach, he finds that a two-component Fermi gas is mechanically stable, although a gas with more than two components is not Heiselberg . A similar conclusion can also be drawn from an examination of the compressibility of a strongly attractive two-component gas Kokkelmans .

In this Letter, we first show heuristically that the maximum inward force arising from attractive interactions in a two-state gas of atomic fermions is limited by quantum mechanics to a value less than the outward Fermi pressure. We then demonstrate this idea experimentally by directly cooling a two-component Fermi gas of 6Li atoms to high degeneracy in an optical trap at magnetic fields near a Feshbach resonance, where strong interactions are observed OharaScience . The highest densities obtained in the experiments far exceed those predicted previously Houbiers  for the onset of mechanical instability. Imaging the cloud both in the trap and after abrupt release, we find no evidence of instability. To quantitatively describe the stability, we present a measurement of the relevant universal many-body parameter $`\beta `$—the mean-field contribution to the chemical potential in units of the local Fermi energy. In Ref. OharaScience , we provided a first estimate of $`\beta `$ assuming a zero temperature gas. Here, we present a method to properly account for the non-zero temperature of the gas. This revised measurement of $`\beta `$ is in good agreement with a prediction of Heiselberg .

The mechanical stability of the gas can be understood through a heuristic discussion of the forces acting on a gas of fermions in a 50-50 mixture of two spin components. The equation of state for a normal, zero-temperature Fermi gas is Stringari

$$ϵ_F(𝐱)+U_{\text{MF}}(𝐱)+U_{\text{trap}}(𝐱)=\mu ,$$
(1)
where $`\mu `$ is the chemical potential, $`ϵ_F(𝐱)`$ is the local Fermi energy, $`U_{\text{MF}}(𝐱)`$ is the mean field contribution to the chemical potential, and $`U_{\text{trap}}(𝐱)`$ is the trap potential. In the local density approximation, $`ϵ_F(𝐱)=\mathrm{}^2k_F^2(𝐱)/(2M)`$, where $`k_F`$ is the local Fermi wavevector, which is related to the density according to $`n(𝐱)=k_F^3(𝐱)/(6\pi ^2)`$ Butts .

For each spin component, there is an effective outward force fermiforce  arising from the Fermi pressure $`P_{\text{Fermi}}=2n(𝐱)ϵ_F(𝐱)/5`$, even in the absence of interatomic interactions. The outward force per unit volume is $`P_{\text{Fermi}}`$. The corresponding force per particle is $`(P_{\text{ Fermi}})/n`$. The local force per particle is then easily shown to be

$$𝐅_{\text{Fermi}}=\mathbf{}ϵ_F(𝐱),$$
(2)
Since the density decreases with distance from the center, $`𝐅_{\text{Fermi}}`$ is outward.

Two body scattering interactions make a contribution $`U_{\text{MF}}`$ to the chemical potential, which is given in the mean field approximation by,

$$U_{\text{MF}}(𝐱)=\frac{4\pi \mathrm{}^2a_{\text{eff}}}{M}n(𝐱),$$
(3)
where $`a_{\text{eff}}`$ is an effective scattering length, which generally depends on the local thermal average of a momentum-dependent scattering amplitude. The potential is attractive when $`a_{\text{eff}}<0`$, and repulsive when $`a_{\text{eff}}>0`$. The local force per particle arising from the mean field contribution is just

$$𝐅_{\text{MF}}=\mathbf{}U_{\text{MF}}(𝐱).$$
(4)
If we assume that $`a_{\text{eff}}`$ is energy-independent and equal to $`a_S`$ for $`a_S<0`$, it is easy to show that the gas becomes unstable for suitably large values of $`a_S`$. In this case, the inward force from the mean-field potential, $`𝐅_{\text{MF}}\mathbf{}n(𝐱)`$ while the outward force from the Fermi pressure, $`𝐅_{\text{Fermi}}\mathbf{}n^{2/3}(𝐱)`$. Then the inward force exceeds the outward when $`|𝐅_{\text{MF}}|>|𝐅_{\text{Fermi}}|`$, i.e., when $`k_F|a_{\text{eff}}|>\pi /2`$. The corresponding density for mechanical instability satisfies $`n>n_0`$ where $`n_0=\pi /(48|a_S|^3)`$. This result matches the previous prediction of Houbiers  which was derived via a more rigorous calculation.

The assumption that $`a_{\text{eff}}`$ is energy-independent, however, is only valid if $`k_F|a_S|1`$. Outside that regime, this assumption violates the quantum mechanical requirement of unitarity. At intermediate densities Heiselberg , where $`a_Sk_F^1n^{1/3}R`$, with $`R`$ the range of the collision potential, two-body scattering is dominant, but the mean field interaction is proportional to a momentum-dependent two-body T-matrix element, which in turn is proportional to the scattering amplitude $`f(k)`$. It is well known that $`f`$ has a magnitude which is limited by unitarity to a maximum of $`1/k`$ Landau . Hence, one expects that in a zero-temperature Fermi gas, one should use an effective scattering length in Eq. 3 with a maximum magnitude on the order of $`a_{\text{eff}}=1/k_F`$. For $`a_{\text{eff}}=1/k_F`$, it is easy to show that the mean field contribution to the chemical potential is OharaScience ,

$$U_{\text{MF}}(𝐱)=\beta ϵ_F(𝐱),$$
(5)
where the universal parameter $`\beta `$ has a maximum value of $`4/(3\pi )=0.42`$ in this heuristic treatment. In this limit, both the inward force, $`𝐅_{\text{MF}}`$ and the outward force, $`𝐅_{\text{Fermi}}`$ are proportional to the same power of the atomic density. The ratio of their magnitudes, however, is given by $`|𝐅_{\text{MF}}|/|𝐅_{\text{Fermi}}|=|\beta |<1`$. Thus, one finds that a two-component Fermi gas is mechanically stable as a result of the quantum mechanical requirement of unitarity. This is in contrast to attractive Bose gases Hulet ; Wieman  and Bose-Fermi mixtures Inguscio  which exhibit dramatic instabilities.

A possible criticism of our heuristic argument is that the two-body relative wavevector $`k`$ lies in the range $`0kk_F`$. Since $`f`$ increases as $`k`$ decreases, $`f(k)`$ should be averaged over $`W(𝐤)`$, the probability distribution for $`𝐤`$, where $`𝑑𝐤W(𝐤)=1`$. For a noninteracting gas at zero temperature, one can show

$$W(𝐤)=\frac{6}{\pi k_F^3}\left[1\frac{3}{2}\frac{k}{k_F}+\frac{1}{2}\left(\frac{k}{k_F}\right)^3\right]\mathrm{\Theta }(k_Fk),$$
(6)
where $`k=|𝐤|`$. Near a Feshbach resonance, where $`|a_S|R`$, we assume that $`f(k)`$ takes the form expected for a zero energy resonance, and that the effective scattering length is determined by the real part of $`f`$, i.e., $`a_{\text{eff}}=a_S/(1+k^2a_S^2)`$, where $`\mathrm{}`$ denotes averaging with Eq. 6. We find that $`|a_{\text{eff}}|`$ has a maximum value of $`1.05/k_F`$ when $`k_F|a_S|=2`$, and decreases slowly for $`k_F|a_S|>1`$, reaching $`0.6/k_F`$ at $`k_F|a_S|=10`$.

Our heuristic estimate, which neglects the effects of the interactions on the free particle wavefunctions, yields a value of $`\beta `$ comparable to that estimated by Heiselberg Heiselberg . In contrast to our heuristic approach, the self consistent many-body approach of Ref. Heiselberg  also predicts that $`\beta `$ is independent of the sign and magnitude of $`k_Fa_S`$ in the intermediate density limit where $`k_F|a_S|1`$. Hence, for strongly interacting fermions, $`\beta `$ is always negative and the effective mean field interaction should always be attractive.

To demonstrate that a strongly attractive, two spin-component Fermi gas of atoms is stable, we employ a 50-50 mixture of the two lowest hyperfine states of 6Li, i.e., the $`|F=1/2,M=\pm 1/2`$ states in the low field basis. This mixture has a predicted broad Feshbach resonance at 860 G Elastic ; Zerocross . A magnetic field of 910 G is applied to produce a very large and negative zero-energy scattering length.

The gas mixture is prepared and rapidly cooled to degeneracy at 910 G by forced evaporation in our ultrastable CO2 laser trap, as described for our previous experiments OharaScience . For our trap, $`\omega _{}=2\pi \times (6625\pm 50\text{Hz})`$, $`\omega _z=2\pi \times (230\pm 20\text{Hz})`$, and $`\overline{\omega }(\omega _{}^2\omega _z)^{1/3}=2\pi \times (2160\pm 65\text{Hz})`$. Following evaporation, the trap is recompressed to full depth over 0.5 s, and allowed to remain for an additional 0.5 s, to ensure thermal equilibrium. The CO2 laser power is then extinguished and the gas is imaged as described in our previous paper OharaScience . After release, the gas expands hydrodynamically, rapidly increasing in the transverse dimension while remaining nearly stationary in the axial direction OharaScience .

We determine the number of atoms by numerically integrating the column density (Table 1). For a typical number, $`N8.0\times 10^4`$ atoms per state, the corresponding Fermi density for the experiments is calculated to be $`n_F4.8\times 10^{13}/\text{cm}^3`$ per state Butts . The density $`n_F`$, can be compared to the density $`n_0=\pi /(48|a_S|^3)`$ predicted for the onset of mechanical instability, assuming a momentum-independent scattering length. Using the best available molecular potentials, which are constrained by measurements of the zero crossing in the s-wave scattering length Zerocross ; Grimm , the zero energy scattering length $`a_S`$ is estimated to be $`10^4a_0`$ at 910 G ($`a_0=0.53\times 10^8`$ cm), within a factor of two Vanessa . Then $`n_0=4.4\times 10^{11}/\text{cm}^3`$, showing that the density $`n_F`$ exceeds $`n_0`$ by one to two orders of magnitude.

Although the molecular potentials are constrained by the measured zero crossing in the scattering length Zerocross ; Grimm , the precise location of the Feshbach resonance still may be incorrect. Since the scattering length has not been directly measured, it is possible that the attractive potential is not as large as expected. However, a simple argument based on our previous experiments shows that $`a_S`$ must be large OharaScience . We observe similar hydrodynamic expansion after release from either full trap depth $`U_0`$ or from a reduced trap depth of $`U_0/100`$. In the latter case, we estimate that $`k_F^1=4000a_0`$. Since the observed hydrodynamic expansion appears to be independent of the trap depth, as expected for a unitarity limited interaction, we conclude that in the shallow trap we must have $`|k_Fa_S|>1`$. Then we must have $`|a_S|>4000a_0`$ and $`n_06.9\times 10^{12}/\text{cm}^3`$.

Thus, we have clearly demonstrated that the gas is mechanically stable in the intermediate density regime, where $`|a_S|`$ is large compared to the interparticle spacing. The remainder of this paper provides a revised estimate of the universal parameter $`\beta `$ which quantitatively determines the stability of the gas. In our previous estimate OharaScience , we assumed a zero temperature gas and calculated $`\beta `$ from an estimate of the release energy. There, we obtained a value of $`\beta `$ substantially smaller in magnitude than predicted Heiselberg . Here, we determine $`\beta `$ from the transverse spatial widths of the expanding gas and properly include both finite temperature and hydrodynamic scaling effects.

The one dimensional transverse spatial profiles after expansion for 0.4-0.8 ms are very well fit by normalized finite-temperature Thomas-Fermi (T-F) distributions which determine the width $`\sigma _x`$ as well as the ratio of the temperature to the Fermi temperature, $`T/T_F`$. We expect that $`T/T_F`$ is approximately constant during the expansion, since we observe hydrodynamic scaling of the transverse radii consistent with an effective potential $`n^{2/3}`$ OharaScience , suggesting an adiabatic process. The T-F shape is not unreasonable despite a potentially large mean field interaction, since the mean field contribution to the chemical potential, Eq. 5 is proportional to the local Fermi energy. In this case, assuming Eq. 1 for a normal degenerate gas at zero temperature Stringari , it is easy to show that the mean field should simply scale the Fermi energy of the trapped cloud without changing the shape from that of a T-F distribution OharaScience . The zero temperature spatial distribution then corresponds to that of a harmonic oscillator potential, with the frequencies scaled so that $`\omega _i^{}=\omega _i/\sqrt{1+\beta }`$, $`i=x,y,z`$. Hydrodynamic expansion then preserves the shape of the trapped atom spatial distribution OharaScience ; Stringari . Assuming that $`T/T_F`$ is small, one expects that a finite temperature spatial distribution for a harmonic potential is a reasonable approximation. A Sommerfeld expansion AshcroftMermin  of the normalized density for $`|x|\sigma _x`$ yields

$$n(x)/N=\frac{16}{5\pi \sigma _x}\left[f_0(x)+5\pi ^2(T/T_F)^2f_2(x)\right],$$
(7)
where $`f_0(x)=g^5(x)`$, $`f_2(x)=g(x)/8g^3(x)/6`$, and $`g(x)=\sqrt{1x^2/\sigma _x^2}`$.

Eq. 7 is used to fit the measured transverse spatial distributions, Fig. 1. We obtain the results given in Table 1 for four trials each at expansion times $`t`$ of 0.6 ms and 0.8 ms. Note that Eq. 7 begins to break down near $`|x|\sigma _x`$ for $`T/T_F>0.15`$. However, an exact treatment using polylogarithm functions yields similar results even for our highest temperatures where $`T/T_F0.18`$.

The value of the parameter $`\beta `$ can be determined from the transverse radii of the trapped cloud, $`\sigma _x(0)=\sqrt{2ϵ_F^{}/M\omega _{}^{}_{}{}^{}2}`$, where $`ϵ_F^{}=\mathrm{}\overline{\omega }^{}(6N)^{1/3}`$ is the Fermi energy including the mean field contribution. Then, $`\sigma _x(0)=(1+\beta )^{1/4}\sigma _{xF}`$ Stringari3 , where $`\sigma _{xF}=\sqrt{2ϵ_F/(M\omega _{}^2)}`$, and $`ϵ_F=\mathrm{}\overline{\omega }(6N)^{1/3}`$ is the Fermi energy in the absence of interactions. As we have pointed out previously, the observed anisotropic expansion can arise from unitarity-limited collisional hydrodynamics or superfluid hydrodynamics OharaScience . In either case, the effective potential is $`n^{2/3}`$, and we can assume $`\sigma _x(t)=\sigma _x(0)b_x(t)`$ OharaScience ; Stringari . Hence,

$$\beta =\left(\frac{\sigma _x(t)}{b_x(t)\sigma _{xF}}\right)^41,$$
(8)
For our trap parameters, $`b_x(0.6\text{ms})=29.74`$ and $`b_x(0.8\text{ms})=39.88`$. The scale factor $`b_x(t)`$ properly includes the spatial anisotropy of the expansion and the correct hydrodynamic scaling. We obtain from Table 1 an average value $`\beta =0.26\pm 0.07`$.

This result is consistent with that obtained from the measured transverse release energy. The release energy can be initially calculated from the zero temperature T-F fits OharaScience  or calculated numerically from second moment of the density. Since $`T/T_F0.15`$ from Table 1, the measured release energy is larger than the zero temperature release energy by a factor $`\eta 1+(2\pi ^2/3)(T/T_F)^2=1.15`$. Reducing the release energy by a factor $`\eta `$ and using the method of Ref. OharaScience , yields an average value of $`\beta `$ consistent with that given above.

The measured average value of $`\beta =0.26\pm 0.07`$ can be compared with predictions of the energy per particle for the conditions of our experiment, where $`k_Fa_S=7.4`$. For our equation of state, the sum of the local kinetic and mean field energy per particle is $`(3/5)(1+\beta )ϵ_F(𝐱)`$ Stringari2 . In Ref. Heiselberg , Heiselberg provides two estimates for the total energy per particle. His Eq. 11, based on the Galitski equations, yields $`\beta _{\text{calc}}=0.54`$ Error . Using the Wigner-Seitz cell approximation, his Eq. 14 provides an alternative estimate, $`\beta _{\text{calc}}=0.33`$. An independent calculation by Steele Steele  using effective field theory yields $`\beta =0.46`$. The theoretical predictions are for zero temperature, and the second is in reasonable agreement with our measurements including the temperature correction to the spatial distributions. However, the predictions do not include the additional temperature dependence of $`\beta `$ which arises from a thermal average. Since this important universal many-body parameter can now be experimentally measured, further refinement of both the experimental measurements and the theory is worthwhile and may permit the observation of the superfluid corrections to $`\beta `$. In addition, mixtures of the three lowest hyperfine states of 6Li may permit observation of predicted mechanical instabilities in a three-component Fermi gas Heiselberg .

This research is supported by the Physics Division of the Army Research Office and the Fundamental Physics in Microgravity Research program of the National Aeronautics and Space Administration, the Chemical Sciences, Geosciences and Biosciences Division of the Office of Basic Energy Sciences, Office of Science, U. S. Department of Energy and by the National Science Foundation.