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., Reference Rivilla, Drozdovskaya, Altwegg, Caselli, Beltrán, Fontani and van der Tak2020; Sousa-Silva et al., Reference Sousa-Silva, Seager, Ranjan, Petkowski, Zhan, Hu and Bains2020; Maciá, Reference Maciá2005) and in metabolic processes of living organisms (Maciá, Reference Maciá2005; Goldford et al., Reference Goldford, Hartman, Smith and Segrè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., Reference Jiménez-Serra, Viti, Quénard and Holdship2018; Chantzos et al., Reference Chantzos, Rivilla, Vasyunin, Redaelli, Bizzocchi, Fontani and Caselli2020). The question of how it became present in Earth’s prebiotic chemistry and biologically available is also an open problem (Goldford et al., Reference Goldford, Hartman, Smith and Segrè2017; Rivilla et al., Reference Rivilla, Drozdovskaya, Altwegg, Caselli, Beltrán, Fontani and van der Tak2020).
Phosphorus nitride (PN) was the first P-bearing molecule detected in space, which occurred in the late eighties in several molecular clouds (Turner & Bally, Reference Turner and Bally1987; Ziurys, Reference Ziurys1987). It was later also observed in dense clouds and in their massive cores (Ziurys, Schmidt and Bernal, Reference Ziurys, Schmidt and Bernal2018). This molecule is one of the most important reservoirs of phosphorus in the gas phase in the ISM (Thorne et al., Reference Thorne, Anicich, Prasad and Huntress1984; Tenenbaum, Woolf and Ziurys, Reference Tenenbaum, Woolf and Ziurys2007; Lefloch et al., Reference Lefloch, Vastel, Viti, Jiménez-Serra, Codella, Podio, Ceccarelli, Mendoza, Lépine and Bachiller2016; Yamaguchi et al., Reference Yamaguchi, Takano, Sakai, Sakai, Liu, Su and Hirano2011) 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., Reference Fontani2019).
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., Reference Lefloch, Vastel, Viti, Jiménez-Serra, Codella, Podio, Ceccarelli, Mendoza, Lépine and Bachiller2016; Tenenbaum, Woolf and Ziurys, Reference Tenenbaum, Woolf and Ziurys2007). Both PN and PO have been identified in massive star-forming regions (Rivilla et al., Reference Rivilla, Fontani, Beltrán, Vasyunin, Caselli, Martn-Pintado and Cesaroni2016) 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., Reference Jiménez-Serra, Viti, Quénard and Holdship2018). 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., Reference Jiménez-Serra, Viti, Quénard and Holdship2018; Rivilla et al., Reference Rivilla, Drozdovskaya, Altwegg, Caselli, Beltrán, Fontani and van der Tak2020; Concepción et al., Reference de la Concepción, Puzzarini, Barone, Jiménez-Serra and Roncero2021).
The most commonly used neutral-neutral reactions in astrochemical databases that explain the formation of PN were proposed by (Millar, Bennett and Herbst, Reference Millar, Bennett and Herbst1987):
\begin{align} {}\textrm{N}(^4\textrm{S}) + \textrm{PO}(^2\Pi) \longrightarrow \textrm{O}(^3\textrm{P}) + \textrm{PN}(^1\Sigma^+) \end{align}
\begin{align} \textrm{N}(^4\textrm{S}) + \textrm{PO}(^2\Pi) \longrightarrow \textrm{P}(^4\textrm{S}) + \textrm{NO}(^2\Pi) \end{align}
\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, Reference Souza, Silva and Galvão2021). However, the rate coefficients employed in astrochemical databases such as UMIST (McElroy et al., Reference McElroy, Walsh, Markwick, Cordiner, Smith and Millar2013) and KIDA (Wakelam et al., Reference Wakelam, Herbst, Loison, Smith, Chandrasekaran, Pavone, Adams, Bacchus-Montabonel, Bergeat and Béroff2012) for reaction (3) was estimated based on the chemical similarity with the reaction N+NH
$\rightarrow$
N
$_2$
+H (Smith, Herbst and Chang, Reference Smith, Herbst and Chang2004), which in turn had its rate constant already studied theoretically and experimentally (Caridade et al., Reference Caridade, Rodrigues, Sousa and Varandas2005; Hack, Wagner and Zasypkin, Reference Hack and Wagner1994; Miller et al., Reference Miller, Branch, Mclean, Chandler, Smooke and Kee1985).
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, Reference Douglas, Gobrecht and Plane2022) 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, Reference Georgievskii and Klippenstein2005) (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:
\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., Reference Wakelam, Herbst, Loison, Smith, Chandrasekaran, Pavone, Adams, Bacchus-Montabonel, Bergeat and Béroff2012; McElroy et al., Reference McElroy, Walsh, Markwick, Cordiner, Smith and Millar2013)
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., Reference Schmidt, Baldridge, Boats, Elbert, Gorgon, Jensen and Koseki1993) and MOLPRO (Werner et al., Reference Werner, Knowles, Knizia, Manby and Schütz2012) packages. For preliminary exploratory work, we chose the density functional theory (DFT) (Kohn & Sham, Reference Kohn and Sham1965) using the aug-cc-pV(T+d)Z (Dunning, Reference Dunning1989; Kendall, Dunning and Harrison, Reference Kendall, Dunning and Harrison1992) basis set and the exchange and correlation functional M06-2X (Zhao & Truhlar, Reference Zhao and Truhlar2008). 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., Reference Szalay, Müller, Gidofalvi, Lischka and Shepard2012) (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., Reference Szalay, Müller, Gidofalvi, Lischka and Shepard2012; Shiozaki, Knizia and Werner, Reference Shiozaki, Knizia and Werner2011) (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, Reference Bartlett1989; Bartlett et al., Reference Bartlett, Watts, Kucharski and Noga1990; Raghavachari et al., Reference Raghavachari, Trucks and Pople1989) 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, Reference Adler, Knizia and Werner2007; Knizia, Adlerand Werner, Reference Knizia, Adler and Werner2009) 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, Reference Lee and Taylor1989; Leininger et al., Reference Leininger, Nielsen, Crawford and Janssen2000; Lee, Reference Lee2003); 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, Reference Bode and Gordon1998; Hanwell et al., Reference Hanwell, Curtis, Lonie, Vandermeersch, Zurek and Hutchison2012).
2.2. Rate coefficients
Capture rate coefficients were calculated using direct variable reaction coordinate transition state theory (Klippenstein, Reference Klippenstein1991, Reference Klippenstein1992) (VRC-TST), as implemented in the distributed computer code VaReCoF (Harding, Georgievskii and Klippenstein, Reference Harding, Georgievskii and Klippenstein2005). 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, Reference Klippenstein2017). 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, Reference Martin and Uzan1998) and Dunning’s AVTZ and AVQZ basis sets.
The MESS software (Georgievskii et al., Reference Georgievskii, Miller, Burke and Klippenstein2013) 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.
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).
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, Reference Walch, Duchovic and Rohlfing1989; Mota & Varandas, Reference Mota and Varandas2007, Reference Mota and Varandas2008; Mota, Galvão, Coura and Varandas, Reference Mota, Galvão, Coura and Varandas2020). 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, Reference Miller, Safron and Herschbach1967).
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, Reference Douglas, Gobrecht and Plane2022).
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, Reference Millar, Bennett and Herbst1987) 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., Reference Yang, Chen, Liu, Liu, Li, Dong, Chen, Yan, Yao and Duan2018; Marchuk et al., Reference Marchuk, Pucher, Karau and Schnick2014; Guo et al., Reference Guo, Yang, Zhu, Yi and Xie2005).
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, Reference Smith, Herbst and Chang2004) 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, Reference Smith, Herbst and Chang2004), 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.
\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, Reference Klippenstein1991, Reference Klippenstein1992) (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}$
.
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, Reference Douglas, Gobrecht and Plane2022) 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, Reference Douglas, Gobrecht and Plane2022) calculated the rate coefficients for both reactions studied in this work, using long-range transition state theory (Georgievskii & Klippenstein, Reference Georgievskii and Klippenstein2005) (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., Reference McElroy, Walsh, Markwick, Cordiner, Smith and Millar2013) and KIDA (Wakelam et al., Reference Wakelam, Herbst, Loison, Smith, Chandrasekaran, Pavone, Adams, Bacchus-Montabonel, Bergeat and Béroff2012) 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.
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., Reference Gomes, Spada, Lefloch and Galvão2022) 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, Reference Douglas, Gobrecht and Plane2022) (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.
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, Reference Douglas, Gobrecht and Plane2022) are orange and brown for the respective mentioned reactions.
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.
Acknowledgements
The authors would like to thank the financial support provided by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) - Finance Code 001, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant 311508-2021-9, Fundação de Amparo à Pesquisa do estado de Minas Gerais (FAPEMIG), and Centro Federal de Educação Tecnológica de Minas Gerais (CEFET-MG). Rede Mineira de Química (RQ-MG) is also acknowledged. We are also thankful for computational resources provided by LNCC and the Santos Dumont supercomputer. Ahren W. Jasper would like to acknowledge the U. S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under Contract Number DE-AC02-06CH11357.
Conflicts of interest
We have no conflicts of Interest to declare.
Data availability
The data underlying this article are available in the article and in its online supplementary material.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/pasa.2023.13.
