Contents

• Excitation conditions in the Orion molecular cloud obtained from observations of ortho- and para-lines of H$_{\sf 2}$
• 1 Introduction
• 2 Observations and data reduction
• 3 Ortho/para ratios and their relationship to v= 1-0 S(0) and S(1) line brightness
• 4 Observational constraints on models
• 5 PDR as a possible source of excitation
• 6 Shocks as a source of H$_{\sf 2}$ excitation
• 7 Comparison of observations with shock models
• 8 Concluding remarks
• References
• 9 Online Material
• Appendix A: Model results for classes A1, A2, B and C and objects 1, 2 and 3

A&A 469, 561-574 (2007)
DOI: 10.1051/0004-6361:20065786

Excitation conditions in the Orion molecular cloud obtained from observations of ortho- and para-lines of H$_{\sf 2}$^,

L. E. Kristensen¹ - T. L. Ravkilde² - D. Field^2, - J. L. Lemaire^1, - G. Pineau des Forêts^3,1

1 - Observatoire de Paris & Université de Cergy-Pontoise, LERMA, UMR 8112 du CNRS, 92195 Meudon Cedex, France
2 - Department of Physics and Astronomy, University of Aarhus, 8000 Aarhus C, Denmark
3 - Institut d'Astrophysique Spatiale, UMR 8617 du CNRS, Université de Paris Sud, 91405 Orsay, France

Received 8 June 2006 / Accepted 6 March 2007

Abstract
Aims. We seek to study excitation mechanisms in the inner region of the Orion Molecular Cloud by comparing observations of ortho- and para-lines of H₂ with theoretical models of slow shocks and photodissociation regions.
Methods. K-band observations of H₂ obtained with the Canada-France-Hawaii 3.6 m telescope using the PUEO adaptive optics system are reported. Data were centered on the Becklin-Neugebauer object northwest of the Trapezium stars. Narrow-band filters were used to isolate emission from the v=1-0 S(1) ortho- and v=1-0 S(0) para-lines at a spatial resolution of 0 $.\!\!^{\prime\prime}$45 ($\sim$200 AU). We are able to combine their intensity to obtain the column densities of rovibrationally excited ortho and para H₂ levels of the molecular gas at high spatial resolution.
Results. The resulting line ratios show variations between 2 and the statistical equilibrium value of 6. We find 4 different classes of emission, characterised by the ratio of the v=1-0 S(1) and S(0) line brightness and the absolute line brightness. Shock models are used to estimate the physical properties of pre-shock density and shock velocity for these 4 classes. We find that the pre-shock density is in the range of 10⁵-10⁷ cm^-3 and shock velocities lie between 10 and 40 km s^-1. Studies of individual objects, using additional constraints of shock velocity and width, allow quite precise physical conditions to be specified in three prominent bow shocks, one with a shock speed of $18\pm2$ km s^-1 and pre-shock density $1\pm0.5\times10^6$ cm^-3 (3$\sigma$) and two with shock speeds of $\sim$$36\pm2$ km s^-1 and pre-shock densities of $7.5\pm2.5\times10^4$ cm^-3.

Key words: ISM: individual objects: OMC1 - ISM: kinematics and dynamics - ISM: molecules - shock waves - ISM: lines and bands

1 Introduction

The aim of this work is to study the nature of shocks in the inner region of the Orion Molecular Cloud (OMC1) using H₂ NIR emission within a radius of $\sim$30 $^{\prime\prime}$ of the BN-IRc2 complex. OMC1, D=460 pc (Bally et al. 2000), is the closest active massive star forming region. Among the wealth of information available reference is only included here to those data and theoretical models which have a direct bearing on H₂ emission in shocks. For general reviews of the Orion region, the reader is referred to O'Dell (2001) and Ferland (2001).

OMC1 and its immediate surroundings form a nursery of OB stars, both exposed, as in the Trapezium, and buried deep within dusty gas (Beuther et al. 2004; Menten & Reid 1995; Gezari et al. 1998). Spatially associated with these massive stars is a cluster of young low mass stars formed in the last 10⁶ years (Hillenbrand 1997). It has long been suggested that outflows from young OB stars may trigger star formation (Elmegreen & Lada 1977). This may serve to explain the general observation that the great majority of low mass stars form in clusters, many with OB associations. How this trigger may operate remains unclear. The present work seeks in a modest way to shed more light on this issue by characterizing shocks associated with OB stars in more detail than has previously been achieved, using both observation and theory.

The gas in OMC1 displays very many flows, varying in velocity between a few  km s^-1 (Gustafsson et al. 2003; Gustafsson 2006; Nissen et al. 2007; Chrysostomou et al. 1997) to several hundred  km s^-1, the latter represented by the well-known fingers and bullets (eg. Allen & Burton 1993). The origin of the flows is threefold arising from outflows from OB stars, which permeate the entire region of OMC1 (O'Dell & Doi 2003; McCaughrean & Mac Low 1997; Allen & Burton 1993; Stone et al. 1995; Doi et al. 2002), local flows arising from outflows from low mass protostellar objects buried within OMC1 (Gustafsson et al. 2003; Nissen et al. 2007), and supersonic turbulence (Gustafsson et al. 2006b,a). Flows in the inner zone generate slow shocks, graphically illustrated in Gustafsson et al. (2003) & Nissen et al. (2007), with a detailed morphology given by K-band IR images of the v=1-0 S(1) line at 2.121 $\mu$m (Vannier et al. 2001; McCaughrean & Mac Low 1997; Gustafsson et al. 2003; Schild et al. 1997; Lacombe et al. 2004; Schultz et al. 1999; Kristensen et al. 2003; Stolovy et al. 1998; Chen et al. 1998). Lacombe et al. (2004), using the NACO adaptive optics system on the VLT, displays a resolution of 30 AU and provides the most detailed morphology available.

In this paper we attempt to find shock speeds and densities within the centre of OMC1 using observational data for ortho- and para- rovibrational lines of H₂. The analysis is performed by confronting data with the latest shock models (Flower & Pineau des Forêts 2003), with recourse to further constraints such as values of radial bulk velocities and shock widths (Gustafsson et al. 2003; Nissen et al. 2007; Lacombe et al. 2004) which were not available in earlier work (e.g. Smith et al. 1997). We consider both the general characteristics of shocks and three specific zones associated with an outflow from a deeply buried massive star located close to BN, source I (Beuther et al. 2004,2006; Greenhill et al. 2004a; Gezari et al. 1998; Menten & Reid 1995; Greenhill et al. 2004b) for which the IR signature of the outflow has recently been discovered (Nissen et al. 2007).

Quite generally, the relative proportions of ortho- and para- populations of H₂ (${\it o/p}$ ratios) are valuable data for establishing the nature of interstellar shocks (Wilgenbus et al. 2000, and references therein). High temperature H atom exchange-reactions provide net conversion of para- into ortho- H₂, tending to create the high temperature limit of ${\it o/p} =3$. Thus in warm post-shock gas the ${\it o/p}$ ratio becomes a measure of the period of time for which the gas is sufficiently hot and dense for para- to ortho- conversion to proceed. This period of time in turn depends on the shock velocity and the pre-shock density. The resulting ${\it o/p}$ ratio, frozen into the cold post-shock gas for upwards of 10⁵ years (Flower et al. 2006), also depends on the initial ${\it o/p}$ ratio. This initial ratio becomes an additional important parameter of the medium (Flower et al. 2006).

Incomplete conversion of para- to ortho-H₂, that is o/p<3has been observed in a variety of different objects ranging from starburst galaxies (e.g. NGC253, Harrison et al. 1998), galactic molecular clouds (Rodríguez-Fernández et al. 2000), star forming regions and Herbig-Haro objects (e.g. Davis et al. 1999; Lefloch et al. 2003; Cabrit et al. 1999; Neufeld et al. 1998). For OMC1, only two observational papers report data related to the ${\it o/p}$ ratio. Rosenthal et al. (2000) takes in the whole NW region (Peak 1) in an aperture of $15\hbox{$^{\prime\prime}$ }\times30\hbox{$^{\prime\prime}$ }$, reporting an averaged value of the ${\it o/p}$ ratio of 3. Smith et al. (1997), in a detailed study of OMC1, used three rotational lines of H₂, employing v=1-0, S(0) and S(2) to obtain a rotational temperature. The S(2) line may be severely absorbed in the atmosphere as Smith et al. note. v=1-0 S(1) brightness was then used to extract an effective ${\it o/p}$ ratio (see Sect. 3). Smith et al. (1997) reported $3.0\pm0.4$ with no strong evidence of variation of the ${\it o/p}$ ratio down to scales of 1 $^{\prime\prime}$. In this connection, as we show in Sect. 4 there is diffuse background of H₂ line emission which permeates at least the northern part of OMC1. The spatially integrated brightness of this weak emission may dominate the total emission in H₂ lines and has a major contribution from the photo dissociation region (PDR) generated by $\theta^1$Ori C, and J-shock contribution. Both of these generate an ${\it o/p}$ratio of 3 under the prevalent conditions.

The present data differ from earlier work in the important respect that images are taken at a spatial resolution that allows identification of features down to 0.45 $^{\prime\prime}$, that is a scale of 200 AU. There are quite extensive zones at high S/N ratio where data are consistent with an ${\it o/p}$ ratio of 3, as earlier work suggests. However our higher resolution and concentration on bright shocked regions of high S/N also shows significant zones where para- to ortho- conversion is incomplete. The current work builds on Kristensen et al. (2003) and extends the observations described there both to include the para-H₂rovibrational line, v=1-0 S(0), as well as a greatly enlarged field of view.

Observations and data reduction are described in detail in Sect. 2. In Sect. 3 we discuss how to estimate quantities related to the ${\it o/p}$ ratio from our observations of v=1-0 S(0) and v=1-0 S(1). The observational constraints to be met by PDR and shock-models are described in Sect. 4. The nature of PDR and shock models is described in Sect. 5 and the grid of shock models investigated and characteristic output is shown in Sect. 6. Fits of models to observations are presented in Sect. 7, showing some degree of success in constraining the physical conditions of specific shocked zones in OMC1.

  
2 Observations and data reduction

Observations were performed using the 3.6 m Canada-France-Hawaii Telescope on Mauna Kea, Hawaii from the 6th to the 8th of December 2000. The observations are centered around the BN-object situated $\sim$70 $^{\prime\prime}$ northwest of the Trapezium stars, that is 0.16 pc at the distance of Orion. The field of view is shown in Fig. 1 where we show continuum subtracted emission from the v=1-0 S(1) H₂ line. The observed field covers the region designated as Peaks 1 and 2 by Beckwith et al. (1978).

Figure 1: The field in OMC1 covered in the present observations. The data consist of three regions marked East, West and North. Data show emission in the H2 v=1-0 S(1) line. The values on the greyscale bar are in units of 10-5 W m-2 sr-1. Axes are labelled in arcseconds from the reference position of TCC0016: 05$^{\rm h}$35$^{\rm m}$14 $.\!\!^{\rm s}$91, -05$^\circ$22$^\prime$39 $.\!\!^{\prime\prime}$31 (J2000).

Observations were performed with the PUEO adaptive optics (AO) system at the Canada-France-Hawaii Telescope (CFHT). PUEO was equipped with the KIR detector ( $1024\times1024$ pixels). At the time of our observations the seeing was rather poor ($\ge$1 $.\!\!^{\prime\prime}$5). The lens set used corresponds to 35 mas/pixel, with a field of view of $36\hbox{$^{\prime\prime}$ }\times 36\hbox{$^{\prime\prime}$ }$. Data were recorded in the ortho-H₂v=1-0 S(1) and para-H₂ v=1-0 S(0) lines.

The region of OMC1 observed consists of three overlapping fields each $36\hbox{$^{\prime\prime}$ }\times 36\hbox{$^{\prime\prime}$ }$, with the entire region centered approximately 5 $^{\prime\prime}$ W and 15 $^{\prime\prime}$ N of TCC0016 (05$^{\rm h}$35$^{\rm m}$14 $.\!\!^{\rm s}$91, -05$^\circ$22$^\prime$39 $.\!\!^{\prime\prime}$31; J2000), which we use as a positional reference throughout, lying itself about 40 $^{\prime\prime}$ N and 15 $^{\prime\prime}$ W of the Trapezium cluster. We designate these three fields as West, East and North (Fig. 1).

Two wavefront reference stars, necessary for functioning of the PUEO-AO system, were used: TCC0016 (m_V=14.0) for fields East and West, and Parenago 1838 (m_V=15.2) for field North. The resolution of the observations is $\sim$0 $.\!\!^{\prime\prime}$45 corresponding to 200 AU for all fields reported, that is, a Strehl ratio between 0.08 and 0.16. Thus correction was rather poor, reflecting the unsatisfactory seeing. It is found that no improvement in spatial resolution could be achieved by deconvolution whilst maintaining acceptable signal-to-noise.

Data reduction to obtain H₂ images is performed so as to take account of any temporal variability of the sky background, spatial variations in the sensitivity of the detector (flat-fielding), differences in the sky brightness at different wavelengths and differing efficiencies of the detection system for the different filters used (see below). Dark counts are subtracted and bad pixels and noise due to cosmic rays removed.

In both sets of data, isolation of spectral lines and observation of the continuum at 2.183 $\mu$m were achieved using narrow-band filters. Line filters were centred on 2.121 $\mu$m for the v=1-0 S(1) line and 2.223 $\mu$m for the S(0) line. The filter width is 0.02 $\mu$m ($\lambda$/100) in each case.

Comparison of line brightness using filters raises the following five issues: (i) emission from other lines than that desired within the filter bandwidth (ii) atmospheric absorption (iii) differential reddening (iv) relative filter transmission and (v) image registration. We consider each of these below.

(i) Within the S(1) filter only additional high v, high Jlines may be present. These lines are negligibly weak in shocks but may be found in PDRs. However PDRs are intrinsically one to two orders of magnitude lower in brightness than the C-type shocks encountered here. Moreover these very high v lines are weak in PDRs (Black & van Dishoeck 1987).

Within the S(0) filter there is also contamination from high v, high J lines. There may also be weak contamination from the v=2-1 S(1) line, which lies 0.024 $\mu$m to longer wavelength than the v=1-0 S(0) line. However transmission through the S(0) filter of the v=2-1 S(1) line is only 3%. We conclude that contamination by other lines is not a problem for either the S(1) or S(0) filters.

The continuum background is subtracted from each filter. The continuum is weak, that is, typically considerably less <10% of either major line brightness. The continuum filter has a line centre 2.183 $\mu$m. Within this filter there are no H₂ lines. There may be some weak contribution from Br$\gamma$ which we neglect here.

(ii) Turning to atmospheric absorption, it is essential that brightness estimates are as free as possible from differential effects between the two lines. The velocity of the gas must be considered in this context. Data obtained (Dec. 2000) on the Canada-France-Hawaii Telescope, using a combination of the PUEO adaptive optics system and Fabry-Perot interferometry ("GriF''; Clénet et al. 2002; Gustafsson et al. 2003; Nissen et al. 2007), as well as extensive data in Chrysostomou et al. (1997), reveal that the region of OMC1 observed contains H₂emission which shows velocity shifts, relative to Earth, of between +60 to -10 km s^-1. Using the atmospheric absorption line atlas of Livingston & Wallace (1991), we find that there is negligible absorption for the v=1-0 S(1) line in all cases, save over a very narrow range of velocities around +30 km s^-1 for which an absorption of 7% is found. For the v=1-0 S(0) line, the situation is similar with a weak absorption feature again of 7% at around +43 km s^-1. GriF data show that the regions studied span the range of velocities which includes these values. Thus differential absorption may introduce systematic errors into estimates of S(1)/S(0) ratio, but of only a few per cent. The effect cannot be accurately determined and we choose to ignore it in the present work.

(iii) Differential reddening is a further consideration. The v=1-0 S(0) line will be less reddened than the v=1-0 S(1) line. The relative magnitude difference between the two is $\sim$ $\left(\lambda_1/\lambda_2\right)^{-1.7}$ (Mathis 1990). If we adopt the extinction law derived by Rosenthal et al. (2000) and an extinction at 2.12 $\mu$m of 1mag, the v=1-0 S(0) line may be overestimated by $\sim$7% compared to the v=1-0 S(1) line. We present results here for data uncorrected for this imprecisely known and spatially variable differential absorption. If included, ratios used below could be reduced by $\sim$7%.

(iv) The relative transmission of the entire system through the v=1-0 S(1) and S(0) filters at the operating temperature of 77 K has been determined directly from an extensive set of observations. This was done by observing between 6 and 10 stars in each region and for each filter. The stars chosen are bright, non-saturated and as free from emission from the nebula as possible. Over the wavelength range of interest, the continuum-emission from these stars is effectively constant. Comparison of the flux from the stars in different filters yield the relative transmission of the entire system. The S(0) filter is found to have a transmission which was found to vary between $38\pm2$% and $43\pm2$% of that of the S(1) filter depending on the region considered. Corresponding figures appropriate to each region were used in calculating line ratios. The continuum-filter shows $70~~\pm~2$% transmission of that of the S(1) filter for each case.

(v) Spatial registration of the images is critical. For each of the three fields, between 6 and 10 stars - those used for filter transmission estimates - were used for image registration in each field. This allowed for registration to better than $\pm$1 pixel between the images in the two lines throughout the entire region. Thus images of the ratio of line brightness could be made without significant loss of spatial resolution.

We also performed an absolute calibration of the emission brightness for the S(1) v=1-0 line using TCC0044 ( m_K' = 10.50), TCC0031 ( m_K' = 9.86) and TCC0011 ( m_K' =10.94) (McCaughrean & Stauffer 1994) as standard stars. In the brightest part of the image, at a position 15 $.\!\!^{\prime\prime}$8 east and 2 $.\!\!^{\prime\prime}$1 south of TCC0016, we find a brightness of $2.4\pm0.4\times10^{-5}$ W m^-2 sr^-1. This figure may be compared with the value of $3.0\pm0.15\times10^{-5}$ W m^-2 sr^-1 quoted in Vannier et al. (2001). The discrepancy arises from the higher spatial resolution of 0 $.\!\!^{\prime\prime}$15 in Vannier et al. (2001) compared with the present data. For the purpose of comparison with shock models, it is more appropriate to use the higher figure in Vannier et al. (2001) which we therefore adopt from hereon. Note that absolute values of brightness may however be underestimated associated with the reddening referred to above (Rosenthal et al. 2000).

  
3 Ortho/para ratios and their relationship to v= 1-0 S(0) and S(1) line brightness

In principle it is necessary to obtain the full set of ortho- and para-lines in order to estimate a meaningful ${\it o/p}$ ratio. A Boltzmann plot of log(column density per sublevel) vs. energy of the level would then show departures from the high temperature equilibrium value of o/p = 3, if such departures exist. However we show below that because of the proximity in energy of the J=2 and J=3 levels in v=1, it is possible to obtain approximate values of an ${\it o/p}$ ratio which are meaningful, using only S(0) and S(1) v=1-0 emission line data. To differentiate these values from the ${\it o/p}$ ratio global to all lines, we name the term derived purely from v=1-0 S(0) and S(1), $\phi _{10}$.

We use the definition of the ${\it o/p}$ ratio found in standard textbooks and used in Wilgenbus et al. (2000); Ramsay et al. (1993); Hora & Latter (1996); Hoban et al. (1991); Chrysostomou et al. (1993); Neufeld et al. (1998). The ${\it o/p}$ ratio at local thermodynamic equilibrium (LTE) at a rotational temperature of $T_{\rm rot}$ is given by

$$o/p({\rm LTE}, T_{\rm rot}) = \frac{\sum_{J\ {\rm odd}} 3(2J+... ... (2J+1) \exp\left( \frac{-E_J}{k_{\rm B} T_{\rm rot}} \right)}$$ (1)

where J is the rotational quantum number and E_J the energy of the rotational state for a given vibrational level v.

For each ortho-level, J, the non-equilibrium ${\it o/p}$ ratio is given by (Wilgenbus et al. 2000):

$$\frac{o/p(J)}{o/p({\rm LTE},T_{\rm rot})} = \frac{N_J}{N_J({\rm LTE},T_{\rm rot})}$$ (2)

where N_J is the value of the column density of the non-equilibrium ortho-line and $N_J({\rm LTE},T_{\rm rot})$ is the expected column density of the ortho-line, had it been observed at LTE at a rotational temperature of $T_{\rm rot}$. Note that we deal throughout with column density. This implicitly ignores any spatial variation in the line of sight, though such variations must of course be present.

Referring to the upper state of the transition v=1-0 S(0), that is v=1, J=2, as i=0 and the upper state of v=1-0 S(1)  that is v=1, J=3, as i=1, one obtains the approximate ${\it o/p}$ ratio, $\phi _{10}$:

$$\phi_{10} = \frac{N_1}{N_0}\frac{g_0}{g_1} \exp\left( \frac{E_1 - E_0}{k_{\rm B} T_{\rm rot}} \right) o/p(LTE, T_{\rm rot})$$ (3)

where g_i [=(2I+1)(2J+1)] is the total multiplicity and $(E_1-E_0)/k_{\rm B}=6947$ K-6474 K=473 K.

The column densities, N_i, in Eqs. (2) and (3) can be obtained from the observed brightness, I_i, using

$$N_i = \frac{4 \pi \lambda_i}{hc} \frac{I_i}{A_i}$$ (4)

where $\lambda$ is the wavelength and A is the Einstein A-value for the corresponding line, $3.47\times10^{-7}$ s^-1 for v=1-0 S(1) and $2.53\times10^{-7}$ s^-1 for v=1-0 S(0) (Wolniewicz et al. 1998). Inserting Eq. (4) in Eq. (3) gives (Harrison et al. 1998):

 

$$\phi_{10} \frac{I_1}{I_0} \frac{\lambda_1}{\lambda_0} \frac{A_0}{A_1} \frac{g_0}{g_1}$$ =

$$\times \exp\left(\frac{E_1-E_0}{k_B T_{\rm rot}} \right)\ o/p({\rm LTE},T_{\rm rot})$$

$$0.4970\ \frac{I_1}{I_0}\ \exp\left(\frac{473~{\rm K}}{T_{\rm rot}} \right)$$ = (5)

where we assumed $T_{\rm rot}>300$ K for which o/p(LTE, $T_{\rm rot})=3$.

Figure 2: A map of the approximate o/p ratio, $\phi _{10}$, calculated from v=1-0 S(0) and S(1) emission, for the field identified in Fig. 1 estimated using Eq. (5). The area in grey represents regions in which emission is below specified signal levels (see text). The colour bar is for $\phi _{10}$. Coordinates are relative to TCC0016 as in Fig. 1. Original images have been smoothed using a boxcar average over $7\times 7$ pixels. The three large squares delineated by grey borders are named East, West and North, as illustrated in Fig. 1. The black rectangles, A1+A2, B and C, delineate regions which have been chosen for special study.

Using ISO-SWS observations, Rosenthal et al. (2000) find that the rotational temperature measured with an aperture of $\sim$15 $^{\prime\prime}$ by 30 $^{\prime\prime}$ is of the order of 3000 K. Le Bourlot et al. (2002) reanalysed the data and found the rotational temperature to be 3300 K. In Kristensen et al. (2003) it was found that the excitation temperature over a small field in region East varied between 2000 and 5000 K. This excitation temperature was calculated from the v=1-0 S(1) and v=2-1 S(1) H₂ lines using high spatial resolution data from the ESO 3.6 m telescope (see also Vannier et al. 2001). Unpublished data recently obtained from the VLT in the same two lines show that the excitation temperature in Peak 1 (NW of BN) is in the interval 2000 K to 4000 K. In the following we have chosen a constant value of $T_{\rm rot}=3500$ K based on all of the above observations.

The systematic errors generated by the energy term in Eq. (5) are small. For example, given that the rotational temperature, $T_{\rm rot}$ is in the interval from 2000 K to 5000 K as suggested by the observations just mentioned, the error introduced by taking a constant value of the rotational temperature in the energy term exp(473 K/ $T_{\rm rot}$) is no greater than $\sim$10%.

We emphasise that $\phi _{10}$ refers only to the ratio in the v=1, J=2 and 3 excited states in that part of the medium in which they are populated, and does not represent the ${\it o/p}$ ratio of all the molecular H₂ present in the medium.

The resulting $\phi _{10}$ map can be seen in Fig. 2. To avoid unacceptable levels of noise in forming this image, all emission in the v=1-0 S(1) and v=1-0 S(0) lines weaker than $8\times10^{-7}$ W m^-2 sr^-1 was excluded. This represents $\sim$2.5% of the maximum in the v=1-0 S(1) line and 9% of the maximum in the weaker v=1-0 S(0) line. Prior to obtaining the ratio the v=1-0 S(1) and v=1-0 S(0), images were smoothed using $7\times 7$ boxcar averaging. This degraded the spatial resolution by $\sim$15%. The map shows surprisingly clear structure in $\phi _{10}$, ranging from $\phi _{10}$ of 1 to 3. In particular, individual clumps of material in region West in Fig. 2 each show structure where $\phi _{10}$ is low (1-1.5) at the centre of emission rising to 3 at the edges.

Structures in $\phi _{10}$ found in our field cannot be caused by differential exinction effects alone as this would require a K-band extinction of more than 6.5 mag. Under these circumstances H₂emission would locally have to be several hundred times greater than that observed, that is, of the order of several times 10^-3 up to 10^-2 W m^-2 sr^-1. This is greatly in excess of any value that shock or PDR models can account for. Moreover data in van Dishoeck et al. (1998) show that H₂ emission lies in part in front of the 9.7 $\mu$m silicon absorption feature. These data indicate that at least some of the H₂ emission is generated in a region relatively unobscured by the main absorbing material.

A comparison may be made between our values of $\phi _{10}$ and the ISO-SWS data reported in Rosenthal et al. (2000). We have performed a weighted average over the aperture of the ISO-observations, using the S(1) brightness as weight. We find that $\phi_{10}=2.5\pm0.3$similar to the value of 3 quoted in Rosenthal et al. (2000).

  
4 Observational constraints on models

In the inner zone of OMC1 studied here, which omits the Orion fingers or bullets to the NW (e.g. Allen & Burton 1993), we may divide the H₂ emission into the following groups, based upon the general characteristics of the emission. The first group consists of blue-shifted emission representing a massive outflow originating between Peaks 1 and 2, in the north-eastern part of region West in Fig. 1. This group of objects is discussed in detail in Nissen et al. (2007). Data obtained with VLT using the NACO adaptive optics system resolve the widths of isolated shocks in this region in a very graphic manner (Lacombe et al. 2004, Sect. 7.2). The second group belongs to Peak 1 and Peak 2 (North and East in Fig. 1). These are especially bright, with overlapping interconnected features and a complex velocity structure (Gustafsson et al. 2003; Nissen et al. 2007). The third group is represented by the faint background emission observed in region North. This does not show small scale spatial structure at our level of sensitivity and spatial resolution. The brightness of this pervasive emission in the v=1-0 S(1) line is $4.0\pm1.3\times10^{-6}$ W m^-2 sr^-1. Brightness in v=1-0 S(0), corresponding to this level of emission in v=1-0 S(1), lies below the noise level. However there remains a good deal of diffuse S(0) emission detectable at around $2.0\pm1.3\times10^{-6}$ W m^-2 sr^-1, noting the brightness ratio of S(1) to S(0) lies between a factor of 2 and 6. This type of emission as characterized by the S(1) line shows no detectable velocity structure (Nissen et al. 2007).

In the following we seek to find a generalized set of shock and PDR models which are consistent with our observations. These observations include both the line brightnesses in v = 1-0 S(0) and S(1) as well as the ratio and also velocities as measured with GriF (Gustafsson et al. 2003; Nissen et al. 2007). In the next section we will also include the width of the bowshock-structures observed in the VLT-NACO data. For a part of the East field, we also have brightness data for the v=2-1 S(1) line (Kristensen et al. 2003).

To put our data in a generalized form, we plot the absolute brightness of the v=1-0 S(1) vs. the line ratio, R₁₀ defined as $I_{v=1-0~{\rm S}(1)}$/ $I_{v=1-0~{\rm S}(0)}$ for the regions A1+A2, B and C whose locations are given in Fig. 2. Results are shown separately for the regions A1+A2, B and C in Fig. 3. Very similar results are obtained with the v=1-0 S(0) data. The choice of location and size of regions A1+A2, B and C as distinctive regions was made partly on the basis of the map of $\phi _{10}$ in Fig. 2 and partly following the results in Nissen et al. (2007). Note that the zone north-east of BN which lies at -15 $^{\prime\prime}$, +17 $^{\prime\prime}$ relative to TCC0016, south-east of A1+A2, has been excluded because of possible artefacts associated with strong continuum emission in this region.

Figure 3: a) Region North: brightness of v=1-0 S(1) vs. the line ratio R10. b) Region West : similarly for the blueshifted clumps in this zone. c) Region East. All data have been rebinned to $500\times 500$ pixels from the original field size of $2000\times 2000$ pixels. The effective pixel size becomes 0 $.\!\!^{\prime\prime}$14 or 3 times better than the resolution. The grey oblongs identify those parts of the data whose characteristics are given in Table 1.

In the A1+A2 region, Fig. 3a, there is a clear tendency for pixels with higher brightness to have higher R₁₀. Two condensations of points located at R₁₀=3.2, $I_{\rm S(1)}=7\times10^{-6}$ W m^-2 sr^-1 and R₁₀=4.2, $I_{\rm S(1)}=1.2\times10^{-5}$ W m^-2 sr^-1 are clearly seen in Fig. 3a. These two classes of points were identified according to the following criterion. The two condensations were first separated by locating the minimum in point density between the two condensations. The contour of this minimum point density was then used around each condensation to form a locus defining each class. These loci are shown in Fig. 3a schematically as oblongs, defining the range of properties which specify points of class A1 and A2.

It is evident that certain regions are associated with either the A1 or A2 classes. That is, the low R₁₀ are found in a restricted zone in the southern and eastern half of the A1+A2 region. Thus the A1 region is specifically that part of the emission. This also turns out to be the more weakly emitting zone. The A2 class of points is restricted to the two high ratio zones in Fig. 2. Figure 3a also shows that a minimum value of brightness is associated with each ratio. This is not an artefact due to a noise level cut-off, which lies at $\sim$ $8\times10^{-7}$ W m^-2 sr^-1, but arises because of the diffuse background. This has a brightness of $\sim$ $4.0\times10^{-6}$ W m^-2 sr^-1 in the S(1) line (see above).

Figure 3b shows data for the blue-shifted clumps in region West. Similar plots restricted to individual blue-shifted clumps show the same structure of higher brightness towards lower values of R₁₀ (Sect. 7.2). Thus here, in contrast to class A1 or A2, positions of data points within the scatter plot are not associated with any particular spatial sub-zone of the chosen region. The loci of points which we call class B is defined by the oblong in Fig. 3b. The criterion here is that we have chosen the subset of data with >65% of the maximum brightness. The reason for this restriction is as follows. In class B, which represents the blue outflow region, much of the data arises from highly localized shocks, some of which take a bow form, judging from the morphology in Lacombe et al. (2004). Data in our chosen subset refer to that brighter emission which lies near the tip or centre of the bow shock. We therefore do not consider the fainter wings of the bow shocks.

Figure 3c, for Region East, shows a different structure, with a central condensation around R₁₀=4.0 and $I_{\rm S(1)}=1.2\times10^{-5}$ W m^-2 sr^-1. We have v=2-1 S(1) data for part of region C (Kristensen et al. 2003) and these yield a diagram of very similar appearance to that shown in Fig. 3c. We define R₁₂ as the line ratio of v=1-0 S(1) to v=2-1 S(1). The oblong, defining points of class C, was obtained as follows. Contours of density were obtained and all data above the half-maximum were included, as schematically outlined by the oblong in Fig. 3c.

 

 

Table 1: Characteristics of the 4 classes of points described in the text and displayed in Fig. 3. Brightness is given in units of 10^-5 W m^-2 sr^-1. The v=2-1 S(1) brightness and R₁₂ are from Kristensen et al. (2003) and radial velocities are from Gustafsson (2006); Nissen et al. (2007). The figures shown as $\pm$ represent the range of values.

Observations	Class A1	Class A2	Class B	Class C
Brightness v=1-0 S(1)	$0.67\pm0.11$	$1.15\pm0.10$	$0.90\pm0.08$	$1.05\pm0.18$
Brightness v=1-0 S(0)	$0.24\pm0.04$	$0.28\pm0.04$	$0.31\pm0.09$	$0.31\pm0.08$
Brightness v=2-1 S(1)				$0.15\pm0.05$
R10	$3.2\pm0.6$	$4.2\pm0.5$	$3.2\pm1.1$	$4.2\pm0.8$
R12				$7\pm2$
$\phi _{10}$	$1.8\pm0.3$	$2.4\pm0.3$	$1.8\pm0.6$	$2.4\pm0.5$
Associated radial velocity / km s-1	11	11	18	8

Our task now is to identify shock models which satisfy the characteristics of data of classes A1, A2, B and C as specified in Fig. 3a-c. These characteristics are listed in Table 1. Each class is defined by a range of characteristic values of absolute brightness and line ratio(s). Also included in table 1 are values of average radial shock velocities taken from GriF data reported in Gustafsson (2006); Nissen et al. (2007).

  
5 PDR as a possible source of excitation

H₂ emission in OMC1 arises from both heating through shocks (e.g. Vannier et al. 2001; Kristensen et al. 2003) and from photon excitation in PDRs (e.g. Sternberg & Dalgarno 1989; Black & van Dishoeck 1987; Störzer & Hollenbach 1999; Black & Dalgarno 1976). We turn first to PDRs.

We now show that the diffuse background of H₂ emission which permeates most of region North (but not region East or West), and to which we have drawn attention in Fig. 3a, may be approximately modelled using results reported from existing PDR codes. In our region $\theta^1$Ori C, an O6 star in the Trapezium located at a projected distance of 0.16 pc from BN, generates a radiation field of 2- $3\times10^5$ times the standard interstellar field (G₀=2- $3\times10^5$). Combined with a high density, for example exceeding $n_{\rm H}>10^5$ cm^-3, collisional events result in a kinetic temperature in a PDR with values greater than 800 K (Sternberg & Dalgarno 1989; Kristensen et al. 2003; Le Petit et al. 2006; Kaufman et al. 1999; Störzer & Hollenbach 1999). The importance of this figure here is that interactions between H and H₂ begin to overcome the activation energy barrier for H atom exchange at these temperatures, scrambling the ortho- and para- populations and creating ${\it o/p} =3$, as mentioned in the introduction (Sternberg & Neufeld 1999).

We use results from the PDR models of both Störzer & Hollenbach (1999) and the "Meudon PDR code'' (Le Petit et al. 2006). We focus upon the weaker background emission without measurable velocity structure because (i) PDRs are unable to reproduce the high brightness of many localized regions (ii) the large bulk motions in the gas, associated with very bright regions, are not characteristic of PDRs. We therefore seek to reproduce a brightness in v=1-0 S(1) of $\sim$ $4.0\times10^{-6}$ W m^-2 sr^-1, with an upper limit of $\sim$ $8\times10^{-7}$ W m^-2 sr^-1 in S(0), the noise level. This implies that R₁₀ must be greater than 5 resulting in a lower limit of $\phi _{10}$ of 2.8 close to the high temperature equilibrium value of the ${\it o/p}$ ratio of 3.

Using the model of Störzer & Hollenbach (1999) with $G_0=2.4\times10^5$, $n_{\rm H}=4.0\times10^6$ cm^-3, including 2.6 km s^-1 of advection, a value of $4.2\times10^{-6}$ W m^-2 sr^-1 arises in the S(1) line. This is in fact the maximum that any models in Störzer & Hollenbach (1999) report and reproduces the observed value of the S(1) background emission seen in region North. The corresponding brightness for the S(0) line is not reported in Störzer & Hollenbach (1999).

Turning to use of the "Meudon PDR code'', we first note this does not include advection. This has the result that the high brightness in v=1-0 S(1) is more difficult to match, at any rate for a simple face-on model. The most extreme conditions explored use $n_{\rm H}=5\times10^6$ cm^-3 and $G_0=5\times10^5$. These yield S(1) brightness of $3.0\times10^{-6}$ W m^-2 sr^-1. The ratio R₁₀ is calculated to be 3.8 and thus S(0) is predicted to be close to the noise level but a little too bright. In this connection, R₁₀ is insensitive to the value of G₀ in the range of high number densities and high G₀ used here.

We conclude that a significant part of the diffuse background in region North is due to the direct action of a PDR generated by $\theta^1$Ori C. We also conclude that the density here is higher than 10⁶ cm^-3 implying that the temperature is >1500 K. Hence changes in the ${\it o/p}$ ratio occur through reactive collisions. The region is of course also subjected to the well-known major outflow from the general area of BN/IRc2. Thus diffuse shocked gas also makes a contribution to the emission (see Sect. 7.1).

  
6 Shocks as a source of H$_{\sf 2}$ excitation

In a shock, H₂ is excited through mechanical heating, at the microscopic level through high temperature H₂-H₂, H-H₂ and He-H₂ collisions (Le Bourlot et al. 1999). As the shock develops, the temperature becomes sufficient that excited vibrational states become significantly populated. Emission is observed in the IR, for example, from J=2 or J=3 states in v=1 to form respectively the S(0) and S(1) lines. We first consider the type of shocks relevant here, that is, whether they are J- or C- type.

6.1 C-type vs. J-type shocks

As we now show it appears very likely that the shocks which give rise to localised bright emission in the central region of OMC1 are magnetic C-type shocks, rather than non-magnetic J-type. First, it has been demonstrated that the region can support substantial magnetic fields (Crutcher et al. 1999; Norris 1984) and the gas is at least weakly ionized. Second, there are numerous features, especially in the central zone (region West in Fig. 1) between Peaks 1 and 2 (regions North and East), which are clearly individual shocks, as imaged at 70 mas resolution (30 AU) using the NACO-VLT adaptive optics system (Lacombe et al. 2004). We return to individual objects in the NACO-VLT field in Sect. 7.2. The component of magnetic flux density transverse to the direction of shock propagation in a C-type shock softens the shock and makes very extensive the region in which high temperatures and accompanying excitation of H₂ are encountered. We find below that it is possible to model observed shock widths of 40-80 AU in dense regions only with C-type shocks.

The occurrence of J-type, non-magnetic shocks (Hollenbach & McKee 1989; Lim et al. 2002) has been discussed in detail in Kristensen et al. (2003). It was shown there, for data in region East, that J-shocks contribute in very restricted areas at the edges of clumps. These zones are not resolved here.

  
6.2 Shock model description and initial conditions

C-type shocks have been investigated by Wilgenbus et al. (2000); Kaufman & Neufeld (1996b); Smith & Brand (1990); Draine et al. (1983); Pineau des Forêts et al. (1988); Kaufman & Neufeld (1996a); Timmermann (1998), whose results were used in Vannier et al. (2001). The most recent version of the models are those of Le Bourlot et al. (2002), Flower et al. (2003) and Flower & Pineau des Forêts (2003). The latter extends the work of Wilgenbus et al. (2000), showing that the critical velocity for C-shocks is lower than previously believed depending on both the initial PAH-abundance and initial density of the shock. PAHs play an important role as they increase the magnetosonic speed (Flower & Pineau des Forêts 2003), given that PAHs may efficiently attach electrons (Field et al. 1999,2004). For a C-shock to propagate, it is necessary that the shock-velocity is below the magnetosonic speed. If this is not satisfied, the shock becomes a J-type shock.

In the following we use the model described in Flower & Pineau des Forêts (2003) and references therein. The most important limitation of the current C-type shock model is that the dynamics is treated in one dimension. Therefore as the postshock gas is compressed it will remain compressed. In reality the compressed gas would diffuse into the non-shocked surrounding medium due to the large pressure gradient. Model results therefore give only an indication of the scale-size which should be associated with a certain set of conditions.

In other respects the treatment of the shock is very detailed, with regard to an inclusive range of chemistry, including etching of grains (May et al. 2000), atomic and molecular excitation and associated cooling. Chemical events, involving 136 species and described by 1040 chemical processes, determine critical parameters such as the degree of ionization in the medium. Initial abundances of the elements included in the model are given in Flower & Pineau des Forêts (2003). Initial species abundances in the shock models are derived from chemical steady state models, where the PAH abundance is set to $n_{\rm PAH}$/ $n_{\rm H}=10^{-6}$. The rate of cosmic ray ionization was taken to be $5\times10^{-17}$ s^-1 per H atom. The model abides by the relationship that transverse magnetic flux density is given by $b \times(n_{\rm H}^{1/2}$) $\mu$G, where $n_{\rm H}=$ number of hydrogen nuclei ($\sim$2 $n_{\rm H_2}$) in units of cm^-3, and b, the magnetic scaling factor, is typically unity (Flower et al. 2003, and references therein). This is in contrast to the models reported for example in Smith et al. (1991) and Smith (1991) which invoke very high magnetic fields.

The critical H atom exchange reaction with H₂ with its associated temperature dependence is treated using results from ab initio quantum mechanical calculations (Le Bourlot et al. 2002, and references therein). The full chemistry and all excitation and cooling processes are integrated in parallel with the magnetohydrodynamic equations. This is necessary because of the strong coupling of the chemistry and physics in the evolution of the shock. A total of 100 rovibrational level populations of H₂ are calculated in parallel with the dynamical and chemical variables, allowing for all radiative transitions and collisional processes which modify level populations. Since the IR radiative transitions are weak, transitions in H₂ are optically thin for line transfer for any relevant column densities.

In all models of J-type shocks we set the magnetic flux density equal to zero. A subset of the J-type shocks included below will dissociate H₂, in which case H₂ reformation on grains becomes an important process. We chose to let the internal energy of the newly formed H₂molecules be proportional to a Boltzmann distribution at a temperature of 17 250 K corresponding to one third of the binding energy of H₂.

  
6.3 Grid of C-type shock models

The parameters at play are (i) shock velocities (ii) the initial ${\it o/p}$ratio (iii) the pre-shock density and (iv) the magnetic flux density. We limit ourselves to slow shocks of between 10 to 50 km s^-1 (Clénet et al. 2002; Gustafsson et al. 2003; Gustafsson 2006; Nissen et al. 2007) lying within the field shown in Fig. 1. These shocks are essentially non-dissociative (save for the relatively unimportant J-type shocks), non-ionizing and do not possess a radiative precursor. The fast shocks of several hundred  km s^-1 (Doi et al. 2002; Lee & Burton 2000, and references therein) associated with the bullets or fingers (Allen & Burton 1993) lie outside our field of view.

The parameters which we wish to establish are pre-shock densities, shock velocities and an indication of the magnetic field. The constraints which we seek to match are outlined in Table 1. In addition we use shock widths, measured in Lacombe et al. (2004) as constraints for individual objects. We first discuss the ${\it o/p}$ ratio of the pre-shocked gas.

  
6.3.1 Values of the initial o/p ratio

The o/p ratio in a dark cloud at T=10 K, left undisturbed for a period in excess of 1 Myr, attains a steady state value close to that of the kinetic temperature of the gas. The conversion process arises through exchange reactions involving H⁺, H₃⁺ and other protonated species (Flower et al. 2006). Under these conditions the o/p ratio would be $\ge$ $2\times10^{-3}$. This is shown in Fig. 4, with an initial ${\it o/p}$ ratio of 3.0 and physical conditions T=10 K, $n_{\rm H}=10^5$ cm^-3, cosmic ray ionization rate $5\times10^{-17}$ s^-1 per H atom and initial degree of ionization $\sim$10^-8, the same conditions as used in the steady state models to describe the pre-shock gas (see Sect. 6.2). Grains are not included in the chemical steady state model (Flower et al. 2005). The timescale for converting the o/p ratio from that of hot gas, that is, a value of 3, to a cold gas value is only weakly dependent on density.

Figure 4: Temporal evolution of the ${\it o/p}$ ratio in cold gas at T=10 K and $n_{\rm H}=10^5$ cm-3 in a chemical steady-state gas-phase model. The initial o/p ratio is set to 3.0.

Due to the high level of activity in Orion, any parcel of gas may have been shocked more than once in its lifetime. If the interval between shocks is greater than a few times $\sim$10⁶ years, then an initial ${\it o/p}$ of <0.1 is relevant for the study of the effects of any subsequent shock. If an earlier shock caused an ${\it o/p}$ ratio of less than 3 to be frozen into the gas - which is a likely event as we show below - then the time interval appropriate to an initial ${\it o/p}<0.1$ may be <10⁶ years. One may imagine a series of shocks impinging on a parcel of gas in the molecular cloud. The first shock will raise the ${\it o/p}$ ratio to 1 (say) from which it will relax before the next shock in the series reaches the parcel. This subsequent shock will further raise the o/p ratio and so on until a value of 3 is achieved. Evidently we need to explore all possible initial o/p ratios. In this connection there are no processes in the shock itself which lend themselves to a reduction in the o/p ratio.

6.3.2 Parameter ranges

For the study of C-type shocks, the range of parameters investigated was as follows:

• preshock density: 10⁴, $5\times10^4$, 10⁵, $5\times10^5$, 10⁶, $5\times10^6$ and 10⁷ cm^-3;
• shock velocity: 10-50 km s^-1, step-size: 1 km s^-1;
• initial ${\it o/p}$ ratio: 0.01, 1.0, 2.0 and 3.0;
• transverse magnetic flux density = $b \times n_{\rm H}^{1/2}$ $\mu$G, with b=1.0 and 5.0.

The detailed results of this grid of models will be published in a forthcoming paper (Kristensen et al. 2007).

A grid of J-type models was also calculated with the same conditions, except that b=0.0. Trial calculations showed that the results of the model were essentially independent of the initial ${\it o/p}$ ratio, the high temperatures leading to complete conversion to ${\it o/p} =3$. Therefore the single value of 3 was used.

  
6.4 Model results

In the following we discuss seven parameters which the grid of models yields. These are the spatially local o/p ratio, the integrated ${\it o/p}$ratio, the rotational temperature within v=1, shock width, v=1-0 S(1) and S(0) line brightness and their ratio, R₁₀. Line brightness is estimated using a face-on geometry and integrating along the length of the shock. The shock width is defined as the length of region for which the gas temperature remains >1000 K.

6.4.1 Local o/p ratios: temperature dependence of para to ortho conversion

The purpose of this section is to demonstrate the interplay between shock velocity, density and associated shock temperatures in para- to ortho- conversion. For this purpose it is instructive to use the local ${\it o/p}$ ratio, see Eq. (1), bearing in mind of course that the observable quantity is an integrated value.

For any pre-shock density we take a range of shock velocities. We then record the set of maximum temperatures which are achieved in each of these shocks. For any maximum temperature we record the maximum spatially local ${\it o/p}$ ratio achieved. This is a measure of the effectiveness of conversion of para to ortho H₂. We repeat this for a set of pre-shock densities. Note that the o/p ratio integrated along the line-of-sight will always be less than any local maximum value (Wilgenbus et al. 2000). Note also that the local maximum once achieved remains unchanged for a period of time of the order of typically 10³-10⁴ years (Fig. 4). The cooling time for the shocks considered here is less than 100 years.

In Fig. 5 we show these data for pre-shock densities of 10⁴ to 10⁷ cm^-3. Clearly, and independently of density, when the maximum kinetic temperature exceeds $\sim$3200 K the conversion of para- to ortho-H₂ is locally complete and the local maximum ${\it o/p}$ ratio will be 3 - although as pointed out above an observed integrated value may be <3. In addition Fig. 5 shows that para- to ortho- conversion sets in weakly at $\sim$800 K and rises sharply above 1300 K.

Figure 5: Maximum local o/p ratios in C-type shocks, with magnetic scaling factor, b equal to 1, as a function of temperature for different preshock densities. Para- to ortho-H2 interconversion begins at $\sim$800 K and is complete at $\sim$3200 K.

In order to demonstrate the behaviour for a single shock we show in Fig. 6 the spatial variation of shock temperature, the local ${\it o/p}$ ratio, the integrated value of $\phi$₁₀, the integrated value of R₁₀ and the integrated ${\it o/p}$ ratio, for a shock velocity of 20 km s^-1, pre-shock density of 10⁶ cm^-3 and initial ${\it o/p}$ of 0.01. It may be readily seen that the local ${\it o/p}$ ratio continues to rise until the temperature falls below about 800 .

To clarify this further, in Fig. 7 we show the maximum kinetic temperature as function of shock velocity for a set of preshock densities. This map of maximum kinetic temperatures, taken in conjunction with Fig. 5, serves to identify combinations of preshock density and velocity for which local para- to ortho conversion will be initiated and for which it will be complete. For example for a 15 km s^-1 shock in pre-shock gas of 10⁶ cm^-3, with initial o/p=0.01, local para- to ortho- conversion is incomplete throughout, achieving a maximum local value of 0.57, since the maximum temperature is only 2100 K. However for a 25 km s^-1 shock with the same pre-shock density, local para- to ortho- conversion may be complete, since the maximum temperature is 3900 K. The integrated value in the latter case however achieves a value of only 1.1.

Figure 6: Comparison of the ${\it o/p}$ ratio, local and integrated as specified, $\phi _{10}$ and R10 as a function of position within a C-type shock with preshock density of $1\times 10^6$ cm-3, shock velocity 20 km s-1, initial ${\it o/p}$ ratio 0.01 and b=1.

Figure 7: Maximum C-type shock temperature as a function of shock velocity for different preshock densities. The shaded area indicates the region in which net interconversion between ortho- and para-H2 takes place ($\sim$800-3200 K). The temperature is independent of initial o/p ratio and b=1.

As one introduces larger model velocities for a given density, the shock velocity will at some stage become greater than the magnetosonic velocity. The shock is then no longer a C-type shock but a J-type shock. This dictates the high velocity limits considered in Fig. 7 and subsequently in Figs. 8-12.

Figure 8: C-type shock widths as a function of shock velocity for different preshock densities (b=1).

Figure 9: Local brightness in v=1-0 S(1), S(0) and v=2-1 S(1) H2 emission lines in the evolution of a C-type shock with initial shock velocity 20 km s-1, preshock density $1\times 10^6$ cm-3, initial o/p ratio 0.01 and the magnetic scaling factor, b, equal to 1. The kinetic temperature of the shock is shown as a solid line. Zero distance is set to be the point at which the kinetic temperature starts to rise.

Figure 10: Brightness of the v=1-0 S(1) line as a function of C-type shock velocity for different preshock densities. The brightness is almost independent of initial ${\it o/p}$ ratio. It is shown here for ${\it o/p}_{\rm ini}=3$.0 and b=1.

Figure 11: Brightness of the v=1-0 S(0) line as a function of C-type shock velocity for different preshock densities and b=1, using the same linestyle as in Fig. 10. Each of the 4 diagrams shows a different value of ${\it o/p}_{\rm ini}$ as labelled.

Figure 12: $\phi _{10}$ as a function of C-type shock velocity for different preshock densities and b=1. Each plot shows a different value of integrated o/p ranging from 0.01 to 3.0. Linestyles specifying pre-shock density are as in the previous Fig. 10.

6.4.2 Calculated integrated o/p ratio

For the integrated o/p ratio, in contrast to the local, there is no clear temperature for the conversion to be complete, that is, we have encountered no temperature beyond which the integrated o/p ratio is 3. If such a temperature exists it must be greater than 10⁴ K. On the other hand there exists a minimum temperature for the onset of para- to ortho-H₂ and this temperature is naturally the same as that for the local o/p ratio (see Fig. 5).

It turns out that models yield roughly the same value for the integrated o/p ratio and $\phi _{10}$ (always integrated) if the preshock gas is assigned an o/p ratio of $\le$2. Under these circumstances the integrated o/p ratio and $\phi _{10}$, which may lie between 0.01 and 3, agree to within 20%.

When computing the integrated o/p ratios obtained from shock models and comparing with observations, we choose 15 K as the kinetic temperature for the post-shock gas as the limit at which the integrated o/p ratio is evaluated.

6.4.3 Shock width

As defined above, the width of a shock is taken as the length of the region over which the gas temperature remains $\ge$1000 K. It is well-known that the width of a C-type shock depends on the strength of the ion-neutral coupling (Draine 1980) and thus on the degree of ionization. For the models presented here the degree of ionization takes on a value of between $\sim$10^-7-10^-8. In Fig. 8 we plot the shock width for different values of the preshock density as a function of shock velocity. It is seen that for high preshock densities and velocities (greater than 10⁵ cm^-3 and 15 km s^-1) the width is almost constant as a function of density. Note that when the shock never achieves a temperature >1000 K, no width is recorded (Fig. 7).

A related point is that shocks measured in different lines will take on different apparent dimensions. Populations of upper states, for example J=3, v=1 for the S(1) line, or J=2, v=1 for the S(0) line are created at different points in the shock. For example if the ${\it o/p}$ ratio is initially 0.01, substantial conversion must take place before S(1) emission can form, even though the temperature may be sufficient to excite v=1. Thus different emission lines probe regions of different extent and different location as seen in Fig. 9. A further example is that emission from rotational states in v=2involves a more restricted zone than that for v=1.

This is illustrated for v=1-0 S(0), S(1) and v=2-1 S(1) in Fig. 9, where the local brightness, rather than integrated, is shown. The v=1-0 S(1) line turns on strongly at a later stage in the shock than the S(0) line. Also both these lines gain more brightness over a larger range than does the v=2-1 S(1) line.

6.5 Predictions of the model

The purpose of this section is to discuss the range of values of brightness of H₂ emission that models yield for the grid of parameters set out earlier. As we have seen, brightness depends on shock velocity and pre-shock density. But in particular the relative brightness of individual ortho- and para-lines is strongly dependent on the initial ${\it o/p}$ ratio as earlier discussion has suggested. Figures 10-12 present a brief summary of the behaviour which is found.

We turn first to Fig. 10. It is evident that the value of S(1) v=1-0 emission brightness alone does not in general specify the characteristics of a shock. Thus a brightness of 10^-6 W m^-2 sr^-1 can arise through pre-shock gas of density 10⁴-10⁶ cm^-3 and shock velocity of 13-42 km s^-1. Additional observational constraints of (say) shock width, shock velocity or the brightness of other lines, are required.

A further point is that for high density systems, an increase in the shock velocity does not necessarily create an increase in emission in the v=1-0 S(1) line. This arises because, for higher shock velocities, lower J and v may become more nearly thermalized with population spread among a greater range of levels.

Figure 11 shows data as for Fig. 10 but for v=1-0 S(0) and for four different initial o/p ratios. There are clear variations in the S(0) brightness which are most pronounced for higher densities. Qualitatively similar behaviour, but significantly different in detail, is found for the v=1-0 S(1) line both with respect to shock velocity and the influence of the initial o/p ratio. This behaviour gives rise to a strong variation in the line ratio, as expressed by $\phi _{10}$, with respect both to shock velocity and initial ${\it o/p}$ ratio. This is shown in Fig. 12. We see that in the case of initial ${\it o/p}=3.0$ the value of $\phi _{10}$ lies close to 2.5 independent of preshock density for densities $\ge$10⁵ cm^-3. The lowest density data are limited to values of shock velocity greater than 20 km s^-1 since below this value there is insignificant excitation of the v=1 populations.

  
7 Comparison of observations with shock models

  
7.1 Physical conditions associated with different classes of data

Our aim is primarily to establish shock velocities and preshock density for all four classes of data defined in Sect. 4. This may be successfully achieved through comparison with a very large number of models taken from the grid described in Sect. 6.3. From the outset we note that there are generally insufficient constraints to exclude anything but a large range of initial o/p values for any of the four classes. The same is true of the magnetic field.

We use a $\chi^2$-method to quantify the best fit models of our observations, calculating $\chi^2 = \sum_{\rm obs} (\frac{X_{\rm obs} - X_{\rm model}}{\sigma_{\rm obs}})^2$ where $X_{\rm obs}$ and $X_{\rm model}$ refer to the observed and modelled quantities, respectively. $\sigma_{\rm obs}$ refers to the uncertainty in the parameter associated with any class, that is, effectively the range of values appropriate to that class. These ranges of values are given in Table 1 for the line brightness. In the case of the velocity, Gustafsson (2006); Nissen et al. (2007) reports only radial velocities. These are effectively minimum velocities and are shown as such in Table 1. The value of $\sigma$ associated with these velocities was the standard deviation of the sample used.

A typical contour plot of confidence intervals, in this case for class A1 data and initial ${\it o/p}=0.01$, with b=1, defining the transverse magnetic flux, can be seen in Fig. 13. Contours of 2, 3, 4 and 5$\sigma$ are shown corresponding to each level of certainty. Similar contour plots were obtained for each value of the initial o/p ratio and of the value of b, for each class. Each contour plot typically covers 200-300 individual shock models. Common to all these contour plots is that they cover a combination of high pre-shock density with low shock velocity to low pre-shock density with high shock velocity. The criterion of fit for each class is taken to be the 3$\sigma$ limit (99.7% confidence). For each value of the initial ${\it o/p}$ ratio, the derived range of values of shock velocity and pre-shock density are shown in the appendix in Table A.1, for both b=1 and b=5. We also show the corresponding range of the post-shock density, the shock width, the integrated o/p ratio and the maximum kinetic temperature, where all values are generated by the shock model.

Figure 13: The confidence intervals of class A1. The contours are given at intervals of $\sigma$ from 2$\sigma$ to 5$\sigma$. The models used to make this plot all have an ${\it o/p}_{\rm ini}$ of 0.01 and b=1.0. See Sect. 7.1 for further description.

There turn out to be rather few general conclusions that may be drawn at this stage from the results in Table A.1 despite the detailed analysis. The underlying reason for this is that we are attempting in the case of data class C, for example, to model all the emission in Region East, which comprises most of Peak 2, in terms of a single set of shock conditions. Nevertheless various general statements may be made which give a useful overview of the characteristics of shocks in the inner part of OMC1. These may be summarised as follows:

• Class C objects, in region East, cannot be modelled with b=5, that is with high magnetic fields for any initial o/p ratio. High magnetic fields are also excluded for classes A2 and B for values of the initial o/p ratio of 0.01. Isolated regions of parameter space may in principle exist where agreement is possible but are not accessed by our model grid.

• Initial ${\it o/p}=0.01$ tends to require higher velocities. This may suggest that higher velocities should be rejected since the required delay between successive shocks to reset the initial ${\it o/p}$ to 0.01 is 10⁷ years, whereas the OMC1 complex is no older than $\sim$10⁶ years (O'Dell 2001; Hillenbrand 1997).

• Higher initial o/p ratio requires higher pre-shock densities leading to higher post-shock densities. This favours high densities in OMC1 clumps, given that low initial o/p seems unlikely for reasons of cloud lifetime.

• The predicted width, in particular the lower limit decreases as higher initial o/p ratios are used. Where widths can be measured, typically values are of the order of 50 to 100 AU (Lacombe et al. 2004). For b=1, this tends to exclude an initial ${\it o/p} =3$ in all classes.

• For higher initial o/p values, higher magnetic fields may also be used to fit the observations. Higher velocities are naturally required because of the cushioning effects of higher fields.

A general conclusion from the above items is that the initial ${\it o/p}$ratio probably lies between 1 and 2. This implies an upper limit of the order of 10⁶ years between successive shocks, consistent with the lifetime of OMC1.

The kinetic gas temperature in OMC1 as measured from for example NH₃, CO or CH₃CCH is $\sim$45-75K (Sweitzer 1978; Churchwell & Hollis 1983; Liszt et al. 1974). At equilibrium the ${\it o/p}$ ratio would be in the range $\sim$0.25-0.9, lower than the initial o/p ratio which we find above. Again this indicates that the gas has probably been shocked previously by jets from protostellar objects in the region or that the PDR generated by massive stars in the region (e.g. $\theta^1$Ori C or BN) have raised the ${\it o/p}$ ratio of the gas.

  
7.2 Individual objects in region West

In region West a group of objects located between 7 $^{\prime\prime}$ to 35 $^{\prime\prime}$ west and -5 $^{\prime\prime}$ to 16 $^{\prime\prime}$ north of our reference, TCC0016, show similar properties regarding the absolute brightness, $\phi _{10}$ and velocity structure (Nissen et al. 2007). For example, the maximum absolute brightness of these objects is $\sim$ $1.0\times10^{-5}$ W m^-2 sr^-1, $\phi _{10}$ is $\sim$1.0-1.5 at the centre of the objects rising to 3 at the edges (see Fig. 14). These objects are of special interest since they are part of the IR counterpart of an outflow identified originally in the radio, originating from a highly obscured massive star (or stars) buried in the depths of OMC1 (source I or n; Shuping et al. 2004; Gustafsson 2006; Nissen et al. 2007; Greenhill et al. 2004c; Menten & Reid 1995).

Figure 14: Map showing $\phi _{10}$ in objects 1, 2 and 3 identified in Fig. 15. The area in grey represents regions in which emission is below specified signal levels (see text). The colour bar is for $\phi _{10}$ and coordinates are in arcsec and relative to TCC0016 (see Fig. 1).

Figure 15: ESO-VLT NACO images of three objects where the bowshocks have been resolved. The greyscale bar is in units of 10-5 W m-2 sr-1 (Lacombe et al. 2004). Coordinates are as in Fig. 1.

We have chosen three objects to model, selected on the basis of their bow shapes. These objects are shown in Fig. 15. Their characteristics are given in Table 2 where widths are obtained from ESO VLT-NACO observations of the region (Lacombe et al. 2004). Note that we now have the additional constraints of shock velocity (but see below) and shock width. In this connection an observed (radial) shock velocity is a few  km s^-1 lower in velocity than the lower limit of the shock speed, since energy is taken into heating in the shock impact and velocity is lost from the impacting material.

In Fig. 16 we show brightness versus R₁₀ for Object 1 (see Fig. 15 for labelling of objects). The oblong identifies the subset of points that we use for comparison with models. Note also the similarity in form with the data in Fig. 3b, which defines this class of objects.

Again we use a $\chi^2$ method to quantify which models fit observations of objects 1, 2 and 3 at the 3$\sigma$ level, using the same grid as earlier. We treat the observed velocity data in the following manner. If the shock velocity in any model is less than the observed radial velocity, then the velocity is included as a constraint in the $\chi^2$ fit. If the velocity is greater than the observed radial velocity, then we do not include this as a constraint. This is in recognition of the fact that the radial velocity is a lower limit to the true velocity. We find below that a fit at 3$\sigma$ is given with a shock velocity essentially equal to the observed radial velocity. The method of analysis adopted ensures that this is not an artefact.

Figure 16: A plot brightness in the v=1-0 S(1) line for object 1 similar to Fig. 3 for object 1, but without spatial rebinning. The oblong encloses those data used for comparison with models.

Figure 17: Confidence intervals for object 2. Model parameters are initial ${\it o/p}=0.01$ and b=1.0.

A contour plot of confidence intervals for object 2, initial o/p in the pre-shock gas =0.01, b=1 can be seen in Fig. 17. Full results are summarized in Table A.2 in the appendix. Because of the extra constraints and our limitation to a single object, we obtain a much narrower range of physical conditions. In fact we can show that object 1 is distinct from objects 2 and 3, reflected in the much lower observed radial velocity.

 

 

Table 2: Characteristics of the 3 objects described in the text and displayed in Fig. 15. The brightness is given in units of 10^-5 W m^-2 sr^-1. The velocities are from Nissen et al. (2007) and the widths from Lacombe et al. (2004).

	Object 1	Object 2	Object 3
Location	-18 $.\!\!^{\prime\prime}$5; +0 $.\!\!^{\prime\prime}$5	-18 $.\!\!^{\prime\prime}$1; -0 $.\!\!^{\prime\prime}$8	-20 $.\!\!^{\prime\prime}$8; -6 $.\!\!^{\prime\prime}$2
Brightness S(1)	$0.94\pm0.10$	$0.87\pm0.10$	$0.65\pm0.07$
Brightness S(0)	$0.34\pm0.08$	$0.29\pm0.08$	$0.34\pm0.08$
R10	$2.8\pm0.5$	$3.2\pm0.7$	$2.8\pm0.5$
$\phi _{10}$	$1.6\pm0.3$	$1.8\pm0.4$	$1.6\pm0.3$
Width / AU	$80\pm30$	$80\pm30$	$40\pm20$
Velocity /  km s-1	$18\pm1$	$37\pm1$	$36\pm1$

The physical conditions in our three objects may be summarised as follows:

• Object 1 requires that the initial o/p ratio be around 1 or above. Moreover the magnetic field cannot be high, that is, b<5. The velocity of the shock lies around 18 km s^-1, suggesting, as do observations, that motion lies approximately in the line-of-sight. It appears difficult to reproduce the observed width, which may be an order of magnitude too low. Whilst strictly the full range of models for object 1 in Table A.2 are of equal validity, on the basis of the width criterion alone perhaps the most satisfactory model overall is that with initial ${\it o/p} = 2$, shock velocity $18\pm2$ km s^-1, pre-shock density $1\pm0.5\times10^6$ cm^-3. At all events, all models at the 3$\sigma$ level show the same pre-shock density, which implies a transverse magnetic flux of 1 mG.

• Objects 2 and 3 may be classed together. Again the derived shock velocities agree with the radial component value, once more suggesting that these shocks are travelling approximately in the line of sight. For b=1, the pre-shock density lies an order of magnitude lower than in Object 1 and is $7.5\pm2.5\times10^4$ cm^-3 with a corresponding transverse magnetic flux of 0.3 mG. Once more the widths are not well reproduced, though in this case they are too large. Higher magnetic fields cannot strictly be excluded but widths are still greater for higher fields. The initial o/p ratio cannot be determined.

8 Concluding remarks

The results presented here show that observations of ortho- and para- lines of H₂ present a useful way of probing the physical conditions in shocked zones. We have introduced the quantity $\phi$₁₀, based on the 2 rovibrational H₂ lines v=1-0 S(0) and S(1), as defined in Eq. (3). A map of $\phi _{10}$, a quantity which we have shown is approximately equal to the true ${\it o/p}$ ratio given a high rotational temperature, demonstrates strong spatial variation, ranging from 1 to the high temperature equilibrium value of 3. Spatially averaged values however are close to 3, in agreement with earlier work.

We have identified 4 classes of objects in OMC1, classified through similar properties with respect to line brightness and values of $\phi _{10}$. This allowed the identification of a diffuse background emission in region North (but not elsewhere) whose presence may be partly attributed to a general PDR arising from the action of $\theta^1$Ori C. The bulk of the work is devoted to the development of a large grid of shock models with a view to identifying the physical conditions associated both with the 4 classes of object and also with specific chosen shocked regions in the field. At the 3$\sigma$ level it was possible to determine a range of shock-models that fit our observations with pre-shock densities ranging from $\sim$10⁵-10⁷ cm^-3 and shock velocities in the range of 10-40 km s^-1. It was found that no J-type shock models fit our observations at the 3$\sigma$ level if we restrict pre-shock densities to <10⁷ cm^-3 for which models are valid.

For individual bow-shocks it was possible to identify relatively precise shock conditions. Working with objects in the massive blue-shifted outflow emerging from between peaks 1 and 2, three objects were examined. A velocity of $\sim$18 km s^-1 and pre-shock density of 10⁶ cm^-3 apply to one such object and a shock velocity of $\sim$36 km s^-1 and pre-shock density of $7.5\times10^4$ cm^-3 apply to the other two. Derived transverse magnetic flux was 1 mG and 0.3 mG respectively. These magnetic fields are similar to those derived from observational data of Norris (1984) and Crutcher et al. (1999).

Observations have recently been made with ESO VLT-NACO with a Fabry-Perot interferometer, scanning the three H₂ rovibrational lines v=1-0 S(0), S(1) and v=2-1 S(1) at a spatial resolution 4-5 times better than reported here. These data will allow us not only to map $\phi _{10}$ at higher spatial resolution and with higher accuracy, but through techniques developed here will also allow us to characterize physical conditions in the emitting objects with greater precision.

Acknowledgements

L.E.K. and J.L.L. would like to acknowledge the support of the PCMI National Program, funded by the CNRS in cooperation with the CEA and IN2P3. T.L.R. and D.F. would like to acknowledge the support of the Aarhus Centre for Atomic Physics (ACAP), funded by the Danish Basic Research Foundation. We also wish to thank the Directors and Staff of CFHT and ESO for making possible the observations reported in this paper. Part of this work was performed using computer resources at the Université de Cergy-Pontoise.

References

 
• Allen, D. A., & Burton, M. G. 1993, Nature, 363, 54
   
• Bally, J., O'Dell, C. R., & McCaughrean, M. J. 2000, AJ, 119, 2919
   
• Beckwith, S., Persson, S. E., Neugebauer, G., & Becklin, E. E. 1978, ApJ, 223, 464
   
• Beuther, H., Zhang, Q., Greenhill, L. J., et al. 2004, ApJ, 616, L31
   
• Beuther, H., Zhang, Q., Reid, M. J., et al. 2006, ApJ, 636, 323
   
• Black, J. H., & Dalgarno, A. 1976, ApJ, 203, 132
   
• Black, J. H., & van Dishoeck, E. F. 1987, ApJ, 322, 412
   
• Cabrit, S., Bontemps, S., Lagage, P. O., et al. 1999, in The Universe as Seen by ISO, ESA SP-427, 449
   
• Chen, H., Bally, J., O'Dell, C. R., et al. 1998, ApJ, 492, L173
   
• Chrysostomou, A., Brand, P. W. J. L., Burton, M. G., & Moorhouse, A. 1993, MNRAS, 265, 329
   
• Chrysostomou, A., Burton, M. G., Axon, D. J., et al. 1997, MNRAS, 289, 605
   
• Churchwell, E., & Hollis, J. M. 1983, ApJ, 272, 591
   
• Clénet, Y., Le Coarer, E., Joncas, G., et al. 2002, PASP, 114, 563
   
• Crutcher, R. M., Troland, T. H., Lazareff, B., Paubert, G., & Kazès, I. 1999, ApJ, 514, L121
   
• Davis, C. J., Smith, M. D., Eislöffel, J., & Davies, J. K. 1999, MNRAS, 308, 539
   
• Doi, T., O'Dell, C. R., & Hartigan, P. 2002, AJ, 124, 445
   
• Draine, B. T. 1980, ApJ, 241, 1021
   
• Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485
   
• Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725
   
• Ferland, G. J. 2001, PASP, 113, 41
   
• Field, D., Lunt, S. L., Hoffmann, S. V., Ziesel, J. P., & Gulley, R. J. 1999, in The Physics and Chemistry of the Interstellar Medium, Proceedings of the 3rd Cologne-Zermatt Symposium, held in Zermatt, September 22-25, 1998, ed. V. Ossenkopf, J. Stutzki, & G. Winnewisser (GCA-Verlag Herdecke), 367
   
• Field, D., Jones, N. C., & Ziesel, J.-P. 2004, Phys. Rev. A, 69, 052716
   
• Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 343, 390
   
• Flower, D. R., Le Bourlot, J., Pineau des Forêts, G., & Cabrit, S. 2003, MNRAS, 341, 70
   
• Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2005, A&A, 436, 933
   
• Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621
   
• Gezari, D. Y., Backman, D. E., & Werner, M. W. 1998, ApJ, 509, 283
   
• Greenhill, L. J., Gezari, D. Y., Danchi, W. C., et al. 2004a, ApJ, 605, L57
   
• Greenhill, L. J., Reid, M. J., Chandler, C. J., Diamond, P. J., & Elitzur, M. 2004b, in IAU Symposium, 155
   
• Greenhill, L. J., Reid, M. J., Chandler, C. J., Diamond, P. J., & Elitzur, M. 2004c, in ed. M. Burton, R. Jayawardhana, & T. Bourke, IAU Symp., 155
   
• Gustafsson, M. 2006, Ph.D. Thesis, University of Aarhus, Denmark
   
• Gustafsson, M., Kristensen, L. E., Clénet, Y., et al. 2003, A&A, 411, 437
   
• Gustafsson, M., Brandenburg, A., Lemaire, J. L., & Field, D. 2006a, A&A, 454, 815
   
• Gustafsson, M., Field, D., Lemaire, J. L., & Pijpers, F. P. 2006b, A&A, 445, 601
   
• Harrison, A., Puxley, P., Russell, A., & Brand, P. 1998, MNRAS, 297, 624
   
• Hillenbrand, L. A. 1997, AJ, 113, 1733
   
• Hoban, S., Reuter, D. C., Mumma, M. J., & Storrs, A. D. 1991, ApJ, 370, 228
   
• Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306
   
• Hora, J. L., & Latter, W. B. 1996, ApJ, 461, 288
   
• Kaufman, M. J., & Neufeld, D. A. 1996a, ApJ, 456, 611
   
• Kaufman, M. J., & Neufeld, D. A. 1996b, ApJ, 456, 250
   
• Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795
   
• Kristensen, L. E., Gustafsson, M., Field, D., et al. 2003, A&A, 412, 727
   
• Kristensen, L. E., Pineau des Forêts, & Lemaire, J. L., & Field, D. 2007, in preparation
   
• Lacombe, F., Gendron, E., Rouan, D., et al. 2004, A&A, 417, L5
   
• Le Bourlot, J., Pineau des Forêts, G., & Flower, D. R. 1999, MNRAS, 305, 802
   
• Le Bourlot, J., Pineau des Forêts, G., Flower, D. R., & Cabrit, S. 2002, MNRAS, 332, 985
   
• Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506
   
• Lee, J.-K., & Burton, M. G. 2000, MNRAS, 315, 11
   
• Lefloch, B., Cernicharo, J., Cabrit, S., et al. 2003, ApJ, 590, L41
   
• Lim, A. J., Raga, A. C., Rawlings, J. M. C., & Williams, D. A. 2002, MNRAS, 335, 817
   
• Liszt, H. S., Wilson, R. W., Penzias, A. A., et al. 1974, ApJ, 190, 557
   
• Livingston, W., & Wallace, L. 1991, An Atlas of the Solar Spectrum in the Infrared from 1850 to 9000 cm^-1, N. S. O. Technical Report No. 91-001, Tech. rep., National Optical Astronomy Observatories
   
• Mathis, J. S. 1990, ARA&A, 28, 37
   
• May, P. W., Pineau des Forêts, G., Flower, D. R., et al. 2000, MNRAS, 318, 809
   
• McCaughrean, M. J., & Mac Low, M. 1997, AJ, 113, 391
   
• McCaughrean, M. J., & Stauffer, J. R. 1994, AJ, 108, 1382
   
• Menten, K. M., & Reid, M. J. 1995, ApJ, 445, L157
   
• Neufeld, D. A., Melnick, G. J., & Harwit, M. 1998, ApJ, 506, L75
   
• Nissen, H. D., Gustafsson, M., Lemaire, J. L., et al. 2007, A&A, 466, 949
   
• Norris, R. P. 1984, MNRAS, 207, 127
   
• O'Dell, C. R. 2001, ARA&A, 39, 99
   
• O'Dell, C. R., & Doi, T. 2003, AJ, 125, 277
   
• Pineau des Forêts, G., Flower, D. R., & Dalgarno, A. 1988, MNRAS, 235, 621
   
• Ramsay, S. K., Chrysostomou, A., Geballe, T. R., Brand, P. W. J. L., & Mountain, M. 1993, MNRAS, 263, 695
   
• Rodríguez-Fernández, N. J., Martín-Pintado, J., Fuente, A., et al. 2000, in ISO Beyond the Peaks: The 2nd ISO Workshop on Analytical Spectroscopy, ESA SP-456, 127
   
• Rosenthal, D., Bertoldi, F., & Drapatz, S. 2000, A&A, 356, 705
   
• Schild, H., Miller, S., & Tennyson, J. 1997, A&A, 318, 608
   
• Schultz, A. S. B., Colgan, S. W. J., Erickson, E. F., et al. 1999, ApJ, 511, 282
   
• Shuping, R. Y., Morris, M., & Bally, J. 2004, AJ, 128, 363
   
• Smith, M. D. 1991, MNRAS, 252, 378
   
• Smith, M. D., & Brand, P. W. J. L. 1990, MNRAS, 242, 495
   
• Smith, M. D., Brand, P. W. J. L., & Moorhouse, A. 1991, MNRAS, 248, 451
   
• Smith, M. D., Davis, C. J., & Lioure, A. 1997, A&A, 327, 1206
   
• Störzer, H., & Hollenbach, D. 1999, ApJ, 515, 669
   
• Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197
   
• Sternberg, A., & Neufeld, D. A. 1999, ApJ, 516, 371
   
• Stolovy, S. R., Burton, M. G., Erickson, E. F., et al. 1998, ApJ, 492, L151
   
• Stone, J. M., Xu, J., & Mundy, L. G. 1995, Nature, 377, 315
   
• Sweitzer, J. S. 1978, ApJ, 225, 116
   
• Timmermann, R. 1998, ApJ, 498, 246
   
• van Dishoeck, E. F., Wright, C. M., Cernicharo, J., et al. 1998, ApJ, 502, L173
   
• Vannier, L., Lemaire, J. L., Field, D., et al. 2001, A&A, 366, 651
   
• Wilgenbus, D., Cabrit, S., Pineau des Forêts, G., & Flower, D. R. 2000, A&A, 356, 1010
   
• Wolniewicz, L., Simbotin, I., & Dalgarno, A. 1998, ApJS, 115, 293
  

 

  
9 Online Material

Appendix A: Model results for classes A1, A2, B and C and objects 1, 2 and 3

In the following we present the results of the models that fit observations at the 3$\sigma$ level. The observational constraints for these models are listed in Tables 1 and 2.

 

 

Table A.1: Best fit models at the 3$\sigma$ limit for class A1, A2, B and C for b=1.0 and b=5.0. If it was not possible to match observations with models at the 3$\sigma$ limit, we have left a horizontal line (--).

	${\it o/p}_{\rm ini}=0.01$	${\it o/p}_{\rm ini}=1.0$	${\it o/p}_{\rm ini}=2.0$	${\it o/p}_{\rm ini}=3.0$	Observations
Class A1, b=1.0
Preshock density / cm-3	$5.0\times10^4$- $4.0\times10^5$	$5.0\times10^4$- $5.0\times10^6$	$1.0\times10^5$- $1.1\times10^7$	$1.0\times10^6$- $1.1\times10^7$
$v_{\rm s}$ / km s-1	22-41	11-40	15-34	10-17
Postshock density / cm-3	$1.6\times10^6$- $6.6\times10^6$	$1.5\times10^6$- $4.0\times10^7$	$2.6\times10^6$- $1.2\times10^8$	$1.3\times10^7$- $8.0\times10^7$
Width / AU	50-210	5-220	2-130	3-20
${\it o/p}_{\rm int}$	1.1-2.3	1.3-2.7	2.4-2.8	3.0
$\phi _{10}$	0.8-2.2	1.1-2.4	2.4-2.5	2.4-2.7	$1.8\pm0.3$
Class A2, b=1.0
Preshock density / cm-3	$5.0\times10^4$- $2.5\times10^5$	$5.4\times10^4$- $6.0\times10^5$	$5.4\times10^4$- $1.3\times10^7$	$1.0\times10^6$- $1.4\times10^7$
$v_{\rm s}$ / km s-1	28-43	21-43	15-43	10-19
Postshock density / cm-3	1.6$\times$106-5.3$\times$106	$1.6\times10^6$- $9.5\times10^6$	$1.6\times10^6$- $1.4\times10^8$	$1.4\times10^7$- $1.0\times10^8$
Width / AU	60-200	30-200	2-200	2-20
${\it o/p}_{\rm int}$	1.5-2.4	1.9--2.8	2.4-2.9	3.0
$\phi _{10}$	1.4-2.3	1.5-2.4	2.4-2.5	2.4-2.8	$2.4\pm0.3$
Class B, b=1.0
Preshock density / cm-3	$1.0\times10^5$- $6.0\times10^5$	$1.0\times10^5$- $4.6\times10^6$	$1.0\times10^5$- $1.2\times10^7$	$1.0\times10^6$- $1.1\times10^7$
$v_{\rm s}$ / km s-1	22-36	12-36	15-36	10-18
Postshock density / cm-3	$2.7\times 10 ^6$- $9.9\times10^6$	$2.7\times 10 ^6$- $4.0\times10^7$	$2.7\times 10 ^6$- $1.3\times10^8$	$1.3\times10^7$- $1.5\times10^8$
Width / AU	30-130	6-130	2-130	2-20
${\it o/p}_{\rm int}$	1.0-1.9	1.4-2.6	2.4-2.9	3.0
$\phi _{10}$	0.9-2.0	1.1-2.4	2.4-2.5	2.4-2.8	$1.8\pm0.6$
Class C, b=1.0
Preshock density / cm-3	2.5$\times$105-4.0$\times$105	$3.0\times10^5$- $4.5\times10^6$	$2.5\times10^5$- $1.5\times10^7$	$1.0\times10^6$- $1.5\times10^7$
$v_{\rm s}$ / km s-1	25-28	12-27	15-28	10-19
Postshock density / cm-3	$5.3\times10^6$- $7.5\times10^6$	$6.1\times10^6$- $4.0\times10^7$	$5.3\times10^6$- $1.7\times10^8$	$1.4\times10^7$- $1.1\times10^8$
Width / AU	40-60	6-60	2-60	2-20
${\it o/p}_{\rm int}$	1.3-1.5	1.4-2.3	2.4-2.7	3.0
$\phi _{10}$	1.3-1.6	1.1-2.2	2.4-2.6	2.4-2.8	$2.4\pm0.5$
Class A1, b=5.0
Preshock density / cm-3	$1.5\times10^5$- $3.5\times10^5$	$2.0\times10^5$- $5.0\times10^6$	$2.0\times10^5$- $5.0\times10^6$	$2.5\times10^5$- $5.0\times10^6$
$v_{\rm s}$ / km s-1	36-43	21-40	21-40	21-38
Postshock density / cm-3	$9.3\times10^5$- $1.8\times10^6$	$1.1\times10^6$- $1.4\times10^7$	$1.1\times10^6$- $1.4\times10^7$	$1.4\times10^6$- $1.4\times10^7$
Width / AU	280-540	30-440	30-440	40-370
${\it o/p}_{\rm int}$	0.6-0.8	1.2-1.8	2.1-2.5	3.0
$\phi _{10}$	1.1-1.9	1.0-2.1	1.7-2.2	2.3-2.4	$1.8\pm0.3$
Class A2, b=5.0
Preshock density / cm-3	--	$2.5\times10^5$- $4.0\times10^5$	$2.5\times10^5$- $4.5\times10^6$	$1.5\times10^5$- $4.5\times10^6$
$v_{\rm s}$ / km s-1	--	36-40	22-40	21-40
Postshock density / cm-3	--	$1.4\times10^6$- $2.0\times10^6$	$1.4\times10^6$- $1.3\times10^7$	$8.6\times10^6$- $1.3\times10^7$
Width / AU	--	250-360	40-360	40-70
${\it o/p}_{\rm int}$	--	1.7-1.8	2.2-2.5	3.0
$\phi _{10}$	--	1.9-2.1	1.8-2.3	2.2-2.3	$2.4\pm0.3$
Class B, b=5.0
Preshock density / cm-3	--	$5.5\times10^5$- $4.5\times10^6$	$5.5\times10^5$- $5.0\times10^6$	$5.5\times10^5$- $5.0\times10^6$
$v_{\rm s}$ / km s-1	--	22-34	22-33	21-33
Postshock density / cm-3	--	$2.6\times10^6$- $1.3\times10^7$	$2.6\times10^6$- $1.5\times10^7$	$2.6\times10^6$- $1.4\times10^7$
Width / AU	--	40-200	40-200	40-200
${\it o/p}_{\rm int}$	--	1.3-1.5	2.2-2.3	3.0
$\phi _{10}$	--	1.1-1.5	1.8-2.1	2.3-2.7	$1.8\pm0.6$
Class C, b=5.0
Preshock density / cm-3	--	--	--	--
$v_{\rm s}$ / km s-1	--	--	--	--
Postshock density / cm-3	--	--	--	--
Width / AU	--	--	--	--
${\it o/p}_{\rm int}$	--	--	--	--
$\phi _{10}$	--	--	--	--	$2.4\pm0.5$

 

 

Table A.2: Best fit C-type shock models at the 3$\sigma$ limit for objects 1, 2 and 3 identified in Table 2 and Fig. 15. If it was not possible to match observations with models at the 3$\sigma$ limit, we have left a horizontal line (--).

	${\it o/p}_{\rm ini}=0$.01	${\it o/p}_{\rm ini}=1.0$	${\it o/p}_{\rm ini}=2.0$	${\it o/p}_{\rm ini}=3.0$
Object 1, b=1.0
Preshock density / cm-3	--	$6.0\times10^5$- $5.0\times10^6$	$6.0\times10^5$- $1.5\times10^6$	$1.0\times10^6$- $1.35\times10^6$
$v_{\rm s}$ / km s-1	--	16-21	16-20	16-18	>18
Postshock density / cm-3	--	$9.5\times10^6$- $6.0\times10^7$	$9.0\times10^6$- $1.8\times10^7$	$1.3\times10^7$- $1.6\times10^7$
Width / AU	--	5-30	20-30	10-20	$80\pm30$
${\it o/p}_{\rm int}$	--	1.6-1.9	2.3-2.4	3.0
$\phi _{10}$	--	1.5-1.6	2.0-2.1	2.7-2.8	$1.6\pm0.3$
Object 2, b=1.0
Preshock density / cm-3	$5.0\times10^4$- $1.0\times10^5$	$5.0\times10^4$- $1.0\times10^5$	$5.0\times10^4$- $1.0\times10^5$	$5.0\times10^4$- $1.0\times10^5$
$v_{\rm s}$ / km s-1	34-41	34-40	35-40	34-40	>37
Postshock density / cm-3	$1.6\times10^6$- $2.6\times10^6$	$1.5\times10^6$- $2.6\times10^6$	$1.5\times10^6$- $2.7\times 10 ^6$	$1.5\times10^6$- $2.6\times10^6$
Width / AU	130-210	130-220	130-220	130-220	$80\pm30$
${\it o/p}_{\rm int}$	1.9-2.3	2.5-2.7	2.8-2.9	3.0
$\phi _{10}$	1.8-2.2	2.3-2.4	2.4-2.5	2.5-2.6	$1.8\pm0.4$
Object 3, b=1.0
Preshock density / cm-3	$7.0\times10^4$- $1.2\times10^5$	$5.0\times10^4$- $1.0\times10^5$	$5.0\times10^4$- $1.0\times10^5$	$5.0\times10^4$- $1.0\times10^5$
$v_{\rm s}$ / km s-1	33-37	34-38	34-38	34-38	>36
Postshock density / cm-3	$2.0\times10^6$- $3.0\times10^6$	$1.4\times10^6$- $2.6\times10^6$	$1.4\times10^6$- $2.6\times10^6$	$1.4\times10^6$- $2.6\times10^6$
Width / AU	120-170	130-230	130-230	130-230	$40\pm20$
${\it o/p}_{\rm int}$	1.9-2.2	2.5-2.7	2.8-2.9	3.0
$\phi _{10}$	1.8-2.0	2.3-2.4	2.4-2.5	2.5-2.6	$1.6\pm0.3$
Object 1, b=5.0
Preshock density / cm-3	--	--	--	--
$v_{\rm s}$ / km s-1	--	--	--	--	>18
Postshock density / cm-3	--	--	--	--
Width / AU	--	--	--	--	$80\pm30$
${\it o/p}_{\rm int}$	--	--	--	--
$\phi _{10}$	--	--	--	--	$1.6\pm0.3$
Object 2, b=5.0
Preshock density / cm-3	$2.5\times10^5$- $4.5\times10^5$	$2.5\times10^5$- $5.0\times10^5$	$2.0\times10^5$- $5.0\times10^5$	$2.0\times10^5$- $5.0\times10^5$
$v_{\rm s}$ / km s-1	36-40	34-40	34-40	34-40	>37
Postshock density / cm-3	$1.4\times10^6$- $2.3\times10^6$	$1.4\times10^6$- $2.4\times10^6$	$1.1\times10^6$- $2.4\times10^6$	$1.1\times10^6$- $2.4\times10^6$
Width / AU	220-360	210-360	210-440	210-440	$80\pm30$
${\it o/p}_{\rm int}$	0.6-0.7	1.5-1.8	2.3-2.5	3.0
$\phi _{10}$	1.2-1.7	1.4-2.1	2.1-2.2	2.3-2.7	$1.8\pm0.4$
Object 3, b=5.0
Preshock density / cm-3	2.5$\times$105-4.5$\times$105	$2.5\times10^5$- $5.5\times10^5$	$2.5\times10^5$- $5.0\times10^5$	$2.5\times10^5$- $4.5\times10^5$
$v_{\rm s}$ / km s-1	35-38	33-38	33-38	34-38	>36
Postshock density / cm-3	$1.4\times10^6$- $2.2\times10^6$	$1.4\times10^6$- $2.6\times10^6$	$1.4\times10^6$- $2.3\times10^6$	$1.4\times10^6$- $2.2\times10^6$
Width / AU	230-370	200-370	210-370	240-370	$40\pm20$
${\it o/p}_{\rm int}$	0.5-0.7	1.4-1.7	2.3-2.5	3.0
$\phi _{10}$	0.8-1.5	1.3-2.0	2.1-2.2	2.4-2.7	$1.6\pm0.3$