Excitation of Molecular Hydrogen in the Orion Bar Photodissociation Region From a Deep Near-Infrared IGRINS Spectrum

TL;DR

This study uses high-resolution near-infrared spectroscopy to analyze molecular hydrogen excitation in the Orion Bar PDR, providing detailed physical conditions and comparing them with previous findings.

Contribution

First detailed rovibrational H_2 emission line analysis of the Orion Bar PDR using IGRINS, constraining temperature and density with model fitting.

Findings

H_2 excitation temperature around 625 K

Hydrogen density estimated at 5x10^3 cm^-3

Results suggest lower density gas than previously reported

Abstract

We present a deep near-infrared spectrum of the Orion Bar Photodissociation Region (PDR) taken with the Immersion Grating INfrared Spectrometer (IGRINS) on the 2.7 m telescope at the McDonald Observatory. IGRINS has high spectral resolution (R~45000) and instantaneous broad wavelength coverage (1.45-2.45 microns), enabling us to detect 87 emission lines from rovibrationally excited molecular hydrogen (H_2) that arise from transitions out of 69 upper rovibration levels of the electronic ground state. These levels cover a large range of rotational and vibrational quantum numbers and excitation energies, making them an excellent probe of the excitation mechanisms of H_2 and physical conditions within the PDR. The Orion Bar PDR is thought to consist of cooler high density clumps or filaments (T=50-250 K, n_H = 10^5 - 10^7 cm^-3) embedded in a warmer lower density medium (T=250-1000 K, n_H=10^4 - 10^5 cm^-3). We fit a grid of simple constant-temperature and constant-density Cloudy models, which recreate the observed H_2 level populations well, to constrain the temperature to a range of 600 to 650 K and the density to n_H = 2.5 x 10^3 to 10^4 cm^-3. The best fit model gives T = 625 K and n_H = 5x10^3 cm^-3. This well constrained warm temperature is consistent with kinetic temperatures found by other studies for the Orion Bar's lower density medium. However, the range of densities well fit by the model grid is marginally lower than those reported by other studies. We could be observing lower density gas than the surrounding medium, or perhaps a density-sensitive parameter in our models is not properly estimated.

Tables

Table 1. Table 1: H 2 lines observed in the Orion Bar

${\lambda }_{\text{vacuum}}$	$\Delta \lambda$	H2 line ID	${\log }_{10}\left(F_i/F_r\right)$	S/N	$v_u$	$J_u$	$E_u/k$	${\log }_{10}\left(A_{ul}\right)$	$\ln \left(\frac{N_u}{g_u}/\frac{N_r}{g_r}\right)$	$N_u/N_m$
($\mu$m)	(${10}^{-6}$ $\mu$m)						(K)	[${\log }_{10}$(s-1)]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
2.406592	0.95	1-0 Q(1)	${1.273}_{-0.001}^{+0.001}$	817.4	1	1	6149	-6.37	${3.982}_{-0.001}^{+0.001}$	0.60
2.413439	1.67	1-0 Q(2)	${0.762}_{-0.001}^{+0.001}$	399.7	1	2	6471	-6.52	${3.743}_{-0.003}^{+0.002}$	0.81
2.223290	0.72	1-0 S(0)	${0.731}_{-0.001}^{+0.001}$	580.0	1	2	6471	-6.60	${3.771}_{-0.002}^{+0.002}$	0.83
2.423730	5.72	1-0 Q(3)	${1.022}_{-0.001}^{+0.001}$	488.6	1	3	6951	-6.56	${2.996}_{-0.002}^{+0.002}$	0.86
2.121834	0.00	1-0 S(1)	${1.193}_{-0.000}^{+0.000}$	1031.9	1	3	6951	-6.46	${3.035}_{-0.001}^{+0.001}$	0.90
2.437489	0.00	1-0 Q(4)	${0.417}_{-0.002}^{+0.002}$	225.6	1	4	7584	-6.58	${2.504}_{-0.004}^{+0.004}$	1.36
2.033758	-0.72	1-0 S(2)	${0.643}_{-0.001}^{+0.001}$	422.1	1	4	7584	-6.40	${2.436}_{-0.002}^{+0.002}$	1.27
1.957559	-2.62	1-0 S(3)	${0.959}_{-0.001}^{+0.001}$	339.8	1	5	8365	-6.38	${1.772}_{-0.003}^{+0.003}$	1.73
1.788050	-21.46	1-0 S(6)	$-{0.091}_{-0.006}^{+0.005}$	78.8	1	8	11521	-6.45	${0.099}_{-0.013}^{+0.013}$	0.60
1.747955	-2.86	1-0 S(7)	${0.285}_{-0.003}^{+0.003}$	132.7	1	9	12817	-6.53	$-{0.096}_{-0.008}^{+0.008}$	0.67
1.714738	-2.26	1-0 S(8)	$-{0.435}_{-0.015}^{+0.015}$	29.2	1	10	14220	-6.63	$-{0.532}_{-0.035}^{+0.034}$	0.40
1.687761	-3.93	1-0 S(9)	$-{0.202}_{-0.005}^{+0.005}$	89.3	1	11	15721	-6.78	$-{0.869}_{-0.011}^{+0.011}$	0.43
1.666475	-1.67	1-0 S(10)	$-{1.029}_{-0.034}^{+0.031}$	13.4	1	12	17311	-6.98	$-{1.306}_{-0.078}^{+0.072}$	0.32
1.650413	0.60	1-0 S(11)	$-{1.021}_{-0.029}^{+0.027}$	15.4	1	13	18979	-7.27	$-{1.789}_{-0.067}^{+0.063}$	0.30
2.355605	-2.38	2-1 S(0)	${0.000}_{-0.004}^{+0.004}$	97.8	2	2	12095	-6.43	${1.769}_{-0.010}^{+0.010}$	1.92
2.247716	1.67	2-1 S(1)	${0.465}_{-0.001}^{+0.001}$	376.5	2	3	12550	-6.30	${1.057}_{-0.003}^{+0.003}$	1.71
2.154216	-1.43	2-1 S(2)	$-{0.012}_{-0.003}^{+0.003}$	159.3	2	4	13150	-6.25	${0.645}_{-0.006}^{+0.006}$	1.30
2.073482	0.00	2-1 S(3)	${0.408}_{-0.001}^{+0.001}$	323.5	2	5	13890	-6.24	${0.244}_{-0.003}^{+0.003}$	1.73
1.679641	9.66	2-0 O(9)	$-{1.421}_{-0.069}^{+0.060}$	6.8	2	7	15763	-7.89	$-{0.692}_{-0.159}^{+0.137}$	1.18
1.522033	-0.60	3-1 O(5)	$-{0.235}_{-0.009}^{+0.009}$	46.6	3	3	17818	-6.70	$-{0.025}_{-0.022}^{+0.021}$	1.18
2.386471	-3.10	3-2 S(1)	$-{0.006}_{-0.006}^{+0.006}$	68.0	3	3	17818	-6.29	$-{0.004}_{-0.015}^{+0.015}$	1.21
1.581171	1.55	3-1 O(6)	$-{0.837}_{-0.029}^{+0.027}$	15.6	3	4	18386	-6.86	$-{0.163}_{-0.066}^{+0.062}$	1.18
2.287045	0.72	3-2 S(2)	$-{0.392}_{-0.006}^{+0.006}$	72.9	3	4	18386	-6.25	$-{0.181}_{-0.014}^{+0.014}$	1.16
2.201399	0.72	3-2 S(3)	$-{0.013}_{-0.003}^{+0.003}$	150.1	3	5	19086	-6.25	$-{0.645}_{-0.007}^{+0.007}$	1.55
2.128015	1.43	3-2 S(4)	$-{0.508}_{-0.006}^{+0.006}$	69.3	3	6	19911	-6.28	$-{0.813}_{-0.015}^{+0.014}$	1.49
2.065584	1.43	3-2 S(5)	$-{0.181}_{-0.004}^{+0.004}$	100.1	3	7	20856	-6.34	$-{1.182}_{-0.010}^{+0.010}$	1.95
1.509865	0.00	4-2 O(3)	${0.000}_{-0.006}^{+0.006}$	70.1	4	1	22079	-6.11	${0.000}_{-0.014}^{+0.014}$	1.00
1.563515	1.07	4-2 O(4)	$-{0.489}_{-0.010}^{+0.010}$	43.8	4	2	22352	-6.29	$-{0.107}_{-0.023}^{+0.023}$	1.00
1.622299	-2.74	4-2 O(5)	$-{0.328}_{-0.009}^{+0.009}$	47.4	4	3	22759	-6.44	$-{0.791}_{-0.021}^{+0.021}$	0.92
1.686462	2.26	4-2 O(6)	$-{1.045}_{-0.036}^{+0.033}$	12.7	4	4	23295	-6.58	$-{1.224}_{-0.082}^{+0.076}$	0.71
1.756281	0.00	4-2 O(7)	$-{0.754}_{-0.014}^{+0.014}$	30.6	4	5	23955	-6.73	$-{1.475}_{-0.033}^{+0.032}$	1.20
2.266764	2.15	4-3 S(4)	$-{1.007}_{-0.044}^{+0.040}$	10.5	4	6	24733	-6.39	$-{1.637}_{-0.100}^{+0.091}$	1.07
2.200974	0.72	4-3 S(5)	$-{0.690}_{-0.012}^{+0.012}$	35.7	4	7	25623	-6.49	$-{1.957}_{-0.028}^{+0.028}$	1.59
2.145873	-0.72	4-3 S(6)	$-{1.322}_{-0.047}^{+0.042}$	9.8	4	8	26616	-6.64	$-{2.124}_{-0.108}^{+0.097}$	1.03
2.099586	5.48	4-2 O(11)	$-{1.647}_{-0.104}^{+0.084}$	4.7	4	9	27706	-7.36	$-{2.449}_{-0.239}^{+0.193}$	1.45
2.100426	4.29	4-3 S(7)	$-{1.192}_{-0.033}^{+0.030}$	13.8	4	9	27706	-6.86	$-{2.542}_{-0.075}^{+0.070}$	1.33
1.549455	-2.62	4-2 Q(11)	$-{0.819}_{-0.023}^{+0.022}$	19.6	4	11	30139	-6.34	$-{3.372}_{-0.052}^{+0.050}$	0.72
1.560736	-1.55	5-3 O(2)	$-{0.509}_{-0.012}^{+0.012}$	35.8	5	0	26606	-5.65	$-{0.020}_{-0.028}^{+0.028}$	1.30
1.613520	1.55	5-3 O(3)	$-{0.211}_{-0.007}^{+0.007}$	64.9	5	1	26735	-5.95	$-{0.787}_{-0.016}^{+0.015}$	0.76
1.671814	8.94	5-3 O(4)	$-{0.671}_{-0.016}^{+0.016}$	27.3	5	2	26992	-6.12	$-{0.846}_{-0.037}^{+0.036}$	0.77
1.515792	-5.01	5-3 Q(4)	$-{0.614}_{-0.022}^{+0.021}$	20.1	5	4	27878	-6.13	$-{1.367}_{-0.051}^{+0.049}$	0.98
1.528648	1.07	5-3 Q(5)	$-{0.320}_{-0.013}^{+0.013}$	34.1	5	5	28498	-6.14	$-{1.954}_{-0.030}^{+0.029}$	1.30
2.057127	-4.77	5-3 O(9)	$-{1.313}_{-0.063}^{+0.055}$	7.4	5	7	30063	-6.81	$-{2.727}_{-0.146}^{+0.127}$	1.27
1.562635	-2.03	5-3 Q(7)	$-{0.552}_{-0.012}^{+0.012}$	36.9	5	7	30063	-6.17	$-{2.728}_{-0.027}^{+0.027}$	1.27
1.608398	-5.36	5-3 Q(9)	$-{0.628}_{-0.015}^{+0.015}$	29.3	5	9	32014	-6.19	$-{3.053}_{-0.035}^{+0.034}$	1.46
1.675032	1.55	6-4 O(2)	$-{0.747}_{-0.015}^{+0.014}$	29.6	6	0	30942	-5.55	$-{0.709}_{-0.034}^{+0.033}$	1.04
1.601534	-4.29	6-4 Q(1)	$-{0.348}_{-0.009}^{+0.009}$	49.0	6	1	31063	-5.85	$-{1.339}_{-0.021}^{+0.020}$	0.74
1.732641	-0.60	6-4 O(3)	$-{0.388}_{-0.012}^{+0.012}$	37.0	6	1	31063	-5.85	$-{1.354}_{-0.027}^{+0.027}$	0.73
1.536891	0.60	6-4 S(0)	$-{0.706}_{-0.019}^{+0.018}$	23.4	6	2	31303	-6.08	$-{1.091}_{-0.044}^{+0.042}$	0.90
1.607390	-1.67	6-4 Q(2)	$-{0.665}_{-0.015}^{+0.014}$	29.6	6	2	31303	-6.00	$-{1.134}_{-0.034}^{+0.033}$	0.86
1.796524	0.60	6-4 O(4)	$-{0.734}_{-0.025}^{+0.023}$	18.2	6	2	31303	-6.01	$-{1.169}_{-0.057}^{+0.054}$	0.83
1.501560	2.03	6-4 S(1)	$-{0.270}_{-0.011}^{+0.010}$	41.3	6	3	31661	-5.94	$-{1.872}_{-0.024}^{+0.024}$	0.95
1.616224	-1.67	6-4 Q(3)	$-{0.468}_{-0.009}^{+0.008}$	50.7	6	3	31661	-6.04	$-{2.031}_{-0.020}^{+0.020}$	0.81
1.628094	-0.60	6-4 Q(4)	$-{0.929}_{-0.026}^{+0.025}$	17.1	6	4	32132	-6.06	$-{2.197}_{-0.060}^{+0.057}$	0.73
2.029684	-3.34	6-4 O(7)	$-{0.998}_{-0.023}^{+0.022}$	19.4	6	5	32711	-6.39	$-{2.659}_{-0.053}^{+0.050}$	1.02
1.661304	4.41	6-4 Q(6)	$-{1.105}_{-0.031}^{+0.029}$	14.3	6	6	33394	-6.08	$-{2.889}_{-0.072}^{+0.067}$	0.85
1.708041	4.53	6-4 Q(8)	$-{1.244}_{-0.044}^{+0.040}$	10.4	6	8	35040	-6.11	$-{3.383}_{-0.101}^{+0.091}$	0.87
1.728799	0.00	7-5 Q(1)	$-{0.506}_{-0.013}^{+0.013}$	33.6	7	1	35057	-5.82	$-{1.702}_{-0.030}^{+0.029}$	0.78
1.735762	4.17	7-5 Q(2)	$-{1.032}_{-0.032}^{+0.030}$	14.1	7	2	35281	-5.97	$-{1.975}_{-0.074}^{+0.069}$	0.51
1.746280	-1.79	7-5 Q(3)	$-{0.670}_{-0.021}^{+0.020}$	21.7	7	3	35613	-6.01	$-{2.488}_{-0.047}^{+0.045}$	0.77
1.620548	1.55	7-5 S(1)	$-{0.539}_{-0.010}^{+0.010}$	44.8	7	3	35613	-5.93	$-{2.438}_{-0.023}^{+0.022}$	0.81
1.760446	-4.05	7-5 Q(4)	$-{1.097}_{-0.037}^{+0.034}$	12.3	7	4	36051	-6.03	$-{2.571}_{-0.085}^{+0.078}$	0.82
2.204989	-0.72	7-5 O(7)	$-{1.182}_{-0.036}^{+0.033}$	12.7	7	5	36588	-6.31	$-{3.188}_{-0.082}^{+0.076}$	0.87
1.561510	-1.55	7-5 S(3)	$-{0.607}_{-0.014}^{+0.014}$	30.8	7	5	36588	-5.86	$-{3.250}_{-0.033}^{+0.032}$	0.82
1.540006	-6.68	7-5 S(4)	$-{1.069}_{-0.048}^{+0.043}$	9.6	7	6	37220	-5.87	$-{3.382}_{-0.111}^{+0.100}$	0.72
1.523623	-0.95	7-5 S(5)	$-{0.717}_{-0.024}^{+0.023}$	18.7	7	7	37941	-5.89	$-{3.758}_{-0.055}^{+0.052}$	1.05
1.512240	2.50	7-5 S(6)	$-{1.251}_{-0.088}^{+0.073}$	5.4	7	8	38743	-5.95	$-{3.902}_{-0.203}^{+0.169}$	0.74
1.979270	-1.31	7-5 Q(11)	$-{1.361}_{-0.066}^{+0.058}$	7.0	7	11	41558	-6.18	$-{4.739}_{-0.153}^{+0.133}$	0.92
2.092904	7.63	7-5 Q(13)	$-{1.576}_{-0.068}^{+0.059}$	6.9	7	13	43693	-6.26	$-{5.154}_{-0.157}^{+0.136}$	0.85
2.041830	2.62	8-6 O(3)	$-{0.684}_{-0.012}^{+0.012}$	35.3	8	1	38708	-5.80	$-{2.005}_{-0.029}^{+0.028}$	0.83
2.210763	2.15	8-6 O(5)	$-{0.993}_{-0.037}^{+0.034}$	12.1	8	3	39219	-6.05	$-{2.888}_{-0.086}^{+0.079}$	0.56
1.763952	-4.05	8-6 S(1)	$-{0.864}_{-0.044}^{+0.040}$	10.4	8	3	39219	-5.97	$-{3.003}_{-0.101}^{+0.092}$	0.50
2.310167	1.67	8-6 O(6)	$-{1.529}_{-0.116}^{+0.091}$	4.3	8	4	39622	-6.17	$-{2.971}_{-0.266}^{+0.210}$	0.54
1.701803	-6.79	8-6 S(3)	$-{1.014}_{-0.033}^{+0.031}$	13.5	8	5	40116	-5.93	$-{3.935}_{-0.077}^{+0.072}$	0.52
1.664584	2.15	8-6 S(5)	$-{1.169}_{-0.038}^{+0.035}$	11.9	8	7	41355	-6.01	$-{4.441}_{-0.088}^{+0.081}$	0.76
2.172715	-1.43	9-7 O(2)	$-{1.280}_{-0.046}^{+0.042}$	9.9	9	0	41903	-5.57	$-{1.645}_{-0.106}^{+0.096}$	0.82
2.073187	-1.43	9-7 Q(1)	$-{1.088}_{-0.039}^{+0.036}$	11.6	9	1	41997	-5.91	$-{2.661}_{-0.090}^{+0.082}$	0.93
2.253724	1.67	9-7 O(3)	$-{0.969}_{-0.029}^{+0.027}$	15.4	9	1	41997	-5.85	$-{2.451}_{-0.067}^{+0.063}$	1.15
2.345581	-8.58	9-7 O(4)	$-{1.406}_{-0.082}^{+0.069}$	5.8	9	2	42185	-5.98	$-{2.518}_{-0.188}^{+0.159}$	0.51
1.987350	-4.05	9-7 S(0)	$-{1.401}_{-0.059}^{+0.052}$	7.9	9	2	42185	-6.18	$-{2.204}_{-0.135}^{+0.119}$	0.70
2.084098	9.06	9-7 Q(2)	$-{1.258}_{-0.076}^{+0.065}$	6.2	9	2	42185	-6.06	$-{2.103}_{-0.175}^{+0.149}$	0.77
2.100664	3.58	9-7 Q(3)	$-{1.237}_{-0.033}^{+0.031}$	13.7	9	3	42462	-6.11	$-{3.385}_{-0.076}^{+0.071}$	0.60
2.151876	3.58	9-7 Q(5)	$-{1.362}_{-0.042}^{+0.038}$	10.8	9	5	43274	-6.16	$-{3.971}_{-0.097}^{+0.088}$	0.69
2.230268	-4.53	9-7 Q(7)	$-{1.558}_{-0.088}^{+0.073}$	5.5	9	7	44392	-6.23	$-{4.545}_{-0.203}^{+0.168}$	0.90
1.548849	0.60	10-7 O(3)	$-{1.054}_{-0.056}^{+0.050}$	8.2	10	1	44903	-5.98	$-{2.722}_{-0.130}^{+0.115}$	0.79
2.176855	1.67	10-8 S(1)	$-{1.388}_{-0.046}^{+0.042}$	9.9	10	3	45317	-6.27	$-{3.314}_{-0.106}^{+0.096}$	0.93
1.648305	-1.67	10-7 O(5)	$-{1.273}_{-0.049}^{+0.044}$	9.3	10	3	45317	-6.29	$-{3.282}_{-0.114}^{+0.102}$	0.97
Columns are as follows. (1) The H2 line vacuum wavelength in $\mu$m calculated from the ground electronic state rovibrational energy levels in Komasa et al. (2011). See $\S$ 3.4 for more details. (2) The observed line centroid wavelength (in the Orion Bar rest frame) minus the expected theoretical line wavelength calculated from the level energies in Komasa et al. (2011) in units of ${10}^{-6}\mu$m ($\S$ 3.4). (3) H2 line rovibrational identifications in spectroscopic notation in the format “W-X Y(Z).” W and X denote the transition’s upper and lower $v$ states. Y denotes the change in $J$, where S is $\Delta J=-2$, Q is $\Delta J=0$, and O is $\Delta J=+2$. Z denotes the upper $J$ state. (4) The base 10 logarithm of the line flux $F_i$ normalized to the 4-2 O(3) reference line flux $F_r$ ($\S$ 3.5). (5) The signal-to-noise ratio for the line flux ($\S$ 3.5). (6) The transition’s upper vibrational state. (7) The transition’s upper rotational state. (8) The energy of the upper state $E_u$ above the ground ($v=0$, $J=0$) divided by the Boltzmann constant $k$ to convert the energies into temperature units ($\S$ 4.3). (9) The base 10 logarithm of the rovibrational radiative transition probability $A_{ul}$ from Wolniewicz et al. (1998), in units of $s^{-1}$ ($\S$ 4.2). (10) The natural logarithm of the column density in a transition’s upper state $N_u$ divided by the quantum degeneracy $g_u$, normalized to $N_r/g_r$ for the reference line 4-2 O(3) ($\S$ 4.2 and $\S$ 4.3). This is the value plotted in the excitation diagram shown in Figure 3. (11) The ratio of the observed column density of the transition’s upper state $N_u$ to the column density predicted by our best-fit model $N_m$ ($\S$ 5.2), as shown in the bottom of Figure 3.

Table 2. Table 2: Best-fit Cloudy Model Parameters

Parameter	Value
Constant Temperature	625 K
Constant Density ($n_H$)	$5\times {10}^3$ cm-3
Turbulence	2 km s-1 a
Abundances	Orionb
Grains	Orionb
Cosmic Ray Flux	Galactic Backgroundb
Incident Radiation Field (O7V star ${\theta }^1$ Ori C)	Kurucz Stellar Atmosphere model, $T_{\mathrm{eff}}=39700$ Ka
No. of Ionizing Photons from ${\theta }^1$ Ori C	$Q\left(\text{H}\right)=8.13\times {10}^{48}$ s-1a
Cloud Face Distance from ${\theta }^1$ Ori C	0.114 pca
Stopping Condition	$A_V=14$ mag
No. of Iterations	10
aParameters from Cloudy models of Orion Bar by Pellegrini et al. (2009) and Shaw et al. (2009).
bStored prescription in Cloudy.

Equations

$$w_{x,v}=(F_{x,v}^{1-0\mathrm{S}(1)}{)}^2.$$

$$w_{x,v}=w_{x,v}/\sum _x\sum _v{w}_{x,v}.$$

$$F=\sum _x\sum _v(w_{x,v}(F_{x,v}-B)).$$

$$N_u=\frac{F_{ul}}{\Delta E_{ul}hc{A}_{ul}},$$

$$g_u^{\mathrm{ortho}}=3(2J_u+1),g_u^{\mathrm{para}}=2J_u+1.$$

$$\left({{N_{u}\over g_{u}}\Big{/}{N_{r}\over g_{r}}}\right)\propto\exp\left(-{E_{u}\over kT}\right),\ \ \ln\left({{N_{u}\over g_{u}}\Big{/}{N_{r}\over g_{r}}}\right)=-{E_{u}\over kT},$$

Full text

(ISM:) photon-dominated region (PDR) – ISM: molecules – ISM: individual objects (Orion Bar) – infrared: ISM – techniques: spectroscopic

1 Introduction

Ultraviolet (UV) radiation regulates the process of star formation and the energetics, ionization state, and chemistry of the interstellar medium (ISM). Photodissociation or Photon-dominated Regions (PDRs) are regions in the ISM at the interface between hot ionized gas and cool molecular gas that are energetically dominated by non-ionizing UV photons. PDRs arise around regions of massive star formation or star death (Tielens & Hollenbach, 1985; Hollenbach & Tielens, 1997, 1999) and make up the bulk of the neutral ISM in star forming galaxies such as the Milky Way. In the most extreme cases, starburst galaxies can have much of their starlight reprocessed and reradiated by PDRs.

The canonical model for a PDR, as presented by Tielens & Hollenbach (1985), is a plane-parallel slab of gas illuminated on one side by stellar UV radiation. The interaction between UV photons and the gas sets up a differentiated structure that can be characterized by the phases of hydrogen, which transition from predominantly ionized (H+), to neutral atomic (H0), and then to molecular (H2). The H+/H0 and H0/H2 interfaces are called the “ionization” and “dissociation” fronts respectively. Extreme-UV (EUV) photons with energies above the H0 ionization potential ($>13.6$ eV, $\lambda<912$ Å) pass through the H+ zone and are absorbed by H0 at the ionization front. Far-UV (FUV) photons with energies just below the Lyman continuum ($11.2$–$13.6$ eV, $912<\lambda<1110$ Å) pass through the H0 zone but are attenuated by dust, elements with lower ionization potentials, and the Lyman and Werner bands of H2. Beyond the dissociation front, the remaining FUV photons are rapidly absorbed as the cloud transitions to cold molecular gas.

The H2 rotational ($J$) and vibrational ($v$), hereafter referred to as “rovibrational,” levels of the ground electronic state (Black & Dalgarno, 1976) can be excited by two main processes: UV excitation (fluorescence) and collisional (thermal) excitation. UV excitation occurs when FUV photons absorbed by H2 excite the molecules to upper electronic states (through the Lyman and Werner bands), from which $\sim 10\%$ of the molecules will dissociate (Field et al., 1966), and the rest decay into bound excited rovibrational levels (Black & Dalgarno, 1976; Black & van Dishoeck, 1987). Collisional excitation occurs via inelastic collisions with other particles in the gas that leave the molecules in the excited rovibrational levels of the ground electronic state. These rovibrationally excited molecules decay via either collisions or a radiative cascade of rovibrational transitions, which have wavelengths ranging from the mid-infrared to the optical. H2 is a homonuclear diatomic molecule lacking a permanent electronic dipole moment, so the rovibrational transitions occur only as electric quadrupole transitions ($\Delta~{}J=~{}0,~{}\pm 2$), which are optically thin under most conditions. Since the line flux from each transition is proportional to the number of molecules in an upper level, observing many lines allows us to calculated the detailed H2 rovibrational level populations. These emission lines have significant diagnostic power to reveal physical conditions within PDRs at the H0/H2 dissociation front where most of the emission arises.

The UV and collisional excitation and de-excitation processes give rise to two limiting cases for the rovibrational level populations: one that is thermal, and the other for pure UV radiative excitation. In dense and/or hot gas (such as in shocks), the rovibrational levels are excited and de-excited by frequent collisions and the level populations approach a thermal distribution. In low-density cool gas exposed to UV radiation, the level populations take on a distinctive non-thermal distribution that does not decline monotonically with increasing excitation energy. Observations of the rovibrational line flux ratios readily distinguish between these two limiting cases. However, in practice, many sources show observed level populations intermediate between UV excited and thermal H2. Two possible reasons are the superposition of spatially unresolved components, or collisional modification of UV excited H2. Sternberg & Dalgarno (1989) and Burton et al. (1990b) show that dense PDRs can exhibit level populations that are modified from the pure UV excited case by collisions. Collisions easily dominate the excitation of H2 into states with low energies above the ground, bringing these states into thermal equilibrium with the gas. States at higher energies are primarily UV excited but collisional de-excitation modifies their populations from the pure UV excited case. This collisional modification of the level populations in UV excited H2 can mimic the spectrum from an unresolved combination of pure UV excited and thermal components.

Dense interstellar PDRs are found in star forming regions where molecular clouds are exposed to UV radiation from newly formed hot massive stars. At a distance of $\sim 400-500$ pc (Schlafly et al., 2014), the Orion Nebula is the nearest example of such a high-mass region, and it serves as an archetype for the more distant star forming regions found elsewhere in the Milky Way and in starburst galaxies. The optically visible part of the Orion Nebula is an H+ (or H II) region where the massive OB-stars that make up the $\theta^{1}$ or Trapezium cluster have ionized the adjacent gas and carved out a blister or cavity shaped region on the surface of the Orion Molecular Cloud (Zuckerman, 1973; Genzel & Stutzki, 1989; O’Dell, 2001; O’Dell et al., 2009). The UV radiation field generated by the Trapezium cluster is fairly well constrained (Ferland et al., 2012), with the O7V star $\theta^{1}$ Ori C contributing most of the UV photons.

The southeastern edge of the blister is a dense ($n~{}\gtrsim~{}10^{5}$ cm*-3*) PDR called the “Orion Bar,” viewed nearly edge on (Tielens et al., 1993; Walmsley et al., 2000; Pellegrini et al., 2009). The H2 emission from its dissociation front is bright and has been well studied. The first observations of the H2 emission (e.g., Hayashi et al. 1985; Gatley & Kaifu 1987) found intermediate rovibrational level populations that they interpreted as a combination of pure UV and shock excitation in the H2. Later studies by Hippelein & Muench (1989), Burton et al. (1990a) and Parmar et al. (1991) found that the H2 line widths in the bar are narrow, with Local Standard of Rest (LSR) radial velocities matching the ambient molecular cloud, inconsistent with shocks, which typically exhibit significant lateral motion (e.g.. such as observed in Orion KL by Oh et al. 2016). These authors suggested instead that the emission arises from collisionally modified UV excited H2. Luhman et al. (1998) came to the same conclusion from their observations of 16 H2 lines in the Bar.

The large spatial scale of the Orion Bar suggests that it is supported in a state of quasi-hydrostatic equilibrium by magnetic pressure that counteracts the radiation pressure from the Trapezium stars (Pellegrini et al., 2009; Shaw et al., 2009). Others, such as Parmar et al. (1991), Goicoechea et al. (2016), and Salgado et al. (2016), argue that the Orion Bar is not in hydrostatic equilibrium but instead represents a slow moving ($<4$ km s*-1*) magnetohydrodynamic shock front of swept-up molecular material supported by a strong compressed magnetic field. Observations of the Orion Bar find that complex molecules in the far-IR, sub-millimeter, and radio (Danby et al., 1988; Simon et al., 1997; Young Owl et al., 2000; Batrla & Wilson, 2003; Lis & Schilke, 2003; Parise et al., 2009; Goicoechea et al., 2011, 2016) trace relatively cool dense gas ($T=50$ to $250$ K, $n_{H}=10^{5}$ to $10^{7}$ cm*-3*). Observations of the collisionally excited pure rotation ($v=0$) lines of H2 (Parmar et al., 1991; Allers et al., 2005; Shaw et al., 2009), ions such as C+ (Tielens et al., 1993; Tauber et al., 1994; Wyrowski et al., 1997), and excited molecules formed in the presence of rovibrationally excited H2 (Nagy et al., 2013) trace warmer lower density gas ($T=250$ to $1000$ K, $n_{H}=10^{4}$ to $10^{5}$ cm*-3*). The emerging consensus is that the Orion Bar PDR does not consist of a single homogenous slab of gas, but instead is composed of cooler dense molecular clumps or filaments embedded in a warmer lower density medium (Burton et al., 1990b; Parmar et al., 1991; Meixner & Tielens, 1993; Andree-Labsch et al., 2017). However, some have argued against the presence of dense clumps near the dissociation front where the H2 emission is strongest (Marconi et al., 1998; Allers et al., 2005). All these observations find that the Orion Bar gas is warmer than models predict, suggesting that an extra heating mechanism, not yet identified, is present. Several candidate heating mechanisms have been proposed including an enhanced flux of cosmic rays trapped by a strong magnetic field (Pellegrini et al., 2007, 2009; Shaw et al., 2009), a larger than expected number of photoelectrons from grains (Allers et al., 2005), X-rays emitted by young stars in the Orion Nebula (Shaw et al., 2009), or collisional de-excitation of formation pumped H2 (Le Bourlot et al., 2012).

In this paper, we use H2 to probe the physical conditions and processes in the zone of the Orion Bar dissociation front. Section 2 describes our deep near-infrared spectrum of the Orion Bar, taken at the location of the peak H2 surface brightness, with the Immersion Grating INfrared Spectrometer (IGRINS). In $\S$ 3, we describe the initial data reduction, wavelength calibration, flux calibration, telluric correction, method for extracting H2 line fluxes, and effects of dust extinction. We discuss how we convert the line fluxes into rovibrational level populations in $\S$ 4. Our spectrum contains a larger number of H2 rovibrational transition emission lines at higher spectral resolution than all previous near-IR studies of the Orion Bar. The lines are all observed simultaneously through the same slit and cover a wide range of H2 upper vibrational states from $v=1$ to $11$ with level energies up to 50,000 K above the ground state ($v=0,~{}J=0$). This gives us an excellent handle on the excitation mechanisms of the H2. In $\S$ 5, we compare the observed H2 rovibrational level populations to those predicted by Cloudy models (Shaw et al., 2005; Ferland et al., 2013), to check whether we can match the observed level populations in the Orion Bar, determine which models provide the best match, and discuss the implications of the temperature and density of the H2 emitting gas derived from the model fits. We present our summary and conclusions in $\S$ 6.

2 Observations

The data were taken with the IGRINS on the 2.7 m Harlan J. Smith Telescope at the McDonald Observatory on the night of 2014 October 24 UT. IGRINS is a near-infrared cross-dispersed echelle spectrometer that uses a silicon immersion grating to achieve high spectral resolution of $R=\lambda/\Delta\lambda\approx 45000$ or 7.5 km s*-1* in two separate H- and K-band channels ($1.45-2.45\ \mu$m; Park et al. 2014). Dark and flat calibration frames were taken with the internal IGRINS calibration unit at the beginning of the night. IGRINS contains a fixed slit that subtends approximately $1^{\prime\prime}\times 15^{\prime\prime}$ on the sky when the instrument is mounted on the 2.7 m telescope at McDonald Observatory. For the Orion Bar observations, IGRINS was rotated at the Cassegrain focus to set the slit position angle (PA) on the sky to 135∘ counterclockwise from north, perpendicular to the dissociation front. Pointing and tracking were performed with the IGRINS slit-viewing camera, which images a $\sim 2^{\prime}\times 3^{\prime}$ field surrounding the slit in the K-band every 10 s. We used the nearby star V1501 Ori at 05h35m15$\fs$55, -05∘24′14$\farcs$0 (J2000) for offslit guiding. The center of the slit was positioned at 05h35m19$\fs$73, -05∘25′26$\farcs$7 (J2000), within the maps from Allers et al. (2005). Figure 1 shows the finder chart and the IGRINS slit position and angle superposed on the Orion Bar. We took three 10 minute exposures on the target and three 10 minute exposures on the sky (30′ south and 30′ west of the target). The sky exposures were used to subtract various backgrounds such as telluric OH emission, H2O emission, thermal emission, bias signal, and dark current. For telluric correction and relative flux calibration, we observed the A0V star HD 34317 with four exposures of 60 s each, nodding the star between two positions along the slit, around the same airmass and time we observed the Orion Bar. We subtracted one nod position from the other to remove sky and background while preserving the signal at each position.

3 Data Reduction and Line Flux Extraction

3.1 Basic Data Reduction and Wavelength Calibration

To reduce the data, we run the data reduction pipeline (IGRINS Pipeline Package [PLP],111IGRINS Pipeline Package (PLP): https://github.com/igrins/plp. Lee 2015). The PLP finds and subtracts the vertical and horizontal medians of the inter-order pixels to remove the readout pattern from each raw frame. Flat frames are combined to create a master flat, which is used to correct pixel to pixel variations on the detector and to trace the aperture for each order in the echellograms. The A0V standard star spectrum is optimally extracted with a weighted sum along the positive and negative traces of the star, which are then summed into a single one-dimensional spectrum. The spectrum of the Orion Bar, which is spatially extended and fills the whole slit, is extracted by subtracting the sky frames from the science frames with no weighting. Cosmic rays are identified and masked by running the Python version of LA-Cosmic (van Dokkum, 2001) on the reduced science frames divided the square of their variance. We use OH emission lines from the sky frames as an initial estimate of the wavelength solution by fitting 2D Chebyshev polynomials to the H- and K-band echellograms. The polynomial solution is refined by cross-correlating the telluric absorption lines in the A0V standard star spectrum with predictions from an atmospheric model generated with the Telfit222Telfit: http://www.as.utexas.edu/~kgulliks/projects.html. code by Gullikson et al. (2014). The final wavelength solution has sub-pixel accuracy with a typical uncertainty of $\pm 0.5$ km s*-1* or $<6\times 10^{-6}$ $\mu$m at any given point in the spectrum. We confirm our final wavelength solution by finding no significant difference between it and solutions derived using a ThAr arc lamp or the OH sky emission lines.

For the rest of this paper, we carry out our calibrations and analysis using our publicly available “plotspec”333Plotspec: https://github.com/kfkaplan/plotspec. python code designed for analyzing the reduced 2D IGRINS spectra of emission line nebulae.

3.2 Telluric Correction and Relative Flux Calibration

Calculation of relative column densities of H2 rovibrational states requires reliable flux ratios for all observed lines across the full wavelength range covered by IGRINS. To obtain a reliable relative spectrophotometric flux calibration, we need to correct for atmospheric (telluric) absorption lines, atmospheric transmission, instrumental throughput, and detector response. A0V stars have a well-known continuum shape, broad H0 absorption lines, and weak metal lines, making them preferred standards for telluric correction and relative flux calibration in the near-IR.

We adopt a technique similar to that used for telluric correction and relative flux calibration in the SpeX data reduction package Spextool444Spextool: http://irtfweb.ifa.hawaii.edu/~spex/. (Vacca et al., 2003). We assume that every A0V star has a continuum shape similar to that of Vega, and modify the model spectrum of Vega vegallpr25.50000resam5 by R. Kurucz555R. Kurucz synthetic stellar spectra: http://kurucz.harvard.edu/stars.html. to match the spectrum of the observed A0V standard star. First, we mask out the H0 absorption lines in the model Vega spectrum and fit a cubic spline to the regions between the absorption lines to define the continuum. Next, we artificially redden this continuum to match the A0V standard HD 34317 using the near-IR extinction law from Rieke & Lebofsky (1985), with $E(B-V)$ calculated from the standard star’s B and V magnitudes. The H0 absorption lines in the spectrum of the standard star are fit by scaling and Gaussian smoothing the H0 lines in the model Vega spectrum and adding them to the artificially reddened continuum to create a synthetic spectrum of the standard star. This synthetic spectrum represents our estimate of the intrinsic spectrum of the standard star. Finally, we divide the IGRINS spectrum of the A0V standard HD 34317 by the synthetic spectrum to derive the counts-to-flux ratio at each wavelength, and simultaneously apply the telluric correction and relative flux calibration by dividing the science spectrum by this ratio.

3.3 Continuum and Residual OH Removal

The faint continuum arises from a combination of starlight scattered off dust grains along with free-free and bound-free emission from the ionized gas. The continuum is subtracted from each order using a two step robust median filtering technique. We start by assuming that the spatial profile of the continuum along the slit has the same shape for all wavelengths in an order, but that the total flux from the continuum may vary with wavelength. For the first step, we estimate the continuum’s spatial profile by finding the median for each row of pixels along the spectral axis. For the second step, we normalize the estimated spatial profile flux at each wavelength to the median of the surrounding $\pm 187$ pixels. The normalized median continuum spatial profile across each order is then subtracted from the order. This technique fits the continuum well, while ignoring narrow features in the spectrum, such as emission and absorption lines, bad pixels, or regions with a large amount of telluric absorption. After continuum subtraction, we splice all orders together into a single 2D spectrum on a single wavelength grid.

Since telluric OH emission lines vary in flux over time, OH line residuals are a possible source of confusion and could in some cases be misidentified as H2 lines. However, OH residuals can be easily identified using the list of Rousselot et al. (2000). Furthermore, their fluxes are roughly uniform in the spatial dimension along the 15″ slit, whereas the H2 lines vary in brightness across the slit (as seen, for example, in Figure 2). To minimize the effect of OH residuals, we apply a first order correction by taking the difference between the first and last sky frames to estimate the variability of the OH lines. We then scale and subtract the difference in sky frames from the science frames, removing most of the flux from OH residuals.

3.4 Line Wavelengths

The spectral resolution of IGRINS is large enough that we can use it to test the rovibrational energy levels used to calculate wavelengths for our line list. To correct for the net LSR, solar, and barycentric velocities, the line vacuum wavelengths were shifted by 8.0 km s*-1*, derived from fitting the centroid offset of the 1-0 S(1) line. The vacuum wavelengths for the H2 lines are calculated from the theoretical ground electronic state rovibrational energy levels given in Komasa et al. (2011), and reported in column 1 of Table LABEL:tab:coldens. We measure the difference between the observed line centroids and the theoretical wavelengths (given as $\Delta\lambda$ in column 2 of Table LABEL:tab:coldens). This difference is well within the wavelength calibration precision of $<6\times 10^{-6}$ $\mu$m for most lines, although 1-0 S(6), 2-0 O(9), 5-3 O(4), 7-5 S(4), 7-5 Q(13), 8-6 S(3), 9-7 O(4), and 9-7 Q(2) show somewhat larger deviations. We observe the same differences in wavelengths in other PDRs, confirming that these deviations are real and not caused by kinematics within the Orion Bar or some other issue. To ensure that all the line centroids are aligned in velocity space, we adjust the wavelengths in our line list to match these small deviations from the theoretical values.

Wavelengths calculated from the theoretical rovibrational energy levels of Komasa et al. (2011) provide much improved agreement between the observed and theoretical line wavelengths over previous values (observed and theoretically derived). For example, wavelengths calculated from the commonly cited rovibrational energy levels in Dabrowski (1984) differ from the observed line centroids by up to $10^{-4}$ $\mu$m which is well in excess of the precision of the IGRINS wavelength calibration and differs by 5-10 times more than the wavelengths calculated from Komasa et al. (2011).

3.5 Line Flux Extraction

We extract line fluxes by interpolating all the H2 rovibrational transition lines we observe in the Orion Bar with IGRINS onto a common position-velocity (PV) grid, on which we use a S/N weighted sum to calculate the flux for each line. PV maps for each line are created by linearly interpolating from wavelength space to a $\pm 100$ km s*-1* velocity grid of 1 km s*-1* wide pixels. H2 lines blended with other lines are removed from consideration. Figure 2 compares the PV diagrams of several of the lines we observed and illustrates our procedure for extracting the flux of each line.

To extract the line fluxes, we use a flux density weighting scheme designed to scale with S/N across each line profile. Figure 2 shows that the lines have similar profiles. We confirm that this is the case for all the H2 lines, by stacking multiple dim lines and comparing the stacked profile to the brightest observed H2 line, 1-0 S(1). We therefore use the bright 1-0 S(1) line as the basis for our weighting scheme, and calculate the weights $w_{x,v}$ by squaring the flux $F_{x,v}^{\rm 1-0\ S(1)}$ found in each pixel in position ($x$) and velocity ($v$) space for the 1-0 S(1) line:

[TABLE]

The weights are then normalized as follows:

[TABLE]

The background $B$ per pixel is determined from the median value of all pixels in the PV diagram that are $\leqslant 0.8\%$ the flux of the brightest pixel. The $0.8\%$ limit was chosen to ensure that no line flux ends up in the background determination. We subtract the background from the flux in each pixel $F_{x,v}$, multiply by the weights $w_{x,v}$, and then sum the result to get the extracted flux $F$:

[TABLE]

Each line extraction is visually inspected to ensure that it is a real feature. Lines that appear to be contaminated by blends, misidentifications such as OH residuals, or noise spikes are rejected. For propagation of the statistical uncertainties, the interpolation and extraction process is repeated for the variance reported by the PLP. Table LABEL:tab:coldens gives the fluxes for all lines with S/N $>3$.

4 Analysis

4.1 Effects of Dust Extinction

The dense molecular gas of the Orion Bar co-exists with copious amounts of dust. If there is enough dust in the foreground of the observed H2 emission, the differential extinction across the H and K bands could be significant enough to affect the line ratios we use to derive the rovibrational level populations. An effective way to measure extinction is to compare the observed to theoretical line flux ratios from pairs of lines arising from the same upper level that are widely separated in wavelength. Two such line pairs in our data with sufficient S/N and widely separated in wavelength are the 3-1 O(5)/3-2 S(1) transitions spanning $\lambda=1.55220-2.238645~{}\mu$m and the 3-1 O(6)/3-2 S(2) transitions spanning $\lambda=1.58115-2.28703~{}\mu$m. Assuming the near-IR extinction law from Rieke & Lebofsky (1985), the observed 3-1 O(5)/3-2 S(1) and 3-1 O(6)/3-2 S(2) line ratios give extinctions of $A_{V}=8.50$ and $8.00$ mag respectively (or $A_{K}=0.99$ and $0.93$ mag). We therefore apply an extinction correction of $A_{V}=8.25$ or $A_{K}=0.96$ to our spectrum before extracting line fluxes. This value of the extinction is consistent with the foreground extinction of $A_{V}\sim 1.3$ mag or $A_{K}\sim 0.15$ mag towards the ionized gas (Weilbacher et al., 2015), allowing for additional extinction between the ionized gas and the region of excited H2. Our value for extinction in the Bar is lower than the values of $A_{K}=2.3\pm 0.8$ mag and $2.6\pm 0.7$ mag for two regions in the Bar $\sim 22$″ NE of the slit measured by Luhman et al. (1998). However, it is possible that the internal extinction is variable depending on the chosen sightline, and that the bright H2 emitting region we targeted is a sightline with low internal extinction.

4.2 Calculating H2 Level Populations

The near-IR H2 lines are optically thin and the line fluxes are linearly proportional to the column density of molecules in the upper states of the transitions. We calculate the column density of H2 in the upper state $N_{u}$ from the following equation,

[TABLE]

where $F_{ul}$ is the flux of the radiative transition from upper ($u$) to lower ($l$) rovibrational states, $\Delta E_{ul}$ is the difference in energy between the states in wavenumbers (cm*-1*), $A_{ul}$ is the transition probability (s*-1*) (we use the values from Wolniewicz et al. 1998, which are the same ones used in Cloudy), $h$ is Planck’s constant, and $c$ is the speed of light.

In the Orion Bar, we measure relative fluxes for 87 lines with S/N $>3$, yielding the relative $N_{u}$ values reported in Table LABEL:tab:coldens. These values are normalized to the population of the $v=4$, $J=1$ level, which is taken to be the reference level $r$, giving $N_{r}$ and $g_{r}$ We selected $v=4$, $J=1$ to be the reference level because it is primarily excited by UV photons and its population is derived from the bright 4-2 O(3) line. In many cases, there are multiple observed transitions arising from the same upper level (e.g. 1-0 S(1) and 1-0 Q(3)), providing independent measurements of $N_{u}$ for those upper states.

4.3 Excitation Diagram of H2 Level Populations

The top panel of Figure 3 shows an excitation (or Boltzmann) diagram for the relative H2 rovibrational level column densities (or level populations) we observe in the Orion Bar. This diagram is a plot of the logarithmic column density of a transition’s upper state $N_{u}$ divided by its quantum degeneracy $g_{u}$ vs. the excitation energy above the ground state ($v=0,J=0$), and is a convenient diagnostic tool for determining excitation mechanisms.

The spin of the two protons in H2 can be either aligned or anti-aligned, forming two distinct spin isomers called ortho-H2 (spins aligned) and para-H2 (spins anti-aligned). Since protons are fermions, the wave function of ortho-H2 can only have odd values of $J_{u}$, while para-H2 can only have even values of $J_{u}$. In collisional equilibrium, the statistical weights for nuclear spin give an ortho-to-para ratio of three. The value of $g_{u}$ depends on the upper rotation state $J_{u}$ and whether the H2 is ortho or para:

[TABLE]

H2 that is primarily excited and de-excited by collisions (e.g. as in gas heated by a shock) has thermal rovibrational level populations. In an isothermal region, the rovibrational level populations follow the Boltzmann distribution:

[TABLE]

where $E_{u}$ is the energy above the ground rovibrational state, $k$ is the Boltzmann constant, and $T$ is the kinetic temperature of the gas; in other words, the level populations follow a linear trend on an excitation diagram with a slope inversely proportional to $T$. If multiple temperature components are present, or there is a temperature gradient, the slope will flatten at higher excitation energies (e.g., Rosenthal et al. 2000). This occurs because hotter gas dominates the excitation of states at the highest energies above ground, while cooler gas dominates the excitation of states at the lowest energies.

UV excitation of H2 is a non-thermal process that leads to populations that do not show a monotonically decreasing trend for all the data-points on an excitation diagram, but instead follow a characteristic “sawtooth” pattern (see Figure 3). The bulk of the H2 in a PDR exists in the pure rotation $v=0$ states, which lie at low enough excitation energies that collisions thermalize their level populations so that they reflect the underlying kinetic temperature of the gas. UV excitation takes a small fraction of the underlying level populations of $J$ at $v=0$ and, in effect, transposes them to higher $v$. The “sawtooth” pattern occurs because because quantum selection rules limit changes in $J$ during radiative transitions, but there are no such restrictions to changes in $v$.

One can fit straight lines to a series of rovibrational states of constant $v$ to derive a “rotation temperature” or across states of constant $J$ to derive a “vibrational temperature,” but one should be careful not to confuse these quantities with the actual kinetic temperature of the gas. Instead, they are shorthand for characterizing the relative level populations. For UV excited H2, the level populations have high vibrational temperatures and lower rotation temperatures. While linear fits of these ladders (trends in constant $v$ or $J$) have been used in past studies of UV excited H2, they are not an ideal description for our information-rich data set, which probes up to high $J$ for many rotation ladders. For example, some of the data-points in the $v=1$ rotation ladder deviate from a linear fit by up to $\sim 5$ orders of magnitude. We therefore forgo the use of rotation or vibration temperatures in favor of comparing the rovibrational level populations we measure in the Orion Bar directly to values predicted by PDR models.

5 Modeling and Interpretation

5.1 Simulating H2 in

the Orion Bar with Cloudy

With IGRINS, we observe 87 NIR H2 emission lines, which arise from 69 independent rovibrational states with excitation energies up to $E_{u}/k=$50,000 K above the ground ($v=0,~{}J=0$) state. Our large dataset allows us to test our understanding of the physics in the Orion Bar by comparing the observed H2 rovibrational level populations to model predictions. For our models, we use version 13.03c of Cloudy666Cloudy: http://nublado.org. (Ferland et al., 2013), a one-dimensional plasma simulation code that solves for the physical conditions of a slab (or sphere) of gas irradiated by a photoionization source. It includes detailed physics for radiative transfer through the gas, and the state of the constituent ions, atoms, molecules, and dust, and predicts the physical conditions of the gas and the emergent spectrum. This version of Cloudy includes a fully self-consistent treatment of H2 including the excited electronic and rovibrational states, radiative and collisional excitation, photodissociation, and reformation on dust grains (Shaw et al., 2005).

Collisions in dense gas, such as the Orion Bar, can modify the H2 rovibrational level populations (Sternberg & Dalgarno, 1989; Burton et al., 1990b). For the Cloudy models, we have replaced the H2–H0 collision rate coefficients from Wrathmall et al. (2007) used by default in Cloudy 13.03c with updated values from Lique (2015). We use the default rates in Cloudy for H2–H2, H2–H+, and H2–He collisions. For collision rate coefficients that have no data (typically high $v$ and $J$), the “g-bar approximation” is used to estimate collision rate coefficients. The g-bar approximation assumes that the rate coefficient for a collisionally induced transition is a function of that transition’s change in energy (van Regemorter, 1962; Shaw et al., 2005). We do not observe significant radial motion in the PV diagrams (see Figure 2) expected from shock heated gas, consistent with the small line widths found by Burton et al. (1990a), Parmar et al. (1991), and Allers et al. (2005). This confirms that we are observing collisionally modified UV exited H2 as opposed to a combination of shocked and low-density UV excited H2. We compare the observed level populations to those for a grid of models with constant temperature and density, and to hydrostatic models of the Orion Bar. In $\S$ 5.2 and 5.3, we present our best-fit model, along with a brief comparison to hydrostatic models to illustrate their advantages and disadvantages.

5.2 Constant Temperature and Density Cloudy Models

5.2.1 Description of the Model Grid

While constant temperature and density models do not properly capture the structure of the full Orion Bar PDR from the ionized zone to the cold molecular regions, such simple models do reproduce the H2 rovibrational level populations within the narrow H2 emitting region. It is possible that the temperature and density are nearly uniform across the narrow observed emitting region, explaining why these models provide good fits. To explore the parameter space, we ran a grid of models with constant temperatures ranging from $T=200$ to $800$ K and constant densities ranging from $n_{H}=6.3\times 10^{2}$ to $10^{5}$ cm*-3*. The gas turbulence and incident radiation field (from the O7V star $\theta^{1}$ Ori C) used in these models are taken from the Orion Bar Cloudy models by Pellegrini et al. (2009) and Shaw et al. (2009). This model grid allows us to explore the effects of different values of temperature and density on the H2 rovibrational level populations. Increasing the density and temperature increases the rate of collisions in the gas. Each rovibrational level has a specific “critical density” for which the rate of collisional de-excitations equals the rate of radiative decays.

5.2.2 Effects of Collisions

Levels of low excitation energy, mainly the pure rotation states ($v=0$), where the majority of the H2 in a PDR lies, have low critical densities so their populations are primarily set by collisional excitation and de-excitation, which brings their populations into thermal equilibrium with the gas. The kinetic temperature(s) of the gas sets the kinetic energy of the collisions, so increasing the temperature raises the populations of the higher $J$ states. For an isothermal region, the Boltzmann distribution (Equation 6) describes these level populations. If the gas is warm enough, collisions can excite some of the molecules to $v=$ 1, 2, and maybe 3.

The populations of levels with high critical densities, typically those with high excitation energies at $v\geqslant 1$, primarily depend on UV excitation and the subsequent radiative cascade set by the transition probabilities ($A_{ul}$) and other physical constants that are mostly invariant to external variables such as the UV radiation field intensity (Black & van Dishoeck, 1987; Sternberg, 1988; Sternberg & Dalgarno, 1989). During UV excitation and the subsequent radiative cascade, quantum selection rules allow transitions with all values of $\Delta v$ but restrict the value of $\Delta J$ to 0 or $\pm 2$. In this sense, UV excitation crudely transposes the distribution of level populations for $J$ at $v=0$ to higher $v$. Since the gas kinetic temperature(s) sets level populations for $J$ at $v=0$, this shift of the $J$ level populations at $v=0$ to higher $v$ via UV excitation sets the relative column density of molecules at high $J$ for all $v$ values. Making the gas warmer increases the relative column density of molecules at high $J$ for all $v$ and compresses each rotation ladder vertically in the logarithmic space of the excitation diagram (Figure 3) while preserving the shape of a given rotation ladder (e.g., making the gas warmer vertically compresses the “bent knee” shape of the $v=1$ rotation ladder).

Increasing the temperature and/or density of the gas increases the rate of collisions, and this has differential effects on the populations of the lower energy $v=0$ and $1$ levels vs. the higher energy $v>1$ levels. As the collision rates increase, the level populations in $v>1$ become increasingly depressed by collisional de-excitation, while the $v=0$ and $1$ levels are less depressed since they are at low enough energy to also be collisionally excited. Increasing the temperature and/or density of the gas increases the rate of collisions, and increases the suppression of the level populations for $v>1$.

Increasing the rate of collisions due to higher temperature and/or density has another effect on the level populations. Collisionally induced transitions do not follow the same route to the ground level as the radiative cascade. Radiative transitions favor low-$J$ and are limited by quantum selection rules ($\Delta J=0$ or $\pm 2$), while collisionally induced transitions favor higher $J$ states and are not constrained by the same quantum rules. This raises the population of the high-$J$ levels in a given rotation ladder and “straightens” the shape of the rotation and vibration ladders as seen on the excitation diagram. This effect occurs simultaneously with, and at high density overwhelms, the compression of the rotation ladders caused by the redistribution of the collisionally excited $v=0$ levels to higher $v$ by UV excitation (e.g., the “bent knee” shape of the $v=1$ rotation ladder gets straightened into a monotonically decreasing trend).

5.2.3 Fitting the Model Grid

With our model grid, we fit our observations of the Orion Bar and pinpoint the gas temperature and density by leveraging the effects of collisions on UV excited H2. The main effect of higher density is to increase the collision rates. Temperature affects both the collision rates and the thermal populations of the $v=0$ ladder from which the relative level populations in $J$ are transposed to higher $v$ via UV excitation. Because of these dual effects of the temperature on the level populations, the model grid provides good leverage in fixing the gas temperature, while the density is less constrained.

We quantify the goodness of fit of the models to the data with a $\chi^{2}$ parameter of the logarithm of the data-to-model ratios $\sum\log_{10}(N_{\rm data}/N_{\rm model})^{2}$. This gives all the data points equal weight regardless of the large dynamic range in the level populations. Figure 4 shows a contour plot of $\sum\log_{10}(N_{\rm data}/N_{\rm model})^{2}$ for the grid of constant density vs. constant temperature models, which are marked as dots. The best-fit model, marked with a star in Figure 4, has $\sum\log_{10}(N_{\rm data}/N_{\rm model})^{2}=2.48$ with parameters of $T=625$ K and $n_{H}=5\times 10^{3}$ cm*-3*. Table 2 shows all input parameters for the best-fit model. The other models in the grid have identical parameters except for temperature and/or density. Figure 3 shows the Orion Bar data with the best-fit model’s predicted level populations, and column 11 in Table LABEL:tab:coldens gives the ratios of the data to the best-fit model. The level populations observed in the data and the best-fit model agree with each other within 0.5 dex.

As expected, we find that only a narrow range of temperatures, between 600 and 650 K, fits the data well. This temperature range is consistent with the warm gas ($T=250$ to $1000$ K) observed via other species thought to coexist with rovibrationally excited H2 in the Orion Bar, including excited pure rotation ($v=0$) lines of H2 (Parmar et al., 1991; Allers et al., 2005; Shaw et al., 2009), ions such as C+ (Tielens et al., 1993; Tauber et al., 1994; Wyrowski et al., 1997), and other excited molecules formed in the presence of excited H2 (Nagy et al., 2013). If the Orion Bar really does consist of a two phase medium with cooler dense clumps embedded in a warmer low-density medium (e.g., Burton et al. 1990b; Parmar et al. 1991; Meixner & Tielens 1993; Andree-Labsch et al. 2017), the H2 emission we observed arises from the warmer low-density gas.

The range of densities that fit the data is, again as expected, less well constrained than the temperature. We get good fits between $n_{H}=2.5\times 10^{3}$ to $10^{4}$ cm*-3*, which are marginally consistent with the values of $n_{H}\geqslant 10^{4}$ cm*-3* reported by nearly all other measurements and estimates from excited H2 and other species in the literature. If we assume pressure equilibrium where $P/k\sim 10^{8}$ cm*-3* K from Goicoechea et al. (2016) and T = 600 to 650 K from our model grid, we get a density of $n_{H}\sim 10^{5}$ cm*-3*. This is at least an order of magnitude greater than the densities of $n_{H}=2.5\times 10^{3}$ to $10^{4}$ cm*-3* best fit by the model grid. It is unclear why the best model fits have lower than expected densities. One possibility is that the UV excited H2 emitting gas at our slit position actually has lower density than the majority of the gas in the Orion Bar, and previous studies of the Orion Bar used species such as the pure rotation lines of H2 that trace the higher density gas. Perhaps we are viewing the lowest density part of the cloud face, where self shielding is lowest and the UV radiation field interacting with the H2 is strongest. Another possibility is that some density-sensitive parameter(s) in the Cloudy models, such as the H2 formation rate or the collisional rate coefficients, are overestimated or underestimated compared to their actual values.

The overall level populations across the different rotation ladders are well matched by the model, but the model over-predicts the observed level populations for high $J$ levels in the $v=1$ ladder. Le et al. (2016) found similar results when fitting models by Draine & Bertoldi (1996) to their IGRINS observations of rovibrationally excited H2 in the NGC 7023 PDR. One possible explanation for this discrepancy is “formation pumping,” where H2 forms on dust grains in excited rovibrational states. The distribution of rovibrational level populations for newly formed H2 assumed in the models might be over-predicting the observed level populations at high $J$. Cloudy assumes the prescription of Takahashi & Uehara (2001) for formation pumping. We ran a separate model grid with the formation pumping prescription of Draine & Bertoldi (1996) and another set of grids with the formation pumping prescription set to thermal (Boltzmann) distributions with temperatures of 1500, 5000, 10,000, and 17,329777The default thermal formation pumping prescription in Cloudy has a temperature of 17,329 K, corresponding to 1.5 eV or one-third of the energy released during the formation of an H2 molecule, as described in Le Bourlot (1991). K. We find that changing the formation pumping prescription in Cloudy does have a large effect on the predicted level populations at high $J$, but these alternate prescriptions do not provide better fits than the default Takahashi & Uehara (2001) prescription. The range of temperatures that best fit our data does not change significantly in any of these grids with alternative formation pumping prescriptions, but the range of best-fit densities approaches $n_{H}=5\times 10^{4}$ cm*-3* for the Boltzmann distribution prescriptions as the temperature is lowered from 17,329 to 1500 K. Since the high $J$ lines are sensitive to the adopted formation pumping prescription, using new formation pumping prescriptions or fine tuning existing prescriptions to fit the high $J$ levels might be an avenue for exploring formation pumping in future studies.

5.3 Hydrostatic Models

We ran a suite of hydrostatic Cloudy models of the Orion Bar, based on the models from Pellegrini et al. (2009) and Shaw et al. (2009). These models were designed to simulate the full structure of the Orion Bar PDR and the H2 emission. We ran these models with varying cosmic-ray fluxes, grain types, magnetic field strengths, temperature floors, and treatments for H2 collisions. While our best-fit constant density and constant temperature model fits the data well ($\sum\log_{10}(N_{\rm data}/N_{\rm model})^{2}=2.48$), all the hydrostatic models provided poorer fits of $\sum\log_{10}(N_{\rm data}/N_{\rm model})^{2}>10$.

The main complication we find is that leaving the g-bar approximation on (as defined in $\S$ 5.1) yields unphysical rovibrational level populations, making it necessary to disable this feature. Disabling the g-bar approximation means omitting some of the H2 physics (Shaw et al., 2005). This could introduce artificial effects between levels with well-known collision rate coefficients (mainly levels with $v\leqslant 3$) and those without, and it is unclear whether the predicted H2 level populations for these models are physically meaningful. Turning the g-bar approximation off has a negligible effect on our constant temperature and density model grid fits. Clearly, there exists an interdependence between the collisional processes for H2 and how the structure of the PDR is calculated in these hydrostatic models, that is less significant for the simpler constant temperature and density models. The hydrostatic model predictions for H2 rovibrational level populations would greatly benefit from well-known collisional rate coefficients for transitions between high $v$ and $J$ states. New and improved collisional data for H2 will ultimately give us a better understanding of PDR physics.

6 Summary and Conclusions

We observed the Orion Bar PDR in a deep pointed observation with IGRINS at the 2.7 m telescope at the McDonald Observatory. The instrument’s high spectral resolution of $R\sim$ 45,000 and broad wavelength coverage of the NIR H and K bands (1.45-2.45 $\mu$m) enables us to detect 87 H2 rovibrational transition emission lines with S/N $>3$. We extract the flux of each line with a robust weighting scheme and calculate the column density of H2 for a total of 69 different rovibrational states, which have excitation energies up to $E_{u}/k=$ 50,000 K above the ground state ($v=0$, $J=0$). The large range in rotational ($J$) levels, vibrational ($v$) levels, and excitation energy covered by the observed transitions allow us to perform a detailed study of the excitation of H2 within the Orion Bar PDR. We compare the observed rovibrational level populations to predictions from one-dimensional Cloudy 13.03c (Ferland et al., 2013) models.

As a result of our analysis, we find the following.

The spectral resolution of IGRINS ($R\sim$ 45,000) is high enough that the wavelengths for the H2 rovibrational transitions calculated from the experimentally determined H2 ground electronic state rovibrational energy levels in Dabrowski (1984) were found to differ from the observed wavelengths by up to $10^{-4}$ $\mu$m. New wavelengths calculated from the theoretical energy levels in Komasa et al. (2011) provide almost an order of magnitude improvement in the agreement between observed and calculated line wavelengths, with the majority of the lines agreeing to within the uncertainty of our wavelength calibration ($<6\times 10^{-6}$ $\mu$m). 2. 2.

The line-of-sight extinction towards the H2 emitting region is $A_{V}=8.25$ or $A_{K}=0.96$, as measured from line pairs arising from common upper states. 3. 3.

Constant temperature and density Cloudy models provide a better fit to the IGRINS H2 data than the hydrostatic models of Shaw et al. (2009) and Pellegrini et al. (2009), which explicitly solve for the structure throughout the PDR and have nearly constant pressure. This could be due to the fact that the $v\geqslant 1$ transitions we observe in the Orion Bar arise from a relatively narrow zone of the overall PDR structure. Another possible explanation for the poorer fit of the hydrostatic models is that this results from disabling the g-bar approximation for the collisional rate coefficients of the high-$v$ levels (van Regemorter, 1962; Shaw et al., 2005), which may omit physical effects important in determining the level populations. 4. 4.

The model grid, combined with the large number of rovibrational levels we probe, constrains the temperature for the observed H2 emitting region to 600 to 650 K, consistent with earlier findings. The best-fit model gives a temperature of 625 K. 5. 5.

The model grid constrains the density to $n_{H}=2.5\times 10^{3}$ to $10^{4}$ cm*-3*, with the best-fit model giving $n_{H}=5\times 10^{3}$ cm*-3*, which is marginally lower than most values in the literature. The reason may be either that this emission arises predominantly in the lower density inter-clump region of a two-component clumpy medium (with which our density is marginally consistent), or that one or more of the assumed parameters in the Cloudy models are sensitive to density and their values are over- or underestimated.

This work used the Immersion Grating INfrared Spectrometer (IGRINS) that was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the US National Science Foundation under grant AST-1229522 to the University of Texas at Austin, and of the Korean GMT Project of KASI. This paper includes data taken at The McDonald Observatory of The University of Texas at Austin. We acknowledge the Cambridge Astronomical Survey Unit and WFCAM Science Archive for making available data that were used for the finder charts in Figure 1. We would like to acknowledge Gary Ferland for helpful discussions on the Cloudy modeling, and Evelyne Roueff for pointing out the theoretical H2 ground electronic state rovibrational energy levels in Komasa et al. (2011) from which we derive improved line wavelengths for the rovibrational transitions ($\S$ 3.4).