The P(⁴S) + NH(³Σ^–) and N(⁴S) + PH(³Σ^–)reactions as sources of interstellar phosphorus nitride

Abstract

Phosphorus nitride (PN) is believed to be one of the major reservoirs of phosphorus in the interstellar medium (ISM). For this reason, understanding which reactions produce PN in space and predicting their rate coefficients is important for modelling the relative abundances of P-bearing species and clarifying the role of phosphorus in astrochemistry. In this work, we explore the potential energy surfaces of the $\textrm{P}(^4\textrm{S}) + \textrm{NH}(^3\Sigma^-)$ and $\textrm{N}(^4\textrm{S}) + \textrm{PH}(^3\Sigma^-)$ reactions and the formation of $\textrm{H}(^2\textrm{S}) + \textrm{PN}(^1\Sigma^+)$ through high accuracy ab initio calculations and the variable reaction coordinate transition state theory (VRC-TST). We found that both reactions proceed without an activation barrier and with similar rate coefficients that can be described by a modified Arrhenius equation ($k(T)=\alpha\!\left( T/300 \right)^{\beta} \exp\!{(\!-\!\gamma/T)})$ with $\alpha=0.93\times 10^{-10}\rm cm^3\,s^{-1}$, $\beta=-0.18$ and $\gamma=0.24\, \rm K$ for the $\textrm{P} + \textrm{NH} \longrightarrow \textrm{H} + \textrm{PN}$ reaction and $\alpha=0.88\times 10^{-10}\rm cm^3\,s^{-1}$, $\beta=-0.18$ and $\gamma=1.01\, \rm K$ for the $\textrm{N} + \textrm{PH} \longrightarrow \textrm{H} + \textrm{PN}$ one. Both reactions are expected to be relevant for modelling PN abundances even in the cold environments of the ISM. Given the abundance of hydrogen in space, we have also predicted rate coefficients for the destruction of PN via H + PN collisions.

1. Introduction

Phosphorus is a crucial and vital element for life in our planet. It is present in essential biological molecules, such as nucleic acids, phospholipids, phosphorylated proteins, nicotinamide adenine dinucleotide (NAD) (Rivilla et al., 2020; Sousa-Silva et al., 2020; Maciá, 2005) and in metabolic processes of living organisms (Maciá, 2005; Goldford et al., 2017). In spite of its crucial role for life as we know it, the chemistry of phosphorus-containing (P-bearing) molecules in the interstellar medium (ISM) is still not very well understood (Jiménez-Serra et al., 2018; Chantzos et al., 2020). The question of how it became present in Earth’s prebiotic chemistry and biologically available is also an open problem (Goldford et al., 2017; Rivilla et al., 2020).

Phosphorus nitride (PN) was the first P-bearing molecule detected in space, which occurred in the late eighties in several molecular clouds (Turner & Bally, 1987; Ziurys, 1987). It was later also observed in dense clouds and in their massive cores (Ziurys, Schmidt and Bernal, 2018). This molecule is one of the most important reservoirs of phosphorus in the gas phase in the ISM (Thorne et al., 1984; Tenenbaum, Woolf and Ziurys, 2007; Lefloch et al., 2016; Yamaguchi et al., 2011) and, therefore, the knowledge of how it is formed and destroyed is very important. However, the lack of observational constraints does not allow for an understanding on even how atomic P gets transformed into molecules (Fontani et al., 2019).

Another important P-bearing molecule is phosphorus monoxide (PO), which represents the first identification of the P-O bond in the ISM (Lefloch et al., 2016; Tenenbaum, Woolf and Ziurys, 2007). Both PN and PO have been identified in massive star-forming regions (Rivilla et al., 2016) and constitute two important molecules in the aim of better understanding the astrochemistry of phosphorus. Apart from high-mass star-forming regions and circumstellar envelopes, these molecules were also observed in giant molecular cloud regions, where energetic phenomena such as UV radiation fields and cosmic rays are present and may influence in their respective abundances (Jiménez-Serra et al., 2018). The PO/PN abundance ratio can change depending on the observed environment. For instance, in the presence of cosmic rays the ratio is larger than 1, while in their absence, it is lower (Jiménez-Serra et al., 2018; Rivilla et al., 2020; Concepción et al., 2021).

The most commonly used neutral-neutral reactions in astrochemical databases that explain the formation of PN were proposed by (Millar, Bennett and Herbst, 1987):

(1) \begin{align} {}\textrm{N}(^4\textrm{S}) + \textrm{PO}(^2\Pi) \longrightarrow \textrm{O}(^3\textrm{P}) + \textrm{PN}(^1\Sigma^+) \end{align}

(2) \begin{align} \textrm{N}(^4\textrm{S}) + \textrm{PO}(^2\Pi) \longrightarrow \textrm{P}(^4\textrm{S}) + \textrm{NO}(^2\Pi) \end{align}

(3) \begin{align} {}\textrm{N}(^4\textrm{S}) + \textrm{PH}(^3\Sigma^-) \longrightarrow \textrm{H}(^2\textrm{S}) + \textrm{PN}(^1\Sigma^+).\end{align}

However, these reactions have not been experimentally probed, and the rate coefficients employed in the models are estimations based on similar reactions. Recently, the rate coefficients employed in the models for reactions (1) and (2) received computational corroboration (Souza, Silva and Galvão, 2021). However, the rate coefficients employed in astrochemical databases such as UMIST (McElroy et al., 2013) and KIDA (Wakelam et al., 2012) for reaction (3) was estimated based on the chemical similarity with the reaction N+NH $\rightarrow$ N $_2$ +H (Smith, Herbst and Chang, 2004), which in turn had its rate constant already studied theoretically and experimentally (Caridade et al., 2005; Hack, Wagner and Zasypkin, 1994; Miller et al., 1985).

It is clear that a better understanding of several other reactions involving PO and PN with abundant and reactive species available in space (such as atomic carbon, nitrogen, and oxygen) is essential to obtain a better description of the abundances of these molecules in different regions of the ISM. In this spirit, a very recent work by (Douglas, Gobrecht and Plane, 2022) estimated the rate coefficients for several reactions involving PO and PN using collision capture rates methodology and with long-range transition state theory (Georgievskii & Klippenstein, 2005) (LR-TST), which is an important contribution to this endeavour.

In this work we aim to enhance our understanding of the reaction of P-bearing molecules by performing accurate calculations on the energies and rate coefficients for reaction (3). We also include in the effort the similar reaction:

(4) \begin{align} {}\textrm{P}(^4\textrm{S}) + \textrm{NH}(^3\Sigma^-) \rightarrow \textrm{H}(^2\textrm{S}) + \textrm{PN}(^1\Sigma^+),\end{align}

which is also potentially relevant in the transformation of atomic phosphorus into molecular species and is not available in the astrochemical databases (Wakelam et al., 2012; McElroy et al., 2013)

2. Methodology

2.1. Structures optimisation and potential energy surfaces calculations

All electronic structure calculations reported here were performed using the GAMESS-US (Schmidt et al., 1993) and MOLPRO (Werner et al., 2012) packages. For preliminary exploratory work, we chose the density functional theory (DFT) (Kohn & Sham, 1965) using the aug-cc-pV(T+d)Z (Dunning, 1989; Kendall, Dunning and Harrison, 1992) basis set and the exchange and correlation functional M06-2X (Zhao & Truhlar, 2008). For simplicity of notation, the aug-cc-pV(X+d)Z basis set and will be hereafter denoted as AVXZ+d. Vibrational analysis was carried out to confirm the minima and transition states (TSs) found.

With the optimised DFT structures in hand, all results were fully reoptimised and refined using two highly accurate ab initio levels. The first uses geometry optimisation and frequencies at the full valence complete active space self-consistent field (Szalay et al., 2012) (CASSCF) level for all stationary structures (11 correlated electrons in 9 active orbitals), followed by a single point energy calculations with the explicitly correlated multireference configuration interaction (Szalay et al., 2012; Shiozaki, Knizia and Werner, 2011) (MRCI-F12) method, including the Davidson correction and also using the AVTZ+d basis set. The second high accuracy method, for benchmark purposes, is the reoptimisation and frequencies calculations at the coupled cluster singles and doubles and perturbative triples (CCSD(T)) level (Bartlett, 1989; Bartlett et al., 1990; Raghavachari et al., 1989) using the AVTZ+d basis set, followed by a single point energy refinement using the explicitly correlated coupled cluster method (CCSD(T)-F12) (Adler, Knizia and Werner, 2007; Knizia, Adlerand Werner, 2009) with the cc-pVQZ-F12 basis set (hereafter, VQZ-F12).

The MRCI-F12/AVTZ+d//CAS/AVTZ+d and CCSD(T)-F12/VQZ-F12//CCSD(T)/AVTZ+d calculations were also used to obtain the magnitude of the dominant configuration in the CASSCF wave functions (C ${_0^2}$ ) and the T $_1$ diagnostic of all stationary points. This was made to assess the multireference character (MR) of each structure, as it is known that a T $_1$ diagnostic higher than 0.02 (Lee & Taylor, 1989; Leininger et al., 2000; Lee, 2003); and a C ${_0^2}$ smaller than 0.90 are indicatives of a significant MR character. The Wxmacmolplt and Avogadro software were used for graphic visualisation and representation of the molecular geometries (Bode & Gordon, 1998; Hanwell et al., 2012).

2.2. Rate coefficients

Capture rate coefficients were calculated using direct variable reaction coordinate transition state theory (Klippenstein, 1991, 1992) (VRC-TST), as implemented in the distributed computer code VaReCoF (Harding, Georgievskii and Klippenstein, 2005). In this approach, variational optimisations are carried out for a series of flexible transition state dividing surfaces suitable for describing barrierless reactions such as fixed centre-of-mass dividing surfaces. Direct VRC-TST has been shown to be generally very accurate for barrierless reactions, often predicting capture rate constants to better than 25% (Klippenstein, 2017). Fragment interaction energies were calculated using CAS(5e,5o)PT2/CBS with the complete basis set (CBS) limit estimated using a two-point formula (Martin & Uzan, 1998) and Dunning’s AVTZ and AVQZ basis sets.

The MESS software (Georgievskii et al., 2013) was employed to calculate the overall reaction rate coefficients, which were proven to be equivalent to the rate of the first step (capture). We have also used this software to calculate the backward reaction rates up to 4889 K, in order to analyse the possibility that the H + PN collision leads to PN destruction in environments of high temperatures of the ISM, such as shock and star-forming regions. Another possible barrierless channel on the N+PH collisions (N+PH $\rightarrow$ P + NH) was also evaluated within this approach to check its influence in the branching ratio of the overall N+PH reactivity.

3. Results

3.1. Potential energy surface

Even though the title reactions can occur either on doublet, quartet or sextet potential energy surfaces (PESs), the ground state products [H( $^2S$ )+PN( $^1\Sigma^+$ )] can only be achieved adiabatically in the doublet state, and therefore we begin the discussion focusing only on the doublet PES. The main stationary structures obtained in this work and their connections are summarised in Figure 1, while Table 1 gathers the geometries, energies and parameters of MR character of all minima and TSs. The Cartesian coordinates and frequencies of all structures can be found in the supplementary information (SI). As seen in Figure 1, the CCSD(T)-F12 and MRCI-F12 energies agree very well, with deviations in the relative energies lying below 6 kJ mol $^{-1}$ , which is about the expected accuracy of these high level methods. The M06-2X results are also in qualitative agreement, being its worse performance on the well depth of the HNP structure, which deviated 24 kJ mol $^{-1}$ from the MRCI-F12 one.

Figure 1. Potential energy diagram for the doublet state of HPN. Energies are given relatively to the H + PN asymptote and are ZPE corrected. The MRCI-F12 energies are given as plain numbers, while the CCSD(T)-F12 and M06-2X ones results given under parenthesis and in boldface, respectively. Atoms are coloured as: hydrogen (white); nitrogen (blue); phosphorus (orange).

Table 1. Properties of the stationary structures of the doublet HPN PES $^a$ .

$^a$ Interatomic distances are given in Å and energies in kJ mol $^{-1}$ relative to the H+PN limit. The energies are ZPE corrected and given at the MRCI-F12/AVTZ+d//CAS/AVTZ+d level.

$^b$ Magnitude of the dominant configuration in the CASSCF wave functions.

$^c$ T $_1$ diagnostic obtained through CCSD(T)/AVTZ+d calculations.

As it can be noted in Table 1, the HPN system does not possesses a high MR character when we analyse the C ${_0^2}$ parameter. The only exception is the transition state termed TS $_{\rm HNP\rightarrow H+PN}$ . However, within the T $_1$ diagnostic both HPN minimum and TS $_{\rm HNP\rightarrow HPN}$ structures possess a considerable MR character. Interestingly, the CCSD(T)-F12 results are nevertheless very close to the MRCI-F12 ones, perhaps due to a cancellation of errors. Since the two high level methods agree very well, we will use in the following discussion only the MRCI-F12/AVTZ+d//CAS/AVTZ+d values for simplicity.

The $\textrm{N} + \textrm{PH} \longrightarrow \textrm{H} + \textrm{PN}$ and $\textrm{P} + \textrm{NH} \longrightarrow \textrm{H} + \textrm{PN}$ reactions are exothermic by 300 and 268 kJ mol $^{-1}$ , respectively. The difference between the two numbers is due to NH having a deeper potential well than PH. Given the isovalency of nitrogen and phosphorus, both reactions are analogous to the $\textrm{N} + \textrm{NH} \longrightarrow \textrm{H} + \textrm{N2}$ one. However the latter is more than two times more exothermic, owing to the strong triple bond present in $\textrm{N}_{2}$ . Another striking difference is the presence of two deep minima in the HPN potential energy surface, while ${\textrm{HN}_{2}}$ only displays a shallow metastable minimum which lies higher in energy than the ${\textrm{H} + \textrm{N}_{2}}$ limit (Walch, Duchovic and Rohlfing, 1989; Mota & Varandas, 2007, 2008; Mota, Galvão, Coura and Varandas, 2020). The presence of such minima will change the dynamics of the phosphorus reactions by allowing energy randomisation among the degrees of freedom during the complex lifetime, which will change the products distributions (Miller, Safron and Herschbach, 1967).

In principle, collisions starting on the quartet PES, which is also attractive, could contribute to the overall reactivity if they could hop to the doublet state after a collision complex is formed, and from there dissociate to the ground state products (recall that the quartet state cannot adiabatically lead to such products). This would be more likely if the quartet and doublet PESs crossed each other along the reaction coordinate, and for this reason, we searched for a such a crossing seam. At the MRCI-F12/AVTZ+d//CAS/AVTZ+d level, we have performed potential energy surface scans for the $\textrm{N} + \textrm{PH} \longrightarrow^4$ HPN and $\textrm{P} + \textrm{NH} \longrightarrow^4$ HNP attacks on the quartet surface by relaxing both the valence angle and diatomic distance along the path, and calculated a single point energy for the doublet state (at the optimised quartet geometries). This would allow us to see if there is a crossing lying on the minimum energy path from the reactants to the quartet minima. However, as it can be noted in Figures S1 and S2 of the SI, we have not found a crossing seam, and found that the quartet minima ( $^4$ HPN and $^4$ HNP) lie much higher in energy than the doublet ones.

Mechanism for reaction (3):

The N+PH collision begins with a barrierless capture towards the formation of the HPN minimum, which lies 387 kJ mol $^{-1}$ below the initial reactants. This structure may go through a barrierless hydrogen loss, or isomerise to the global minimum HNP with a barrier of 89 kJ mol $^{-1}$ (2.4 kJ mol $^{-1}$ relative to the final products). Assuming a long lived HPN structure and sufficient randomisation among the vibrational degrees of freedom, both possibilities will occur with similar probability. If isomerisation occurs, the HPN structure can either suffer a hydrogen loss, amounting to an indirect mechanism such as $\textrm{N} + \textrm{PH} \longrightarrow \textrm{HPN} \longrightarrow \textrm{PNH} \longrightarrow \textrm{H} + \textrm{PN}$ , or alternatively yield P+NH. The latter process will lead to different and less exoergic products than formation of PN, and should occur with a much lower probability (its rate coefficients are also given in following sections). Recall that the destruction of the PH molecule via N+PH towards H + PN is available in the astrochemical databases, and our work supports the assumption that this reaction is fast at low temperatures, since it does not involve any step requiring energies higher than that of the reactants.

It is worth mentioning that we found a negligible barrier for hydrogen loss from HPN above the H+PN channel of 2.3 kJ mol $^{-1}$ (Figure S3 of the SI file), lying below the accuracy of our calculations and is not given in Figure 1. The presence of such barrier would not change the overall rate coefficients, as it lies well below reactants, and thus is of no consequence for the present mechanism.

Mechanism for reaction (4):

The collisions between P+NH also occur without an initial activation barrier and leads to the HNP global minimum. Similarly to the previous one, after HNP formation, it may go directly to the PN products by hydrogen loss, but this time with a barrier of 187 kJ mol $^{-1}$ (8.4 kJ mol $^{-1}$ relative to the final products), or either isomerise with a barrier of 181 kJ mol $^{-1}$ . Note that we could not find the TS $_{\rm HNP\rightarrow H+PN}$ within the CCSD(T) approach, perhaps due to its MR character. The destruction of the NH molecule by atomic phosphorus may be an important source of interstellar PN, and this reaction is not available in most astrochemical database yet. As it will be discussed later, these rate coefficients have been very recently estimated (Douglas, Gobrecht and Plane, 2022).

H+PN:

Destruction of PN by collision with atomic hydrogen can only lead to endothermic products, and for this reason it may only occur in the hotter and more extreme environments in space such as shocks and circumstellar envelopes. This corroborates with the assumption made by (Millar, Bennett and Herbst, 1987) that the depletion of the PN molecule in its singlet ground state in cold environments of the ISM can only occur through ion-molecules reactions. However, the collision towards the HPN minimum is likely barrierless, and may be achieved even in colder regions. The formation of this transient species may play a significant role in the inelastic collisions H + PN( $\upsilon$ ,j) $\rightarrow$ H + PN( $\upsilon$ ’,j’), which will thermalise the PN molecule by the highly abundant hydrogen atoms. Interestingly, in spite of being an exotic molecule in the gas-phase, HPN (phosphorus nitride imide) is being recently studied in materials sciences as nanotubes, and also studied as a possible clean source of energy (Yang et al., 2018; Marchuk et al., 2014; Guo et al., 2005).

3.2. Reaction rates

The thermal rate coefficient for the $\textrm{N}(^4\textrm{S}) + \textrm{PH}(^3\Sigma^-)$ reaction found in the astrochemical databases was estimated (Smith, Herbst and Chang, 2004) mainly based on the chemical similarity with the reaction $\textrm{N}(^4\textrm{S}) + \textrm{NH}(^3\Sigma^-)$ . Using the modified Arrhenius formula (Equation 5), the rate coefficient reported by (Smith, Herbst and Chang, 2004), at the temperature range of 10–300 K, is of $\alpha = 5.00\times 10^{-11}\rm cm^3\,s^{-1}$ ( $\beta$ and $\gamma = 0.00\, \rm K$ ). This relatively high rate coefficient that has been estimated indicates that this reaction is possibly relevant in the cold environments of the ISM. However, accurate theoretical calculations or experimental data are relevant for a confirmation of the assumptions made and for a better understanding of the role of this reaction in modelling the PN abundance in the ISM. Furthermore the P+NH reaction studied here is not incorporated in the astrochemical databases.

(5) \begin{align}k(T)=\alpha\!\left( T/300 \right)^{\beta} \exp\!{(\!-\!\gamma/T).}\end{align}

As mentioned in the methodology section, the rates were calculated using direct variable reaction coordinate transition state theory (Klippenstein, 1991, 1992) (VRC-TST). Figure 2 shows the rate coefficients for both $\textrm{N} + \textrm{PH} \longrightarrow \textrm{H} + \textrm{PN}$ and $\textrm{P} + \textrm{NH} \longrightarrow \textrm{H} + \textrm{PN}$ reactions. Our results indicate that both have similar rate coefficients, with the $\textrm{P} + \textrm{NH}$ slightly higher, possibly due to a more attractive interaction between reactants due to the higher polarisability of the P atom. A fit of our VRC-TST results to the modified Arrhenius equation (Equation (5)) yields $\alpha = 0.88\times 10^{-10}\rm cm^3\,s^{-1}$ , $\beta=-0.18$ and $\gamma=1.01\, \rm K$ for $\textrm{N} + \textrm{PH}$ and $\alpha = 0.93\times 10^{-10}\rm cm^3\,s^{-1}$ , $\beta=-0.18$ and $\gamma=0.24\, \rm K$ for $\textrm{P} + \textrm{NH}$ .

Figure 2. Rate coefficients as a function of temperature. The points refer to the calculated values at fixed temperatures, while the curves are fits to the modified Arrhenius expressions. Solid lines correspond to the $\textrm{N} + \textrm{PH} \longrightarrow \textrm{H} + \textrm{PN}$ reaction, while dashed lines correspond to the $\textrm{P} + \textrm{NH} \longrightarrow \textrm{H} + \textrm{PN}$ . Our VRC-TST results are shown in black, while those calculated by (Douglas, Gobrecht and Plane, 2022) are brown. The values currently in use in the major astrochemical databases for $\textrm{N} + \textrm{PH}$ reaction is in orange.

A very recent work by (Douglas, Gobrecht and Plane, 2022) calculated the rate coefficients for both reactions studied in this work, using long-range transition state theory (Georgievskii & Klippenstein, 2005) (LR-TST). Their reported rates are also given in Figure 2 for comparison. Differently from our VRC-TST results, their values show that the rate coefficient for reaction P+NH is about three times higher than N+PH, with no temperature dependence for the former. Nevertheless, their results for N+PH agree well with the present ones. Since the VRC-TST method is more reliable, we trust that our results should be more accurate.

While the results already in use in the UMIST (McElroy et al., 2013) and KIDA (Wakelam et al., 2012) coincide well with our calculated ones for the N+PH (see Figure 2) reaction below 100 K, for high temperatures our results are about twice as high. It must be stressed here once again that the P+NH reaction is not incorporated in the databases, but it may also contribute for a better modelling the abundances of PN in the ISM.

We have also checked the possibility of $\textrm{N} + \textrm{PH} \longrightarrow \textrm{P} + \textrm{NH}$ (which is also slighly exothermic and barrierless) and its reverse, and the thermal rate coefficients are gathered in Figure 3. However our results show that the rate coefficient for this reaction at the temperature of 328 K $(0.57\times 10^{-12}\rm cm^3\,s^{-1}$ , see Figure 3) is 157 $\times$ lower than the $\textrm{N} + \textrm{PH} \longrightarrow \textrm{H} + \textrm{PN}$ one, and therefore should not be of relevance and contribute to the branching ratio. The reverse reaction is slightly endothermic and displays a typical Arrhenius behaviour of a reaction that shows an activation energy.

Figure 3. Rate coefficients as a function of temperature for the $\textrm{P} + \textrm{NH} \longrightarrow \textrm{N} + \textrm{PH}$ (black) and $\textrm{N} + \textrm{PH} \longrightarrow \textrm{P} + \textrm{NH}$ (blue) reactions. The points refer to the calculated values at fixed temperatures, while the curves are fits to the modified Arrhenius expressions.

3.3. Destruction of PN by H atoms in high temperature regions

PN is also present in high temperature environments, such as star forming regions and shock ones. A recent work performed by our group (Gomes et al., 2022) showed for the first time that in such high temperature environments, atomic nitrogen plays an important role in the destruction, and therefore, the respective abundance of PN molecules. That being noticed, it would be interesting to see at what temperature this molecule can be destructed by H atoms, which is the most abundant atomic species in space. However, as it can be noted in Figure 4, the rate coefficients are very low, even at 5 000 K, making this destruction route a not particularly important one. We have compared these results with the ones reported by (Douglas, Gobrecht and Plane, 2022) (see Figure 4), and apart the expected differences arising from the different methods of calculation, we arrive at the same qualitatively conclusions. Table 2 gathers the rate coefficients as a function of temperature for all reactions studied in this work and discussed in this and in the previous sections, using a modified Arrhenius equation fit.

Figure 4. Rate coefficients as a function of temperature for the destruction of PN by H atoms. The points refer to the calculated values at fixed temperatures, while the curves are fits to the modified Arrhenius expressions. Solid lines correspond to the $\textrm{H} + \textrm{PN} \longrightarrow \textrm{N} + \textrm{PH}$ reaction, while dashed lines correspond to the $\textrm{H} + \textrm{PN} \longrightarrow \textrm{P} + \textrm{NH}$ Our results using MESS are shown in black, while those calculated by (Douglas, Gobrecht and Plane, 2022) are orange and brown for the respective mentioned reactions.

Table 2. Rate coefficients as a function of temperature for all reactions studied in this work using a modified Arrhenius equation fit.

4. Conclusions

The elucidation of PN formation in the ISM, and the general astrochemical behaviour of phosphorus is very important for the understanding of how this element became so important in our planet. In this work, we performed MRCI-F12/AVTZ+d//CAS/AVTZ+d calculations, in order to study the potential energy surfaces of two reactions that lead to the important phosphorus bearing molecule PN. We also calculated their rate coefficients, using direct variable reaction coordinate transition state theory (VRC-TST). We found that both reactions proceed without an activation barrier being relevant for PN formation in the interstellar medium. The P+NH reaction is not incorporated in the databases and, therefore, should be included for better modelling the abundances of PN.

Regarding the kinetics, collision $\textrm{N}(^4\textrm{S}) + \textrm{PH}(^3\Sigma^-)$ possesses a rate coefficient of $\alpha=0.88\times 10^{-10}\rm cm^3\,s^{-1}$ , $\beta=-0.18$ and $\gamma=1.01\, \rm K$ , while $\textrm{P}(^4\textrm{S}) + \textrm{NH}(^3\Sigma^-)$ the rate coefficient is $\alpha=0.93\times 10^{-10}\rm cm^3\,s^{-1}$ , $\beta=-0.18$ and $\gamma=0.24\, \rm K$ . Therefore, both reactions contribute equally for products formation. The other barrierless channel (N+PH $\rightarrow$ P+NH) possesses a much lower rate constant, and can be considered negligible for the total branching ratio. The results for the destruction of PN by hydrogen atoms support previous assumptions that this can only happen in cold environments of the ISM through ion-molecule reactions.