# Collisional excitation of interstellar PN by H 2: New interaction potential and scattering calculations

Benjamin Desrousseaux, Ernesto Quintas-Sánchez, Richard Dawes, Sarantos Marinakis, François Lique

► To cite this version:

Benjamin Desrousseaux, Ernesto Quintas-Sánchez, Richard Dawes, Sarantos Marinakis, François Lique. Collisional excitation of interstellar PN by H 2: New interaction potential and scattering calculations. The Journal of Chemical Physics, 2021, 154 (3), pp.034304. 10.1063/5.0039145. hal-03115612

HAL Id: hal-03115612

https://hal.science/hal-03115612v1

Submitted on 18 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Collisional excitation of interstellar PN by H₂: new interaction potential and scattering calculations

Benjamin Desrousseaux,¹ Ernesto Quintas-Sánchez,² Richard Dawes,² Sarantos Marinakis,³,⁴ and François Lique¹, a)

¹) LOMC, UMR 6294, CNRS-Université du Havre, 25 rue Philippe Lebon, BP 1123, 76063 Le Havre Cedex, France

²) Department of Chemistry, Missouri University of Science and Technology, Rolla, Missouri 65409, United States

³) School of Health, Sport and Bioscience, University of East London, Stratford Campus, Water Lane, London E15 4LZ, UK

⁴) Department of Chemistry, University of Patras, Patras, GR-26504, Greece

(Dated: 29 December 2020)

Rotational excitation of interstellar PN molecules induced by collisions with H₂ is investigated. We present the first *ab initio* four-dimensional potential energy surface (PES) for the PN–H₂ van der Waals system. The PES was obtained using an explicitly-correlated coupled cluster approach with single, double, and perturbative triple excitations [CCSD(T)-F12b]. The method of interpolating moving least squares was used to construct an analytical PES from these data. The equilibrium structure of the complex was found to be linear, with the H₂ aligned at the N end of the PN molecule, at an intermolecular separation of 4.2 Å. The corresponding well-depth is 224.3 cm⁻¹. The dissociation energies were found to be 40.19 cm⁻¹ and 75.05 cm⁻¹ for complexes of PN with *ortho*-H₂ and *para*-H₂, respectively. Integral cross-sections for rotational excitation in PN–H₂ collisions were calculated using the new PES, and were found to be strongly dependent on the rotational level of the H₂ molecule. These new collisional data will be crucial to improve the estimation of PN abundance in the interstellar medium from observational spectra.

## I. INTRODUCTION

Phosphorus (P) is known to be an essential element for life as we know it. Indeed, it has a central role in biochemical processes since it is a key element in forming deoxyribonucleic acid (DNA) and cell membranes.¹ P-bearing compounds have also been proposed as key catalysts and chemical buffers for the formation of nucleotides.² However, the way it arrived on Earth is still not fully understood, especially since the abundance of phosphorus in space is low³ (P/H ~ 3 × 10⁻⁷). Hence, there has been an intense interest in searching for P-bearing compounds in astrophysical environments and especially in star-forming regions. So far, only a few simple P-bearing species (PN, PO, CP, HCP, CCP, and PH₃) have been identified, mostly, towards the circumstellar envelopes of evolved stars.⁴

Among the interstellar P-bearing molecules, PN is perhaps the most interesting one. The first interstellar detection of this molecule was reported more than 30 years ago by Turner & Bally⁵ and Ziurys⁶ in Ori(KL), W51M, and SgrB2 sources, by observations of the $j_1 = 2 \to 1$, $j_1 = 3 \to 2$, $j_1 = 5 \to 4$, and $j_1 = 6 \to 5$ rotational transitions ($j_1$ designates the PN rotational level). This was the first example of a phosphorus-containing species seen in interstellar molecular clouds. Since then, the PN molecule has been observed in both star-forming regions⁷⁻⁹ and circumstellar gas.¹⁰,¹¹ Until the recent detection of PO,¹² PN remained the only P-bearing molecule detected in interstellar molecular clouds. With a few exceptions in the last five years (e.g. Rivilla et al.¹³ and references therein), the number of PN observations has been relatively

low compared to other small nitrogen containing species such as CN or NO.

The chemistry of interstellar P-bearing molecules (including PN) has been the object of several studies.¹⁴,¹⁵ However, the gas phase chemistry of phosphorus is still not well understood, partly because some of these molecules are very toxic and difficult to study in a laboratory. In particular, there is still an active debate on whether PN or PO is more abundant in various astrophysical objects.¹³,¹⁶ Studying the amount of PN in astrophysical objects can provide information for other P-bearing species and help in better constraining the phosphorus chemistry in the universe. An accurate determination of PN abundance in the interstellar medium (ISM) from spectral line data requires collisional rate coefficients with the most abundant interstellar species—usually, the H₂ molecule—since collisional processes contribute, in competition with the radiative ones, to the molecular excitation. Without these rates, only approximate estimates of the molecular abundance are possible assuming local thermodynamic equilibrium (LTE), which is generally not a good approximation.¹⁷

Excitation studies of the PN molecule started more than a decade ago. The first study on collisions of PN with He (as a surrogate for H₂) was performed by Tobola *et al.*¹⁸ using quantum scattering calculations. These authors provided collisional rate coefficients between the first 31 rotational energy levels of PN for temperatures up to 300 K. These data were used in astrochemical models in order to derive the PN abundance from the observational spectra.¹⁰,¹¹,¹³,¹⁹ However, it is now well established¹⁷ that rate coefficients for collisions with H₂ are generally different from those for collisions with He, and that strong abundance determination errors can result from the use of uncertain collisional data.²⁰,²¹

More recently, Najar *et al.*²² studied the rotational exci-

<footnote>a) Electronic mail: francois.lique@univ-lehavre.fr</footnote>
tation of PN by collisions with H$_2$, and provided rate coefficients for the rotational (de-)excitation of PN by para-H$_2$. They found significant differences with the PN-He data of Toboła *et al.*,¹⁸ confirming that actual H$_2$ rate coefficients are needed. However, those authors computed and used a potential energy surface (PES) averaged over the H$_2$ rotation, resulting in a simplified dynamical calculation only valid for the collisional excitation by H$_2$ ($j_2 = 0$) ($j_2$ designates the H$_2$ rotational level) at low temperatures. Therefore, the validity of the approximation used by Najar *et al.*²² also needs to be checked. Moreover, as PN is also observed in astrophysical environments where the *ortho*-to-*para*-H$_2$ ratio is expected to be large, it seems also important to provide new rate coefficients for the rotational excitation of PN by H$_2$ ($j_2 > 0$). Such collisional data will provide the astronomical community with the necessary tools to interpret PN emission and to make the most of the observations.

The aim of this study is to obtain inelastic cross-sections for collisions between PN and H$_2$ using quantum scattering calculations. The computation of such data usually takes place within the Born-Oppenheimer approximation for the separation of electronic and nuclear motions. Scattering cross-sections are thus obtained by solving for the motion of the nuclei on an electronic PES, which is independent of the masses and spins of the nuclei.

In this paper, we present the calculation of a new four-dimensional (4D) PES for the ground electronic state of the PN-H$_2$ collisional system. Then, the dissociation energy of the complex with *ortho*-H$_2$ and *para*-H$_2$ is provided. Such data can be useful for future spectroscopical studies. Finally, using this new PES, rotational excitation of PN with both *para*-H$_2$ and *ortho*-H$_2$ is studied. The paper is organized as follows: Section II describes the calculation of the PES. Section III presents the bound states calculations. Section IV contains a concise description of the collisional excitation calculations and of the scattering results. Concluding remarks are presented in Section V.

## II. PN-H₂ POTENTIAL ENERGY SURFACE

The coordinates used to define the 4D PN-H$_2$ PES: R, θ₁, θ₂, and φ, are depicted in Figure 1. $\vec{R}$ is the vector between the centers of mass of the two fragments, $\vec{r}_1$ and $\vec{r}_2$ are vectors aligned with the bond axes of the fragments ($\vec{r}_1$ for PN and $\vec{r}_2$ for H$_2$). Coordinate R is the length of vector $\vec{R}$, while coordinates θ₁ and θ₂ represent the angles between $\vec{R}$ and the vectors $\vec{r}_1$ and $\vec{r}_2$, respectively. The fourth coordinate is the dihedral (out of plane) torsional angle, labeled φ, which is the angle between the vectors $\vec{R} \times \vec{r}_1$ and $\vec{R} \times \vec{r}_2$.

### A. Electronic structure calculations

In this study, both monomers are approximated as rigid, held at their ground state vibrationally-averaged structures. It is a good approximation in this application to consider only the inter-monomer coordinates because their frequencies

<figure><img src="image_2.png" /><figcaption>FIG. 1. Coordinates used to describe the PN-H$_2$ interaction. R: center-of-mass separation (length of $\vec{R}$); θ₁ and θ₂: angles between $\vec{R}$ and the vectors $\vec{r}_1$ and $\vec{r}_2$; and φ: torsional angle.</figcaption></figure>

are all much less than those of the intra-monomer coordinates. The bond distance for H$_2$ was fixed at 0.7666511 Å, the vibrationally-averaged bond distance for *para*-H$_2$, while the $^{31}$P$^{14}$N molecule was held rigid at 1.49350 Å, consistent with the experimental rotational constant of 23495.20419(41) MHz.²³ Masses of 1.00782503223, 14.00307400443, and 30.97376199842 amu were used for $^{1}$H, $^{14}$N, and $^{31}$P respectively.

All *ab initio* calculations were performed using the MOLPRO electronic structure code package.²⁴ The final high-level PES was constructed using *ab initio* data computed with explicitly-correlated coupled-cluster theory.²⁵ The complete basis set limit was estimated by extrapolating calculations at the CCSD(T)-F12b/VTZ-F12 and CCSD(T)-F12b/VQZ-F12 levels, using the $l^{-3}$ formula.²⁶ Exploiting the system's symmetry, energies were only needed in the reduced angular range: 0 < θ₁ < 180°, 0 < θ₂ < 90°, and 0 < φ < 180°. As discussed below, a lower-level guide surface (to avoid computing high level data at inaccessible geometries) was also constructed; this was done using *ab initio* data at the CCSD(T)-F12b/VDZ-F12 level of theory.

### B. Analytical representation

As we have done in the past for other van der Waals (vdW) linear dimers,²⁷⁻³⁶ the 4D PES analytical representation was constructed using an automated interpolating moving least-squares methodology, which has been recently released as a software package under the name AUTOSURF.³⁷ As usual,³⁸ a local fit was expanded about each data point, and the final potential is obtained as the normalized weighted sum of the local fits. The procedure has been described in detail elsewhere.³²,³⁷,³⁹

The shortest intermonomer center-of-mass distance considered is R = 2.4 Å, with the additional restriction of a maximum repulsive energy of 6 kcal/mol (~2100 cm⁻¹) relative to the separated monomers asymptote. The fitted range of the PES extends to R = 22.0 Å. For each of the local fits, a fitting basis of 301 functions is used. As has been pre-
<figure><img src="image_3.png" /><figcaption>FIG. 2. R-optimized contour plot of the PES as a function of the extended angles $\tilde{\theta}_1$ and $\tilde{\theta}_2$. For each pair of angles, the energy (given in cm$^{-1}$) is minimized with respect to the center-of-mass distance R. The position of each stationary point—and the corresponding molecular configuration—is also highlighted. The dashed line represents the disrotatory energy path connecting equivalent minima. See the text for details.</figcaption></figure>

viously described,³⁸ our default fitting basis is a product of radial functions and angular functions (associated Legendre polynomials with a cosine-based torsional function). Note that due to the use of many local expansions, each with their own coefficients, the overall flexibility of the basis is much greater than what would be implied by these parameters if used as a single expansion. To guide the placement of the high-level data—and avoid computing and discarding computationally expensive *ab initio* energies in highly repulsive regions—a lower-level guide PES was constructed using 2176 symmetry-unique points, distributed using a Sobol sequence⁴⁰ biased to sample the short range region more densely. For the high-level PES, the global root-mean-squared fitting error was 0.16 cm$^{-1}$ (0.05 cm$^{-1}$ for energy-regions below the asympote) and the total number of automatically generated symmetry-unique points needed to reach that target was 3759. The analytical representation of the PES is included as Supporting Information.

### C. Characterization of the PES

Figure 2 shows a 2D representation of the PES (denoted R-optimized) for planar configurations, as a function of the extended angles $\tilde{\theta}_1$ and $\tilde{\theta}_2$. The position of stationary points and the corresponding molecular configurations are also high-

lighted in the figure. The extended-angle coordinates have been described elsewhere.²⁸ For planar geometries ($\varphi = 0°$ for quadrants II and IV, and $\varphi = 180°$ for quadrants I and III), the plot describes the complete ranges of $\tilde{\theta}_1$ and $\tilde{\theta}_2$, relaxing the intermonomer distance coordinate R for each pair of angles. This type of plot provides unique insight into the isomers in the system, since for many systems—those (such as this one) without non-planar minima—the plot will include all isomers and the relevant isomerization paths between them, making it easy to visualize planar motions during which $\varphi$ changes from 0 to 180°. There is only one physical well in the PES, which appears twice in the extended angles plot, reflecting the symmetry of the H₂ molecule.

The PES is characterized by a single isomer (labeled

<table><caption>TABLE I. Geometric parameters and potential energy for the global minimum (GM) and transition structures (TS) of the PN-H₂ complex. Units are Angströms, degrees, and cm$^{-1}$.</caption><thead><tr><th></th><th>GM</th><th>TS1</th><th>TS2</th></tr></thead><tbody><tr><td>R</td><td>4.218</td><td>4.265</td><td>4.078</td></tr><tr><td>θ<sub>1</sub></td><td>0.0</td><td>180.0</td><td>58.3</td></tr><tr><td>θ<sub>2</sub></td><td>0.0</td><td>90.0</td><td>58.2</td></tr><tr><td>φ</td><td>—</td><td>—</td><td>0.0</td></tr><tr><td>V</td><td>-224.3</td><td>-75.0</td><td>-23.0</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 3. The energy along various radial cuts defined by different relative orientations of the monomers is plotted. Lines represent the fitted PES, points represent *ab initio* calculations (not used in the fit).</figcaption></figure>

GM) corresponding to a linear structure—with H$_2$ aligned at the N end of the PN molecule—with a well-depth ($D_e$) of 224.3 cm$^{-1}$. The energies and geometric parameters of the global minimum and transition structures are reported in Table I.

In constructing their 2D PES, Najar *et al.*²² employed standard coupled cluster with a triple-zeta basis set, CCSD(T)/AVTZ, but also added mid-bond functions and corrected for basis set superposition error. They report a well-depth of 219.8 cm$^{-1}$ at that level (slightly shallower than our PES), but also list a number of other benchmarks including a well-depth value of 224.0 cm$^{-1}$ for CCSD(T)/CBS (very close to our result of 224.3 cm$^{-1}$ obtained with an explicitly-correlated method). The electronic structure energy values reported by Najar *et al.* are quite similar to ours for a variety of configurations, such that it seems that the main difference between the two studies is that of their reduced dimensionality and associated averaging over orientations.

Various motions of the monomers can move the system from the global minimum to an equivalent structure on the PES. The minimum energy path between the minima, highlighted in Figure 2 (dashed line), represents a disrotatory (or geared) motion that moves the system along a planar low-energy path, passing through the transition structure 1 (TS1). A disrotatory path between minima is common in complexes formed by two linear monomers with planar minima.²⁸,²⁹,³² It is easily identified as it follows an angle of roughly 45 degrees from lower left to upper right in the extended angles plot (cf. Figure 2). Following the path, exiting the plot on the right side, the path would reenter the plot from the left (at the same value of $\theta_2$), thus connecting with the next well up in the plot. Also in Figure 2, a second planar path is seen, moving almost vertically between minima in the plot, passing through a much higher-energy transition structure (TS2). This

movement effectively rotates the H$_2$ molecule 180° while the PN molecule relaxes slightly and then back again to $\theta_1 = 0°$ without flipping. In this system, given only one physical well, unlike numerous other van der Waals systems we have studied, no tunneling splittings are anticipated in the rovibrational levels.

Figure 3 shows 1D cuts of the potential as a function of *R*, upon approach, holding fixed the relative orientations of the monomers. The variation in those cuts gives some indication of the anisotropy of the interactions. It is worth highlighting that all of the *ab initio* energies plotted in Figure 3 were not used in the PES-construction process, and thus serve to confirm the overall quality of our fitted potential. In order to better appreciate the topography of the PES near the global minimum, and the anisotropy of the interaction with respect to the PN and H$_2$ rotations, additional 2D plots of the potential were made and are shown in Figure 4. The plots in Figure 4 explore the region around the minimum by presenting both the ($R, \theta_1$) and ($R, \theta_2$) planes that cross it.

A subroutine of the PES is provided as supplementary material.

### III. PN–H$_2$ DISSOCIATION ENERGY

The PN–H$_2$ dissociation energy was computed with a close coupling method using the BOUND program.⁴¹ The coupled equations were solved using the log-derivative method of Manolopoulos.⁴² The propagator step size was set to 0.05 Å in order to converge the bound states.

The two molecules were considered as rigid-rotors, adopting the data compiled by Herzberg & Huber:⁴³ we used molecular constants $B_{e1} = 0.7864854 \text{ cm}^{-1}$, $\alpha_{e1} = 0.0055364 \text{ cm}^{-1}$, and $D_{e1} = 1.091 \times 10^{-6} \text{ cm}^{-1}$ for the PN
<figure><img src="image_5.png" /><figcaption>FIG. 4. 2D cuts of the PES exploring the region around the global minimum for planar geometries ($\varphi = 0^\circ$), holding the angle $\theta_2$ (upper panel) and $\theta_1$ (lower panel) fixed. Energies are in cm$^{-1}$. Red contours represent positive energies (intervals of 60 cm$^{-1}$), the black contour is the zero, and negative energies are represented by blue dashed contours (intervals of 30 cm$^{-1}$).</figcaption></figure>

molecule and $B_{e2} = 60.853$ cm$^{-1}$, $\alpha_{e2} = 3.062$ cm$^{-1}$, and $D_{e2} = 4.71 \times 10^{-2}$ cm$^{-1}$ for the H₂ molecule. The basis describing the rotation of PN and H₂ molecules included the first 23 and 6 rotational states, respectively.

The expansion of the potential $V(R, \theta_1, \theta_2, \varphi)$ was performed over angular functions for each R distance as described by Green:⁴⁴

$$V(R, \theta_1, \theta_2, \varphi) = \sum_{l_1, l_2, l} v_{l_1, l_2, l}(R) A_{l_1, l_2, l}(\theta_1, \theta_2, \varphi) \quad (1)$$

where $A_{l_1, l_2, l}(\theta_1, \theta_2, \varphi)$ is constructed from coupled spherical functions $\tilde{Y}_{l_i, m_i}(\theta_i, \varphi)$ and the rotational angular momenta of

PN and H₂. In our case, the potential was expanded including $0 \le l_1 \le 18$ for the PN molecule, and $0 \le l_2 \le 6$ for the H₂ molecule.

Our calculated ground state dissociation energies ($D_0$) are equal to 40.19 cm$^{-1}$ and 75.05 cm$^{-1}$ for complexes with *ortho-* and *para*-H₂, respectively. The $D_e$ and $D_0$ values obtained in the present work suggest that the zero-point energy associated with the intermolecular vibrations and internal rotations accounts for ca. 2/3 of the well depth.

To the best of our knowledge, no measurements of the binding energy $D_0$ of PN–H₂ have been reported. Such measurements would be useful for validation of the PES, complementing the spectroscopic level structure. Indeed, scattering calculations depend on the shape and the well-depth of the PES and the dissociation energy is influenced by both. The PES provided here will be useful for any future study of the rovibrational bound states.

## IV. SCATTERING CALCULATIONS

We computed integral cross-sections for the rotational excitation of PN induced by both *para*-H₂($j_2 = 0$) and *ortho*-H₂($j_2 = 1$) using the quantum close-coupling approach implemented in the MOLSCAT non-reactive scattering code.⁴⁵ The hyperfine structure of the PN molecules, due to the non-zero nuclear spin of the nitrogen atom, is neglected in these calculations since the hyperfine structure is not resolved in the astronomical spectra. It can however be considered from the present data using recoupling techniques as described by Lanza & Lique.⁴⁶ The expansion of the potential $V(R, \theta_1, \theta_2, \varphi)$ over angular functions was performed as described in Section III for the calculation of the PN–H₂ dissociation energy. As for dissociation energy calculations, the close-coupled equations of inelastic scattering were solved using the log-derivative propagator of Manolopoulos,⁴² starting from $R_{min} = 2.12$ Å to $R_{max} = 21.7$ Å.

Inelastic cross-sections were obtained between the 11 first rotational levels of PN ($0 \le j_1 \le 10$). The integration parameters were chosen to ensure convergence of these cross-sections for collisional energies ranging from 1 cm$^{-1}$ up to 150 cm$^{-1}$. At each energy, the maximum value of the total angular momentum needed to obtain cross-sections converged to better than $1 \times 10^{-3}$ Å$^2$ was automatically set by the MOLSCAT program.

The size of the rotational basis was studied carefully in order to ensure convergence of the cross-sections. The first 19 rotational states ($0 \le j_1 \le 18$) of PN were included in the basis set. Two sets of calculations were performed: calculations have been done either considering only the ground rotational states of *para-* and *ortho*-H₂ ($j_2 = 0$ and $j_2 = 1$, respectively) or expanding the basis to the two first rotational levels ($j_2 = 0, 2$ and $j_2 = 1, 3$).

In Figure 5, a comparison between those two sets of calculations is presented for 4 different values of the total collisional energy. As can be seen, considering only the first rotational state of the H₂ molecule in the rotational basis leads to inaccuracies higher than a factor of 2 for the cross-section
<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison between cross-sections (in units of Å²) for the collisional (de)excitation of PN by H₂ obtained including only the ground rotational level of H₂ ($j_2 = 0$) or including the two first rotational levels ($j_2 = 0, 2$) in the rotational basis. The two dashed lines delimit the region where the cross-sections differ by less than a factor of 2.</figcaption></figure>

values at all the considered energies. It should be noted that differences higher than a factor of 4 can be observed at low collisional energy (20 cm⁻¹), where the cross-sections exhibit strong shape and Feshbach resonances. It is then essential to add at least one closed H₂ level in the basis set and we also found (not shown here) that including the $j_2 = 4$ level of H₂ in the basis set does not affect the magnitude of the cross-sections. Similar conclusions were also found for the ortho-H₂ basis so that we have retained the first two rotational states of para- and ortho-H₂ in the rotational basis set.

Collision energy variations of the cross-sections for the rotational de-excitation of PN($j_1 = 7$) induced by both para- and ortho-H₂ are illustrated in Figure 6 for selected transitions with $\Delta j_1 = 1$ and $\Delta j_1 = 2$.

For collisions with para-H₂, the cross-sections exhibit strong resonances at low collisional energy. Resonances with lower amplitudes are also observed for collisions with ortho-H₂. This relates to the creation of quasi-bound states within the deep vdW well of the PN-H₂ complex before its dissociation, as already discussed by Najar *et al.* 22

For transitions with $\Delta j_1 = 1$, the values of the cross-sections for collisions with ortho-H₂ are significantly larger than those for collisions with para-H₂, the differences being typically of a factor of 4 at low energy and up to almost an order of magnitude at highest energies. However, this trend is no longer observed for transitions with $\Delta j_1 = 2$, cross-sections for collisions with ortho-H₂ being even lower than those for collisions with para-H₂ at high energy. Actually, we found that the ortho-H₂ cross-sections are much larger than the para-H₂ only for odd $\Delta j_1$ transitions (see Figure 7).

The $j_2$ dependence of the cross-sections is due to the large

<table><caption>TABLE II. Comparison between cross-sections previously obtained by Najar *et al.* 22 and those obtained in this work for a total energy of 100 cm⁻¹.</caption><thead><tr><th rowspan="2">j<sub>1</sub> → j'<sub>1</sub></th><th colspan="2">Najar et al.<sup>22</sup></th><th colspan="2">This work</th></tr><tr><th>CC</th><th>CC, H<sub>2</sub>(j = 0,2)</th><th>CC, H<sub>2</sub>(j = 0)</th></tr></thead><tbody><tr><td>1 → 0</td><td>1.86</td><td>2.74</td><td>1.24</td></tr><tr><td>2 → 0</td><td>6.16</td><td>6.16</td><td>6.32</td></tr><tr><td>2 → 1</td><td>2.69</td><td>3.96</td><td>1.75</td></tr><tr><td>6 → 4</td><td>13.09</td><td>12.83</td><td>13.68</td></tr><tr><td>10 → 8</td><td>33.70</td><td>31.15</td><td>31.55</td></tr><tr><td>10 → 7</td><td>2.29</td><td>3.42</td><td>1.75</td></tr><tr><td>10 → 6</td><td>4.56</td><td>3.35</td><td>5.69</td></tr></tbody></table>

anisotropy of the PES with respect of H₂ rotations and can be explained by the magnitude of the radial coefficients $v_{l_1 l_2 l}(R)$ of Eq. 1. The radial coefficients that mainly contribute to cross-sections with $j_2 \to j_2'$ transitions are those with $l_2$ in the $|j_2 - j_2'| < l_2 < |j_2 + j_2'|$ range. Hence, for collisions with para-H₂ ($j_2 = 0$), only the $l_2 = 0$ term contribute whereas for collisions with ortho-H₂($j_2 = 1$), the $l_2 = 0$ and 2 terms contribute. The radial coefficients with $l_2 = 2$ are not negligible compared to the one with $l_2 = 0$ and explain the difference for collisions with H₂($j_2 = 0$) and H₂($j_2 = 1$). The differences between para- and ortho-H₂ cross-sections are more marked for transitions with odd $\Delta j_1$ since $v_{l_1 l_2 l}(R)$ terms with odd $l_1$ and $l_2 = 2$ are stronger than those with even $l_1$ and $l_2 = 2$.

In Figure 7, we present de-excitation cross-sections out of PN($j_1 = 7$) induced by collisions with para- and ortho-H₂ at a collision energy of 100 cm⁻¹. The propensity rules depend on the quantum state of the H₂ molecule. Collisions with H₂($j_2 = 0$) show a remarkably strong propensity for even $\Delta j_1$ transitions. The explanation of such behavior could be found in the features of the PES. As also found by Najar *et al.*, 22 the averaged PES over the rotational wave function of H₂($j_2 = 0$) (i.e. a spherical average of the PES) exhibits an almost near-homonuclear symmetry, leading to propensity rules in favor of even $\Delta j_1$. For collisions with H₂($j_2 = 1$), a simple propensity in favor of $\Delta j_1 = 1$ is observed and the cross-sections follow an energy-gap law behavior with $\Delta j_1 = 1 > \Delta j_1 = 2 > \Delta j_1 = 3$ etc...

In Table II, we compare, at a total energy of 100 cm⁻¹, the cross-sections previously obtained by Najar *et al.* 22 and those obtained in this work. The cross-sections are reported for a total energy of 100 cm⁻¹, avoiding bias by shifted resonances at lower energy.

When comparing results by Najar *et al.* 22 to our data (CC, H₂($j = 0,2$)), differences higher than 30 % are generally observed with the exception of cross-section values for transitions involving $\Delta j_1 = 2$.

Two main causes can be considered for the differences observed: (i) different levels of electronic structure theory were employed to generate the PESs, (ii) Najar *et al.* 22 averaged their PES over the angular motion of the H₂ molecule, neglecting the anisotropy effects with respect to the H₂ rotation.

As already discussed above, considering a structure-less H₂ molecule for scattering calculations of this system leads to inaccuracies that can be higher than a factor of 2. In Table II, we
<figure><img src="image_7.png" /><figcaption>FIG. 6. Collision energy variation of the cross-sections for the rotational de-excitation of PN induced by collisions with *para*- (left panels) and *ortho*-H$_2$ (right panels) for selected transitions with $\Delta j_1 = 1$ (upper panels) and $\Delta j_1 = 2$ (lower panels).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Rotational de-excitation cross-sections of PN($j_1 = 7$) to PN($j'_1$) induced by collisions with *para*-H$_2$ (circles) and *ortho*-H$_2$ (squares) for a collision energy of 100 cm$^{-1}$.</figcaption></figure>

also report our cross-sections obtained when including only the H$_2$ ($j_2 = 0$) level in the rotational basis set. The comparison of these data with those of Najar *et al.*$^{22}$ will probe the direct impact of the two PESs since the two scattering calculations are then equivalent. Differences between 20 and up to more than 50 % can be seen in the general case. Again, transitions with $\Delta j_1 = 2$ present a better agreement with differences lower than 10 %.

From this comparison, given the similarity found between the employed electronic structure methods (discussed above), one can conclude that the differences between the two sets of data are primarily due to the effects of the H$_2$ rotational structure neglected by Najar *et al.*$^{22}$ Therefore, we recommend the use of our new collisional data in astrophysical applications.

## V. CONCLUSIONS

We have computed the first 4D PES for the PN–H$_2$ vdW system at the CCSD(T)-F12b level of theory, together with an extrapolation to the complete basis set limit. The global
minimum ($\Delta E = -224.3 \text{ cm}^{-1}$) corresponds to a linear geometry with the H₂ molecule approaching the N-atom of the PN molecule. The new PES, in addition to accounting for the full H₂ rotation, is slightly deeper than the one of Najar et al.<fnref n="22" /> The calculated dissociation energy of the *ortho*-H₂ ($D_0 = 75.05 \text{ cm}^{-1}$) and *para*-H₂ ($D_0 = 40.19 \text{ cm}^{-1}$) complexes were also determined. The PES permits calculation of rovibrational bound states suitable for experimental comparisons.

Scattering calculations were performed using this new high-level PES. The results can be summarized as follows:

* The inelastic cross-sections depend significantly on the rotational level of the colliding H₂ molecule. The cross-sections for collisions with *para*-H₂($j_2 = 0$) are globally smaller than the ones for collisions with H₂($j_2 = 1$) even though the differences strongly depend on the $\Delta j_1$ of the rotational transitions.

* For *para*-H₂ collisions, a strong propensity for even $\Delta j_1$ appears whereas, for *ortho*-H₂, cross-sections follow an energy-gap law behavior with $\Delta j_1 = 1 > \Delta j_1 = 2 > \Delta j_1 = 3$ etc...

* We found that obtaining accurate collisional cross-sections requires the inclusion of H₂ closed channels.

These results can be used to improve the modeling of PN abundance in the ISM and hence, improve our knowledge of phosphorus chemistry in space. This may be a key to understanding the origin of life on Earth. Finally, the present calculations will be extended to higher rotational levels and higher collision energies to better cover the astrophysical needs, especially for the interpretation of observations in warm regions of the ISM.

## SUPPLEMENTARY MATERIAL

A Fortran subroutine of the potential energy surface is available as the supplementary material.

## ACKNOWLEDGMENTS

R.D. and E.Q.S. are supported by the U.S. Department of Energy, Office of Basic Energy Sciences (Award DE-SC0019740). F.L. and S.M. acknowledge the support from COST Action CM1401 "Our Astrochemical History". F.L. acknowledges the Institut Universitaire de France, and the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

<footnote n="1">¹M. A. Pasek and D. S. Lauretta, "Aqueous corrosion of phosphide minerals from iron meteorites: A highly reactive source of prebiotic phosphorus on the surface of the early earth," *Astrobiology* **5**, 515–535 (2005).</footnote>

<footnote n="2">²M. W. Powner, B. Gerland, and J. D. Sutherland, "Synthesis of activated pyrimidine ribonucleotides in prebiotically plausible conditions," *Nature* **459**, 239–242 (2009).</footnote>

<footnote n="3">³M. Asplund, N. Grevesse, A. J. Sauval, and P. Scott, "The Chemical Composition of the Sun," *Ann. Rev. Astron. Astrophys.* **47**, 481–522 (2009).</footnote>

<footnote n="4">⁴B. A. McGuire, "2018 Census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disk, and Exoplanetary Molecules," *ApJS* **239**, 17 (2018).</footnote>

<footnote n="5">⁵B. E. Turner and J. Bally, "Detection of Interstellar PN: The First Identified Phosphorus Compound in the Interstellar Medium," *ApJL* **321**, L75 (1987).</footnote>

<footnote n="6">⁶L. M. Ziurys, "Detection of Interstellar PN: The First Phosphorus-bearing Species Observed in Molecular Clouds," *ApJL* **321**, L81 (1987).</footnote>

<footnote n="7">⁷B. E. Turner, T. Tsuji, J. Bally, M. Guelin, and J. Cernicharo, "Phosphorus in the Dense Interstellar Medium," *ApJ*. **365**, 569 (1990).</footnote>

<footnote n="8">⁸P. Schilke, T. D. Groesbeck, G. A. Blake, Phillips, and T. G., "A Line Survey of Orion KL from 325 to 360 GHz," *ApJS* **108**, 301–337 (1997).</footnote>

<footnote n="9">⁹C. Mininni, F. Fontani, V. M. Rivilla, M. T. Beltrán, P. Caselli, and A. Vasyunin, "On the origin of phosphorus nitride in star-forming regions," *MNRAS* **476**, L39–L44 (2018).</footnote>

<footnote n="10">¹⁰S. N. Milam, D. T. Halfen, E. D. Tenenbaum, A. J. Apponi, N. J. Woolf, and L. M. Ziurys, "Constraining Phosphorus Chemistry in Carbon- and Oxygen-Rich Circumstellar Envelopes: Observations of PN, HCP, and CP," *ApJ*. **684**, 618–625 (2008).</footnote>

<footnote n="11">¹¹L. M. Ziurys, D. R. Schmidt, and J. J. Bernal, "New Circumstellar Sources of PO and PN: The Increasing Role of Phosphorus Chemistry in Oxygen-rich Stars," *ApJ*. **856**, 169 (2018).</footnote>

<footnote n="12">¹²B. Lefloch, C. Vastel, S. Viti, I. Jiménez-Serra, C. Codella, L. Podio, C. Ceccarelli, E. Mendoza, J. R. D. Lepine, and R. Bachiller, "Phosphorus-bearing molecules in solar-type star-forming regions: first PO detection," *MNRAS* **462**, 3937–3944 (2016).</footnote>

<footnote n="13">¹³V. M. Rivilla, M. N. Drozdovskaya, K. Altwegg, P. Caselli, M. T. Beltrán, F. Fontani, F. F. S. van der Tak, R. Cesaroni, A. Vasyunin, M. Rumin, F. Lique, S. Marinakis, L. Testi, and the ROSINA team, "ALMA and ROSINA detections of phosphorus-bearing molecules: the interstellar thread between star-forming regions and comets," *MNRAS* **492**, 1180–1198 (2020).</footnote>

<footnote n="14">¹⁴S. B. Charnley and T. J. Millar, "The Chemistry of Phosphorus in Hot Molecular Cores," *MNRAS* **270**, 570 (1994).</footnote>

<footnote n="15">¹⁵I. Jiménez-Serra, S. Viti, D. Quénard, and J. Holdship, "The Chemistry of Phosphorus-bearing Molecules under Energetic Phenomena," *ApJ*. **862**, 128 (2018).</footnote>

<footnote n="16">¹⁶I. Jiménez-Serra, S. Viti, D. Quenard, and J. Holdship, "The chemistry of phosphorus-bearing molecules under energetic phenomena," *ApJ* **862**, 128 (2018).</footnote>

<footnote n="17">¹⁷E. Roueff and F. Lique, "Molecular Excitation in the Interstellar Medium: Recent Advances in Collisional, Radiative, and Chemical Processes," *Chem. Rev.* **113**, 8906–8938 (2013).</footnote>

<footnote n="18">¹⁸R. Tobola, J. Kłos, F. Lique, G. Chalasi'nski, and M. Alexander, "Rotational excitation and de-excitation of PN molecules by He atoms," *A&A* **468**, 1123–1127 (2007).</footnote>

<footnote n="19">¹⁹F. Fontani, V. M. Rivilla, F. F. S. van der Tak, C. Mininni, M. T. Beltrán, and P. Caselli, "Origin of the PN molecule in star-forming regions: the enlarged sample," *MNRAS* **489**, 4530–4542 (2019).</footnote>

<footnote n="20">²⁰M. Lanza, Y. Kalugina, L. Wiesenfeld, A. Faure, and F. Lique, "New insights on the HCl abundance in the interstellar medium," *MNRAS* **443**, 3351–3358 (2014).</footnote>

<footnote n="21">²¹F. Lique, F. Daniel, L. Pagani, and N. Feautrier, "Hyperfine excitation of N₂H⁺ by H₂: towards a revision of N₂H⁺ abundance in cold molecular clouds," *MNRAS* **446**, 1245–1251 (2015).</footnote>

<footnote n="22">²²F. Najar, M. Naouai, H. E. Hanini, and N. Jaidane, "Rotationally inelastic scattering of PN by *para*-H₂($j=0$) at low/moderate temperature," *Monthly Notices of the Royal Astronomical Society* **472**, 2919–2925 (2017).</footnote>

<footnote n="23">²³G. Cazzoli, L. Cludi, and C. Puzzarini, "Microwave spectrum of P¹⁴N and P¹⁵N: Spectroscopic constants and molecular structure," *J. Mol. Spectrosc.* **780–781**, 260–267 (2006).</footnote>

<footnote n="24">²⁴H. Werner, P. Knowles, G. Knizia, F. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, et al., *MOLPRO, version 2019.1, a package of ab initio programs* (University College Cardiff Consultants</footnote>
Ltd.: Cardiff, UK, 2019).

<footnote n="25">²⁵H.-J. Werner, G. Knizia, and F. R. Manby, "Explicitly correlated coupled cluster methods with pair-specific geminals," Mol. Phys. **109**, 407–417 (2011).</footnote>

<footnote n="26">²⁶D. Feller, K. A. Peterson, and T. D. Crawford, "Sources of error in electronic structure calculations on small chemical systems," J. Chem. Phys. **124**, 054107 (2006).</footnote>

<footnote n="27">²⁷J. Brown, X.-G. Wang, R. Dawes, and T. Carrington Jr, "Computational study of the rovibrational spectrum of (OCS)₂," J. Chem. Phys. **136**, 134306 (2012).</footnote>

<footnote n="28">²⁸R. Dawes, X. G. Wang, and T. Carrington, "CO dimer: New potential energy surface and rovibrational calculations," J. Phys. Chem. A **117**, 7612–7630 (2013).</footnote>

<footnote n="29">²⁹X.-G. Wang, T. Carrington Jr, and R. Dawes, "Computational study of the rovibrational spectrum of (CO₂)₂," J. Mol. Spectrosc. **330**, 179–187 (2016).</footnote>

<footnote n="30">³⁰E. Castro-Juárez, X.-G. Wang, T. Carrington Jr, E. Quintas-Sánchez, and R. Dawes, "Computational study of the ro-vibrational spectrum of CO–CO₂," J. Chem. Phys. **151**, 084307 (2019).</footnote>

<footnote n="31">³¹J. Brown, X.-G. Wang, T. Carrington Jr, G. S. Grubbs, and R. Dawes, "Computational study of the rovibrational spectrum of CO₂–CS₂," J. Chem. Phys. **140**, 114303 (2014).</footnote>

<footnote n="32">³²R. Dawes, X.-G. Wang, A. W. Jasper, and T. Carrington Jr, "Nitrous oxide dimer: A new potential energy surface and rovibrational spectrum of the nonpolar isomer," J. Chem. Phys. **133**, 134304 (2010).</footnote>

<footnote n="33">³³G. Donoghue, X.-G. Wang, R. Dawes, and T. Carrington Jr, "Computational study of the rovibrational spectra of CO₂–C₂H₂ and CO₂–C₂D₂," J. Mol. Spectrosc. **330**, 170–178 (2016).</footnote>

<footnote n="34">³⁴A. Barclay, A. McKellar, N. Moazzen-Ahmadi, R. Dawes, X.-G. Wang, and T. Carrington, "Infrared spectrum and intermolecular potential energy surface of the CO–O₂ dimer," Phys. Chem. Chem. Phys. **20**, 14431–14440 (2018).</footnote>

<footnote n="35">³⁵B. Desrousseaux, E. Quintas-Sánchez, R. Dawes, and F. Lique, "Collisional Excitation of CF⁺ by H₂: Potential Energy Surface and Rotational Cross Sections," J. Phys. Chem. A **123**, 9637–9643 (2019).</footnote>

<footnote n="36">³⁶C. T. Bop, F. A. Batista-Romero, A. Faure, E. Quintas-Sánchez, R. Dawes, and F. Lique, "Isomerism effects in the collisional excitation of cyanoacetylene by molecular hydrogen," ACS Earth and Space Chemistry **3**, 1151–1157 (2019).</footnote>

<footnote n="37">³⁷E. Quintas-Sánchez and R. Dawes, "Autosurf: A freely available program to construct potential energy surfaces," J. Chem. Inf. Model. **59**, 262–271 (2019).</footnote>

<footnote n="38">³⁸R. Dawes and E. Quintas-Sánchez, "The Construction of Ab Initio?Based Potential Energy Surfaces," in *Reviews in Computational Chemistry* vol. 31 (Wiley, 2018) Chap. 5, pp. 199–264.</footnote>

<footnote n="39">³⁹M. Majumder, S. A. Ndengue, and R. Dawes, "Automated construction of potential energy surfaces," Mol. Phys. **114**, 1–18 (2016).</footnote>

<footnote n="40">⁴⁰I. M. Sobol, "Uniformly distributed sequences with an additional uniform property," USSR Comput. Math. Math. Phys. **16**, 236–242 (1976).</footnote>

<footnote n="41">⁴¹(1993), J. M. Hutson, BOUND computer code, version 5 (1993), distributed by Collaborative Computational Project No. 6 of the Science and Engineering Research Council (UK).</footnote>

<footnote n="42">⁴²D. E. Manolopoulos, J. Chem. Phys. **85**, 6425 (1986).</footnote>

<footnote n="43">⁴³G. Herzberg and K. P. Huber, *Molecular Spectra and Molecular Structure* (Van Nostrand Reinhold, New York London, 1979).</footnote>

<footnote n="44">⁴⁴S. Green, "Rotational excitation in H₂-H₂ collisions: Close-coupling calculations," J. Chem. Phys. **62**, 2271–2277 (1975).</footnote>

<footnote n="45">⁴⁵Hutson, J. M. and Green, S., "Molscat computer code, version 14 (1994), distributed by collaborative computational project no. 6 of the engineering and physical sciences research council (uk)," (1994).</footnote>

<footnote n="46">⁴⁶M. Lanza and F. Lique, "Hyperfine excitation of linear molecules by para- and ortho-H₂: Application to the HCl-H₂ system," J. Chem. Phys. **141**, 164321 (2014).</footnote>