Revisiting the Formation of HeH⁺ in the Planetary Nebula NGC 7027

Bound molecular complexes involving helium are extremely rare, but early laboratory studies of electron impact ionization of a H₂/He gas mixture suggested the presence of the helium hydride cation HeH⁺ (Hogness & Lunn 1925). Immediately after the spectroscopic identification of HeH⁺ in the laboratory (Dabrowski & Herzberg 1978), astronomers postulated its presence in the interstellar medium and developed chemical schemes for its formation (Black 1978; Flower & Roueff 1979; Roberge & Dalgarno 1982). A prime environment for the presence of HeH⁺ was predicted to be near the Strömgren radius of a planetary nebula (PN) with its formation dominated by the radiative association process,

$$\begin{eqnarray}&&{\mathrm{He}}^{+}+{\rm{H}}\to {\mathrm{HeH}}^{+}+{\rm{h}}\nu .\end{eqnarray} \tag{ 1 }$$

At about the same time, HeH⁺ was proposed to play a role in the formation of ${{\rm{H}}}_{2}^{+}$ and H₂ during the epoch of recombination in the early universe through the following sequence of reactions (Black 1978; Lepp & Shull 1984):

$$\begin{eqnarray}&&{{\rm{H}}}^{+}+\mathrm{He}\to {\mathrm{HeH}}^{+}+{\rm{h}}\nu ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\mathrm{HeH}}^{+}+{\rm{H}}\to {{\rm{H}}}_{2}^{+}+\mathrm{He},\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{{\rm{H}}}_{2}^{+}+{\rm{H}}\to {{\rm{H}}}_{2}+{{\rm{H}}}^{+}\end{eqnarray} \tag{ 4 }$$

(for detailed models see Galli & Palla 1998; Stancil et al. 1998). Over the next four decades, astronomers diligently searched for astrophysical signatures of HeH⁺ in PNs (Moorhead et al. 1988; Liu et al. 1997; Dinerstein & Geballe 2001), supernova ejecta (Miller et al. 1992), and high-redshift quasars (Zinchenko et al. 2011). While hints of its presence were tantalizing, only upper limits to its abundance could be obtained, leaving the existence of HeH⁺ in astrophysical sources in question. Finally, Güsten et al. (2019) made an unambiguous detection of the $J=1\to 0$ rotational transition line of HeH⁺ in the PN NGC 7027. This was enabled by the development of a high spectral-resolution detector, termed the German Receiver for Astronomy at Terahertz Frequencies (GREAT), and its subsequent deployment on the Stratospheric Observatory for Infrared Astronomy. In particular, this detection was made possible by the ability of GREAT to separate the HeH⁺ emission line from the hyperfine components of a nearby blended CH feature, which had plagued past efforts. This detection was a important advance with a number of significant implications for astrochemistry. However, modeling the intensity of the HeH⁺ feature proved to be problematic. A comprehensive nonequilibrium chemical and level population model predicted the HeH⁺ line intensity to be a factor of four smaller than the observation. Güsten et al. (2019) argued that this discrepancy could be remedied by scaling the literature values of the rate coefficient for process (1) by a constant value.

More recently, Neufeld et al. (2020) observed the 1−0 P(1) and P(2) rovibrational lines of HeH⁺ in NGC 7027. These emission features appear in the near-infrared and were obtained with the iSHELL spectrograph on NASA's Infrared Telescope Facility, confirming the discovery of Güsten et al. (2019).

In this work, we demonstrate that the prior rate coefficients computed by other authors are robust and that an arbitrary scaling is unjustified. Further, using more modern theoretical approaches and molecular data, we provide updated rate coefficients for the radiative association process (1) with an estimated uncertainty. We conclude by discussing possible resolutions to the observed and predicted HeH⁺ line intensity, though a detailed model of NGC 7027 is beyond the scope of this work.

Following the work on HeH⁺ photodissociation by Miyake et al. (2011), the ab initio potential energies for the ground X ${}^{1}{{\rm{\Sigma }}}^{+}$ and excited A ${}^{1}{{\rm{\Sigma }}}^{+}$ electronic states are taken from Bishop & Cheung (1979) and Kraemer et al. (1995), respectively. Extensions of the ab initio potentials to both short- and long-range were performed as described in Miyake et al. (2011). The A ${}^{1}{{\rm{\Sigma }}}^{+}\,\leftarrow$ X ${}^{1}{{\rm{\Sigma }}}^{+}$ transition dipole moment (TDM) function was adopted from Kraemer et al. (1995) with short- and long-range extrapolations as given in Miyake et al. (2011). The accuracies of the potential energies are discussed in the respective manuscripts, while the photodissociation cross section computed by Miyake et al. (2011) was found to be consistent with the only available experimental value, suggesting that the TDM is reliable. For the ground electronic state, we computed 162 bound rovibrational states consistent with the calculations of Zygelman et al. (1998).

In this work, we utilized two different theoretical approaches to compute the cross sections and rate coefficients for the radiative association process (1). In the first method, we used the standard two-state quantum-mechanical Fermi golden rule (FGR) approach which should be valid in the zero-density limit (ZDL), henceforth referred to as FGR-ZDL. The FGR-ZDL cross section is given by

$$\begin{eqnarray}&&\sigma (E)=\displaystyle \sum _{{J}^{{\prime} }}\displaystyle \sum _{v^{\prime\prime} }{\sigma }_{{J}^{{\prime} }}(v^{\prime\prime} ,E),\end{eqnarray} \tag{ 5 }$$

where the partial cross section is given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{J}^{{\prime} }}(v^{\prime\prime} ,E) & = & \displaystyle \frac{64}{3}\displaystyle \frac{{\pi }^{5}}{{c}^{3}}\displaystyle \frac{{\nu }^{3}}{{k}^{2}}p\left[{J}^{{\prime} }{M}_{v^{\prime\prime} ,{J}^{{\prime} }-1;k,{J}^{{\prime} }}^{2}\right.\\ & & +\,\left.({J}^{{\prime} }+1){M}_{v^{\prime\prime} ,{J}^{{\prime} }+1;k,{J}^{{\prime} }}^{2}\right],\end{array}\end{eqnarray} \tag{ 6 }$$

and ν is the photon frequency, k the initial wavevector, and p the approach probability factor in the initial electronic state, which is 1/4 in this case.⁴ J' is the initial rotational quantum number, while v'', J'' is the vibrational and rotational quantum numbers for the final product HeH⁺, and M is the electric TDM element between the initial and final states.

In the second approach, a kinetic model was adopted which accounts for both direct and indirect radiative association processes (Forrey 2013, 2015). This method is applied to two limiting cases in which the gas is taken to be in local thermodynamic equilibrium (LTE) and in non-LTE (NLTE) in the ZDL; these cases are expected to be valid at high and low densities, respectively, with typical astrophysical conditions lying between the two limits. Radiative association results in the two limits will be referred to as LTE and NLTE-ZDL, respectively. It has been shown in our earlier work that LTE and NLTE-ZDL rate coefficients merge at relatively high temperatures, while the NLTE-ZDL and FGR-ZDL results are nearly identical for all temperatures (see, for example, Cairnie et al. 2017).

Using the LTE, NLTE-ZDL, and FGR-ZDL methods, we computed total cross sections and rate coefficients for process (1). Figure 1 displays the resulting rate coefficients as a function of kinetic temperature.⁵ The current LTE and NLTE-ZDL results are seen to converge for temperatures above ∼300 K. Our FGR-ZDL results are not shown as they are indistinguishable from the NLTE-ZDL rate coefficients on the scale of the figure. It should also be noted that two independent FGR-ZDL codes were utilized in the calculations and gave consistent results.

Zoom In Zoom Out Reset image size

The current rate coefficients are compared to prior calculations using the FGR-ZDL method (Zygelman & Dalgarno 1990; Kraemer et al. 1995), a quantum decay method (Zygelman et al. 1989), and approximate results (Sando et al. 1971) as reported in Roberge & Dalgarno (1982). The results of Kraemer et al. (1995), which adopted two different sets of ab initio potentials, are in excellent agreement with the current rate coefficients for T < 2000 K. The results of Roberge & Dalgarno (1982), Zygelman et al. (1989), and Zygelman & Dalgarno (1990) are slightly larger than the current FGR-ZDL and NLTE-ZDL calculations, but follow the same temperature dependence. This difference is likely related to the adoption of different TDM functions and potential energies, while the difference at high temperatures with the Kraemer et al. (1995) results may be due to convergence errors in the latter (e.g., the maximum kinetic energy in their cross sections occurs near 10,000 cm⁻¹, which would result in an underestimation of the rate coefficient for ∼10,000 K). Given the fact that the current results agree above 300 K and that they were obtained with three independent theoretical approaches and four independent calculations, we argue that the NLTE-ZDL and FGR-ZDL results are robust with an uncertainty of less than 5%. The LTE radiative association rate coefficients are also robust, and serve as upper limits corresponding to high density conditions. Numerical rate coefficients are presented in Table 1.

Recently, Vranckx et al. (2013) studied the reverse process of photodissociation,

$$\begin{eqnarray}&&{\mathrm{HeH}}^{+}+{\rm{h}}\nu \to {\mathrm{He}}^{+}+{\rm{H}},\end{eqnarray} \tag{ 7 }$$

using a time-dependent approach. Using these results, they estimated the radiative association cross section (1). However, their photodissociation calculations neglected the rotational dependence of the cross sections (i.e., only J'' = 0 was considered). As a consequence, they multiplied each v'' cross section by the total number of rotational states for that vibrational state. Güsten et al. (2019) presented a temperature-independent rate coefficient, based on the above cross section, which we show in Figure 1. This estimate is about a factor of 1.6 smaller than the current detailed calculations, which conversely sum over transitions to all v'', J'' levels as indicated by Equation (5). Hence, the approximation proposed in Vranckx et al. (2013) is not very reliable and, in fact, unnecessary.

Güsten et al. (2019) performed a chemical model of NGC 7027 near the Strömgren radius where the HeH⁺ abundance is predicted to peak. Their model considers reactions (1) and (3), as well as the dissociation recombination process

$$\begin{eqnarray}&&{\mathrm{HeH}}^{+}+e\to \mathrm{He}+{\rm{H}}.\end{eqnarray} \tag{ 8 }$$

They adopted updated rates for the latter two processes, but found that use of the radiative association rate coefficient from Vranckx et al. (2013) resulted in a predicted abundance of a factor of ∼4 too small. It should be noted that they also performed a NLTE spectral model of the HeH⁺ rotational emission assuming excitation by electron collisions. Their NLTE spectra model underpredicted the HeH⁺ $J=1\to 0$ rotational line intensity. To obtain agreement with observation, Güsten et al. (2019) scaled the rate coefficient deduced from Vranckx et al. (2013) by a factor of 4.3 as shown in Figure 1. They justified this arbitrary scaling by arguing that the radiative association rate coefficients in the literature vary widely. If one removes the approximate value deduced from the cross sections of Vranckx et al. (2013) and the highest temperature points from Kraemer et al. (1995), which we believe are too low due to convergence issues, then there is a rather tight agreement between all previous calculations over the 1000–10,000 K range. The radiative association rate coefficient is robustly determined from prior data to be (2.0 ± 0.5) × 10⁻¹⁶ cm³ s⁻¹ in the relevant temperature range. The current set of calculations significantly reduce this uncertainty to (2.14 ± 0.1) × 10⁻¹⁶ cm³ s⁻¹ at 5000 K (near the local minimum). We therefore argue that this arbitrary scaling of 4.3 (or now 2.8 if the current result is adopted) is unjustified, being at odds with statements given in Güsten et al. (2019). In fact, the uncertainties in the other two key reactions are probably larger. Further, scaling of rate coefficients to match observed intensities more likely hides neglected physics in the model.

For the case of dissociative recombination, a new storage ring measurement (Novotný et al. 2019) found a slight increase in the DR rate coefficient by a factor of 1.4 ± 0.3 at 10,000 K. Unfortunately, this would have the affect of reducing the modeled HeH⁺ abundance.

Likewise, the NGC 7027 model of Güsten et al. (2019) adopted the reactive scattering rate coefficients of Bovino et al. (2011), but comparison of their cross sections to earlier measurements suggest the scattering calculations are likely too large by a factor of ∼1.6 near 10⁴ K. Such calculations involving a system of three atoms/ions are far more difficult to perform than two-body radiative association. Further, the time-dependent approach used by Bovino et al. (2011) cannot distinguish the product channels as it incorporates a flux-absorbing imaginary potential. Instead, their result is a sum over reactive and inelastic channels, a so-called depletion rate. Considering these reaction uncertainties suggests that the HeH⁺ modeled abundance may be increased. On balance, a discrepancy factor between ∼3 and 4 remains. Similar conclusions were also reach by Neufeld et al. (2020), who updated the model in Güsten et al. (2019).

One mechanism not considered in the current model is collisional excitation of HeH⁺ due to H collisions; this process which might be relevant as the H/e⁻ fractional abundance ratio is ∼10 just outside the Strömgren radius, where the abundance of HeH⁺ peaks. Unfortunately, we are unaware of any rotational inelastic collisional rate coefficients of HeH⁺ due to heavy particles. If one assumes an upper limit due to the HeH⁺ + H depletion rate of 1.2 × 10⁻⁹ cm³ s⁻¹ compared to the electron rate coefficient of 2.7 × 10⁻⁷ cm³ s⁻¹ at 10,000 K, the contribution to the line intensity due to H collisions is about ∼5% of that due to electrons. The omission of this collisional process is therefore unlikely to reduce the aforementioned discrepancy.

The consideration of a rotationally and/or vibrationally resolved chemistry may resolve some of the discrepancy. However, the sparsity of quantum-state-resolved chemical rates has generally limited such studies. Agúndez et al. (2010) considered the affect of reactions involving vibrationally excited H₂ in a number of environments, including NGC 7027, while Walker et al. (2018) studied a full rovibrationally resolved chemistry for H₂ and ${{\rm{H}}}_{2}^{+}$ in the early universe. The former do find a significant enhancement in the abundance of CH⁺ in their NGC 7027 model, but it is formed by the endothermic process

$$\begin{eqnarray}&&{{\rm{H}}}_{2}+{{\rm{C}}}^{+}\to {\mathrm{CH}}^{+}+{\rm{H}},\end{eqnarray} \tag{ 9 }$$

while reaction (3) is extremely exothermic. While we are unaware of rovibrationally resolved reaction data for process (3), the rate coefficients are unlikely to be very sensitive to internal excitation for kinetic temperatures of 5000–10,000 K. In fact, the recent DR measurements by Novotný et al. (2019) find the rate coefficients for reaction (8) to be nearly independent of rotational excitation for T > 800 K.

Finally, infrared (IR) pumping of excited HeH⁺ rovibrational levels due to thermal dust emission may be important. NGC 7027 has a prominent dust peak near 30 μm. IR radiative pumping may enhance the J = 1 population and should be considered in the excitation model. For the case of ${{\rm{H}}}_{2}^{+}$ in the recombination era, Walker et al. (2018) find that its level populations are driven to a thermal distribution characterized by the temperature of the cosmic background radiation field.

In summary, the first detection of HeH⁺ in an astronomical source by Güsten et al. (2019) is an important advance for molecular astrophysics. It will revitalize interest in helium chemistry. However, we show in this study that of all the reactions relevant to HeH⁺ chemistry in NGC 7027, the radiative association of He⁺ and H is probably the most accurately known, while the recent measurements of Novotný et al. (2019) have put our knowledge of HeH⁺ DR on a much sounder footing. Through accurate calculations performed in this work, we have narrowed the uncertainty in the radiative association rate coefficient further, and argue that arbitrary scaling of its rate coefficient to bring observations and models into agreement is not justified, and in fact, may hide important neglected physics. While detailed astrophysical modeling of the HeH⁺ emission from NGC 7027 is beyond the scope of this work, we propose possible improvements to the chemistry, many of which we are explicitly exploring in a robust investigation of helium reactions in the early universe (E. Courtney et al. 2020, in preparation).

This work was partially supported by NSF grant No. PHY-1806180 and NASA grant NNX15AI61G. ITAMP is supported by NSF grant No. PHYS-1521560 to Harvard University and the Smithsonian Astrophysical Observatory. We thank David Neufeld for helpful comments.