© The Authors 2025

1 Introduction

Phosphorus is present in several molecules that are essential for life as we know it, although it is still a mystery how it became biologically available on Earth, possibly from an inorganic mineral or even extraterrestrial sources (Pasek & Lauretta 2005; Schwartz 2006). When combined with four oxygen atoms, it yields the phosphate ion (${\mathrm{P}\mathrm{O}}_4^{3-}$), which is highly soluble in water and whose salts are plentiful on Earth and were readily available throughout the development of life (Hunter 2012).

There is now biochemical evidence that indicates that phosphate may indeed have played an important role in the chemical origins of life. Several feasible environments capable of producing phosphates are currently conceivable, providing numerous potential solutions to the issue of phosphate availability (Fernández-García et al. 2017). Phosphorus, as well, can form P−N (phosphoramidate), P−S (phosphorothioate) and P−C (phosphonate) linkages, all of which are present in nature. For instance, phosphorylation of amino acids such as histidine, lysine and arginine to generate high energy phosphoramidate bonds may have been crucial under prebiotic circumstances (Hunter 2012).

The prebiotic chemistry (including P-bearing molecules) on our planet may have been enriched by the bombardment of extraterrestrial bodies such as comets and asteroids (Pasek & Lauretta 2005; Altwegg et al. 2016; Fernández-Ruz et al. 2023). For this reason, the investigation of chemical processes taking place on different astronomical environments such as molecular clouds and star forming-regions is of great importance to our understanding of the primordial chemistry of phosphorus on Earth.

The detection of most P-bearing species in space has been restricted so far to the circumstellar envelopes of evolved stars (Guélin et al. 1990; Agúndez et al. 2007; Halfen et al. 2008; Milam et al. 2008; Agúndez et al. 2014a; Agúndez et al. 2014b) in our Galaxy and planetary atmospheres (Bregman et al. 1975; Larson et al. 1977). However, special attention has been given to phosphorus monoxide (PO) and phosphorus nitride (PN), which were the only ones also detected towards star-forming regions (Fontani et al. 2016; Rivilla et al. 2016; Mininni et al. 2018; Bergner et al. 2019; Fontani et al. 2019; Rivilla et al. 2020; Bernal et al. 2021), in addition to the interstellar medium (ISM) (Turner & Bally 1987; Ziurys 1987), circumstellar envelopes of evolved stars (Tenenbaum et al. 2007; Milam et al. 2008; Ziurys et al. 2018), protostellar shocks (Yamaguchi et al. 2011; Lefloch et al. 2016) and Galactic Centre molecular clouds (Rivilla et al. 2018).

It has already been discussed that PO and PN may constitute the main reservoir of phosphorus in the gas phase in the ISM (Thorne et al. 1984; Tenenbaum et al. 2007; De Beck et al. 2013; Lefloch et al. 2016). In low-mass star-forming regions, PN was the first P-bearing molecule to be detected (Yamaguchi et al. 2011), while PO was first detected five years later towards the same solar-type star-forming region (L1157-B1) (Lefloch et al. 2016). The abundances of these species in such environments seem to be comparable, with PO being systematically more abundant than PN (Rivilla et al. 2016; Lefloch et al. 2016; Rivilla et al. 2018; Bergner et al. 2019). However, the small number of detection makes it difficult to precisely elucidate the mechanisms of formation of these molecules.

More recently, the presence of PO and PN in the gas phase has also been reported for the first time at the edge of our galaxy (Koelemay et al. 2023). The characteristic rotational lines of these P-bearing species were observed in the dense cloud WB89-621 (located 22.6 kpc from the Galactic Centre) and their abundances are comparable to values near the Solar System (Bernal et al. 2021; Koelemay et al. 2023). Specifically, the phosphorus nitride molecule has been detected for the first time in an extra-galactic environment towards two giant molecular clouds of NGC 253 (Haasler et al. 2022). The authors report a derived abundance ratio between PN and the shock-tracer SiO that follows the same trend already found towards Galactic sources (Rivilla et al. 2018; Haasler et al. 2022).

Noticeably, the [PO]/[PN] ratio has consistently been found to be >1, independently of the source, with values up to ∼6 (Haasler et al. 2022; Fernández-Ruz et al. 2023). However, the available astrochemical models most frequently fail to reproduce this trend, usually predicting [PO]/[PN] ratios <1 for various different physical conditions (Jiménez-Serra et al. 2018; Chantzos et al. 2020; Sil et al. 2021; Fernández-Ruz et al. 2023). Various factors could contribute to such disparities, including the inadequacy of the phosphorus chemical network and the considerable uncertainties associated with reaction rate coefficients (Fernández-Ruz et al. 2023). Currently, there is compelling evidence indicating that PO emerges in the gaseous phase prior to PN and then vanishes, whereas PN lingers on for a certain time (Gomes et al. 2023b). The precise mechanisms accountable for their formation and decomposition are yet to be better delineated and their study may solve the inconsistency between astrochemical models and observations.

A computational approach can help in understanding the chemistry of phosphorus in the gas phase since experimental studies involving these P-bearing compounds are scarce in the literature. García de la Concepción et al. (2021) obtained rate coefficients and branching ratios for the reactions P + OH → PO + H and P + H₂O → PO + H₂ using ab initio transition state theory (TST) and employing high-level quantum chemical calculations. The authors conclude that the formation of PO, in the gas phase, from OH and P in their ground state is highly favourable, but should not occur if P is in its first excited electronic state. As for the reaction between P and H₂O, the predicted trend is the opposite. More recently, the P(⁴S) + O₂(³Σ⁻) and P(⁴S) + O₂(¹Δ) reactions were also investigated by Gomes et al. (2022) and García de la Concepción et al. (2024) from a theoretical point of view, and the mechanisms towards PO formation could be elucidated. It has been concluded that multireference methods are needed to properly describe the system and that the P(⁴S) + O₂(³Σ⁻) reaction dominates the formation of PO when OH is not sufficiently available (for instance, in regions where shock speeds are below 40 km s⁻¹). Structural, molecular and spectroscopic properties of the HPO system were studied using the coupled cluster method by Francisco (2003) and Puzzarini (2006), and isomerisation and dissociation energies were obtained, showing that the two processes may compete with each other.

Plane et al. (2021) estimated a series of rate coefficients for several reactions of the H_xPO_y system from electronic structure calculations and subsequent application of statistical rate theory. In the following year, Douglas et al. (2022) gathered an extensive network of reactions involving phosphorus with estimates for their rate coefficients obtained from electronic structure calculations at the Complete Basis Set (CBS-QB3) level (Montgomery et al. 2000).

Our research group has also approached the NPO system using high-level electronic structure methods (Souza et al. 2021) and quantitatively confirmed that the depletion of PO by N atoms occurs fast, with a branching ratio favourable to the O + PN channel relative to P + NO. The HPN system was also addressed (Gomes et al. 2023a) using ab initio calculations and the variable reaction coordinate transition state theory (VRC-TST) and it was predicted that the reactions P(⁴S) + NH(³Σ⁻) and N(⁴S) + PH(³Σ⁻) occur without an activation barrier and are relevant for PN formation in the interstellar medium.

Astrochemical models, however, appear to lack important PN destruction pathways (Fernández-Ruz et al. 2023). This may be one reason for the failure to correctly predict the [PO]/[PN] ratio. Previous studies from our research group indicate that PN destruction by H atoms should not occur (Gomes et al. 2023a). However, PN destruction by N atoms may occur in high temperature environments, such as in shock regions (Gomes et al. 2023b). Jiménez-Serra et al. (2018) drew attention to the fact that the reactions N + CP → PN + C and P + CN → PN + C, proposed by Agúndez et al. (2007), have very uncertain rate coefficients and that the results of different models are little affected by the inclusion of these reactions. In the work of Fernández-Ruz et al. (2023), the parameters associated with the analogous reaction N + CN → N₂ + C were used in the rate coefficients of the N + CP and P + CN reactions, assuming that the similarity between N and P should lead to similar chemical behaviour.

Therefore, in the present work, we focus on the computational study of the CPN system using the highest-level electronic structure methods available combined with TST. The main goals are to elucidate the viable reaction pathways involving simple atom-diatom collisions within this system (C³P) + PN(¹ Σ⁺), N(⁴S) + CP(² Σ⁺) and P(⁴S) + CN(²Σ⁺)), providing the respective rate coefficients and evaluating a possible route for the destruction of PN in the interstellar medium. We summarise other relevant reactions from the literature, and include them in the Kinetic Database for Astrochemistry (KIDA) (Wakelam et al. 2012) network to investigate the most important reactions that drive the abundances of PO and PN in the ISM.

2 Methodology

2.1 Electronic structure calculations

All electronic structure calculations reported here were performed using the MOLPRO (Werner et al. 2015) package. For the initial exploration of the potential energy surface (PES) of the PCN system, a density functional theory (DFT) (Kohn & Sham 1965) approach was employed using the aug-cc-pV(T+d)Z (Dunning 1989; Kendall et al. 1992) basis set (henceforth referred to as AV(T + d) Z) with the M06-2X (Zhao & Truhlar 2008) exchange and correlation functional. At this stage, stationary points were found and vibrational analysis was performed to confirm whether the obtained structures corresponded to energy minima or to first order saddle points on the PES. Intrinsic reaction coordinate (IRC) calculations were performed for all transition states (TSs) found in order to ensure their connection to the correct energy minima structures. All the obtained energies were zero-point energy (ZPE) corrected. In order to improve the accuracy of the energy values, further single-point energy calculations at the explicitly correlated coupled cluster (CCSD(T)-F12/cc-pVTZ-F12) (Adler et al. 2007; Knizia et al. 2009) level were performed over the geometries optimised at the DFT level and we name this methodology CCSD(T)-F12/VTZ-F12//M06-2X/AV(T+d)Z + ZPE(M06-2X/AV(T+d)Z). We have not employed scaling factors in the vibrational frequencies.

In a second stage, in order to account for a possible multiconfigurational character of the PCN system, all the previously obtained results were fully re-optimised within the ab initio full valence complete active space self-consistent field (CASSCF) (Szalay et al. 2012) level with AV(T + d) Z basis set. In this approach, the 1 s orbitals of C and N, and the 1s, 2s and 2p orbitals of P were considered as core orbitals, leading to an active space of 12 orbitals with 14 electrons. Vibrational analysis was also carried out and ZPEs were obtained for all stationary points on the CASSCF PES. This time, the accuracy of the energy values was enhanced by performing singlepoint energy calculations over all structures with the explicitly correlated multireference configuration interaction (MRCI-F12) (Szalay et al. 2012; Shiozaki et al. 2011) method, including the relaxed-reference Davidson correction to approximate quadruple excitations and also using the AV(T + d) Z basis set. The results reported for this methodology are referred here as MRCIF12/AV(T+d)Z//CAS/AV(T+d)Z + ZPE(CAS/AV(T+d)Z).

The methodologies employed here could also provide the magnitude of the dominant configuration in the CI wavefunctions (${\mathrm{C}}_0^2$) and the T₁ diagnostic of all stationary points, which allows the assessment of the multireference character (MR) of the system (Lee & Taylor 1989; Leininger et al. 2000; Lee 2003). Lastly, the Wxmacmolplt and Avogadro software were used for graphical visualisation and representation of the molecular geometries (Bode & Gordon 1998; Hanwell et al. 2012).

2.2 Rate coefficients

In order to assess the effects of the reactions studied in this work on the relative abundances of the involved species, their rate coefficients must be calculated so that the chemical processes taking place in an interstellar environment can be simulated. We have obtained rigid rotor and harmonic oscillator rate coefficients using standard TST and including the Eckart tunnelling correction (Eckart 1930) employing the Master Equation System Solver (MESS) software package (Georgievskii et al. 2013). Phase space theory (PST) (Pechukas & Light 1965; Chesnavich 1986), as implemented in MESS, is used for the barrierless capture steps. For such cases, we have fitted the MRCI-F12/AV(T+d)Z energies obtained along the minimum energy path as a function of the internuclear separation (r) to the −C₆/r⁶ function to obtain C₆, which is then used in the PST calculation.

2.3 Astrochemical modelling

We utilised the Nautilus gas-grain code to simulate chemical abundances under interstellar conditions, focusing on the gas-phase and grain-surface processes described in the three-phase model by Ruaud et al. (2016). The gas-phase chemical network is based on the comprehensive kida.uva. 2014 database (Wakelam et al. 2015), which compiles extensive reaction pathways and kinetic data. Although the kida.uva. 2024 update has introduced new species and reaction channels, it does not modify phosphorus chemistry (Wakelam et al. 2024). In this work, we expanded the chemical network within Nautilus to include updated reaction pathways for phosphorus-containing molecules. This enabled the prediction of abundances for key astrochemical species such as PO and PN, whose formation mechanisms in extreme environments remain poorly understood. The study also assessed the significance of individual reaction pathways in phosphorus chemistry. Cosmic ray and X-ray ionisation rates, crucial drivers of interstellar chemical reactions, were systematically varied in the simulations. The default cosmic ray ionisation rate (ζ_CR) is typically set to 1.3 × 10⁻¹⁷ s⁻¹ (Dalgarno 2006), but can reach 10⁻¹⁵ s⁻¹ in regions like the central molecular zone of Sgr B2 (Goto et al. 2014). Similarly, X-ray ionisation rates (ζ_X) influence the chemistry of small neutrals and molecular ions, particularly in protostellar envelopes (Stäuber et al. 2005; Notsu et al. 2021). We tested ζ_CR values from 1.3 × 10⁻¹⁷ to 1.3 × 10⁻¹⁶ s⁻¹ and ζ_x values from 1 × 10⁻¹⁸ to 1 × 10⁻¹⁷ s⁻¹ to evaluate their effects on phosphorus chemistry.

Simulations were conducted to examine phosphorus chemistry under two distinct astrophysical conditions: diffuse and translucent clouds and dense clouds (Corby et al. 2018; Phan et al. 2018; Wolfire et al. 2022). In the diffuse and translucent model, the H₂ gas density was set to 10³ cm⁻³, with a temperature of 70 K and a visual extinction (A_V) of less than 5 mag. The simulation covered a timescale of 1 Myr, during which cosmic ray and X-ray ionisation rates were varied to assess their impact on chemical processes. For the dense molecular cloud model, the gas density was increased to 10⁴ cm⁻³, while the temperature was lowered to 30 K and A_V was fixed at 10 mag. Simulations extended from 1 × 10⁶ yr up to 1 × 10⁸ yr to explore long-term chemical evolution. Ionisation rates were held constant at ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹ and ζ_x = 1 × 10⁻¹⁸ s⁻¹. In both models, abundances were calculated as functions of time. The non-linear nature of the chemical network and the dynamic interplay of reactions prevented the system from reaching full equilibrium, influencing the steady-state abundances (Roueff & Le Bourlot 2020; Dufour et al. 2023). In the dense cloud model, extended simulations provided insights into the stability and evolutionary patterns of phosphorus-bearing species. Individual reaction contributions were evaluated to quantify their roles in shaping chemical pathways and the long-term persistence of these molecules in interstellar environments.

3 Results

3.1 Reaction products

As we expect the reactants (C(³P) + PN(¹Σ⁺)) to be mainly available in their electronic ground state, we have only considered electronic states with triplet spin multiplicities. The C_s symmetry of the system was taken into consideration and the first electronic level of each irreducible representation (A^′ and A^′′) was calculated. The gas-phase chemical reactions investigated in this work, with the respective estimated energies at MRCI-F12/AV(T+d)Z//CAS/AV(T+d)Z + ZPE(CAS/AV(T+d)Z) level, are listed below: $$\mathrm{C}\left({}^3\mathrm{P}\right)+\mathrm{P}\mathrm{N}\left({}^1{\mathrm{\Sigma }}^{+}\right)⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)-140\mathrm{k}\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$$(1) $$⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2\mathrm{\Pi }\right)-33\mathrm{k}\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$$(2) $$⟶\mathrm{N}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{P}\left({}^2{\mathrm{\Sigma }}^{+}\right)83\mathrm{k}\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$$(3) $$\mathrm{N}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{P}\left({}^2{\mathrm{\Sigma }}^{+}\right)⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)-223\mathrm{k}\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}.$$(4)

The only exoergic products obtained from the C + PN reactants are P + CN. Note that the formation of excited state cyano radical CN(²Π) is also exoergic. The cyano radical has been widely observed in several astronomical environments (McKellar 1940; Adams 1941; Jefferts et al. 1970; Federman et al. 1984; Henkel et al. 1988; Fray et al. 2005; Han et al. 2015; Paron et al. 2021) and one should expect it to be often involved in high energy chemistry. Indeed, from the reactions studied here, CN is not expected to be destroyed by collisions with P (no exoergic products available).

Our calculations also indicate that carbon monophosphide (CP), which has already been detected in space (Guélin et al. 1990), may be promptly destroyed by collisions with neutral N atoms in the highly exoergic N + CP ⟶ P + CN reaction (Eq. (4)). Another possible pathway for CP destruction by N atoms is N + CP ⟶ C + PN (Eq. (3), the reverse of Eq. (3)) with an energy release of −83 kJ mol⁻¹. Our computational results on this reaction corroborates the conclusion of Chantzos et al. (2020) that this reaction may be an important route of PN formation. In addition, Chantzos et al. (2020) also predict that the main formation pathway for PN in the early stage of evolution of diffuse and translucent clouds is P + CN ⟶ C + PN, while according to our calculations this reaction (Eq. (1), the reverse of Eq. (1)) is estimated to be significantly endoergic. Therefore, it should not contribute to PN formation in cold environments such as the ones modelled in the work of Chantzos et al. (2020).

In fact, Chantzos et al. (2020) included the N + CP ⟶ C + PN and P + CN ⟶ C + PN reactions in their chemical network of phosphorus as proposed by Agúndez et al. (2007), which employed a rate coefficient equal to that of the nitrogen analogous reaction (N + CN ⟶ C + N₂). The value of k = 3 × 10⁻¹⁰ cm³ s⁻¹ was recommended over the range 300–2500 K and, at low temperatures, the measurements are reported to scatter over four orders of magnitude (Baulch et al. 1992). Moreover, the N + CN ⟶ C + N₂ reaction has also been studied computationally (Moskaleva & Lin 2001) and is predicted to be highly exothermic (with an energy release of −205 kJ mol⁻¹), presenting an energy landscape that is very different from that of the true chemical process investigated here. The use of the nitrogen analogous reaction is thus inadequate to properly describe the chemical evolution of these P-bearing species in cold regions of the ISM.

3.2 Potential energy surface and reaction mechanisms

Along with the thermodynamic analysis presented in the previous section, an investigation of the kinetic viability of the reactions is necessary. The existence of small potential energy barriers, for instance, can turn the formation of thermodynamically favourable products negligible in cold parts of the ISM.

To elucidate the possible reaction pathways of the CPN system, we have explored its triplet A^′ and A^′′ potential energy surfaces and obtained four relevant potential energy minima and four TSs. Table 1 gathers their geometries and energies, while Fig. 1 shows a potential energy diagram summarising our results and the connections between the stationary structures. The Cartesian coordinates and frequencies of all structures can be found in the supplementary information file.

Starting with the C + PN collision, the ground state for both product channels (P + CN and N + CP) can only be achieved adiabatically in the ³A^′′ state (see Fig. 1). The CN radical in its first excited electronic state can also be obtained adiabatically, but on the ³A^′ surface. Despite being an exoergic process, it involves surpassing an energy barrier of 42 kJ mol⁻¹ in the entrance channel for the C(³P) + PN(¹Σ⁺) reaction, which should inhibit the formation of excited CN(A²Π) at low temperatures.

It can be seen from Fig. 1 that the relative energies obtained at CCSD(T)-F12 and MRCI-F12 levels are in good agreement, deviating from each other by a maximum of 6 kJ mol⁻¹, except for the TS_{CPN → CNP} and TS_{CPN → N + CP} cases. This may be attributed to a high multiconfigurational character of these structures (the values for the T₁ diagnostic is 0.056 and 0.080, respectively). The M06-2X results agree qualitatively with the MRCIF12 ones, with relative energy deviations mostly within 10%. The largest deviation observed in the DFT approach was for the CNP (³A^′) minimum, reaching about 22% (an absolute deviation of 30 kJ mol⁻¹). From now on the discussion will be based on the more trustworthy MRCI-F12/AV (T+d)Z//CAS/AV (T+d)Z + ZPE(CAS/AV(T+d)Z) results.

We first address the possible routes of PN destruction by a carbon atom. If the incoming C atom attacks the phosphorus end of the PN molecule, one can notice in Fig. 1 that this can lead to a bent CPN intermediate through a barrierless path on the ³A^′′ PES. Subsequently, this minimum can either isomerise to a linear CNP intermediate or dissociate to N(⁴S) + CP(²Σ⁺). This abstraction mechanism C + PN ⟶ CPN ⟶ N + CP requires, nevertheless, the surpassing of a 84 kJ mol⁻¹ reaction barrier (relative to C + PN) that separates the CPN minimum from the products. Besides the high energy barrier, the products themselves lie 83 kJ mol⁻¹ above the reactants. Therefore, this mechanism can be considered closed for typical ISM conditions. Furthermore, the isomerisation of the CPN minimum to the CNP one occurs through a conventional TS that lies 40 kJ mol⁻¹ above CPN, but is submerged relatively to reactants. The CNP isomer can also undergo a second isomerisation to yield the linear NCP intermediate, which corresponds to the global minimum of this triatomic system. This second isomerisation also occurs through a conventional TS that is also submerged relatively to reactants. Direct dissociation from the CNP and NCP minima can lead to P(⁴S) + CN(²Σ⁺), without an exit barrier. The C+ PN ⟶ P(⁴ S) + CN(²Σ⁺) reaction is thus a viable PN destruction route (exoergic and without an activation energy). It is open for whatever temperature, and may be important to model the abundance of PN in the ISM.

The P(⁴S) + CN(²Σ⁺) products can also be obtained through a simpler mechanism. For example, by the initial approach of the carbon atom to the diatomic PN via the N atom to form the CNP intermediate, which can then lose the P atom in an overall barrierless process (C+PN ⟶ CNP ⟶ P(⁴S) + CN(²Σ⁺)).

In summary, the ³A^′′ PES allows for the possibility of several mechanisms of PN destruction by a carbon atom, all of which are listed below: $$\mathrm{R}⟶\mathrm{C}\mathrm{N}\mathrm{P}⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)$$(5) $$\mathrm{R}⟶\mathrm{C}\mathrm{N}\mathrm{P}⟶\mathrm{N}\mathrm{C}\mathrm{P}⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)$$(6) $$\mathrm{R}⟶\mathrm{C}\mathrm{P}\mathrm{N}⟶\mathrm{C}\mathrm{N}\mathrm{P}⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)$$(7) $$\mathrm{R}⟶\mathrm{C}\mathrm{P}\mathrm{N}⟶\mathrm{C}\mathrm{N}\mathrm{P}⟶\mathrm{N}\mathrm{C}\mathrm{P}⟶\mathrm{P}\left({}^4\mathrm{S}\right)+\mathrm{C}\mathrm{N}\left({}^2{\mathrm{\Sigma }}^{+}\right)$$(8)

where R means reactants (C(³P)+PN(¹Σ⁺)).

When it comes to N + CP collisions, one can see in Fig. 1 that it has two exoergic product channels that may compete: C(³P) + PN(¹Σ⁺) and P(⁴S) + CN(²Σ⁺). The abstraction of a carbon atom by the incoming nitrogen atom to produce P(⁴S) + CN(²Σ⁺) can easily take place through a direct mechanism that involves the formation of the global minimum of this triatomic system. This N + CP ⟶ NCP ⟶ P(⁴S) + CN(²Σ⁺) pathway is barrierless and should be open regardless of the collision energy. The abstraction of the P atom by the incoming N atom to yield C + PN, in turn, presents a small energy barrier in the entrance channel. This barrier of just 1 kJ mol⁻¹ associated with the N + CP ⟶ CPN ⟶ C + PN reaction path lies below the accuracy of our quantum chemical calculations, but this reaction can also occur without an activation energy through the N + CP ⟶ NCP ⟶ CNP ⟶ C + PN mechanism.

Fig. 1 Potential energy diagram for the 3A′ (red) and 3A′′ (black) electronic states of the CPN system. Energies are given in kJ mol−1 relative to the C+PN limit and are ZPE corrected. The MRCI-F12 energies are given in boldface, CCSD(T)-F12 results are given in italics and M06-2X ones are given under parenthesis. Atoms are coloured as: carbon (grey); phosphorus (orange); nitrogen (blue).

3.3 Reaction rates

Among the reactions addressed in this work, only processes P + CN ⟶ C + PN and N + CP ⟶ C + PN have rate constants available in the UMIST Database for Astrochemistry (UDfA) (McElroy et al. 2013), both having the same value of k = 3 × 10⁻¹⁰ cm³ s⁻¹ and none of them are currently present in KIDA (Wakelam et al. 2012). However, the former reaction is shown here to be endothermic, and the present value available in UDfA may lead to wrong conclusions on astrochemical models. As discussed in Section 3.1, this was estimated based on the analogous nitrogen reaction (N + CN → N₂ + C) assuming similar chemical behaviour between the involved species. However, the PN molecule is not as thermodynamically stable as N₂, and we have shown that these chemical processes are actually quite different and present very distinct energy landscapes. The P + CN ⟶ C + PN reaction for instance, is significantly endoergic and should not present such a high rate coefficient, nor should it be independent of temperature. As for the N + CP ⟶ C + PN reaction, our calculations predict the existence of a small entrance barrier for the most straightforward mechanism, contrary to what occurs for the analogous nitrogen reaction (Moskaleva & Lin 2001). Therefore, despite being exoergic, one should expect the N + CP ⟶ C + PN reaction to present a lower rate constant than the currently assigned one.

We have used our accurate computations at the MRCI-F12/AV(T+d)Z//CAS/AV(T+d)Z + ZPE(CAS/AV(T+d)Z) level on the CPN system to calculate the rate coefficients of reactions in Eqs. (1), (3), (4) and their reverses as a function of temperature using the MESS package. The rate coefficients obtained for all six reactions for several different temperature values have been fitted to the modified Arrhenius equation: $$k(T)=\alpha (T/300{)}^{\beta }\exp (-\gamma /T).$$(9)

and the fitting parameters (α, β and γ constants) are compiled in Table 3. One can also find in this table results from previous studies of our research group as well as the ones proposed by Douglas et al. (2022) and García de la Concepción et al. (2021). In this way, we offer a wide set of reaction rates for modelling the abundances of P-bearing molecules. We also refer the reader to the supplementary information file where one can find more details about the determination of the rate coefficients of the PO + N ⟶ PN + O and PO + N ⟶ NO + P reactions.

The calculated rate coefficients with the respective fitted curves for the exoergic reactions are given in Fig. 2. One can notice that these reactions show a very small temperature dependence, as is usually the case for barrierless exoergic processes. As expected, our predictions that the branching ratio between N + CP ⟶ P + CN and N + CP ⟶ C + PN largely favours the former, although only the latter is available in UDfA. For typical ISM conditions, the difference is about two orders of magnitude.

The obtained rate coefficients for the reaction in Eq. (1), which may have an important impact on the abundance of PN in the ISM as the only reaction able to remove this species through neutral-neutral collision, are shown as blue circles in Fig. 2. They are essentially constant over a wide range of temperatures with a mean value of k = 1.35 × 10⁻¹⁰ cm³ s⁻¹, which is consistent with what is expected from a barrierless exoergic reaction. Nonetheless, a rather unusual increase in the reaction rate coefficients as the temperature decreases can be observed in the low temperature regime (<200 K). This can be explained due to the electronic degeneracy of the carbon atom in the ³P state when we consider the splitting of the ground energy level due to spin-orbit coupling. This interaction gives rise to three finely separated energy levels with total angular momentum quantum number J = 0, 1 and 2 and electronic degeneracy number g = 1, 3 and 5, respectively. As temperature decreases, the states with lower degeneracy (and also lower energy) become more populated, increasing the values of the rate coefficients due to the higher statistical weight factors associated with these states.

The reactions in Eqs. (1) and – (1) have also been approached by Douglas et al. (2022), who essentially set the rate coefficient of the direct process equal to the capture rate calculated using long-range transition state theory (Georgievskii & Klippenstein 2005). The rate constant of the reverse process was then calculated by detailed balance. The authors justify this procedure by considering that the reaction in Eq. (1) does not have a complex PES and they seem to treat only the CNP intermediate as relevant to those reactions, neglecting other possible multistep pathways. It should still be a good approximation and, therefore, they are also shown in Figs. 2 and 3 as dashed lines for comparison. There is clearly good agreement between our results and those obtained by Douglas et al. (2022) for these reactions.

The results for the endoergic reactions are given in Fig. 3 in an Arrhenius plot. The activation energies associated with these processes are related to the slope of the curves. As expected, the rate coefficients of these reactions (Eqs. – (1), (3) and – (4)) present a strong temperature dependence. By analysing the results shown in Fig. 3, it becomes evident that the value assigned in UDfA to the rate constant of the P + CN → C + PN reaction differs drastically to the calculated ones. According to our calculations, only at 1000 K the rate coefficient for the P + CN → C + PN process reaches k = 1.0 × 10⁻¹⁷ cm³ s⁻¹, seven orders of magnitude lower than the temperature independent value based on the analogous N reaction (Baulch et al. 1992; Agúndez et al. 2007). While high-temperature values are relevant for studies of exoplanetary atmospheres (Al-Refaie et al. 2022; Fleury et al. 2025), our focus here is on the relatively cold regions of the ISM. The highest rate coefficient in Fig. 3 is the one of the reaction in Eq. (3), and it only reaches the order of magnitude of 10⁻¹⁰ cm³ s⁻¹ at 7000 K. Hence, we should not expect these endoergic reactions to play a relevant role in the chemical evolution of P-bearing species in cold astronomical environments.

Fig. 2 Rate coefficients as a function of temperature obtained for the exoergic reactions. The points refer to the calculated values at fixed temperatures and the curves correspond to a fit to our results using the modified Arrhenius formula. The rate coefficient proposed by Douglas et al. (2022) for the C + PN ⟶ P + CN reaction is shown as a dashed line for comparison.

Fig. 3 Rate coefficients as a function of inverse temperature obtained for the endoergic reactions. The points refer to the calculated values at fixed temperatures and the curves correspond to a fit to our results using the modified Arrhenius formula. The rate coefficient proposed by Douglas et al. (2022) for the P + CN ⟶ C + PN reaction is shown as a dashed line for comparison.

4 Astrochemical modelling of P-bearing species

Evaluating the behaviour of computed reactions in astrochemical models is essential for understanding their impact on the abundances of phosphorus-bearing species. The results are obtained from gas-grain modelling using the Nautilus code, which calculates the initial abundances provided in Table 2 and incorporates the chemical reactions listed in Table 3 into the KIDA network database (see Sect. 2.3). The primary objective is to analyse the temporal evolution of these reactions and predict the abundances of P, PO, PN, and PH (phosphinidene) under conditions that approximate the transition from diffuse to dense cloud environments.

4.1 Diffuse and translucent cloud

In the life cycle of the ISM, diffuse molecular clouds are low density regions primarily composed of molecular hydrogen (H₂) and dust, within which translucent clouds are embedded. These clouds play a crucial role in the early stages of star formation, as they provide the necessary environment for the collapse of matter under gravitational forces, eventually leading to the birth of stars. To investigate the formation of P-bearing species in a diffuse molecular environment, we modelled the system under physical conditions with a gas density of 10³ cm⁻³, a standard gas-to-dust ratio mass ratio of 100:1 (Tricco et al. 2017), and a temperature of 70 K (Corby et al. 2018). The initial elemental composition of the gas, including both atoms and ions, is presented in Table 2.

By examining the PO and PN abundances under varying cosmic ray and X-ray ionisation rates, we observed corresponding changes in their ratios. The results are exhibited in Fig. 4, which shows the resultant PO/PN ratios as a function of time for the different models: a), b) and c). It is observed that early in the simulation (t < 10² yr), the impact of the ionisation rates on the ratios is less pronounced. However, as time evolves, the differences in PO/PN ratios among the models become significantly more pronounced. In the final stages (t > 10⁵ yr), the models show PO/PN ratios constrained within an interval of roughly 15 to 60 as the cosmic ray ionisation rates increase from 1.3 × 10⁻¹⁷ s⁻¹ to 1.3 × 10⁻¹⁶ s⁻¹. Higher ionisation rates and energetic phenomena might affect PO chemistry, influencing the abundance of PO relative to PN. This result is consistent with previous studies that combined quantum chemical and kinetic calculations of PO, also evaluating velocity effects within the context of chemical models in shocked regions (García de la Concepción et al. 2021 and references therein). In comparison to other cases, high PO/PN ratios have been observed in comet 67P/Churyumov-Gerasimenko, Rivilla et al. (2020) estimated that PO is significantly more abundant than PN, with PO/PN ratios exceeding 10.

We calculated the abundances of phosphorus-bearing species as a function of time, maintaining the gas temperature and density conditions of the diffuse and translucent cloud described previously. To account for the significance of ionisation rates, we utilised model B to compute the abundances, since it represents an intermediate case with ζ_CR = 6.3 × 10⁻¹⁷ s⁻¹ and ζ_X = 5 × 10⁻¹⁸ s⁻¹. The results are shown in Fig. 5, which displays the abundances of P, PO, PN, and the radical PH, over a timescale of 10⁶ yr. The analysis yields key findings, such as the elemental phosphorus remains the most abundant species throughout the simulation. Furthermore, PO consistently exhibits higher abundances than PN. Notably, the predicted abundances of P, PO, PN, and PH (phosphinidene) reach a quasisteady state after approximately 10⁵ years–significantly earlier than the typical chemical age of 10⁶ years, adopted for diffuse molecular clouds (e.g. Albertsson et al. 2014). In this work, we include the results for PH, represented by a dashed line in Fig. 5. Compounds such as phosphinidenes are unstable and highly reactive, presenting a significant challenge for detection. Nevertheless, they play a crucial role in phosphorus chemistry.

Specifically, PH contributes to the formation of both PO and PN through neutral reactions: PH + O ⟶ PO + H and PH + N ⟶ PN + H (Table 3).

Fig. 4 Ratios of PO/PN as a function of time in the diffuse and translucent cloud model: a) ζCR = 1.3 × 10−16 s−1 and ζX = 1 × 10−17 s−1; b) ζCR = 6.3 × 10−17 s−1 and ζX = 5 × 10−18 s−1; and c) ζCR = 1.3 × 10−17 s−1 and ζX = 1 × 10−18 s−1. The dashed line signifies the threshold at which PO/PN equals 1.

Fig. 5 Evolution of the abundances of P, PO, PN and the radical PH (dashed line) in the diffuse and translucent cloud model, adopting a H2 gas density and temperature of 103 cm−3 and 70 K, respectively, with ζCR = 6.3 × 10−17 s−1 and ζX = 5 × 10−18 s−1.

4.2 Dense molecular cloud

A second model of phosphorus chemistry was computed under dense cloud conditions. In this context, the final abundances from the diffuse and translucent model become the initial abundances for the dense cloud. Additionally, considering the radiative processes related to gas cooling, we computed the abundances based on a H₂ density of 10⁴ cm⁻³ and a temperature of 30 K. It is important to note that molecular clouds provide effective shielding against ultraviolet and X-ray photons. Consequently, for the dense cloud simulation we adopted A_V = 10 mag and utilised the lowest ionisation rates tested in this study, as discussed in Model C, ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹ and ζ_X = 1 × 10⁻¹⁸ s⁻¹. Previous research has highlighted the complexities involved in accurately simulating the transitions from low-density, ionised ISM gas to dense, neutral cloud environments (e.g. Phan et al. 2018). Thus, we aim to investigate the phosphorus network with a focus on the following objectives: examining the evolution of abundances, identifying key reaction contributions, analysing implications for steady-state conditions, and the dynamic chemical balance over timescales ranging from 1 × 10⁶ to 1 × 10⁸ yr.

The results of this model are presented in Fig. 6, with the principal contributions of key reactions analysed in Fig. 7. The temporal evolution of the abundances of P, PO, PN, and PH are presented analogously to Fig. 5; however, the analysis is conducted over an extended timescale ranging from 10⁶ to 10⁸ years. This hypothetical prolonged period aims to discuss abundance stabilisation and the dynamics of the chemical network. As shown in Fig. 6, the abundance curves of PO and PN exhibit similar trends, with values fluctuating between ∼ 5 × 10⁻¹¹ and 5 × 10⁻⁸. Except during the initial phase, the model shows that PO maintains a higher abundance than PN throughout the evolution. Consequently, in the final stage, the PO abundance reaches approximately 3.3 × 10⁻¹⁰, whereas the PN abundance, closely resembling that of PH, stabilises around 1.3 × 10⁻¹⁰. Figure 6 also includes the abundance curve of PH, which appears decoupled from the behaviour of PO and PN. The formation of PH is primarily driven by the recombination of molecular ions with electrons, as in reactions like HPO⁺ + e⁻ ⟶ O + PH and HPN⁺ + e⁻ ⟶ N + PH. In contrast, its destruction pathways contribute to the synthesis of PO and PN, as seen in the diffuse and translucent model.

We have identified key mechanisms influencing the chemistry of PN and PO. In Fig. 7(a) and (b), we present the percentage contributions of the principal formation and destruction reactions for PN and PO, respectively, as a function of time. In the case of PN (Fig. 7a), it is important to note that its primary formation pathway is attributed to the PO + N reaction, which initially contributes slightly less than 100% and gradually decreases to approximately 50%. Furthermore, the reaction of PH with N contributes to PN formation (e.g. Gomes et al. 2023a), though to a lesser degree. Initially, this pathway contributes for only a minor fraction of the total PN production but gradually increases, reaching approximately 50% in the final stage of the model. Regarding the destruction mechanisms, PN is primarily destroyed through two pathways involving reactions with H⁺ and C, with the latter being one of the major reactions explored in this work. The curves associated with these mechanisms exhibit a degree of symmetry, which may arise from the complex interactions inherent in the computations performed by the code. This symmetry indicates a certain predictability and alternation; for instance, at t ≈ 1.5 × 10⁷ yr it can be observed that the contribution of H⁺ to PN destruction diminishes as the contribution from C becomes more significant.

In the analysis of major reactions influencing PO abundance (Fig. 7b), the primary formation mechanism is the reaction between P and OH (e.g. García de la Concepción et al. 2021), which initially accounts for ∼ 88% of the PO production. However, as the chemical network evolves, its contribution gradually decreases towards the final stage of the model. In contrast, the PH + O reaction follows an opposite trends, since it begins with a minimal contribution but steadily increases. As a result, the contribution profiles of these two pathways shift over time, intersecting at approximately t = 6 × 10⁶ yr. To a lesser extent, the reaction HPO⁺ + e⁻ also contributes to the PO production, though their impact is significantly lower than that of the P + OH and PO + H pathways. Regarding the major destruction pathways of PO, two key reactions contribute to its depletion: PO + N and PO + H⁺. The relative significance of these mechanisms evolves over time. Initially, the PO + N reaction is the dominant process; however, its contribution gradually declines, allowing the PO + H⁺ reaction to become the primary destruction pathway in the later stages.

Our findings indicate that the computed kinetic equations and percentage contributions reveal a pattern that is both somewhat predictable and alternates in behaviour. The formation routes (P + OH and PH + O) appear to establish a dynamic equilibrium with the destruction mechanisms (PO + N and PO + H⁺), suggesting that the steady-state abundance of PO is regulated by the competition between these pathways, which in turn has an effect on the chemistry and abundance of PN. In future work, we will explore this equilibrium in greater detail, assessing its implications for the predictability of chemical abundances under interstellar conditions.

Fig. 6 Similar to Fig. 5, depicting a stage from 106 yr to 108 yr for a dense cloud model. The parameters are set to a temperature of 30 K, a H2 density of 104 cm−3 and adopting ζCR = 1.3 × 10−17 s−1 and ζX = 1 × 10−18 s−1.

Fig. 7 Major reactions for a) PN and b) PO. The contributions of each reaction are represented as percentages over time. Solid and dashed lines differentiate production and destruction mechanisms, respectively. The time is presented on a logarithmic scale to highlight contributions beyond 106 years, with a horizontal grey line at 50% for assessing balance throughout the duration of the long-term trends in the chemical network.

5 PO/PN ratios and observational implications

PO and PN have been identified across a range of astronomical environments, including star-forming regions, protostellar objects, shock regions, and giant molecular clouds. Generally, these molecules have been observed in the interstellar medium and within the circumstellar envelopes of evolved stars. Among the two species, PO is widely recognised as being both more abundant and more prevalent than PN. Studies in the literature consistently indicate that PO is typically more abundant than PN by factors of approximately 1.4 to 3, regardless of the specific physical characteristics of the observed sources (Fernández-Ruz et al. 2023). In modelling a diffuse and translucent cloud, we consistently found that PO is more abundant than PN, with both species displaying a dependence on factors such as visual extinction and influenced by factors like X-ray and cosmic ray ionisation rates. Specifically, these findings suggest that the PO/PN ratios vary over time and increase as the cosmic ray ionisation rate rises from ζ_CR = 1.3 × 10⁻¹⁷ to 1.3 × 10⁻¹⁶ s⁻¹. However, under conditions of low ionisation rates, the PO/PN ratio may drop below 1. This fluctuation highlights the sensitivity of phosphorus chemistry to ionisation mechanisms, emphasising its critical role in modelling the relative abundances of PO and PN over time.

In a study of diffuse and translucent clouds, Chantzos et al. (2020) conducted observations along the line of sight towards B0355+508, focusing on the analysis of P-bearing species. For PO and PN, the derived upper abundance limits were determined to be <5.0 × 10⁻¹⁰ and <4.9 × 10⁻¹¹, respectively. Based on our theoretical predictions (Fig. 5), we estimate PO and PN abundances reaching ∼ 1 × 10⁻¹⁰ and 2 × 10⁻¹², respectively, which are consistent with these upper limits. In a study of P-bearing species on a galactic scale, Koelemay et al. (2023) determined the abundances of PO and PN in a star-forming region in the outer Galaxy to be ∼ 2.0 × 10⁻¹¹ and 3.0 × 10⁻¹², respectively. They also estimated a fit for the PO and PN abundances as a function of distance from the Galactic Centre, suggesting a slight gradient–a decrease in abundance by factors of approximately 1.5 and 2.3 for PO and PN, respectively, between 8.5 and 22.6 kpc.

In a model simulating the conditions of a denser cloud (n ≃ 10⁴ cm⁻³), adopting the standard cosmic ray ionisation rate (ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹), we predicted the abundances of P, PO, PN, and PH, along with the major reaction pathways involving these species. In this model, the final abundances of PO and PN were estimated to be ∼ 3 × 10⁻¹⁰ and 1 × 10⁻¹⁰, respectively. Through a detailed examination of the chemical reactions involving both species, we identified that the primary formation pathway of PN occurs via the reaction of PO with N. This finding is further supported by recognising that the primary mechanism for PO destruction is its reaction with N, which leads to PN formation (see Fig. 7). In an early work, Thorne et al. (1984), using laboratory results and kinetic modelling, highlighted the significance of PO as a pivotal species in phosphorus chemistry, serving as a primary reservoir of phosphorus with the potential for other species, such as PN and CP. Consequently, PO is a molecule that, once synthesised, accumulates in sufficient quantities to support the formation of additional species, such as PN, without being considerably depleted.

Bernal et al. (2021) reported relative abundances of PO ∼ 1.6 × 10⁻¹⁰ and PN ∼ 6.1 × 10⁻¹¹ in the Orion Plateau, yielding a PO/PN ratio of approximately 3. They discussed the presence of PO in high-mass star-forming regions, with the potential for its occurrence in low-mass sources as well. In the Orion-KL region, the clear association between PO and PN suggests that these two molecules are likely formed together through interconnected reactions, as it has been found in the present work. In the context of solar-type protostellar objects, the shock region L1157-B1 is characterised by active and complex chemistry. The detection and formation of various molecules, including HCN, HNC, SiS, HC₃N, PO and PN, have been investigated in detail, with a particular focus on shock-induced gas-phase chemistry (Lefloch et al. 2016; Mendoza et al. 2018; Mota et al. 2021; Lefloch et al. 2021). In a recent study, Lefloch et al. (2024) reported that the PO/PN abundance ratios range from 1 to 5. They argued that the PO/PN ratios observed throughout B1 indicate that the phosphorous chemistry is time-dependent. It has also been suggested that neither species is sputtered from icy mantles; rather, PO forms in the gas phase, followed by the subsequent synthesis of PN (see also Lefloch et al. 2016). García de la Concepción et al. (2021) carried out quantum chemical computations to determine the reaction rates and branching ratios for the reactions P + OH → PO + H and P + H₂O → PO + H₂. Additionally, they conducted chemical modelling in a shocked region, predicting PO/PN ratios exceeding 1, ranging from 3 to 9. Wurmser & Bergner (2022) conducted a search for PN and PO lines towards a sample of solar-type protostars with well-characterised outflows. They identified a relatively narrow range of PO/PN ratios, extending from 0.6 to 2.2. Their analysis highlighted the critical role of outflows and shocked gas in driving gas-phase phosphorus chemistry in these types of sources. Consistent with findings in the literature, gas-phase phosphorus reactions account for the observed abundances of species such as PO and PN. However, follow-up studies are essential to investigate the chemical dynamics between PO and PN in greater detail, as these species play a fundamental role in complex processes related to prebiotic chemistry and astrobiology.

6 Conclusions

Current astrochemical models lack important PN destruction routes (Fernández-Ruz et al. 2023). Furthermore, reactions such as N + CP → PN + C and P + CN → PN + C have very uncertain rate constants available in the literature (Jiménez-Serra et al. 2018), and the rate coefficient of the nitrogen analogous reaction (N + CN → N₂ + C) is often used to describe the chemical evolution of these P-bearing molecules in the interstellar medium (Fernández-Ruz et al. 2023). In this work we tackled these issues by performing high-accuracy theoretical calculations (MRCI-F12/AV(T+d)Z//CAS/AV(T+d)Z and CCSD(T)-F12/VTZ-F12//M06-2X/AV(T+d)Z) and employing transition state theory to study the potential energy landscape of the CPN system and to elucidate the kinetics of the C(³P) + PN(¹Σ⁺) destruction of phosphorus nitride. The outcomes of collisions involving N(⁴S) + CP(²Σ⁺) and P(⁴S) + CN(²Σ⁺) species were also investigated. We have then integrated a revised chemical network for P-bearing molecules into the default KIDA network and ran the Nautilus gas-grain Astrochemical code to investigate the temporal evolution of these reactions within diffuse and dense interstellar cloud conditions.

We have found several mechanisms for PN destruction by collisions with a carbon atom that can lead to P(⁴S) + CN(²Σ⁺) formation without an activation barrier. This is the only product channel that is open for the C(³P) + PN(¹Σ⁺) reaction for whatever collision energy. Regarding the kinetics, the temperature-dependent rate coefficients for the reactions addressed in this work were fitted to the modified Arrhenius equation: k(T) = α(T/300)^β exp (−γ/T). The C(³P) + PN(¹Σ⁺) → P(⁴S) + CN(²Σ⁺) reaction has a rate coefficient with parameters α = 1.29 × 10⁻¹⁰ cm³ s⁻¹, β = 0.025 and γ = −10.0 K, consisting in a relevant PN destruction route that should be important to properly model the abundances of P-bearing species in the ISM.

The N(⁴S) + CP(²Σ⁺) formation processes were found to be considerably endoergic and mechanisms towards CP formation should not take place unless very high temperatures are involved. The destruction of CP, however, is expected to readily occur by collisions with nitrogen atoms even at low temperatures. The obtained rate coefficient for the N(⁴S) + CP(²Σ⁺) → P(⁴S) + CN(²Σ⁺) reaction has α = 1.74 × 10⁻¹⁰ cm³ s⁻¹, β = 0.172 and γ = 0 and is higher than the one obtained for the N(⁴S) + CP(²Σ⁺) → C(³P) + PN(¹Σ⁺), which has α = 9.86 × 10⁻¹² cm³ s⁻¹, β = 0.634 and γ = 0. To the best of our knowledge, both rate constants are being reported for the first time. We have also concluded that the energy landscape of the CPN system is very different from that of NCN, and the use of the rate coefficient of the nitrogen analogous reaction is not adequate to describe the evolution of reactions within the CPN system.

We have presented two chemical models to predict the abundances of phosphorus-bearing species, specifically P, PO, PN, and PH, under varying astrophysical conditions. In the diffuse cloud scenario, the influence of ionising agents on phosphorus chemistry was investigated through PO/PN abundance ratios. Our results indicate that models with higher X-ray and cosmic ray ionisation rates strongly favour PO formation, leading to PO/PN ratios significantly greater than 1. In contrast, models with lower ionisation rates can result in slightly higher PN abundances, with PO/PN approaching or slightly below 1 around t = 10³ yr. In the final phase of the model (t ≈ 1 × 10⁶ yr), PO/PN range from 15 to 60 depending on the ionisation rates. A diffuse cloud model was computed for an intermediate ionisation rate scenario, assuming ζ_CR = 6.3 × 10⁻¹⁷ s⁻¹ and ζ_X = 5 × 10⁻¹⁸ s⁻¹. As a result, quasi-stable abundances of PO and PN were obtained for t > 10⁵ yr, with values of ≈ 1 × 10⁻¹⁰ and ≈ 2 × 10⁻¹², respectively.

For a subsequent second model adopting dense cloud conditions, the chemical network was computed over an extended hypothetical time-scale ranging from 10⁶ to 10⁸ yr, with ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹ and ζ_X = 1 × 10⁻¹⁸ s⁻¹. The results indicate that PO and PN abundances slightly fluctuate over time, with values ranging from 5 × 10⁻¹¹ to 5 × 10⁻⁸. The balance of the chemical network and its influence on steady-state abundances were analysed by examining the temporal percentage contributions of key reactions. PO is consistently the most abundant species compared to PN and PH under these conditions. A detailed analysis of the key formation and destruction pathways revealed that PO is primary produced via the P + OH and PH + O reactions, whose relative contributions vary over time. With a lower contribution (up to 20%), the HPO⁺ + e⁻ reaction also plays a role in the PO formation. Regarding PO destruction mechanisms, the most significant channels are PO + N and PO + H⁺. An analysis of the major formation and destruction pathways of PN was carried out. The primary formation of PN occurs via the reaction of PO with N, contributing more than 50%. This underscores the critical role of PO as a reservoir species that accumulates sufficiently to drive subsequent PN formation. Nevertheless, PO remains as the most abundant P-bearing molecule. The second most important mechanism, the PH + N reaction, contributes to the remaining percentage. Finally, regarding PN destruction, we find that its depletion is primarily driven by reactions with H⁺ and C. These findings align with earlier studies, confirming the interplay of phosphorus chemistry under varying ionisation environments.

In future work, we intend to further enhance the accuracy of the obtained rate coefficients for the approached chemical reactions by performing quasi-classical trajectory simulations. It involves, however, the calculation of a huge set of ab initio energies for sufficiently many nuclear arrangements so that an analytic expression can be fitted to such data and used to model the interaction potential. Despite being a very timeconsuming process, a reliable analytical model offers the advantage of providing fast, continuous and differentiable representations of the PES (Rocha & Varandas 2020). A robust scheme, called the Combined-Hyperbolic-Inverse-Power-Representation (CHIPR) method, has been proposed by Varandas (Varandas 2013a,b), and the code of a general program for analytically representing PESs of any triatomic molecules with CHIPR has already been provided (Rocha & Varandas 2020). This method should be used as the next step in our study of the CPN system.

Data availability

The data underlying this article are available in the article and in its online supplementary material at https://doi.org/10.5281/zenodo.15002910

Acknowledgements

We would like to thank the anonymous reviewer for their insightful comments. 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, grants 311508-2021-9 and 405524-2021-8), Fundação de Amparo à Pesquisa e Inovação do Espírito Santo (FAPES), Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG, grants APQ-03705-23 and RED-00045-23), Centro Federal de Educação Tecnológica de Minas Gerais (CEFET-MG) and Universidade Federal do Espírito Santo (UFES). E.M. acknowledges support under the grant “María Zambrano” from the UHU funded by the Spanish Ministry of Universities and the “European Union NextGenerationEU.” This project has also received funding from the European Union’s Horizon 2020 research and innovation program under Marie Sklodowska-Curie grant agreement No. 872081, and grant PID2022-136228NBC21 (M.C.) funded by MCIN/AEI/10.13039/501100011033, and, as appropriate, by “ERDF A way of making Europe”, the “European Union”, or the “European Union NextGenerationEU/PRTR”. This work is also supported by the Consejería de Transformación Económica, Industria, Conocimiento y Universidades, Junta de Andalucía and European Regional Development Fund (ERDF 2021-2027) under the project EPIT1462023 (M.C. and E.M.).

References

All Tables

All Figures

	Fig. 1 Potential energy diagram for the 3A′ (red) and 3A′′ (black) electronic states of the CPN system. Energies are given in kJ mol−1 relative to the C+PN limit and are ZPE corrected. The MRCI-F12 energies are given in boldface, CCSD(T)-F12 results are given in italics and M06-2X ones are given under parenthesis. Atoms are coloured as: carbon (grey); phosphorus (orange); nitrogen (blue).
In the text

	Fig. 2 Rate coefficients as a function of temperature obtained for the exoergic reactions. The points refer to the calculated values at fixed temperatures and the curves correspond to a fit to our results using the modified Arrhenius formula. The rate coefficient proposed by Douglas et al. (2022) for the C + PN ⟶ P + CN reaction is shown as a dashed line for comparison.
In the text

	Fig. 3 Rate coefficients as a function of inverse temperature obtained for the endoergic reactions. The points refer to the calculated values at fixed temperatures and the curves correspond to a fit to our results using the modified Arrhenius formula. The rate coefficient proposed by Douglas et al. (2022) for the P + CN ⟶ C + PN reaction is shown as a dashed line for comparison.
In the text

	Fig. 4 Ratios of PO/PN as a function of time in the diffuse and translucent cloud model: a) ζCR = 1.3 × 10−16 s−1 and ζX = 1 × 10−17 s−1; b) ζCR = 6.3 × 10−17 s−1 and ζX = 5 × 10−18 s−1; and c) ζCR = 1.3 × 10−17 s−1 and ζX = 1 × 10−18 s−1. The dashed line signifies the threshold at which PO/PN equals 1.
In the text

	Fig. 5 Evolution of the abundances of P, PO, PN and the radical PH (dashed line) in the diffuse and translucent cloud model, adopting a H2 gas density and temperature of 103 cm−3 and 70 K, respectively, with ζCR = 6.3 × 10−17 s−1 and ζX = 5 × 10−18 s−1.
In the text

	Fig. 6 Similar to Fig. 5, depicting a stage from 106 yr to 108 yr for a dense cloud model. The parameters are set to a temperature of 30 K, a H2 density of 104 cm−3 and adopting ζCR = 1.3 × 10−17 s−1 and ζX = 1 × 10−18 s−1.
In the text

Fig. 7 Major reactions for a) PN and b) PO. The contributions of each reaction are represented as percentages over time. Solid and dashed lines differentiate production and destruction mechanisms, respectively. The time is presented on a logarithmic scale to highlight contributions beyond 106 years, with a horizontal grey line at 50% for assessing balance throughout the duration of the long-term trends in the chemical network.

In the text