# Many-body Theory Calculations of Positron Binding to Halogenated
Hydrocarbons

J. P. Cassidy<EMAIL_ADDRESS>School of Mathematics and Physics, Queen’s
University Belfast, University Road, Belfast BT7 1NN, United Kingdom  J.
Hofierka  School of Mathematics and Physics, Queen’s University Belfast,
University Road, Belfast BT7 1NN, United Kingdom  B. Cunningham  School of
Mathematics and Physics, Queen’s University Belfast, University Road, Belfast
BT7 1NN, United Kingdom  C. M. Rawlins  School of Mathematics and Physics,
Queen’s University Belfast, University Road, Belfast BT7 1NN, United Kingdom
C. H. Patterson  School of Physics, Trinity College Dublin, Dublin 2, Ireland
D. G. Green<EMAIL_ADDRESS>School of Mathematics and Physics, Queen’s
University Belfast, University Road, Belfast BT7 1NN, United Kingdom

###### Abstract

Positron binding energies in halogenated hydrocarbons are calculated _ab
initio_ using many-body theory. For chlorinated molecules, including planars
for which the interaction is highly anisotropic, very good to excellent
agreement with experiment and recent DFT-based model-potential calculations is
found. Predictions for fluorinated and brominated molecules are presented. The
comparative effect of fluorination, chlorination and bromination is elucidated
by identifying trends within molecular families including dihaloethylenes and
halomethanes based on global molecular properties (dipole moment,
polarizability, ionization energy). It is shown that relative to brominated
and chlorinated molecules, fluorinated molecules generate a less attractive
positron-molecule potential due to larger ionization energies and smaller
density of molecular orbitals close to the HOMO, resulting in very weak, or in
most cases loss of, positron binding. Overall, however, it is shown that the
global molecular properties are not universal predictors of binding energies,
exemplified by consideration of CH3Cl vs.  _cis._ -C2H2F2: despite the latter
having a larger dipole moment, lower ionization energy and similar
polarizability its binding energy is significantly smaller (25 meV vs. 3 meV,
respectively), owing to the important contribution of multiple molecular
orbitals to, and the anisotropy of, the positron-molecule correlation
potential.

Trap-based positron beams have enabled resonant-annihilation-based
measurements of positron binding energies for around 90 molecules [1, 2, 3, 4,
5, 6, 7, 8, 9, 10]. Whilst the corresponding theory of positron capture into
vibrational Feshbach resonances is well established [11, 12, 2], accurate
calculations of positron binding energies have been realised only relatively
recently (see e.g., [13, 14, 15, 16, 17, 18, 7, 19, 20, 21]). Attempts have
been made to relate the observed binding energies to the global molecular
properties including the dipole moment $\mu$, isotropic polarizability
$\alpha$ and ionization potential $I$ [3, 5, 20], but no such accurate
universal formula has yet been found. Recently we developed an _ab initio_
many-body theory (MBT) approach that quantified the role of strong many-body
correlations, and beyond the interplay of the global properties, highlighted
the importance of individual molecular orbital contributions to the positron-
molecule potential, e.g., the enhancement of binding due to $\pi$ bonds [21,
19] that was also deduced from experiment [3, 9, 10] 111It also successfully
described positron annihilation on small molecules [38].

Also recently, a model-polarization-potential method [16] was used to
calculate binding in chlorinated hydrocarbons, in a joint theory-experimental
study [7]. Although good agreement was found with experiment for many of the
molecules considered, for planar molecules the calculations substantially
overestimated the measured binding energies, with the suggestion that this was
due to the model assuming an isotropic long-range positron-molecule
interaction 222Specifically, their model constructed the positron-molecule
potential as a sum of positron-atom potentials, using a hybrid polarizability
of an atom in a certain chemical environment. At long range the IPP takes the
isotropic asymptotic form $-\alpha/2r^{4}$; while this is true for spherical-
top molecules, the true asymptotic form of the polarization potential is
generally anisotropic, and is given by $-1/2r^{6}\sum x_{i}x_{j}\alpha_{ij}$,
where the $x_{i}$ are the Cartesian coordinates and $\alpha_{ij}$ are the
Cartesian components of the polarizability tensor [16, 19]. By contrast, DFT-
model calculations for planar chloroethylenes [19] accounted for the
anisotropy approximately and found better overall agreement with experiment.
The method relied on an adjustable gradient parameter $\beta$, whose value the
authors of Ref. [19] were able to chose to replicate the binding energies of
dichloroethylenes to within around 10 meV, but this value led to
underestimated binding energies for tri- and tetrachloroethylene, at worst by
30 meV. The anisotropy of the positron-molecule potential, not captured by the
global molecular properties, is thus important, and _ab initio_ calculations
are demanded for fundamental understanding and description of the body of
experimental data.

The purpose of this Letter is twofold. First, we apply our many-body theory
approach 333We use the fixed-nuclei approximation: compared to the
correlations, vibrational effects have been found to have a relatively small
effect on the binding energy [14, 39, 40, 41, 42, 21] to study positron
binding in the chlorinated hydrocarbons considered in the recent model
calculations [7, 19] and experiment [7], accounting for the positron-molecule
correlations and anisotropic potential _ab initio_. We find very good
(excellent in cases) agreement with experiment and DFT-based model
calculations, including for the planar molecules. Secondly, we go beyond the
previous chlorinated studies [7, 19] and make predictions for fluorinated and
brominated molecules, and elucidate the comparative effects of fluorination,
chlorination and bromination. We find that compared to their brominated and
chlorinated counterparts, fluorinated molecules generate a successively less
attractive positron-molecule potential resulting in very weak or loss of
binding. We identify trends in binding based on global molecular properties
($\alpha$, $\mu$ and $I$) for families including the sequences of
cis/(Z)-dihaloethylenes C2H2Br2 $\to$ C2H2BrCl $\to$ C2H2Cl2 $\to$ C2H2ClF
$\to$ C2H2F2 444Here the descriptor $(Z)$ refers to $(E)$/$(Z)$ isomerism,
where $(Z)$ means that the highest priority groups, (i.e. the halogen atoms)
are on the same side of the C-C double bond., and halomethanes. However, we
find the global properties to be poor universal indicators of binding
energies, exemplified by CH3Cl and cis-C2H2F2 which have similar $\alpha$,
$\mu$ and $I$ but significantly different positron binding energies (25 meV vs
3 meV). We explain this and the overall results, and provide further
fundamental insight by considering the individual MO contributions to the
positron-molecule correlation potential, showing that e.g., the decrease (or
loss of) binding for bromination$\to$chlorination$\to$fluorination is due to
successively higher molecular orbital ionization energies and smaller density
of states close to the HOMO.

Theoretical approach.—A detailed description of our MBT approach is given in
[21]. Briefly, we solve the Dyson equation [26, 27]
$(\hat{H}_{0}+\hat{\Sigma}_{\varepsilon})\psi_{\varepsilon}(\mathbf{r})=\varepsilon\psi_{\varepsilon}(\mathbf{r})$
self-consistently for the positron wave function
$\psi_{\varepsilon}(\mathbf{r})$ with energy $\varepsilon$. Here $\hat{H}_{0}$
is the zeroth-order Hamiltonian of the positron in the static (Hartree-Fock)
field of the molecule and $\hat{\Sigma}_{\varepsilon}$ is the positron self
energy (an energy-dependent, non-local correlation potential) [28]. We
calculate it using a diagrammatic expansion in electron-electron and electron-
positron interactions, see Fig. 1 of [21], involving three main diagram
classes: the GW diagram, which describes polarization, screening of the
electron-positron Coulomb interaction, and electron-hole interactions; the
virtual-positronium (vPs) formation ladder series, which describes the
temporary tunnelling of an electron to the positron, denoted
$\Sigma^{\Gamma}$; and the positron-hole repulsion ladder series, denoted
$\Sigma^{\Lambda}$ The significant enhancement and enabling of binding due to
these correlations were delineated in [21]. Here we quote results only for our
most sophisticated self-energy $\Sigma^{GW+\Gamma+\Lambda}$ 555We use screened
Coulomb interactions in the ladder diagrams and MO energies calculated in the
random phase approximation. The $\varepsilon_{b}$ in this approximation are
found to be within 4 meV of those calculated using bare Coulomb interactions
and HF MO energies for chlorinated and brominated molecules (maximum relative
error of $7.4\%$), and 0.1 meV for fluorinated molecules.. We expand the
electron and positron wave functions in Gaussian basis sets, using aug-cc-pVXZ
bases (X=T,Q) [30] on atomic centres as well as additional hydrogen aug-cc-
pVXZ bases on “ghost” centres 1${\rm{\AA}}$ from the molecule to resolve
regions of maximum positron density. For all of the molecules considered, we
placed 5 ghosts around each halogen atom in the molecule in the shape of a
square-pyramidal cap, with each ghost 1Å from the halogen (see Supplemental
Material “SM”). We also use diffuse even-tempered positron bases of the form
$10s9p8d7f6g$, with exponents $\zeta_{0}\times\beta^{k-1}$
($\zeta_{0}=0.00001$– $0.006$ and $\beta=2$– $3$), ensuring the positron is
described well at large distances $r\sim 1/\kappa$, where
$\kappa=\sqrt{2\varepsilon_{b}}$. For molecules with $>2$ chlorines the
positron wave function is delocalized (Fig. 1), and we found that accurate
description of the vPs contribution requires a prohibitively large basis set
666The vPs contribution to the self energy requires diagonalization of dense
matrices of size $(N_{+}N_{-})^{2}$, where $N_{+}$ ($N_{-}$) is the number of
virtual positron (electron) states used. For example, our C2HCl3 calculations
have $N_{+}N_{-}=280800$, corresponding to $630$ GB of memory. With other
memory costs considered, our approach regularly requires more than 1 TB of
memory. (for our current computing resources), and our _ab initio_
calculations are not converged, though are lower bounds. Thus we also
performed MBT-based model calculations approximating $\Sigma\approx
g\Sigma^{(2)}+\Sigma^{(\Lambda)}$, using the second-order self-energy scaled
to approximate the virtual-Ps contribution as introduced and justified in
[21]: _ab initio_ calculations give $g$ in the range 1.4 to 1.5 for the HOMOs
[see [21] and also Fig. 2 (d).] This approach still calculates the anisotropic
polarization potential _ab initio_ , but is much less computationally
expensive.

Chlorinated molecules: comparison with experiment and model calculations.—Our
calculated positron binding energies $\varepsilon_{b}$ for the chlorinated
hydrocarbons considered in the recent isotropic-polarization-potential (IPP)
[7] and DFT model calculations [19] and experiment [7], and our predictions
for their fluorinated and (select) brominated counterparts are presented in
Table 1. Figure 1 summarizes this for the chlorinated molecules, and also
presents the calculated bound-state positron Dyson orbitals for chlorinated
and select chloro-fluorinated molecules, showing that the positron localizes
around the halogens. Overall, very good agreement is found between the _ab
initio_ MBT calculations and experiment. For CH3Cl, our calculated
$\varepsilon_{b}=25$ meV is in excellent agreement with both experiment and
the IPP model calculations. We find excellent agreement with experiment for
CH2Cl2, and for cis-C2H2Cl2 (for which both the IPP and DFT models
substantially overestimate) and trans-C2H2Cl2, and reasonable agreement for
vinylidene chloride C2H2Cl2. Overall, our ab initio results are in good
agreement with the DFT-based calculations [19] (including vinyl chloride, for
which there is no measurement). The results of the MBT-based model
calculation, which importantly augment our unconverged _ab initio_ results for
the molecules with $>2$ chlorines, are presented in the final column of Table
1. The model calculations with $g\sim 1.5$ generally give excellent agreement
with experiment (with the exception of ethylene).

Table 1:  Calculated MBT positron binding energies (meV) for halogenated hydrocarbons compared with experiment and model-potential calculations. For calculations denoted ‘$<0$’ binding was not observed. Where $\varepsilon_{b}<1$ meV, we quote values to 1 decimal place. Molecules marked ‘*’ are those for which we believe our ab initio calculations to be unconverged and we recommend the model-MBT result (final column, see text). Also shown are calculated HF dipole moments, isotropic dipole polarizabilities (calculated at the $GW$@BSE level) and ionization energies (calculated at the $GW$@RPA level and used in the energy denominators of the self-energy analytic expressions [21]).  |  |  |  |  | Present _ab initio_ MBT |  | Model-potential calculations
---|---|---|---|---|---|---|---
Molecule | Formula |  $\mu$ (D) |  $\alpha$ (a.u.) |  $I$ (eV) | $\Sigma^{GW+\Gamma+\Lambda}$ | Exp. [7] | IPP777Model-polarization-potential calculations of Swann and Gribakin, assuming isotropic asymptotic interaction [7]. | DFT 888DFT is the density-functional theory using the Perdew-Burke-Ernzerhof exchange functional result from Suzuki _et al._ [19]. | Present MBT-based model999Using a scaled self-energy $\Sigma=g\Sigma^{(2)}+\Sigma^{(\Lambda)}$ with $g$ ranging from 1.4 to 1.5 to account for vPs formation [21].
Methane | CH4 | 0 | 13.83 | 14.18 | $<0$ | – | $<0$ | – | $<0$
Chloromethane | CH3Cl | 2.15 | 27.80 | 11.78 | 25 | $26\pm 6$ | 29, 26 | – | 8–23
Dichloromethane | CH2Cl2 | 1.83 | 40.61 | 11.93 | 27 |  $32\pm$4 | 34, 30 | – | 15–31
Trichloromethane∗ | CHCl3 | 1.19 | 53.34 | 11.95 | 25∗ |  $37\pm$3 | 40, 34 | – | 16–37
Tetrachloromethane∗ | CCl4 | 0 | 64.14 | 12.02 | 29∗ |  $55\pm$10 | 55, 47 | – | 22–50
Ethylene | C2H4 | 0 | 24.40 | 10.75 | 1 | $20\pm 10$ | 5 | – | $<0$
Vinyl chloride | C2H3Cl | 1.68 | 38.67 | 10.57 | 27 | – | 54, 50 | 27 | 8–28
Vinylidene chloride | C2H2Cl2 | 1.62 | 51.04 | 10.50 | 41 |  $30\pm$5 | 79, 72 | 25 | 13–30
cis-1,2-dichloroethylene | C2H2Cl2 | 2.13 | 51.18 | 10.34 | 63 |  $66\pm$10 | 107, 99 | 80 | 43–75
trans-1,2-dichloroethylene | C2H2Cl2 | 0 | 52.79 | 10.27 | 10 |  $14\pm$3 | 29, 25 | 10 | 2–12
Trichloroethylene∗ | C2HCl3 | 1.01 | 64.84 | 10.16 | 35∗ |  $50\pm$10 | 84, 75 | 64 | 23–51
Tetrachloroethylene | C2Cl4 | 0 | 87.02 | 9.46 | – |  $57\pm$6 | 103, 92 | 54 | 34–70
1-chloro-1-fluoroethylene | C2H2ClF | 1.49 | 38.46 | 10.71 | 5 | – | – | – | 2–10
(Z)-chlorofluoroethylene | C2H2ClF | 2.37 | 38.52 | 10.53 | 32 | – | – | – | 22–39
Fluoromethane | CH3F | 1.94 | 15.56 | 13.99 | 0.3 | 0.3101010Molecule is VFR active, but $\varepsilon_{b}$ is too small to measure [32]. 0.3 meV was derived from the Z${}_{\text{eff}}$ fit of the VFR-based annihilation spectrum [33]. | – | – | 0.2–0.6
Difluoromethane | CH2F2 | 2.09 | 16.15 | 13.70 | 0.2 | – | – | – | 0.1–0.3
Trifluoromethane | CHF3 | 1.75 | 16.66 | 15.17 | $<0$ | – | – | – | –
Tetrafluoromethane | CF4 | 0 | 17.13 | 16.26 | $<0$ |  0– 111111CF4 is not VFR active [32]. | – | – | –
Vinyl fluoride | C2H3F | 1.47 | 26.11 | 10.92 | 0.3 | – | – | – | 0–0.6
Vinylidene fluoride | C2H2F2 | 1.30 | 26.22 | 10.88 | $<0$ | – | – | – | –
cis-1,2-difluoroethylene | C2H2F2 | 2.49 | 26.50 | 10.73 | 3 | – | – | – | 1–7
trans-1,2-difluoroethylene | C2H2F2 | 0 | 26.25 | 10.68 | $<0$ | – | – | – | –
Trifluoroethylene | C2HF3 | 1.37 | 26.43 | 10.75 | $<0$ | – | – | – | –
Bromomethane | CH3Br | 2.18 | 34.75 | 10.93 | 56 | 40121212From Ref. [32], where the uncertainty in the $Z_{\text{eff}}$ peak positions from which $\varepsilon_{b}$ was measured was reported to be between 10 and 15 meV. | – | – | 23–41
cis-1,2-dibromoethylene | C2H2Br2 | 1.97 | 64.67 | 10.09 | 109 | – | – | – | 58–108
(Z)-bromochloroethylene | C2H2BrCl | 2.04 | 57.72 | 10.20 | 80 | – | – | – | 42–87

Figure 1: Top: calculated positron binding energies compared with experiment
for chlorinated molecules: present MBT (red circles, and striped circles for
molecules difficult to converge _ab initio_); MBT-based model calculations
using $\Sigma=g\Sigma^{(2)}+\Sigma^{\Lambda}$ with $g=1.5$ (red squares);
(isotropic) polarization potential model calculation of [7] (black triangles;
two for each molecule reflecting two choices of cut-off parameter); DFT-model
calculation of [19] (blue squares). Vertical error bars are plus-minus the
maximum difference of our calculations using screened Coulomb interactions and
$GW$@RPA MO energies vs bare Coulomb interactions and HF MO energies [21].
Bottom: positron (Dyson) wave function at 80% maximum for chlorinated and
fluorinated molecules with $\varepsilon_{b}\geq 1$ meV; a) Chloromethane; b)
Dichloromethane; c) Trichloromethane; d) Tetrachloromethane (at 93%); e) Vinyl
chloride; f) Vinylidene chloride; g) cis-1,2-dichloroethylene; h)
trans-1,2-dichloroethylene (at 90%); i) Trichloroethylene; j)
cis-1,2-difluoroethylene; k) 1-chloro-1-fluoroethylene; l)
(Z)-chlorofluoroethylene. Figure 2: Dependence of positron binding energies on
global molecular properties and individual MOs. (a)-(c): calculated
$\varepsilon_{b}$ vs. calculated polarizabilities, dipole moments and
ionization energy for the brominated (orange), bromochlorinated (black),
chlorinated (red), chlorofluorinated (green) and fluorinated (blue) molecules;
symbols denote molecular families: squares are cis-dihaloethylenes C2H2XY,
triangles are halomethanes CH3X (X,Y = Br, Cl, F) and circles are isomers of
C2H2Cl2. Dashed lines are guides; (d) the positron-molecule correlation
strength parameters $\mathcal{S}_{n}^{\Gamma}$ (circles) and
$\mathcal{S}_{n}^{2+\Gamma}$ (squares), and the ratio
$g_{n}\equiv\mathcal{S}_{n}^{2+\Gamma}/\mathcal{S}_{n}^{2}$ (crosses) for each
MO $n$ against the MO HF ionization energies (vertical lines between panels)
for the cis-dihaloethylenes sequence [colours as in (a)-(c)]. (e) the
corresponding cumulative $\mathcal{S}^{2+\Gamma}$ obtained by summing from the
HOMO to the core orbitals. (f) the cumulative strength
$\mathcal{S}^{2+\Gamma}$ for CH3Cl (red; asterisks denote double degeneracy)
and cis-C2H2F2 (blue). (g) the calculated unenhanced ($\gamma_{i}=1$) and
enhanced contact densities for molecules with a $\Sigma^{GW+\Gamma+\Lambda}$
bound state. Colours and symbols as in (a)-(c), diamonds are remaining
molecules from Table 1.

Fluorinated molecules: predictions.—Compared to the chlorinated molecules, in
the fluorinated counterparts we find (see Table 1) that positron binding is
either lost or greatly reduced (as explained in the next section). We predict
bound states for fluoromethane, difluoroethylene, vinyl fluoride (a few tenths
of a meV each) and cis-1,2-difluoroethylene ($\varepsilon_{b}\sim 3$ meV).
Although fluoromethane is known to be VFR active, $\varepsilon_{b}$ was found
to be too small to measure [32]. However, our prediction of a weak bound-state
for fluoromethane of $\sim 0.3$ meV is in agreement with that derived from the
Z${}_{\text{eff}}$ fit of the annihilation spectrum of CH3F, which until now
had not been corroborated with any theoretical calculations [33]. This
contradicts a recent machine-learning-based prediction that fluoromethane does
not bind a positron [20]. Our prediction of a bound state for CH2F2 with
$\varepsilon_{b}=0.2$ meV concurs with the $0.4$ meV prediction by an earlier
empirical model [3]. Our lack of binding in CF4 is consistent with experiment;
this molecule is known to not be VFR active [32]. We also considered
1-chloro-1-fluoroethylene and (Z)-chlorofluoroethylene, and report binding
energies of 5 meV and 32 meV. These values lie between the fully chlorinated
and fluorinated binding energies (see below).

_Comparative effect of fluorination, chlorination and bromination; the role of
MO energies and density of states.—_ Figure 2 (a)-(c) show the calculated
$\varepsilon_{b}$ as a function of the global molecular properties $\alpha$,
$\mu$, and $I$ for the dihaloethylenes (cis/Z-C2H2XY and the isomers of
C2H2Cl2) and halomethanes CH3X, where X,Y$=$ F, Cl, or Br. These present three
distinct cases. Across the cis-dihaloethylenes $I$ and $\mu$ vary weakly, and
the increase in $\varepsilon_{b}$ going from X,Y$=$ F2 $\to$ ClF $\to$ … $\to$
Br2 follows an increase in $\alpha$: in a given family a more polarizable
target is more attractive to the positron. Across the halomethanes, $\mu$ is
almost constant, and the increase in $\varepsilon_{b}$ going from F to Br
follows both an increasing $\alpha$ and decreasing $I$ (the less tightly bound
electrons are more susceptible to perturbation from the positron). For the
isomers of C2H2Cl2, $\alpha$ and $I$ are vary weakly, and the decrease in
$\varepsilon_{b}$ from cis-C2H2Cl2 to vinylidene chloride to the non-polar
trans-C2H2Cl2 is due to successively decreasing $\mu$. These three distinct
cases highlight that the global molecular properties can explain trends in
$\varepsilon_{b}$ for families of molecules, but they are not reliable
universal predictors of binding energies, as exemplified by considering CH3Cl
and cis-C2H2F2. These have very similar $\alpha$, but whilst cis-C2H2F2 has a
larger $\mu$ and lower $I$, it has a lower binding energy (3 meV vs. 25 meV).
To explain this, and the reduction or lack of binding in fluorinated molecules
in general, we consider the individual molecular orbital contributions to the
correlation potential. We do so via the strength parameter
$\mathcal{S}=-\sum_{\nu}\bra{\nu}\hat{\Sigma}_{\varepsilon}\ket{\nu}/\varepsilon_{\nu}$
[34, 21], where $\nu$ is an excited positron Hartree-Fock (HF) orbital of
energy $\varepsilon_{\nu}$, with the self energy taken as
$\Sigma\approx\Sigma^{(2+\Gamma)}$, i.e., the sum of the bare polarization
$\Sigma^{(2)}$ and the virtual-Ps $\Sigma^{(\Gamma)}$ diagrams. Figure 2 (d)
shows $\mathcal{S}^{(\Gamma)}$, $\mathcal{S}^{(2+\Gamma)}$ and the ratio
$g=\mathcal{S}^{(2+\Gamma)}/\mathcal{S}^{(2)}$ for individual MOs as a
function of the MO energy for the sequence of cis-dihaloethylenes: the
strength parameters mainly decrease with increasing MO ionization energy
because more tightly bound orbitals are more difficult for the positron to
perturb [21]. Additionally, Fig. 2 (e) shows the cumulative
$\mathcal{S}^{(2+\Gamma)}$ obtained by summing from the HOMO to the core
orbitals. Moving from C2H2Br2 through to C2H2F2 sees both the total
$\mathcal{S}^{(2+\Gamma)}$ and the density of states near the ionization
energy decrease: e.g., in C2H2F2 there is a $\sim 5$ eV gap between the HOMO
and the HOMO$-1$, while this gap is approximately half as wide for C2H2Cl2 and
C2H2ClF and half as wide again for C2H2Br2. Further, the contributions to the
cumulative $\mathcal{S}^{(2+\Gamma)}$ below the HOMO for C2H2F2 are smaller
than those for the other three molecules as the MOs have larger $I$. In
general the transition from Br to Cl to F either shifts all the energy states
to more negative energies, or at least drives the sub-HOMO energies further
from the HOMO energy, inhibiting the molecule’s ability to bind the positron
(SM Fig. S1 shows MO energies of all molecules considered). We now consider
CH3Cl and cis-C2H2F2 [red triangle and blue square in Fig. 2 (a)-(c)]. Figure
2 (f) shows their cumulative $\mathcal{S}^{(2+\Gamma)}$ strength parameter. We
see that although CH3Cl has a larger $I$, its HOMO is doubly degenerate, and
contributes relatively more to the strength than the singly-degenerate HOMO of
CH2F2 (a second doubly degenerate state of $\pi$ character also contributes
strongly at $\sim 17$ eV for CH3Cl). Thus, in spite of CH3Cl having a smaller
dipole moment (which governs the strength of the static potential [35]), its
larger correlation potential (which contributes to binding non-linearly; see
Extended Data Fig. 3 of [21]) ultimately results in stronger binding.

Annihilation and contact densities.— The positron bound state annihilation
rate $\Gamma{\rm[ns^{-1}]}=50.47\,\delta_{ep}{\rm[a.u.]}$, where
$\delta_{ep}=\sum_{i=1}^{N_{e}/2}\gamma_{i}\int\lvert\varphi_{i}(\mathbf{r})\rvert^{2}\lvert\psi_{\varepsilon}(\mathbf{r})\rvert^{2}d^{3}\mathbf{r}$
is the electron-positron contact density. Here $\varphi_{i}$ is the $i$-th
electron MO, $\psi$ is the positron bound-state wavefunction (see e.g., Fig.
1) renormalized to
$a=\mathopen{}\mathclose{{}\left(1-\partial\varepsilon/\partial
E|_{\varepsilon_{b}}}\right)^{-1}<1$ [21] and $\gamma_{i}\geq 1$ are vertex
enhancement factors that account for short-range electron-positron attractions
[36, 37]. We found that they followed
$\delta_{ep}=(F/2\pi)\sqrt{2\varepsilon_{b}}$ [12] with $F=0.67$, remarkably
close to $F\approx 0.66$ for atoms [12], see Fig. 2 (g).

Summary.—Many-body theory calculations of positron binding to chlorinated
hydrocarbons were found to be in good to excellent agreement with experiment
and recent model-potential-based DFT calculations. Additionally, new
predictions elucidated the comparative effects of fluorination, chlorination
and bromination: trends within molecular families based on the global
molecular properties $\mu$, $\alpha$ and $I$ were identified, as was the
importance of describing the positron-molecule potential anisotropy, and
accounting for the energies and density of electron states (at least near the
HOMO). We suggest that any accurate universal formula for positron binding
energies should thus include these molecular properties. As well as providing
fundamental insight, our results provide benchmarks and can inform other
computational approaches to the positron-molecule and many-electron problems.

###### Acknowledgements.

_Acknowledgements_.—We thank Sarah Gregg, Gleb Gribakin and Andrew Swann for
useful discussions. This work was funded by the European Research Council
grant 804383 ‘ANTI-ATOM’ and a DfE Northern Ireland postgraduate research
studentship (JPC), and used the NI HPC Service and the ARCHER2 UK National
Supercomputing Service.

## References

  * Gilbert  _et al._ [2002] S. J. Gilbert, L. D. Barnes, J. P. Sullivan, and C. M. Surko, Vibrational-Resonance Enhancement of Positron Annihilation in Molecules, Phys. Rev. Lett. 88, 043201 (2002).
  * Gribakin  _et al._ [2010] G. F. Gribakin, J. A. Young, and C. M. Surko, Positron-molecule interactions: Resonant attachment, annihilation, and bound states, Rev. Mod. Phys. 82, 2557 (2010).
  * Danielson  _et al._ [2009] J. R. Danielson, J. A. Young, and C. M. Surko, Dependence of positron-molecule binding energies on molecular properties, J. Phys. B 42, 235203 (2009).
  * Danielson  _et al._ [2010] J. R. Danielson, J. J. Gosselin, and C. M. Surko, Dipole Enhancement of Positron Binding to Molecules, Phys. Rev. Lett. 104, 233201 (2010).
  * Danielson  _et al._ [2012a] J. R. Danielson, A. C. L. Jones, J. J. Gosselin, M. R. Natisin, and C. M. Surko, Interplay between permanent dipole moments and polarizability in positron-molecule binding, Phys. Rev. A 85, 022709 (2012a).
  * Danielson  _et al._ [2012b] J. R. Danielson, A. C. L. Jones, M. R. Natisin, and C. M. Surko, Comparisons of Positron and Electron Binding to Molecules, Phys. Rev. Lett. 109, 113201 (2012b).
  * Swann  _et al._ [2021] A. R. Swann, G. F. Gribakin, J. R. Danielson, S. Ghosh, M. R. Natisin, and C. M. Surko, Effect of chlorination on positron binding to hydrocarbons: experiment and theory, Phys. Rev. A 104, 012813 (2021).
  * Danielson  _et al._ [2021] J. R. Danielson, S. Ghosh, and C. M. Surko, Influence of geometry on positron binding to molecules, J. Phys. B 54, 225201 (2021).
  * Ghosh  _et al._ [2022] S. Ghosh, J. R. Danielson, and C. M. Surko, Resonant annihilation and positron bound states in benzene, Phys. Rev. Lett. 129, 123401 (2022).
  * Danielson  _et al._ [2022] J. R. Danielson, S. Ghosh, and C. M. Surko, Enhancement of positron binding energy in molecules containing $\pi$ bonds, Phys. Rev. A 106, 032811 (2022).
  * Gribakin [2000] G. F. Gribakin, Mechanisms of positron annihilation on molecules, Phys. Rev. A 61, 022720 (2000).
  * Gribakin [2001] G. F. Gribakin, Theory of positron annihilation on molecules, in  _New Directions in Antimatter Chemistry and Physics_, edited by C. M. Surko and F. A. Gianturco (Springer Netherlands, Dordrecht, 2001) pp. 413–435.
  * Tachikawa  _et al._ [2011] M. Tachikawa, Y. Kita, and R. J. Buenker, Bound states of the positron with nitrile species with a configuration interaction multi-component molecular orbital approach, Phys. Chem. Chem. Phys. 13, 2701 (2011).
  * Tachikawa [2014] M. Tachikawa, Positron-attachment to acetonitrile, acetaldehyde, and acetone molecules: Vibrational enhancement of positron affinities with configuration interaction level of multi-component molecular orbital approach, J. Phys.: Conf. Ser. 488, 012053 (2014).
  * Romero  _et al._ [2014a] J. Romero, J. A. Charry, R. Flores-Moreno, M. T. d. N. Varella, and A. Reyes, Calculation of positron binding energies using the generalized any particle propagator theory, J. Chem. Phys. 141, 114103 (2014a).
  * Swann and Gribakin [2018] A. R. Swann and G. F. Gribakin, Calculations of positron binding and annihilation in polyatomic molecules, J. Chem. Phys. 149, 244305 (2018).
  * Swann and Gribakin [2019] A. R. Swann and G. F. Gribakin, Positron binding and annihilation in alkane molecules, Phys. Rev. Lett. 123, 113402 (2019).
  * Swann and Gribakin [2020] A. R. Swann and G. F. Gribakin, Model-potential calculations of positron binding, scattering, and annihilation for atoms and small molecules using a gaussian basis, Phys. Rev. A 101, 022702 (2020).
  * Suzuki  _et al._ [2020] H. Suzuki, T. Otomo, R. Iida, Y. Sugiura, T. Takayanagi, and M. Tachikawa, Positron binding in chloroethenes: Modeling positron-electron correlation-polarization potentials for molecular calculations, Phys. Rev. A 102, 052830 (2020).
  * Amaral and Mohallem [2020] P. H. R. Amaral and J. R. Mohallem, Machine-learning predictions of positron binding to molecules, Phys. Rev. A 102, 052808 (2020).
  * Hofierka  _et al._ [2022] J. Hofierka, B. Cunningham, C. M. Rawlins, C. H. Patterson, and D. G. Green, Many-body theory of positron binding to polyatomic molecules, Nature 606, 688 (2022).
  * Note [1] It also successfully described positron annihilation on small molecules [38].
  * Note [2] Specifically, their model constructed the positron-molecule potential as a sum of positron-atom potentials, using a hybrid polarizability of an atom in a certain chemical environment. At long range the IPP takes the isotropic asymptotic form $-\alpha/2r^{4}$; while this is true for spherical-top molecules, the true asymptotic form of the polarization potential is generally anisotropic, and is given by $-1/2r^{6}\sumop\slimits@x_{i}x_{j}\alpha_{ij}$, where the $x_{i}$ are the Cartesian coordinates and $\alpha_{ij}$ are the Cartesian components of the polarizability tensor [16, 19].
  * Note [3] We use the fixed-nuclei approximation: compared to the correlations, vibrational effects have been found to have a relatively small effect on the binding energy [14, 39, 40, 41, 42, 21].
  * Note [4] Here the descriptor $(Z)$ refers to $(E)$/$(Z)$ isomerism, where $(Z)$ means that the highest priority groups, (i.e. the halogen atoms) are on the same side of the C-C double bond.
  * Dickhoff and Neck [2008] W. H. Dickhoff and D. V. Neck,  _Many-body Theory Exposed! - Propagator Description of Quantum Mechanics in Many-Body Systems \- 2nd ed._ (World Scientific, Singapore, 2008).
  * Fetter and Walecka [2003] A. L. Fetter and J. D. Walecka,  _Quantum theory of many-particle systems_ (Dover, New York, 2003).
  * Bell and Squires [1959] J. S. Bell and E. J. Squires, A formal optical model, Phys. Rev. Lett. 3, 96 (1959).
  * Note [5] We use screened Coulomb interactions in the ladder diagrams and MO energies calculated in the random phase approximation. The $\varepsilon_{b}$ in this approximation are found to be within 4 meV of those calculated using bare Coulomb interactions and HF MO energies for chlorinated and brominated molecules (maximum relative error of $7.4\%$), and 0.1 meV for fluorinated molecules.
  * Kendall  _et al._ [1992] R. A. Kendall, T. H. Dunning Jr, and R. J. Harrison, Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions, J. Chem. Phys. 96, 6796 (1992).
  * Note [6] The vPs contribution to the self energy requires diagonalization of dense matrices of size $(N_{+}N_{-})^{2}$, where $N_{+}$ ($N_{-}$) is the number of virtual positron (electron) states used. For example, our C2HCl3 calculations have $N_{+}N_{-}=280800$, corresponding to $630$ GB of memory. With other memory costs considered, our approach regularly requires more than 1 TB of memory.
  * Young and Surko [2008] J. A. Young and C. M. Surko, Feshbach-resonance-mediated positron annihilation in small molecules, Phys. Rev. A 78, 032702 (2008).
  * Gribakin and R. [2006] G. F. Gribakin and C. M. R. Lee, Positron annihilation in molecules by capture into vibrational feshbach resonances of infrared-active modes, Phys. Rev. Lett. 97, 193201 (2006).
  * Dzuba and Gribakin [1994] V. A. Dzuba and G. F. Gribakin, Correlation-potential method for negative ions and electron scattering, Phys. Rev. A 49, 2483 (1994).
  * Gribakin and Swann [2015] G. F. Gribakin and A. R. Swann, Effect of dipole polarizability on positron binding by strongly polar molecules, J. Phys. B 48, 215101 (2015).
  * Green and Gribakin [2015] D. G. Green and G. F. Gribakin, $\gamma$ spectra and enhancement factors for positron annihilation with core electrons, Phys. Rev. Lett. 114, 093201 (2015).
  * Green and Gribakin [2018] D. G. Green and G. F. Gribakin, Enhancement factors for positron annihilation on valence and core orbitals of noble-gas atoms, Concepts, Methods and Applications of Quantum Systems in Chemistry and Physics, Prog. Theor. Chem. and Phys. 31, 243 (2018).
  * Rawlins  _et al._ [2023] C. M. Rawlins, J. Hofierka, B. Cunningham, C. H. Patterson, and D. G. Green, Many-body theory calculations of positron scattering and annihilation in ${\mathrm{H}}_{2}$, ${\mathrm{N}}_{2}$, and ${\mathrm{CH}}_{4}$, Phys. Rev. Lett. 130, 263001 (2023).
  * Gianturco  _et al._ [2006] F. A. Gianturco, J. Franz, R. J. Buenker, H.-P. Liebermann, L. c. v. Pichl, J.-M. Rost, M. Tachikawa, and M. Kimura, Positron binding to alkali-metal hydrides: The role of molecular vibrations, Phys. Rev. A 73, 022705 (2006).
  * Romero  _et al._ [2014b] J. Romero, J. A. Charry, R. Flores-Moreno, M. Varella, and A. Reyes, Calculation of positron binding energies using the generalized any particle propagator theory, J. Chem. Phys. 141, 114103 (2014b).
  * Buenker and Libermann [2007] R. J. Buenker and H. Libermann, Role of the electric dipole moment in positron binding to the ground and excited states of the BeO molecule, J. Chem. Phys. 126, 104305 (2007).
  * Buenker and Liebermann [2008] R. J. Buenker and H.-P. Liebermann, Configuration interaction calculations of positron binding to molecular oxides and hydrides and its effect on spectroscopic constants, Nucl. Instrum. and Meth. B 266, 483 (2008).

in 1,2,3,4 See pages \x, of pos-bin-halogenated_prl_supplement.pdf


# JWST/MIRI detection of suprathermal OH rotational emissions: probing the
dissociation of the water by Lyman alpha photons near the protostar HOPS 370

David A. Neufeld William H. Miller III Department of Physics and Astronomy,
The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD, US  P.
Manoj Department of Astronomy and Astrophysics, Tata Institute of Fundamental
Research, Mumbai 400005, India  Himanshu Tyagi Department of Astronomy and
Astrophysics, Tata Institute of Fundamental Research, Mumbai 400005, India
Mayank Narang Academia Sinica Institute of Astronomy & Astrophysics, Taipei
10617, Taiwan  Department of Astronomy and Astrophysics, Tata Institute of
Fundamental Research, Mumbai 400005, India  Dan M. Watson University of
Rochester, Rochester, NY, US  S. Thomas Megeath University of Toledo, Toledo,
OH, US  Ewine F. Van Dishoeck Leiden Observatory, Leiden University, P.O. Box
9513, 2300 RA Leiden, NL  Max-Planck Institut für Extraterrestrische Physik,
Garching bei München, DE  Robert A. Gutermuth University of Massachusetts
Amherst, Amherst, MA 01003, US  Thomas Stanke Max-Planck Institut für
Extraterrestrische Physik, Garching bei München, DE  Yao-Lun Yang Star and
Planet Formation Laboratory, RIKEN Cluster for Pioneering Research, Wako,
Saitama 351-0198, Japan  Adam E. Rubinstein University of Rochester,
Rochester, NY, US  Guillem Anglada Instituto de Astrofísica de Andalucía,
CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, ES  Henrik Beuther Max
Planck Institute for Astronomy, Heidelberg, Baden Wuerttemberg, DE  Alessio
Caratti o Garatti INAF-Osservatorio Astronomico di Capodimonte, IT  Neal J.
Evans II Department of Astronomy, The University of Texas at Austin, Austin,
TX 78712, US  Samuel Federman University of Toledo, Toledo, OH, US  William J.
Fischer Space Telescope Science Institute, 3700 San Martin Drive, Baltimore,
MD 21218, US  Joel Green Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218, US  Pamela Klaassen United Kingdom Astronomy
Technology Centre, Edinburgh, GB  Leslie W. Looney Department of Astronomy,
University of Illinois, 1002 West Green St, Urbana, IL 61801, USA  National
Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903 US
Mayra Osorio Instituto de Astrofísica de Andalucía, CSIC, Glorieta de la
Astronomía s/n, E-18008 Granada, ES  Pooneh Nazari European Southern
Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany  John J.
Tobin National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville,
VA 22903 US  Lukasz Tychoniec Leiden Observatory, Leiden University, P.O. Box
9513, 2300 RA Leiden, NL  Scott Wolk Center for Astrophysics Harvard &
Smithsonian, Cambridge, MA, US

###### Abstract

Using the MIRI/MRS spectrometer on JWST, we have detected pure rotational,
suprathermal OH emissions from the vicinity of the intermediate-mass protostar
HOPS 370 (OMC2/FIR3). These emissions are observed from shocked knots in a
jet/outflow, and originate in states of rotational quantum number as high as
46 that possess excitation energies as large as $E_{U}/k=4.65\times 10^{4}$ K.
The relative strengths of the observed OH lines provide a  powerful diagnostic
of the ultraviolet radiation field in a heavily-extinguished region
($A_{V}\sim 10-20$) where direct UV observations are impossible. To high
precision, the OH line strengths are consistent with a picture in which the
suprathermal OH states are populated following the photodissociation of water
in its $\tilde{B}-X$ band by ultraviolet radiation produced by fast ($\sim
80\,\rm km\,s^{-1}$) shocks along the jet. The observed dominance of emission
from symmetric ($A^{\prime}$) OH states over that from antisymmetric
($A^{\prime\prime}$) states provides a distinctive signature of this
particular population mechanism. Moreover, the variation of intensity with
rotational quantum number suggests specifically that Ly$\alpha$ radiation is
responsible for the photodissociation of water, an alternative model with
photodissociation by a 104 K blackbody being disfavored at a high level of
significance. Using measurements of the Br$\alpha$ flux to estimate the
Ly$\alpha$ production rate, we find that $\sim{\color[rgb]{0,0,0}4}\%$ of the
Ly$\alpha$ photons are absorbed by water. Combined with direct measurements of
water emissions in the $\nu_{2}=1-0$ band, the OH observations promise to
provide key constraints on future models for the diffusion of Ly$\alpha$
photons in the vicinity of a shock front.

Protostars (1304), Molecular gas (1073), Infrared astronomy (786)

##  1 Introduction

In the realm of molecular astrophysics, one of the most remarkable results
obtained by Spitzer was the detection of highly suprathermal OH rotational
emissions. The observed transitions, detected with the Short-Hi module of the
Infrared Spectrometer (IRS) toward the Herbig-Haro object HH 211 (Tappe et al.
2008), originate in pure rotational states with rotational quantum numbers $N$
as high as 34 and energies as high as $E/k=2.8\times 10^{4}$ K. They are
naturally explained as the “prompt emission” that follows the
photodissociation of water via the $\tilde{B}-X$ band (also known as the
“second absorption band”) by radiation in the 114 – 134 nm wavelength range;
this spectral region includes the strong Ly$\alpha$ line emitted by fast
interstellar shocks. This explanation is supported by both laboratory and
theoretical studies of water photodissociation through the $\tilde{B}-X$ band,
which indicate that OH states as high as $N=47$ can be populated (Harich et
al. 2000, van Harrevelt & van Hemert 2000). Suprathermal OH emissions
resulting from the photodissociation of water were subsequently observed in
protostellar disks with Spitzer: the protostellar disk of DG Tau, in
particular, has been the subject of a detailed analysis by Carr & Najita
(2014).

Spitzer could not perform high-spectral resolution observations shortwards of
10$\mu$m, the short wavelength cutoff of the Short-Hi module of the IRS, and
at shorter wavelengths the Short-Lo module on Spitzer/IRS provided a spectral
resolving power, $\lambda/\Delta\lambda$, of only $\sim 60$, which was
insufficient to detect suprathermal OH emissions. By contrast, the MIRI MRS
spectrometer on JWST provides coverage down to the OH band-head at 9.13 $\mu$m
(and below), yielding spectra with a spectral resolving power $\sim 3000$.
This unique capability opens up the possibility of detecting suprathermal OH
emission in the 9 – 10 $\mu$m range, a possibility that has been realized in
observations of the Orion Bar reported very recently (Zannese et al. 2023),
providing a powerful test of model predictions for the spectrum of the OH
prompt emission (e.g. Tabone et al. 2021, hereafter T21).

In this Letter, we discuss JWST/MIRI observations of suprathermal OH emissions
in the vicinity of the protostar HOPS 370. HOPS 370, a.k.a. OMC2/FIR3, is an
intermediate-mass Class 0/I protostar (Furlan et al. 2016). It is located
north of the Orion Nebula in the OMC2 region of the integral shaped filament
at an estimated distance of 392 pc (Kounkel et al. 2018, Tobin et al. 2020,
hereafter T20). Its central protostellar mass, determined from Keplerian
motions, is $2.5\,M_{\odot}$, and its bolometric luminosity is
$314\,L_{\odot}$ (T20). Extensive observations of HOPS 370 have been carried
out with multiple observatories – including Herschel, SOFIA, VLA, ALMA, and
now JWST – and together reveal an actively accreting protostar with a bipolar
jet/outflow that is orthogonal to a rotating disk of estimated mass 0.05 – 0.1
$M_{\odot}$. It powers a large outflow traced in millimeter and far-IR lines,
which suggests that it is in a state of rapid accretion (T20; Manoj et al.
2013; González-García et al. 2016; Sato et al. 2023). This outflow consists of
both a wide-angle wind and a collimated jet, the latter containing shocks that
are also seen in non-thermal radio emission (Osorio et al. 2017). The
orientation of the disk, with an estimated radius of 100 au, indicates that
this source is observed at a high inclination angle of $\sim 72^{\arcdeg}$
(T20; Federman et al. 2023, and references therein). Luminous shocked knots in
the northern outflow lobe are characterized by strong emissions from a variety
of molecules and atomic ions detected in our observations, including H2, H2O,
CO, OH, Fe+, and Ne+.

In Section 2, we discuss the MIRI and NIRSpec observations carried out toward
HOPS 370 and the methods used to reduce the data. The resultant spectra and
spectral line maps are presented in Section 3, with particular emphasis on the
suprathermal OH emissions from the shocked knots. The origin of those
emissions is discussed in Section 4, in the context of a model in which water
is photodissociated by shock-produced Ly$\alpha$ radiation. A brief summary
follows in Section 5.

##  2 Observations and data reduction

The observations of HOPS 370 were performed as part of the Cycle 1 medium GO
program “Investigating Protostellar Accretion (IPA),” (PID 1802, Megeath et
al. 2021), which carried out NIRSpec and MIRI IFU observations toward five
protostars spanning five orders of magnitude in luminosity (see Federman et
al. 2023). A set of 2 x 2 mosaics was obtained with NIRSpec using the
G395M/F290LP disperser-filter combination, which provides coverage of the 2.87
– 5.10 $\mu$m spectral region at a spectral resolving power
$\lambda/\Delta\lambda\sim 1000$, and with all channels of the MIRI/MRS to
provide complete mid-infrared coverage from 4.9 to 27.9 $\mu$m at spectral
resolving power that ranged from 1500 to 4000 (Jones et al. 2023). The
mosaicking was performed with a 10$\%$ overlap and a 4-point dither pattern.
The total observing time was about 7.5 hours, including overheads. Further
details of the observing strategy have been presented by Narang et al. (2023)
and Federman et al. (2023).

For the reduction of NIRSpec IFU data, we utilized JWST pipeline version 1.9.5
and the JWST Calibration References Data System (CRDS) context version
jwst$\\_$1069.pmap. In our analysis, we identified hot pixels not captured by
the JWST outlier detection step by applying a custom outlier detection
algorithm specific to NIRSpec observations. More information on the NIRSpec
data reduction and the custom flagging routine can be found in Federman et al.
(2023).

The MIRI MRS data reduction utilized JWST pipeline version 1.12.5 along with
the JWST CRDS context version jwst$\\_$1179.pmap. We used the standard Stage 1
JWST pipeline Detector1Pipeline to reduce the MIRI MRS data starting from
uncal data.

In the subsequent Stage 2 (Spec2Pipeline), we performed pixel-by-pixel
background subtraction using dedicated background observations. This process
effectively removed all identified bad pixels, resulting in background-
subtracted cal products. However, we observed extended H2 S(1) and H2 S(2)
emissions in the dedicated background observations, which led to reduced flux
for these lines in the final data. Consequently, we repeated the Spec2Pipeline
without background subtraction. In this case, we encountered hot pixels in the
detector data, which we removed using the VIP package (Gomez Gonzalez et al.,
2017; Christiaens et al., 2023). Furthermore, we performed residual fringe
correction during Stage 2 for both scenarios, with and without background
subtraction.

In Stage 3 (Spec3Pipeline), the CubeBuildStep was set to band mode, generating
separate FITS files for each channel and band. We also generated data cubes
without dedicated background subtraction, with the outlier rejection function
turned off in these cases.

We measured and applied an astrometric offset calibration to the NIRSpec and
MIRI IFU data to improve feature alignment and link the coordinates to the
Gaia DR3 standard. The offset measurement process and listed offsets applied
with uncertainties are presented in Federman et al. (2023).

Additional data reduction tasks were performed using a suite of Python scripts
we developed to (1) extract spectra within a circular region of any specified
position and radius; (2) fit and subtract a continuum from the extracted
spectra; (3a) fit Gaussian lines to continuum-subtracted spectra obtained from
task (2) above; or (3b) fit Gaussian lines with a first-order baseline at each
IFU position and for each spectral line we targeted, thereby enabling us to
generate spectral line maps.

The second of these tasks (continuum fitting) was accomplished using a
procedure that lacked any knowledge of the wavelengths of expected spectral
lines. This “zero-knowledge” feature avoids the risk of artificially creating
spectral lines where lines are expected. Here, for each spectral channel, we
fit a third-order polynomial to the fluxes measured within a 17-channel window
centered on that spectral channel (i.e. with 8 spectral channels on either
side of the central one). The fit was optimized to achieve the best fit to any
10 of the 17 spectral channels in the window, and the continuum flux value for
the central channel was then assigned in accordance with that fit. For
spectral regions where lines cover less than 7/17 $\sim 40\%$ of the spectral
samples, this procedure yields a reliable separation of the continuum
(including instrumental baseline ripples) from the lines. For the third task,
Gaussian fitting, we used the Levenberg–Marquardt algorithm; here, the line
centroid and width were allowed to vary over a narrow range and the line
intensity was allowed to vary freely, as were the continuum level and slope
for task (3b).

##  3 Results

The IFU data acquired toward HOPS 370 are extraordinarily rich, revealing
literally hundreds of spectral lines with a signal-to-noise ratio adequate for
mapping. These data have and will be presented and discussed in series of
papers, some already published (Federman et al. 2023; Rubinstein et al. 2023;
Nazari et al. 2024; Brunken et al. 2024) and some in preparation. Here, we
focus on the suprathermal OH lines and a small set of ancillary lines that are
directly relevant to their interpretation.

###  3.1 Spectral line maps

In Figure 1, we present maps of several spectral lines: a strong well-isolated
water line within the $\nu_{2}=1-0$ vibrational band; the average of nine pure
rotational lines of OH, with upper states with $N_{U}$ between 34 and
43111Here, we excluded the OH $N=37-36$ line, which lies very close to the
much stronger S(3) line of H2; the Br$\alpha$ line at 4.05 $\mu$m, which
traces the Ly$\alpha$ radiation responsible for the photodissociation of water
to produce suprathermal OH emissions; the $v=0-0$ S(3) line of H2 at 9.66
$\mu$m, one of eight pure rotational lines detected with MIRI/MRS that may be
used to estimate the extinction toward the source; the [Fe II] 5.34 $\mu$m
line, a transition recently shown by Narang et al. (2023) to be an excellent
tracer of collimated jets in another IPA target source, IRAS 16253-2429; and
the [Ne II] 12.81 $\mu$m line, a signature of fast, ionizing shocks. The maps
are masked in the vicinity of a bright continuum source in the south, MIPS
2301 (Megeath et al. 2012), where the line fits are unreliable. The RA and Dec
offsets are given in arcsec relative to the ALMA source position (T20, green
star): RA = 83.865142 deg, Dec = –5.159561 deg (J2000).

Red circles near the lower right of each panel indicate the half power beam
width (HPBW) at the relevant wavelength, as determined by the linear fit given
by Law et al. (2023). All these emissions peak roughly 0.8′′ north of the ALMA
source position, near the location of the shocked knots identified by Federman
et al. (2023). The maps presented in Figure 1 exhibit a remarkable dynamic
range: they are shown with a logarithmic stretch extending down to 0.1$\%$ of
the peak intensity. In units of $\rm 10^{-4}\,erg\,cm^{-2}\,s^{-1}sr^{-1}$,
the peak velocity-integrated line intensities are 73, 1.14, 12.7, 36, 67, and
60 respectively for the H2O, OH, Br$\alpha$, H2, [Fe II] and [Ne II] lines.
For all the mapped lines other than the suprathermal OH emissions, a line is
securely detected in every spaxel right up to the edges of the mapped region.
While the [Fe II] and [Ne II] fine structure emissions primarily trace a
collimated bipolar jet, the H2 emissions are much less strongly collimated
(Federman et al. 2023) and appear to trace a wide-angle wind. The Br$\alpha$,
OH and H2O emissions show an intermediate degree of collimation. Velocity
shifts, although smaller than the instrumental linewidths, are clearly
detected and indicate that the northern jet is tilted towards us and the
southern jet away from us. They will be the subject of a future study.

###  3.2 OH suprathermal emission spectra

In Figure 2, we present the 8.8 – 13.4 $\mu$m spectra obtained toward the
shocked knots within the circular region indicated by the white circle in
Figure 1. This aperture has a radius of 0.8$\arcsec$ and is centered at a
projected distance of 316 au from the protostar (ALMA position) on the OH
emission peak at offset $(\Delta\alpha{\rm
cos}\delta,\Delta\delta)=(+0.1\arcsec,+0.8\arcsec)$. The spectral region shown
in Figure 2 covers 24 securely-detected lines of OH, originating in states
with $N_{U}$ ranging from 23 to 46, along with 5 fine structure lines of [Ni
II], [Co II], [Cl I] and [Ne II], and two pure rotational lines of H2, S(2)
and S(3).

Spectral line maps obtained toward HOPS 370, shown with a logarithmic stretch.
RA and Dec offsets are given in arcsec relative to the ALMA source position
(green star). The white circle demarks a 0.8′′ radius region centered on the
shocked knot. The maps are masked near a bright continuum source in the south
where the line fits are unreliable. The red circles shown the beam size
(HPBW).

8.8 – 13.4 $\mu$m spectra obtained toward the shocked knots. From top to
bottom: Band 2B spectrum with continuum fit in blue; continuum-subtracted Band
2B, 2C, and 3A spectra. Red numbers above the OH lines indicate the value of
$N_{U}$.

From top to bottom, separate panels show the observed MIRI Band 2B (Channel 2
sub-band B) spectrum, with the continuum fit (blue) obtained using the
procedure for task (2) described in Section 2 above; and the continuum-
subtracted Band 2B, 2C, and 3A spectra. Colored vertical lines at the bottom
of the lower three panels show the positions of suprathermal OH lines (Brooke
et al. 2016), following the color coding indicated in the second panel from
the top. Each rotational transition $N\rightarrow N-1$ is split into a quartet
of lines by the combined effects of lambda-doubling and spin-orbit coupling,
and the four lines are spectrally-resolvable by MIRI/MRS except at the highest
values of $N$.

Two of the four transitions connect so-called $A^{\prime}$ states, which are
symmetric with respect to reflection about the plane of rotation of the
molecule, and two connect antisymmetric $A^{\prime\prime}$ states. (A further
hyperfine splitting associated with the nuclear spin of H cannot be resolved
spectrally with MIRI/MRS for any of the observed transitions.) The observed
emission is completely dominated by intraladder transitions involving
symmetric $A^{\prime}$ states of OH (i.e. the lower $e$ lambda doublets of the
${}^{2}\Pi_{3/2}$ ladder and the lower $f$ lambda doublets of the
${}^{2}\Pi_{1/2}$ ladder, shown with green and red lines.)

In Figure 3, we show zoomed spectra of the suprathermal OH lines (yellow
histogram). Here, the black lines show Gaussian fits to each line. These were
obtained with the central wavelengths allowed to vary over a narrow range but
with the wavelength separation of the two components fixed at the laboratory
value and the flux ratio of the two components fixed at unity. The velocity
scale is referenced to the average wavelength of the two components. At the
spectral resolution of MIRI, the separation of the ${}^{2}\Pi_{3/2}(e)$ and
${}^{2}\Pi_{1/2}(f)$ transitions is unresolved for the highest-$N_{U}$ lines
detected and fully-resolved for $N_{U}\leq 23$. Line positions are marked with
vertical lines for each component of the OH quartet, with the same color-
coding as in Figure 2: only the A′ states (red and green) are detected. Some
transitions with $N_{U}<20$ are detected, but most lie in spectral regions
where flux measurements are unreliable due to instrumental fringing. They are
not plotted here, and their fluxes are not used in the analysis presented in
Section 4.1 below.

Spectra of suprathermal OH lines observed toward the shocked knots. Yellow:
observed spectrum. Line positions are marked with vertical lines for each
component of the OH quartet, with the same color-coding as in Figure 2. Black
histogram: Gaussian fit to the A′ components (see text).

###  3.3 $H_{2}$ rotational diagram and inferred extinction

While the H2 emissions from HOPS 370 will be discussed in detail in a future
publication, their present relevance is simply in providing a valuable
extinction estimate. Their usefulness for this purpose arises because the S(3)
line lies close to a local maximum in the extinction curve – associated with
the silicate absorption feature – and therefore provides excellent leverage on
the line-of-sight extinction. Using the intensities of the S(1) through S(8)
pure rotational lines of H2, measured with MIRI/MRS, we constructed the
rotational diagram shown in Figure 4. Here, we convolved the S(2) – S(8) MIRI
maps with 2D-Gaussian kernels of the widths needed to degrade the spatial
resolution to a common value for all lines. We then obtained average
intensities for each line within the circular aperture indicated by the white
circle in Figure 1.

Following Neufeld et al. (2006), for example, we fit the rotational diagram
with the sum of two components each in local thermodynamic equilibrium (LTE):
a warm component at temperature $T_{w}$, with an aperture-averaged column
density, $N_{w}$; and a hot component at temperature $T_{h}$, with an
aperture-averaged column density, $N_{h}$. These components were allowed to
have separate ortho-to-para ratios, $\rm OPR_{w}$ and $\rm OPR_{h}$, yielding
six free parameters to describe the rotational state of H2. The line-of-sight
extinction was treated as a seventh free parameter that was adjusted, along
with the other six, to optimize the fit (red and blue dashed curves). The
best-fit values are indicated on Figure 4, and are typical of other
protostellar outflows observed with Spitzer (e.g. Neufeld et al. 2006). 222The
positive curvature of the rotational diagram, which we account for with a
simple two-component model, suggests that a range of gas temperatures is
present (although the temperature distribution need not be bimodal, e.g.
Neufeld et al. 2009). The value of $\rm OPR_{w}$, lying significantly below
the value of 3 expected in LTE at temperature $T_{w}$, is suggestive of
transient heating in a shock wave; here, the OPR retains a fossil record of
its cooler preshock state, there having been insufficient time for it to reach
equilibrium (Neufeld et al. 2006, and references therein).

H2 rotational diagram obtained toward the shocked knots. Black points: no
reddening correction. Blue and red points: reddening correction applied. Blue
and red dashed lines: best fits to the rotational diagram for ortho and
para-H2.

Following Narang et al. (2024), we adopted the KPv5 extinction curve
(Pontoppidan et al. 2024) presented by Chapman et al. (2009), who cited a 2009
unpublished study by K. Pontoppidan for its origin and found that it provided
the best fit to the mid-IR extinction and ice features observed in the Spitzer
c2d program. The black points indicate the column densities in each rotational
state inferred without any extinction correction, while the red and blue
points show the values inferred from the extinction-corrected line fluxes. The
best-fit extinction optical depth at 9.7$\mu$m is
$\tau_{9.7}={\color[rgb]{0,0,0}1.84}$, and other fitting parameters are
specified in Figure 4.

To evaluate the sensitivity of our conclusions to our choice of extinction law
and aperture size, we have also analysed the H2 rotational emissions within an
aperture of radius 0.4$\arcsec$ instead of 0.8$\arcsec$ and for two additional
mid-IR extinction laws that have appeared in the literature. The results are
discussed in Appendix A, both as they pertain to the H2 analysis discussed
above and to the OH analysis discussed below. The relative OH lines fluxes
favor the KPv5 extinction curve over the alternative extinction laws
considered in Appendix A, but the primary conclusions of our study are similar
regardless of which mid-IR extinction law or aperture size we adopt.

##  4 Discussion

###  4.1 Relative strengths of the suprathermal OH emission lines

The high signal-to-noise ratio achieved in our observations of suprathermal OH
emissions facilitates a demanding test of theoretical models for their origin.
In Figure 5, we plot the mean photon intensity, extinction-corrected with the
KPv5 extinction curve, as a function of $N_{U}$ (red crosses with bars showing
the statistical errors). Here, we excluded the OH transitions with $N_{U}=37$,
28, and 25, which lie very close to the H2 S(3), [Cl I] 11.33 $\mu$m, and H2
S(2) lines, respectively.

The results are in excellent agreement with the predictions presented in
Appendix D of T21, which lists the number of OH line photons expected for each
value of $N_{U}$ divided by the number of water photodissociations in the
${\tilde{B}}-X$ band. These calculations, which rest upon theoretical
calculations of the photodissociation dynamics (van Harrevelt & van Hemert
2000) and on experimental measurements (Harich et al. 2000), were presented by
T21 for four different radiation fields. Those expected following water
photodissociation by Ly$\alpha$ radiation are shown by the blue curve. There
is only one free parameter in this comparison: an overall vertical scaling
that is proportional to the photodissociation rate within the beam. If every
available UV photon led to a water photodissociation via the ${\tilde{B}}-X$
band, the required UV photon intensity would be ${\color[rgb]{0,0,0}I_{\rm
UV}=\rm 1.67}\times 10^{9}\rm\,photons\,cm^{-2}\,s^{-1}sr^{-1}.$ With 24
observed line intensities, the number of degrees of freedom here was $N_{\rm
dof}=23$.

Our analysis here is closely-related to that of T21. The only difference is
that we present the minimum possible photon intensity, $I_{\rm UV}$, that
would account for the absolute intensities of the observed OH emissions if
every UV photon were absorbed locally by water. This photon intensity is a
factor 4$\pi$ smaller than the quantity $\Phi$ introduced by T21 and referred
to there as the column density of H2O photodissociated per second.

OH photon intensity as a function of $N_{U}$. Red points: observed
intensities. Also shown are the Tabone et al. (2021) predictions for H2O
photodissociation by Ly$\alpha$ radiation (blue) and by a 104 K blackbody
(red).

The observational intensities clearly show systematic errors that are not
fully captured by the statistical error bars. Assuming (1) that the predicted
curve for Ly$\alpha$ photodissociation (blue) represents the true behavior,
(2) that the statistical and systematic errors both have Gaussian
distributions with dispersions that may be added in quadrature, and (3) that
the fractional systematic error has the same r.m.s., $\epsilon$, for all
lines, we adjusted $\epsilon$ to achieve a reduced $\chi^{2}$ of unity for the
best-fit scaling. The result was $\epsilon$ = 0.105.

While the blue curve provides an excellent fit to the dependence of the line
strengths on $N_{U}$, one aspect of the T21 predictions is in conflict with
the observations. Whereas T21 predict roughly equal populations in the
symmetric (${}^{2}\Pi_{3/2}(e)$ and ${}^{2}\Pi_{1/2}(f)$) and antisymmetric
states (${}^{2}\Pi_{1/2}(e)$ and ${}^{2}\Pi_{3/2}(f)$) of OH, the observations
indicate that the symmetric $A^{\prime}$ states are strongly favored; indeed,
the antisymmetric $A^{\prime\prime}$ states are not detected and are at least
a factor $\sim 10$ less populated than the $A^{\prime}$ states. Regardless of
the relative rates at which the symmetric and antisymmetric states are
populated, the predictions presented in Appendix D of T21 are expected to
apply to the total emission in all four $N_{U}\rightarrow N_{U}-1$ transitions
(T21).

This behavior, also noted in the recent paper of Zannese et al. (2023), is in
fact entirely consistent with a recent theoretical study of the
photodissociation process by Zhou et al. (2015), which indicates a population
ratio $A^{\prime}/A^{\prime\prime}\sim 40$ at the Ly$\alpha$ photon energy
(their Figure 7). The astrophysical data thus provide a clear confirmation of
the molecular physics. A less-pronounced difference ($\sim$ factor 2) between
the line fluxes for the $A^{\prime}$ and $A^{\prime\prime}$ transitions had
previously been measured by Carr & Najita (2014) in Spitzer observations of
the protostellar disk in DG Tau. These authors discussed the effect in detail,
with reference to two possible origins for OH: photodissociation of H2O in the
$\tilde{B}-X$ band, and chemical pumping following formation via reaction of
O(${}^{1}D$) with H2. The larger difference observed in HOPS 370 may indicate
that chemical pumping is relatively less important in this source, at least
for the $N_{U}\geq 20$ transitions discussed here.

The $N_{U}$-dependence of the OH line intensities provides information about
the ultraviolet radiation field. The red curve in Figure 5 shows the
predictions given by T21 (2021) for a blackbody radiation field at 104 K
instead of a Ly$\alpha$ radiation field. These tend to overpredict the fluxes
for $N_{U}<30$ relative to those for $N_{U}>40$. For $\epsilon$ = 0.105, the
minimum reduced $\chi^{2}$ for this case is $\chi^{2}_{\rm
red}={\color[rgb]{0,0,0}2.98}$, implying that the blackbody radiation field is
disfavored at the $[N_{dof}(\chi^{2}_{\rm
red}-1)]^{1/2}\sigma={\color[rgb]{0,0,0}6.7}\sigma$ significance level.

As we did for the H2 emissions discussed in Section 3.3, we have also analysed
the OH emissions within an aperture of radius 0.4$\arcsec$ instead of
0.8$\arcsec$ and for two additional mid-IR extinction laws that have appeared
in the literature. The results are discussed in Appendix A.

###  4.2 Fraction, $f_{w}$, of Ly$\alpha$ photons absorbed by water

The number of Ly$\alpha$ photons available to photodissociate water may be
estimated from the Br$\alpha$ flux observed within the circular aperture
centered on the shocked knots. In this analysis, we assume a geometry in which
the Br$\alpha$ and OH emissions are generated within the shocked knots and
viewed directly rather than as a result of scattering, which is favored due
their location in distinct shock knots. After degrading the resolution of the
Br$\alpha$ map to have the same HPBW as the OH lines, we obtain a value of
$\rm 1.4\times 10^{-14}\,erg\,cm^{-2}\,s^{-1}$ for the Br$\alpha$ flux. If we
apply an extinction correction assuming the value of $\tau_{9.7}$ obtained in
Section 3.3 above, this corresponds to an intrinsic flux of $\rm 3.3\times
10^{-14}\,erg\,cm^{-2}\,s^{-1}.$ As discussed in Appendix B, shock models
appropriate for this source predict typical Ly$\alpha$/Br$\alpha$ luminosity
ratios of $\sim{\color[rgb]{0,0,0}900}$, a factor of several larger than the
Case B recombination ratio because collisional excitation preferentially
enhances Ly$\alpha$.

This would imply a Ly$\alpha$ flux within the aperture of
$\rm{\color[rgb]{0,0,0}3.0}\times 10^{-11}\,erg\,cm^{-2}\,s^{-1},$ or
equivalently 1.8 $\rm\,photons\,cm^{-2}\,s^{-1}$. This corresponds to an
aperture-averaged intensity of $\rm{\color[rgb]{0,0,0}3.8}\times 10^{10}$
$\rm\,photons\,cm^{-2}\,s^{-1}\,sr^{-1}$, a factor of
$\sim{\color[rgb]{0,0,0}23}$ times as large as the minimum intensity of UV
photons needed to account for the OH line fluxes (Section 4.1 above).
Therefore, only a fraction $f_{w}={\color[rgb]{0,0,0}4.3}\%$ of the available
Ly$\alpha$ photons would need to be absorbed by water to explain the OH
emission.

###  4.3 Interpretation of $f_{w}$

Ly$\alpha$ photons are unlikely to travel far without being absorbed by dust
or water. For Ly$\alpha$ radiation, we obtain a grain absorption cross-section
per H nucleus of $\sigma_{\rm abs}({\rm Ly}\alpha)=1.9\times 10^{-21}\,\rm
cm^{2}$, adopting the wavelength dependence and albedo given by KPv5; here,
the overall scaling was chosen to match the average $N_{H}/A_{J}$ ratio of
$5.6\times\rm 10^{21}cm^{-2}\,mag^{-1}$ determined by Vuong et al. 2003 from
X-ray absorption observations in several nearby dense clouds333This scaling is
also consistent with the standard $N_{H}/A_{V}$ ratio in diffuse molecular
clouds (Bohlin et al. 1978), but yields a $N_{H}/A_{K}$ ratio $\sim 40\%$
smaller than the average values determined toward diskless pre-main-sequence
stars in Serpens and Orion by Winston et al. (2010) and Pillitteri et al.
(2013). The water photodissociation cross-section for Ly$\alpha$ is $\rm
1.53\times 10^{-17}\,\rm cm^{2}$ (Heays et al. 2017, and references therein),
and thus the ratio of the water absorption rate to the grain absorption rate
for Ly$\alpha$ photons is $R=8\times 10^{3}\,x(\rm H_{2}O)$, where $x({\rm
H_{2}O})=n({\rm H_{2}O})/n_{\rm H}$ is the water abundance relative to H
nuclei. The corresponding fraction of Ly$\alpha$ photons absorbed by water is
$f_{w}=R/(1+R).$ The estimate of $R$ given above is critically dependent on
the (poorly known) properties of grains in the outflow. Indeed, it assumes
that grains are present in protostellar outflows – as suggested by Cacciapuoti
et al. (2024) and references therein – and moreover that their properties are
similar to those in the dense interstellar medium. The water abundance
required to explain a given value of $f_{w}$ scales linearly with the adopted
value of $\sigma_{\rm abs}({\rm Ly}\alpha)$.

If $f_{w}={\color[rgb]{0,0,0}0.043}$ as determined in Section 4.2, and given
the grain absorption cross-section assumed above, required water abundance is
${\color[rgb]{0,0,0}5}\times 10^{-6}$, amounting to only $\sim 1\%$ of the
gas-phase oxygen abundance 444As noted in Appendix B, our estimate of the
required water abundance is quite strongly dependent upon the adopted grain
properties: an alternative and widely-used grain model in the literature
yields a value of only ${\color[rgb]{0,0,0}1}\times 10^{-6}$.. This is the
average value, ${\bar{x}}(\rm H_{2}O)$, encountered by the Ly$\alpha$ photons
as they suffer repeated scatterings with H atoms and execute a random walk
prior to their eventual absorption. In the region of Ly$\alpha$ production,
the gas is warm ($T\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}6000$ K)
and/or ionized and the water abundance will be extremely small. But if the
photons escape the region where they are produced without being absorbed by
dust, then the water abundance could plausibly exceed $10^{-4}$ if all oxygen
nuclei were driven into gaseous water and $R$ could exceed unity. In this
scenario, the average water abundance is less meaningful, and the quantity
$f_{w}$ might primarily reflect the probability that a Ly$\alpha$ photon
escapes the warm region where it is originally generated and enters a region
where water is abundant. The transfer of Ly$\alpha$ radiation is a complex
process (e.g. Neufeld 1990) that can be profoundly affected by velocity shifts
associated with shock waves (Neufeld & McKee 1989). We defer a detailed
treatment of this process to a future study.

###  4.4 Lower limit on the water abundance from H2O $\nu_{2}=1-0$ emissions

The rovibrational water line map shown in the upper left panel of Figure 1
shows just one of several dozen emission lines detected in the H2O
$\nu_{2}=1-0$ band, which collectively have a total equivalent width of
$\sim{\color[rgb]{0,0,0}0.120}$ $\mu$m. Figure 6 shows the 5.8 – 7.0 spectral
$\mu$m region that is dominated by these emission lines. Unless the density is
extraordinarily high ($n_{\rm H}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{9}\,\rm cm^{-3}$), these lines
are too strong to be produced by collisional excitation. Colored symbols in
Figure 6 show the line positions, with stars denoting transitions of ortho-
water and crosses denoting those of para-water. A color code (top left)
indicates the minimum energies, $E_{\rm min}$, of the $v=0$ states that must
be pumped radiatively to excite each transition. We note here that $E_{\rm
min}$ may be smaller than the energy, $E_{L}$, of the lower state of the
observed rovibrational transition, since radiative pumping via a given
transition may be followed by radiative decay in a different transition of
longer wavelength. Roughly $90\%$ of the water emission emerges in transitions
that can be pumped radiatively out of the lowest 9 rotational states of water
(i.e. those with $J\leq 2$ and $E/k<160$ K). This behavior suggests a low
rotational temperature within H2O $v=0$ state, most likely because the states
are subthermally populated, and supports the hypothesis of radiative pumping.

Although a full treatment of the H2O $\nu_{2}=1-0$ emissions is beyond the
scope of the present study, we may obtain a lower limit on the mean water
abundance, $x({\rm H_{2}O})$, by assuming that the observed lines are
radiatively pumped by radiation from the protostar and that the observed
continuum is radiation from the protostar that has been scattered by dust. The
equivalent width of the water lines is then

$W_{\rm H_{2}O}=\sum{F_{\rm H_{2}O}\over
F_{c}}\leq\sum{\int\sigma_{\lambda}({\rm H_{2}O})d\lambda\over\sigma_{\rm
sca}}\biggl{(}{n_{l}({\rm H_{2}O})\over n_{\rm H}}\biggr{)}$ $None$

where the sum is taken over all lines in the $\nu_{2}=1-0$ band, $W_{\rm
H_{2}O}={\color[rgb]{0,0,0}0.120}\,\mu$m is the total equivalent width (summed
over all ${\rm H_{2}O}$ lines), $F_{\rm H_{2}O}$ is the wavelength-integrated
line flux for a given line, $F_{c}$ is the continuum flux at the line
wavelength, $\sigma_{\rm sca}$ is the grain scattering cross-section per H
nucleus, $n_{l}({\rm H_{2}O})$ is the number density of water molecules in the
lower state, and $\sigma_{\lambda}({\rm H_{2}O})$ is the ${\rm H_{2}O}$ cross-
section for a given rovibrational line (per molecule in the lower state),
which has an integral over wavelength $\lambda$ given by
$A_{ul}\lambda^{4}/(8\pi c)$, where $A_{ul}$ is the spontaneous radiative
rate. The equality in equation (1) applies only if the pumping lines are
optically-thin.

Average 5.8 – 7.0 $\mu$m spectrum obtained toward the shocked knots, with
rovibrational transitions of ortho- and para-water marked with stars and
crosses.

Because the KPv5 grain model suggests that $\sigma_{\rm sca}$ varies only
slowly over the band, we may take $\sigma_{\rm sca}$ as a constant and
approximate the sum of $A_{ul}\lambda^{4}n_{l}({\rm H_{2}O})$ as
$A_{band}\bar{\lambda}^{4}n({\rm H_{2}O})$, where $A_{\rm band}=24\,\rm
s^{-1}$ is the total spontaneous radiation rate for the band and
$\bar{\lambda}=6.3\,\mu$m is the average wavelength. Given the grain
scattering cross-section per H nucleus at $6.3\,\mu$m implied by KPv5,
$\sigma_{\rm sca}(6.3)=1.2\times 10^{-23}\,\rm cm^{-2}$, we then obtain555A
similar argument can be used to determine the CO abundance from observations
(Rubinstein et al. 2023) of the CO $v=1-0$ band: here, we obtain a total
equivalent width of 0.140 $\mu$m which implies a minimum CO abundance of
$1.1\times 10^{-4}$ relative to H nuclei.

$x({\rm H_{2}O})={n({\rm H_{2}O})\over n_{\rm H}}\geq{8\pi c\sigma_{s}W_{\rm
H_{2}O}\over A_{\rm band}\bar{\lambda}^{4}}=3\times 10^{-5}.$ $None$

Our water abundance of $3\times 10^{-5}$ is a lower limit, and – depending on
the water linewidths – would likely increase if optical depth effects are
included. But even this minimum estimate is a factor 6 larger than the value
needed to yield the inferred value of $f_{w}$ (Section 4.3 above). This
supports a picture in which most Ly$\alpha$ photons are absorbed by dust in a
warm and/or ionized zone very close to where they are created in a fast shock
and only a minority escape to the region of significant water abundance that
is probed by the water rovibrational emissions. Like the estimate of
${\bar{x}}({\rm H_{2}O)}$ derived in section 4.3 above, this independent
estimate of the minimum water abundance, derived from $W_{\rm H_{2}O}$, is
also dependent on the grain properties in the outflow; it scales linearly with
the value adopted for $\sigma_{\rm sca}(6.3)$. If grains were depleted in the
outflow (while maintaining the ratio of $\sigma_{\rm sca}(6.3)$ to
$\sigma_{\rm abs}({\rm Ly}\alpha)$), both water abundance estimates would
decrease proportionally. Future detailed analyses of the H2O $\nu_{2}=1-0$
spectrum and how it varies spatially will be needed to discriminate between
the various mechanisms that release water from ices into the gas-phase: these
include thermal desorption, sputtering in shocks, and UV photodesorption.

##  5 Summary

We have presented a study of OH in an outflow jet from the HOPS 370 protostar
observed with MIRI and NIRSpec as part of the IPA program.

1\. We have detected pure rotational, suprathermal OH emissions from the
vicinity of the intermediate-mass protostar HOPS 370 (OMC2/FIR3). These
emissions are observed from shocked knots in a jet/outflow, and originate in
states of rotational quantum number as high as 46 that possess excitation
energies as large as $E_{U}/k=4.65\times 10^{4}$ K. Only symmetric
$A^{\prime}$ states of OH are observed.

2\. The relative OH line strengths are entirely consistent with a picture in
which the suprathermal OH states are populated following the photodissociation
of water in its $\tilde{B}-X$ band by Ly$\alpha$ radiation produced locally be
a fast, ionizing shock. Photodissociation by a blackbody radiation field at
104 K is found to provide a significantly worse fit to the relative OH line
strengths.

3\. Using measurements of the Br$\alpha$ flux to estimate the Ly$\alpha$
production rate in shocked gas near HOPS 370, we find that
$\sim{\color[rgb]{0,0,0}4}\%$ of the Ly$\alpha$ photons are absorbed by water.

4\. The fraction of Ly$\alpha$ photons absorbed by water implies a mean water
abundance (relative to H nuclei) in the range ${\bar{x}}(\rm H_{2}O)\sim
1-5\times 10^{-6}$, the derived value depending upon the adopted grain
properties (Appendix B). This estimate is proportional to the grain absorption
cross-section assumed at the Ly$\alpha$ wavelength (121.6 nm), and represents
the average abundance within the region where Ly$\alpha$ photons scatter prior
to being absorbed by dust or water.

5\. Assuming that the H2O $\nu_{2}$ band emissions observed from HOPS 370 are
radiatively-pumped and that the continuum is scattered light, we obtain a
minimum water abundance in the range $x_{\rm min}=0.7-3\times 10^{-5}$, the
derived value depending upon the adopted grain properties (Appendix A). This
minimum value is proportional to the grain scattering cross-section assumed at
$6\,\mu$m, and would be exceeded if the pumping lines are optically-thick. It
is a factor of several larger than ${\bar{x}}(\rm H_{2}O)$, suggesting that
most Ly$\alpha$ photons are absorbed by dust in a warm and/or ionized zone
very close to where they are created in a fast shock and that only a minority
escape to the region of significant water abundance that emits the water
rovibrational emissions we observe.

6\. Suprathermal OH emissions promise to help elucidate the processes whereby
Lyman $\alpha$ radiation first escapes from fast shocks and then enters nearby
water-rich surroundings where water has been released from grain mantles by
radiative heating or slower non-dissociative shocks, or produced in the gas-
phase by neutral-neutral reactions that are slow at low temperatures but rapid
at the elevated temperatures attained behind shock fronts. Detailed models,
beyond the scope of this Letter, will be needed to understand the transfer of
Ly$\alpha$ radiation and to fully model the rovibrational water emissions
observed from HOPS 370.

We thank the referee, Benoît Tabone, for a very detailed and helpful review
containing multiple suggestions that improved this paper. This work is based
on observations made with the NASA/ESA/CSA James Webb Space Telescope. The
data were obtained from the Mikulski Archive for Space Telescopes at the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127
for JWST. These observations are associated with program #1802. D.A.N. was
supported by grant SOF08-0038 from USRA. P.M. and H.T. acknowledge support of
the Department of Atomic Energy, Government of India, under Project
Identification No. RTI 4002. Support for S.F., A.E.R., S.T.M., R.G., W.F.,
J.G., J.J.T. and D.W. in program #1802 was provided by NASA through a grant
from the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-03127. A.C.G. has been supported by PRIN-MUR 2022 20228JPA3A
“The path to star and planet formation in the JWST era (PATH)” and by INAF-GoG
2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)”.
G.A. and M.O., acknowledge financial support from grants PID2020-114461GB-I00
and CEX2021-001131- S, funded by MCIN/AEI/10.13039/501100011033. Y.- L.Y.
acknowledges support from Grant-in-Aid from the Ministry of Education,
Culture, Sports, Science, and Technology of Japan (20H05845, 20H05844,
22K20389), and a pioneering project in RIKEN (Evolution of Matter in the
Universe). W.R.M.R. is grateful for support from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 101019751 MOLDISK). All the data presented in
this article were obtained from the Mikulski Archive for Space Telescopes
(MAST) at the Space Telescope Science Institute. The specific observations
analyzed can be accessed via https://doi.org/10.17909/3kky-t040 (catalog
10.17909/3kky-t040)

##  Appendix A Dependence on adopted extinction curve and aperture size

We have evaluated the sensitivity of our conclusions to our choice of
extinction law and aperture size. Results are presented in Table 1 for six
cases. We considered three different extinction curves – denoted KP (KPv5), WD
(Weingartner and Draine, 2001, with the modifications described in Draine
2003); and MM (McClure 2009) – and two different aperture radii (0.8$\arcsec$
and 0.4$\arcsec$, denoted by 1 and 2). The standard model, adopted in the main
text, is KP1 (KPv5 extinction curve with the 0.8$\arcsec$ radius aperture).

Different rows in Table 1 show the values obtained for key parameters for all
six cases. The first row below the horizontal line lists the optical depth at
9.7$\,\mu$m, $\tau_{9.7}$, derived from our fit to the H2 lines. It ranges
from 1.53 to 3.12, with the MM extinction law (which is significantly less
dominated by the silicate peak) requiring the largest $\tau_{9.7}$ and the WD
extinction law requiring the smallest, but shows very little variation with
the adopted aperture size. The second row lists, $I_{UV}$, the required UV
photon intensity if every available UV photon led to a water photodissociation
via the ${\tilde{B}}-X$ band. In determining $I_{UV}$, the OH lines were
extinction-corrected using the specified extinction curve and the
corresponding value of $\tau_{9.7}$, and Ly$\alpha$ was assumed to be
responsible for water photodissociation. As discussed in Section 4.1, we
assumed equal fractional systematic errors, $\epsilon$, for each OH flux
measurement and adjusted $\epsilon$ to yield a reduced $\chi^{2}$ of unity
when comparing case KP1 and KP2 fluxes with the predictions for Ly$\alpha$
photodissociation. The third row lists the reduced $\chi^{2}$ obtained for
each of the six cases. The fourth row lists the corresponding values obtained
for a 104 K blackbody radiation field instead of Ly$\alpha$. The fifth and
sixth rows indicate the significance with which each case is disfavored
relative to KP1 or KP2 with photodissociation by Ly$\alpha$. The values
plotted here indicate (1) that the KP extinction law yields a significantly
better fit to the data than either WD or MM; (2) for any extinction curve,
assuming photodissociation by Ly$\alpha$ radiation yields a significantly
better fit to the data than does assuming photodissociation by a $10^{4}$
blackbody.

The seventh row in Table 1 lists the extinction-corrected (ec) Br$\alpha$
intensity, $I(Br\alpha)_{\rm ec}$ for each of the six cases. The next row
lists the corresponding values of $f_{w}$, the fraction of Ly$\alpha$ photons
absorbed by dust. The quantity $f_{w}$ was computed for an assumed
Ly$\alpha$/Br$\alpha$ ratio of 900 (see Appendix B) and is proportional to
$I_{UV}/I(Br\alpha)_{\rm ec}$. Although $I_{UV}$ varies over a factor of 4 – 5
as different extinction curves are adopted, $f_{w}$ varies by less than a
factor 2. This is because there is some degree of cancellation between the
effects on $I_{UV}$ and $I(Br\alpha)_{\rm ec}$. The MM extinction curve, for
example, requires a significantly larger $\tau_{9.7}$ leading to a
significantly larger $I_{UV}$: but because the opacity is less strongly peaked
at silicate feature, $I(Br\alpha)_{\rm ec}$ also increases by a factor 2.

The conclusions of our sensitivity analysis are (1) similar results are
obtained for both aperture sizes; (2) as the origin of the suprathermal OH
emissions, water photodissociation by Ly$\alpha$ radiation is robustly favored
over photodissociation by a $10^{4}$ K blackbody; (3) the KPv5 extinction is
favored over the other two extinction laws considered here; and (4) the
parameter $f_{w}$ lies in the range 2.6 - 4.3 $\%$ for any of the extinction
laws we considered.

Rows 9 - 13 are relevant to the water abundances discussed in Sections 4.3 and
4.4. Here, not only is the relative extinction at different wavelengths of
relevance, but so too is the grain albedo and the ratio of extinction to
column density, $N_{\rm H}$. Both these parameters are available for the KP
and WD extinction laws but not for MM. Row 9 lists the grain absorption cross-
section per H nucleus, $\sigma_{\rm abs}({\rm Ly}\alpha)$, given by each
extinction law. The average water abundance, ${\bar{x}}(\rm H_{2}O)$, needed
to account for $f_{w}$ is listed in Row 10 (see Section 4.3). Row 11 lists the
grain scattering cross-section, $\sigma_{\rm sca}(6.3)$, needed for the
analysis presented in Section 4.4, while Row 12 lists the total $\nu_{2}=1-0$
equivalent width, $W_{\rm H2O}$, for each aperture. The resultant minimum
water abundances, $x_{\rm min}(\rm H_{2}O)$ (see Section 4.3), are given in
Row 13.

The results presented here for ${\bar{x}}(\rm H_{2}O)$ and $x_{\rm min}(\rm
H_{2}O)$ are clearly more sensitive to uncertainties in the dust properties.
Nevertheless, the conclusion that ${\bar{x}}(\rm H_{2}O)<x_{\rm min}(\rm
H_{2}O)$ appears to be robust, with the implications for Ly$\alpha$ radiative
transfer discussed in Section 4.4.

Table 1: Results for different extinction laws and aperture sizes

Row | Quantity | KP1 | WD1 | MM1 | KP2 | WD2 | MM2
---|---|---|---|---|---|---|---
|  | 0.8$\arcsec$ | 0.8$\arcsec$ | 0.8$\arcsec$ | 0.4$\arcsec$ | 0.4$\arcsec$ | 0.4$\arcsec$
1 | $\tau_{9.7}$ | 1.84 | 1.53 | 3.12 | 1.85 | 1.53 | 3.05
2 |  $I_{UV}$ (109 photons cm-2 s-1 sr-1) | 1.67 | 1.10 | 5.11 | 3.75 | 2.42 | 1.07
3 | Reduced $\chi^{2}$ (w.r.t. Ly$\alpha$ model) | 1.00 | 1.63 | 1.55 | 1.00 | 2.65 | 2.33
4 | Reduced $\chi^{2}$ (w.r.t. 104 K BB model) | 2.98 | 6.13 | 4.60 | 3.74 | 8.40 | 7.30
5 | Significancea (for Ly$\alpha$ model) | 0.0$\sigma$ | 4.4$\sigma$ | 3.6$\sigma$ | 0.0$\sigma$ | 6.2$\sigma$ | 5.5$\sigma$
6 | Significancea (for 104 K BB model) | 6.7$\sigma$ | 10.9$\sigma$ | 9.1$\sigma$ | 7.9$\sigma$ | 13.0$\sigma$ | 12.0$\sigma$
7 |  $I(Br\alpha)_{\rm ec}$ (10-3 erg cm-2 s-1 sr-1) | 0.71 | 0.72 | 3.64 | 1.49 | 1.52 | 7.22
8 | $f_{w}$ | 0.043 | 0.028 | 0.025 | 0.046 | 0.029 | 0.027
9 |  $\sigma_{\rm abs}({\rm Ly}\alpha)$ ($10^{-21}\,\rm cm^{2}$) | 1.73 | 0.70 | N/A | 1.73 | 0.70 | N/A
10 | ${\bar{x}}(\rm H_{2}O)$ | $5.1\times 10^{-6}$ | $1.3\times 10^{-6}$ | N/A | $5.1\times 10^{-6}$ | $1.3\times 10^{-6}$ | N/A
11 |  $\sigma_{\rm sca}(6.3)$ ($10^{-23}\,\rm cm^{2}$) | 1.26 | 0.29 | N/A | 1.26 | 0.29 | N/A
12 |  $W_{\rm H2O}$ ($\mu m$) | 0.120 | 0.120 | 0.120 | 0.168 | 0.168 | 0.168
13 | $x_{\rm min}(\rm H_{2}O)$ | $3.0\times 10^{-5}$ | $7.0\times 10^{-6}$ | N/A | $4.2\times 10^{-5}$ | $1.0\times 10^{-5}$ | N/A

aafootnotetext: Significance with which a given model and extinction curve is
disfavored relative to the model with Ly$\alpha$ photodissociation and the KP
extinction curve.

##  Appendix B Shock model predictions for Ly$\alpha$/Br$\alpha$

We have used publicly-available the MAPPINGS V shock model (Sutherland &
Dopita 2017; Sutherland et al. 2018) to estimate the Ly$\alpha$/Br$\alpha$
luminosity ratio within the shocked region where suprathermal OH emissions
were detected. The upper states of these lines may be populated both following
recombination of H+ and by direct colllisional excitation of neutral hydrogen
from the ground state. The Ly$\alpha$/Br$\alpha$ ratio can significantly
exceed the Case B recombination value $\sim 300$, particularly for lower
velocity shocks where collisional excitation is most important, so the use of
shock model predictions is important here.

We ran a grid of models with preshock densities, $n_{0}$, spanning the range
$10^{-1}$ to $10^{5.5}$ H nuclei per cm-3 in steps of 0.1 dex; and with shock
velocities, $v_{s}$, spanning the range 30 to 220 km s-1 in steps of 5 km s-1.
The preshock ionization state was determined self-consistently. The preshock
magnetic field was taken as $0.5\,(n_{0}/\rm cm^{-3})^{1/2}\,\mu G$, and
undepleted solar abundances were adopted. The collision strengths for [Fe II]
fine structure transitions were updated to the values given in the recent
study of Tayal & Zatsarinny (2018).

Figure 7 shows contours of the predicted Ly$\alpha$/Br$\alpha$ luminosity
ratio in the plane of $v_{s}$ and log${}_{10}n_{0}$. Here, the red, cyan and
blue contours show where the Ly$\alpha$/Br$\alpha$ ratio is predicted to be
1000, 1500, and 2000, and black contours show intermediate values spaced by
100. The [Ne III] 15.6 $\mu$m to [Ne II] 12.8 $\mu$m flux ratio is an
excellent tracer of shock velocity. The observed, extinction-corrected value
of 0.026 is obtained for shock parameters lying along the locus marked with
the green solid curve. The green band indicates the region where the predicted
value lies within a factor 1.5 of that observed, with the dotted/dashed
boundaries applying to larger/smaller line ratios. As a probe of the preshock
density, we have considered the [Fe II] 17.9 $\mu$m to [Fe II] 5.3 $\mu$m flux
ratio, which has an observed extinction-corrected value of 4.17. The magenta
curves and magenta band represent analogous results for the [Fe II] line
ratio. The constraint on density is less tight, as indicated by the width of
the magenta band, and must be considered less reliable because recent
independent estimates of the collision strengths show significant differences.
Nevertheless, the intersection of the green and magenta solid lines suggest
that Ly$\alpha$/Br$\alpha$ ratio of 900 is consistent with these diagnostic
line ratios. For any preshock density in the range 102 to 105 cm-3, the flux
[Ne III] 15.6 $\mu$m to [Ne II] 12.8 $\mu$m flux ratio alone suggests a
Ly$\alpha$/Br$\alpha$ ratio in range 700 – 1050.

Contours of the predicted Ly$\alpha$/Br$\alpha$ luminosity ratio in the plane
of $v_{s}$ and log${}_{10}n_{0}$. Red, cyan and blue contours:
Ly$\alpha$/Br$\alpha$ = 1000, 1500, and 2000. Green band: region where the
predicted [Ne III] 15.6 $\mu$m to [Ne II] 12.8 $\mu$m flux ratio lies within a
factor 1.5 of that observed. Magenta band: region where the predicted [Fe II]
17.94 $\mu$m to [Fe II] 5.34 $\mu$m flux ratio lies within a factor 1.5 of
that observed.

## References

  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132. doi:10.1086/156357
  * Brooke et al. (2016) Brooke, J. S. A., Bernath, P. F., Western, C. M., et al. 2016, J. Quant. Spec. Radiat. Transf., 168, 142. doi:10.1016/j.jqsrt.2015.07.021
  * Cacciapuoti et al. (2024) Cacciapuoti, L., Testi, L., Podio, L., et al. 2024, ApJ, 961, 90. doi:10.3847/1538-4357/ad0f17
  * Carr & Najita (2014) Carr, J. S. & Najita, J. R. 2014, ApJ, 788, 66. doi:10.1088/0004-637X/788/1/66
  * Chapman et al. (2009) Chapman, N. L., Mundy, L. G., Lai, S.-P., et al. 2009, ApJ, 690, 496. doi:10.1088/0004-637X/690/1/496
  * Christiaens et al. (2023) Christiaens, V., Gonzalez, C., Farkas, R., et al. 2023, The Journal of Open Source Software, 8, 4774. doi:10.21105/joss.04774
  * Federman et al. (2023) Federman, S., Megeath, S. T., Rubinstein, A. E., et al. 2023, arXiv:2310.03803. doi:10.48550/arXiv.2310.03803
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5. doi:10.3847/0067-0049/224/1/5
  * Gomez Gonzalez et al. (2017) Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7. doi:10.3847/1538-3881/aa73d7
  * González-García et al. (2016) González-García, B., Manoj, P., Watson, D. M., et al. 2016, A&A, 596, A26. doi:10.1051/0004-6361/201527186
  * Harich et al. (2000) Harich, S. A., Hwang, D. W. H., Yang, X., et al. 2000, J. Chem. Phys., 113, 10073. doi:10.1063/1.1322059
  * Heays et al. (2017) Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105. doi:10.1051/0004-6361/201628742
  * Jones et al. (2023) Jones, O. C., Álvarez-Márquez, J., Sloan, G. C., et al. 2023, MNRAS, 523, 2519. doi:10.1093/mnras/stad1609
  * Law et al. (2023) Law, D. R., E. Morrison, J., Argyriou, I., et al. 2023, AJ, 166, 45. doi:10.3847/1538-3881/acdddc
  * Manoj et al. (2013) Manoj, P., Watson, D. M., Neufeld, D. A., et al. 2013, ApJ, 763, 83. doi:10.1088/0004-637X/763/2/83
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192. doi:10.1088/0004-6256/144/6/192
  * Megeath et al. (2021) Megeath, T., Anglada, G., Atnagulov, P., et al. 2021, JWST Proposal. Cycle 1, 1802
  * Narang et al. (2023) Narang, M., Manoj, P., Tyagi, H., et al. 2023, arXiv:2310.14061. doi:10.48550/arXiv.2310.14061 (N23)
  * Neufeld & McKee (1988) Neufeld, D. A. & McKee, C. F. 1988, ApJ, 331, L87. doi:10.1086/185241
  * Neufeld (1990) Neufeld, D. A. 1990, ApJ, 350, 216. doi:10.1086/168375
  * Neufeld et al. (2006) Neufeld, D. A., Melnick, G. J., Sonnentrucker, P., et al. 2006, ApJ, 649, 816. doi:10.1086/506604
  * Neufeld et al. (2009) Neufeld, D. A., Nisini, B., Giannini, T., et al. 2009, ApJ, 706, 170. doi:10.1088/0004-637X/706/1/170
  * Osorio et al. (2017) Osorio, M., Díaz-Rodríguez, A. K., Anglada, G., et al. 2017, ApJ, 840, 36. doi:10.3847/1538-4357/aa6975
  * Pillitteri et al. (2013) Pillitteri, I., Wolk, S. J., Megeath, S. T., et al. 2013, ApJ, 768, 99. doi:10.1088/0004-637X/768/2/99
  * Pontoppidan et al. (2024) Pontoppidan, K. et al. 2024, RNAAS, in preparation
  * Rubinstein et al. (2023) Rubinstein, A. E., Tyagi, H., Nazari, P., et al. 2023, arXiv:2312.07807. doi:10.48550/arXiv.2312.07807
  * Tabone et al. (2021) Tabone, B., van Hemert, M. C., van Dishoeck, E. F., et al. 2021, A&A, 650, A192. doi:10.1051/0004-6361/202039549
  * Sato et al. (2023) Sato, A., Takahashi, S., Ishii, S., et al. 2023, ApJ, 944, 92. doi:10.3847/1538-4357/aca7c9
  * Sutherland & Dopita (2017) Sutherland, R. S. & Dopita, M. A. 2017, ApJS, 229, 34. doi:10.3847/1538-4365/aa6541
  * Sutherland et al. (2018) Sutherland, R., Dopita, M., Binette, L., et al. 2018, Astrophysics Source Code Library. ascl:1807.005
  * Tappe et al. (2008) Tappe, A., Lada, C. J., Black, J. H., et al. 2008, ApJ, 680, L117. doi:10.1086/589998
  * Tayal & Zatsarinny (2018) Tayal, S. S. & Zatsarinny, O. 2018, Phys. Rev. A, 98, 012706. doi:10.1103/PhysRevA.98.012706
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Reynolds, N., et al. 2020, ApJ, 905, 162. doi:10.3847/1538-4357/abc5bf (T20)
  * van Harrevelt & van Hemert (2000) van Harrevelt, R. & van Hemert, M. C. 2000, J. Chem. Phys., 112, 5787. doi:10.1063/1.481154
  * Vuong et al. (2003) Vuong, M. H., Montmerle, T., Grosso, N., et al. 2003, A&A, 408, 581. doi:10.1051/0004-6361:20030942
  * Winston et al. (2010) Winston, E., Megeath, S. T., Wolk, S. J., et al. 2010, AJ, 140, 266. doi:10.1088/0004-6256/140/1/266
  * Zannese et al. (2023) Zannese, M., Tabone, B., Habart, E., et al. 2023, arXiv:2312.14056. doi:10.48550/arXiv.2312.14056
  * Zhou et al. (2015) Zhou, L., Xie, D., & Guo, H. 2015, J. Chem. Phys., 142, 124317. doi:10.1063/1.4915536