Reaction dynamics of P(⁴S) + O₂(X³Σ⁻_g)  →  O(³P) + PO(X²Π) on a global CHIPR potential energy surface of PO₂(X²A₁): implications for atmospheric modelling

The reaction dynamics of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) are thought to be important in atmospheric and interstellar chemistry. Based on the state-of-the-art ab initio energy points, we analytically constructed a global potential energy surface (PES) for the ground-state PO₂(X²A₁) using the combined-hyperbolic-inverse-power-representation (CHIPR) method. A total of 6471 energy points were computed by the multireference configuration interaction method with the Davidson correction and aug-cc-pV5Z basis set. The analytical CHIPR PES reproduces ab initio energies accurately with a root-mean-square deviation of 91.5 cm⁻¹ (or 0.262 kcal mol⁻¹). The strongly bound valence region of the PES has complicated topographical features with multiple potential wells and barriers. The attributes of the important intermediates are carefully validated with our geometry optimization results, as well as previous experimental and computational results. Finally, the reaction probability, integral cross sections, and rate constants for P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) are calculated using the quasi-classical trajectory and time-dependent wave packet methods. The trends of probability and integral cross section versus the collision energy can be divided into three stages, which are governed by the entrance barriers or exothermicity of the reaction. The rate constant demonstrates strong Arrhenius linear behaviour at relatively low temperatures but deviates from this pattern at high temperatures. The calculated cross sections and rate constants are helpful for modelling the phosphorus chemistry in atmospheric and interstellar media.

1 Introduction

Phosphorus (P) is one of the essential biogenic elements found in all known organisms. Its compounds are abundant in biological systems and greatly contribute to basic biological functions (Maciá, 2005). Much of the P on Earth's surface is locked up in mineral phosphate forms with fairly poor bioavailability, while the P compounds with low oxidation are generally more reactive and accessible for potential prebiotic chemistry (Todd, 2022) but unstable under terrestrial redox conditions (Pasek, 2008).

One source of P compounds with low oxidation on Earth's surface is interstellar dust particles (IDPs), which account for 99 % of the total amount of incoming extra-terrestrial material each year (Plane et al., 2018), containing P with about 8 % abundance (Lodders, 2003). The ablation of IDPs in the upper atmosphere of terrestrial planets delivers PO and P (Carrillo-Sánchez et al., 2020), a part of which might undergo a series of chemical processes to form bioavailable H₃PO₃ before reaching the Earth's surface. The corresponding P chemistry networks (Douglas et al., 2019, 2020; Plane et al., 2021) were predicted by electronic structure theory and given by

$$\begin{array}{}\text{(R1)}& {\displaystyle }& {\displaystyle }\mathrm{P}\left({}^{\mathrm{4}}\mathrm{S}\right)+{\mathrm{O}}_{\mathrm{2}}\to \mathrm{O}+\mathrm{PO},\text{(R2)}& {\displaystyle }& {\displaystyle }\mathrm{PO}+{\mathrm{O}}_{\mathrm{2}}\to \mathrm{O}+{\mathrm{PO}}_{\mathrm{2}},\text{(R3)}& {\displaystyle }& {\displaystyle }{\mathrm{PO}}_{\mathrm{2}}+\mathrm{H}(+M)\to {\mathrm{HPO}}_{\mathrm{2}},\text{(R4)}& {\displaystyle }& {\displaystyle }{\mathrm{HPO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}(+M)\to {\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{3}}.\end{array}$$

This suggests that meteor-ablated P is likely oxidized by the Reaction (R1) first. Hence, the dynamics study of Reaction (R1) will contribute to a deeper understanding of the chemical evolution of P and PO in the Earth's atmosphere.

Reaction R1 is also important in astrochemistry. For example, PO has been widely observed in the interstellar medium (ISM) (Tenenbaum et al., 2007; De Beck et al., 2013; Ziurys et al., 2018; Lefloch et al., 2016; Rivilla et al., 2020) and is considered to be the main reservoir for gas-phase P in the ISM (Ziurys et al., 2018; Rivilla et al., 2020). There is some evidence to suggest that PO was present in cometary ice before the birth of the Sun (Rivilla et al., 2020), and the comets possibly provided a major source of prebiotic phosphorus compounds (Maciá et al., 1997). Hence, investigation of the formation of PO is helpful for modelling its abundance in the ISM.

Given its important role in interstellar and atmospheric chemistry, Reaction (R1) was investigated in several experimental and theoretical studies (Douglas et al., 2019; Husain and Norris, 1977; Husain and Slater, 1978; Clyne and Ono, 1982; Henshaw et al., 1987; Gomes et al., 2022). Previous experiments presented its rate constants at a specified temperature of about 300 K (Husain and Norris, 1977; Husain and Slater, 1978; Clyne and Ono, 1982; Henshaw et al., 1987), and their values diverged by an order of magnitude. A recent experiment has determined the rate constants of P(⁴S) + O₂ → O + PO at temperatures ranging from ∼200 to 750 K (Douglas et al., 2019). In the same work, the rate constants of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) were also computed using the Rice–Ramsperger–Kassel–Markus (RRKM) theory with the molecular geometry optimization at the B3LYP/aug-cc-pVQZ (AVQZ) level. Subsequently, the state-of-the-art ab initio method was used to predict the minimum energy path (MEP) of this reaction (Gomes et al., 2022), and the rate constants were computed by employing the standard transition-state theory (TST). These two theoretical rate constants reproduced the experimental results to some extent, and the difference between them came mainly from the predicted barrier heights. The limitation of the work by Douglas et al. (2019) is that the B3LYP method is not good at predicting the barrier heights (Zhao and Truhlar, 2008; Peverati and Truhlar, 2012). Also, the statistical RRKM and TST theories may fail due to unincluded non-statistical dynamics effects (Carpenter, 1998; Thomas et al., 2008), so a dynamics study for this reaction is required. Moreover, their analytical forms of the rate constants were obtained by fitting predictive data lower than 1000 K, so the fitted rate constants at high temperatures may not be accurate. Accurately modelling the rate constants of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) in a wide temperature range is desired, because the temperatures of ablation of IDPs could reach more than 2500 K.

The reactants P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and products O(³P) + PO(X²Π) are connected to the lowest doublet and quartet states of PO₂, in which the doublet state (X²A₁) plays a major role in this reaction (Gomes et al., 2022). The molecular geometries and vibration frequencies of PO₂(X²A₁) have been well studied by several theoretical (Lohr, 1984; Kabbadj and Liévin, 1989; Jarrett-Sprague and Hillier, 1990; Cai et al., 1996; Francisco, 2002; Xianyi et al., 2008; Liang et al., 2013; Xu et al., 1996; Bauschlicher, 1999) and experimental (Cordes and Witschel, 1965; Davies and Thrush, 1968; Drowart et al., 1972; Verma and McCarthy, 1983; Kawaguchi et al., 1985; Hamilton, 1987; Qian et al., 1995; Lei et al., 2001; Lawson et al., 2011) methods. The recent measurement reported R_PO=2.771 a₀ and θ_OPO=135.3^∘ with the C_2v symmetry (Kawaguchi et al., 1985), along with the vibration frequencies of 1075.4, 397.3 and 1327.54 cm⁻¹ for the symmetrical stretching (ω₁) (Lei et al., 2001), bending (ω₂) (Lei et al., 2001), and antisymmetric stretching (ω₃) (Lawson et al., 2011), respectively. It is worth noting that the potential energy surface (PES) can yield physical insights into the reaction path, energy transfer, and structure of intermediates. The analytical form of a global PES modelled by ab initio energies can accurately predict the barrier height, and it is the first step toward molecular simulations, such as the reactive or non-reactive collisions and photodissociation within the system (Conway et al., 2020; Caridade et al., 2013; Schmidt et al., 2013). Therefore, a dynamics study carried out on such a PES is reliable, and the information including the reaction probability and integral cross section (ICS) can be well predicted. The first analytical PES of PO₂(X²A₁) was constructed using the B3P86/6-311$++$G(3df, 3pd) energy points (Zeng and Zhao, 2012), but the dissociation energy was well beyond the experimental value, and the intermediates were not all predicted. For performing a high-precision dynamics study of the title reaction, it is necessary to develop a global PES of PO₂(X²A₁), in which the potential energies should be calculated with the state-of-the-art ab initio method, and the reaction path should be well reproduced and carefully validated.

This work aims to establish a global PES for the ground-state PO₂(X²A₁) so as to present a dynamics study on P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π). The state-of-the-art ab initio method was applied to calculate the potential energies of PO₂(X²A₁). The analytical PES was then generated using the combined-hyperbolic-inverse-power-representation (CHIPR) method (Varandas, 2013; Rocha and Varandas, 2020, 2021). The rate constants for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) reaction at temperatures ranging below 5000 K were obtained using the quasi-classical trajectory (QCT) (Peslherbe et al., 1999; Li et al., 2014) and time-dependent wave packet (TDWP) (Zhang and Zhang, 1993, 1994) methods, and they were compared with available experimental and theoretical results. Moreover, the state-specified reaction probability and ICSs were provided in order to get a deeper understanding of this reaction.

2 Ab initio calculations

All ab initio calculations of the ground-state PO₂(X²A₁) were carried out using the MOLPRO 2015 software package (Werner et al., 2020; Eckert et al., 1997) with the C_s(A^′) symmetry point group. The Dunning-type aug-cc-pV5Z (AV5Z) basis set (Dunning et al., 2001; Martin and Uzan, 1998; Woon and Dunning, 1993) was applied. The calculation processes are as follows. Firstly, the Hartree–Fock (HF) method was used to obtain the single-configuration wave function of PO₂(X²A₁). The relevant multi-configuration wave functions were generated by the full-valence complete active space self-consistent field (CASSCF) method (Knowles and Werner, 1985) based on the HF wave function. Finally, the dynamics correlation energies were considered by the internally contracted multireference configuration-interaction method including the Davidson correction (MRCI(Q)) (Knowles and Werner, 1988; Werner and Knowles, 1988), in which the CASSCF wave functions were set as reference. In the CASSCF and MRCI(Q) calculations, 12 active molecular orbitals ($\mathrm{9}{A}^{\prime }+\mathrm{3}{A}^{\prime \prime }$) involved 17 valence shell electrons, and the remaining 14 inner shell electrons were closed into 7 core orbitals ($\mathrm{6}{A}^{\prime }+\mathrm{1}{A}^{\prime \prime }$). A total of 6471 ab initio energy points were generated from two grids defined in Jacobi coordinates and six additional grids defined around the important intermediates. For example, the R_A-BC, r_BC, and γ_A-BC for the A-BC channel of an ABC molecule are defined as the distance of the A atom relative to the centre of mass of BC, the bond distance of BC, and the angle between R_A-BC and r_BC, respectively. In the P-O₂ channel, the grids were defined by $\mathrm{2.0}\le r_{\mathrm{OO}}/a_{\mathrm{0}}\le \mathrm{5.3}$, $\mathrm{0}\le R_{\mathrm{P}\text{-}\mathrm{OO}}/a_{\mathrm{0}}\le \mathrm{15.0}$, and $\mathrm{0}\le {\mathit{\gamma }}_{\mathrm{P}\text{-}\mathrm{OO}}/\mathrm{deg}\le \mathrm{90}$. In the O-PO channel, the ranges were defined by $\mathrm{2.4}\le r_{\mathrm{PO}}/a_{\mathrm{0}}\le \mathrm{4}$, $\mathrm{2.0}\le R_{\mathrm{O}\text{-}\mathrm{PO}}/a_{\mathrm{0}}\le \mathrm{15.0}$, and $\mathrm{0}\le {\mathit{\gamma }}_{\mathrm{O}\text{-}\mathrm{PO}}/\mathrm{deg}\le \mathrm{180}$. The additional grids around the equilibrium geometry, local minimum, and transition states were constructed to be dense enough according to the geometry optimization (OPTG) (Eckert et al., 1997) results obtained from MOLPRO 2015 (Werner et al., 2020; Eckert et al., 1997).

3 The CHIPR potential energy surface

According to the spin-spatial Wigner–Witmer correlation, the dissociation scheme of the ground-state PO₂(X²A₁) can be described by

$$\begin{array}{}\text{(R5)}& {\displaystyle }{\mathrm{PO}}_{\mathrm{2}}\left({\mathrm{X}}^{\mathrm{2}}{\mathrm{A}}_{\mathrm{1}}{/}^{\mathrm{2}}{\mathrm{A}}^{\prime }\right)& {\displaystyle }\to \mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right)+\mathrm{PO}\left({\mathrm{X}}^{\mathrm{2}}\mathrm{\Pi }\right),\text{(R6)}& {\displaystyle }& {\displaystyle }\to \mathrm{P}\left({}^{\mathrm{4}}\mathrm{S}\right)+{\mathrm{O}}_{\mathrm{2}}\left({\mathrm{X}}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}\right),\text{(R7)}& {\displaystyle }& {\displaystyle }\to \mathrm{P}\left({}^{\mathrm{4}}\mathrm{S}\right)+\mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right)+\mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right).\end{array}$$

The ground-state PO₂(X²A₁) dissociates adiabatically into O(³P) + PO(X²Π), P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$), and P(⁴S) + O(³P) + O(³P). Assuming the energy of infinitely separated P(⁴S) + O(³P) + O(³P) atoms to be the zero point, the global adiabatic CHIPR PES of the ground-state PO₂(X²A₁) has the following many-body expansion (MBE) (Murrell et al., 1984; Varandas and Murrell, 1977) form:

$$\begin{array}{}\text{(1)}& V\left(R\right)=V_{{\mathrm{O}}_{\mathrm{2}}}^{\left(\mathrm{2}\right)}\left(R_{\mathrm{1}}\right)+V_{\mathrm{PO}}^{\left(\mathrm{2}\right)}\left(R_{\mathrm{2}}\right)+V_{\mathrm{PO}}^{\left(\mathrm{2}\right)}\left(R_{\mathrm{3}}\right)+V_{{\mathrm{PO}}_{\mathrm{2}}}^{\left(\mathrm{3}\right)}\left(R\right),\end{array}$$

where V⁽²⁾ represents the two-body fragments represented by the diatomic potential energy curves (PECs) of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π). V⁽³⁾ is the three-body fragment. In the CHIPR method (Rocha and Varandas, 2020, 2021; Varandas, 2013), two-body fragments are given by

$$\begin{array}{}\text{(2)}& V^{\left(\mathrm{2}\right)}\left(R\right)={\displaystyle \frac{Z_{\mathrm{A}}{Z}_{\mathrm{B}}}{R}}\sum _{k=\mathrm{1}}^L{C}_k{y}^k,\end{array}$$

where Z_A and Z_B are the nuclear charges of atoms A and B, respectively. C_k represents expansion coefficients of an Lth-order polynomial, and y is the basis set consisting of the linear combination of R-dependent functions (see Eq. 4 shown below). For the AB₂-type species like PO₂, the CHIPR three-body fragment can be simplified to the following permutation symmetric form (Rocha and Varandas, 2020, 2021; Varandas, 2013):

$$\begin{array}{}\text{(3)}& V^{\left(\mathrm{3}\right)}\left(R\right)=\sum _{i,j,k=\mathrm{0}}^L{C}_{i,j,k}\left[y_{\mathrm{1}}^i\left(y_{\mathrm{2}}^j{y}_{\mathrm{3}}^k+y_{\mathrm{2}}^k{y}_{\mathrm{3}}^j\right)\right],\end{array}$$

where $C_{i,j,k}$ represents expansion coefficients for an Lth-order polynomial, and y_p represents basis sets of coordinates (p=1, 2, 3) for the reference geometry, which are expressed in terms of the Mth-order distributed-origin contracted basis set (Rocha and Varandas, 2020, 2021):

$$\begin{array}{}\text{(4)}& y_{\mathrm{p}}=\sum _{\mathit{\alpha }=\mathrm{1}}^{M-\mathrm{1}}{c}_{\mathit{\alpha }}{\mathit{\varphi }}_{\mathrm{p},\mathit{\alpha }}+c_{\mathrm{M}}{\mathit{\varphi }}_{\mathrm{p},\mathrm{M}},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\mathit{\varphi }}_{\mathrm{p},\mathit{\alpha }}=\mathrm{sech}\left\{{\mathit{\gamma }}_{\mathrm{p},\mathit{\alpha }}\left[R_{\mathrm{p}}-\mathit{\zeta }{\left(R_{\mathrm{p}}^{\mathrm{ref}}\right)}^{\mathit{\alpha }-\mathrm{1}}\right]\right\}\end{array}$$

and

$$\begin{array}{}\text{(6)}& {\mathit{\varphi }}_{\mathrm{p},\mathrm{M}}={\left[{\displaystyle \frac{\tanh \left(\mathrm{0.2}{R}_{\mathrm{p}}\right)}{R_{\mathrm{p}}}}\right]}^{\mathrm{6}}\mathrm{sech}\left\{{\mathit{\gamma }}_{\mathrm{p},\mathrm{M}}\left[R_{\mathrm{p}}-\mathit{\zeta }{\left(R_{\mathrm{p}}^{\mathrm{ref}}\right)}^{M-\mathrm{1}}\right]\right\}\end{array}$$

are primitive bases. γ_p,α represents non-linear parameters, and $R_{\mathrm{ref}}^{\mathrm{p}}$ represents the origin. The fitting process was carried out in the CHIPR 4.0 program (Rocha and Varandas, 2021). For more detailed descriptions of the CHIPR method, we refer the reader to the related manuals (Rocha and Varandas, 2020, 2021). In recent years, there has been a lot of global PESs constructed based on this method, such as H₃ (Varandas, 2013), HO₂ (Varandas, 2013), C₃ (Rocha and Varandas, 2019a), C₃H (Rocha and Varandas, 2019b), PH₂ (Chen et al., 2022), NH₂ (Li et al., 2022), Si₂C (Li et al., 2023), and SiC₂ (Rocha et al., 2022).

3.1 Two-body fragment

Based on Eq. (1), the PECs of the ground states O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π) were fitted to the CHIPR form of Eq. (2) to obtain the two-body fragments of the MBE potential of PO₂(X²A₁). The ab initio potential energies were calculated at the MRCI(Q) level of theory with the AV5Z basis set. For O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$), D_2h symmetry was chosen with the consideration of 10 active molecular orbitals (3 A_g + 1 B_3u + 1 B_2u + 3 B_1u + 1 B_2g +1 B_3g) and 2 closed orbitals (1 A_g + 1 B_1u), while 10 active molecular orbitals (6 A₁ + 2 B₁ + 2 B₂) and 6 closed orbitals (4 A₁ + 1 B₁ + 1 B₂) of C_2v symmetry were applied for PO(X²Π). A total of 43 and 34 ab initio energy points were obtained for O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π), respectively.

During the fitting, the 4th-order contracted bases (Eq. 4) and 8th-order polynomial expression (Eq. 2) were applied for O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π). Figure 1a presents the final CHIPR PECs of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π) with the root-mean-square deviations (RMSDs) of 2.05 and 1.85 cm⁻¹, respectively. Figure 1b shows the deviations between ab initio energy points and the corresponding CHIPR energies, which are within ±6 cm⁻¹. The obtained CHIPR PECs of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π) reproduce well with the ab initio energy points and exhibit a smooth strongly bound valence region and reasonable dissociation behaviours. As shown in Table 1, the spectroscopic constants of our CHIPR PECs agree well with those of previous theoretical and experimental results. Hence, the CHIPR PECs of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π) are reliable and can be used as two-body fragments to construct a global PES of PO₂(X²A₁). The analytic functions of the O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π) PECs are collected in the ready-to-use Fortran code in the Supplement.

Figure 1 (a) The PECs of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π). The unit of potential energy is the atomic unit (hartree, E_h). Solid lines are the CHIPR PECs. Solid circles are ab initio energies. (b) The deviations between ab initio energy points and the corresponding CHIPR energies.

Table 1 Spectroscopic constants of the CHIPR PECs of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) and PO(X²Π), along with previous theoretical and experimental results.

^a The equilibrium geometry in unit of a₀. ^b The dissociation energy in unit of eV. ^c The units of ω_e, ω_eχ_e, α_e, and B_e are cm⁻¹. ^d The CHIPR PECs. ^e Results at the icMRCI(Q)/CBS(56) + CV + DK level (Liu et al., 2014, 2017). ^f Results at the CCSDT/AVQZ level (Sordo, 2001). ^g Ref. (Huber and Herzberg, 1979). ^h Results at the MRCI(Q)-r/aug-cc-wCV5Z level (Prajapat et al., 2017).

3.2 Three-body fragment

For the CHIPR PES of PO₂(X²A₁), the 4th-order contracted bases (Eq. 4) and 12th-order polynomial expression (Eq. 3) were used to fit the three-body fragment. The first trial of the fitted PES presented the complex topographical features around the reaction path, including the important intermediates of 1 global minimum (GM), 1 local minimum (LM), and 4 transition states (TSs). Then, the OPTG was carried out for these configurations, and the additional energy grids were calculated based on the OPTG results. During the subsequent fitting processes, the weights for the additional grids were set to be 5, and those for the remaining energy points were set to be 1. The final PES was constructed from 6471 ab initio energy points, covering an energy range up to 500 kcal mol⁻¹ above the GM. Table 2 lists the RMSDs between the analytical CHIPR energies and ab initio energies in the indicated energy range above the GM and those for the additional energy grids. The GM, LM, and TSs of PO₂(X²A₁) were well reproduced by the CHIPR method with RMSDs lower than 25 cm⁻¹. The total RMSD is 91.5 cm⁻¹ (or 0.262 kcal mol⁻¹), implying the high fitting accuracy and reliability of the CHIPR PES of PO₂(X²A₁). The whole PO₂(X²A₁) PES is collected in the ready-to-use Fortran code in the Supplement.

Table 2 The root-mean-square deviations (RMSDs) in the indicated energy range above the GM and those for the additional energy grids.

^a The number of energy points in the corresponding range. ^b The maximum deviation in the corresponding range, cm⁻¹. ^c The RMSD for the corresponding range, cm⁻¹. ^d The number of energy points with a deviation larger than the RMSD. ^e The indicated energy range above the GM, kcal mol⁻¹. ^f The additional energy grids.

4 Features of CHIPR PES

Figures 2–4 display the topographical features and the relevant stationary points of the CHIPR PES for the ground-state PO₂(X²A₁). Table 3 compares the main attributes of the stationary points for the CHIPR PES with other theoretical and experimental results, including the interatomic distances for OO (R₁) and PO (R₂ and R₃), the bond angle (θ) between R₂ and R₃, the vibration frequencies (symmetrical stretching ω₁, bending ω₂, and antisymmetric stretching ω₃), and the potential energies (E) relative to the P(⁴S) + O(³P) + O(³P) asymptote. In particular, our OPTG results for GM, LM, TS1, TS2, TS3, and TS4 are also collected in Table 3, which are calculated at the MRCI(Q) level as implemented in MOLPRO 2015 (Werner et al., 2020). As shown in Table 3, the attributes of all the stationary points reproduced by the CHIPR method are very similar to the OPTG results.

Figure 2 (a) The contour plot for bond stretching in the OPO bending configuration with θ fixed at 135.1^∘. Contours are equally spaced by 0.02 E_h, starting from −0.4 E_h. (b) The contour plot for the insertion of the P atom into the O₂ fragment with the insertion angle of 50^∘. Contours are equally spaced by 0.011 E_h, starting from −0.38 E_h. (c) The contour plot for P moving around one of the O atoms with R₁ fixed at 2.7 a₀. Contours are equally spaced by 0.006 E_h, starting from −0.4 E_h. (d) The contour plot for bond stretching of OOP colinear configuration. Contours are equally spaced by 0.01 E_h, starting from −0.4 E_h.

Figure 3 (a) The contour plot for P moving around the centre of mass of O₂ with R₁ fixed at 2.285 a₀. Contours are equally spaced by 0.011 E_h, starting from −0.22 E_h. Dashed areas are contours equally spaced by 0.0005 E_h, starting from −0.192 E_h. (b) The contour plot for O moving around the centre of mass of PO with R₂ fixed at 2.788 a₀. Contours are equally spaced by 0.013 E_h, starting from −0.41 E_h. Dashed areas are contours equally spaced by 0.0001 E_h, starting from −0.226 E_h.

Table 3 Attributes of the global minimum (GM), local minimum (LM), transition state (TS) and second-order saddle (SP) of PO₂(X²A₁) CHIPR PES.

^a Results from the geometry optimization (OPTG). ^b Results at the CCSD(T)/AVQZ level (Francisco, 2002). ^c Result at the CCSD(T)/(CBS + SO + DK + CV + ZPE) level (Bauschlicher, 1999). ^d Results at the CCSD(T)/AV5Z level (Liang et al., 2013). ^e Results at the B3LYP/AVQZ level (Douglas et al., 2019). ^f Results at the MRCI(Q)/AVTZ + d level (Gomes et al., 2022). ^g equilibrium geometry deduced from the observed rotational constants (Kawaguchi et al., 1985). ^h The vibrational frequencies from the laser-induced fluorescence spectrum (Hamilton, 1987). ⁱ Origin of the v₃ fundamental band obtained from the infrared absorption spectrum (Qian et al., 1995). ^j The experimental atomization energy D₀(PO₂) (Drowart et al., 1972). ^k Vibration frequencies from fluorescence emission and laser fluorescence excitation spectra (Lei et al., 2001). ^l Term value of the v₃ fundamental band obtained from laser absorption spectrum (Lawson et al., 2011). ^m Results at the SCF/3-21G^* level (Kabbadj and Liévin, 1989). ⁿ Results at the SCF/6-31G^* level (Lohr, 1984).

Figure 2a shows the topographical features of the PO₂(X²A₁) PES for the PO bond stretching from θ fixed at the equilibrium angle of 135.1^∘. As shown in Fig. 2a, the GM of the CHIPR PES has a C_2v symmetric configuration with bond lengths of R₁=5.142 a₀ and $R_{\mathrm{2}}=R_{\mathrm{3}}=\mathrm{2.788}$ a₀, which can be seen in Figs. 3b and 4 as well. The vibrational frequencies are computed to be 1137.6, 464.7, and 1460.0 cm⁻¹ for ω₁, ω₂, and ω₃, respectively. According to Table 3, these attributes of GM agree well with previous theoretical and experimental results, excepted that ω₃ is slightly higher than experimental results and those calculated by the CCSD(T) method. Figure 2a presents the channel of an O atom dissociated from PO₂(X²A₁), and the corresponding dissociation energy D_e(O-PO) of the CHIPR PES is 5.305 eV, which is similar to the experiment value of 5.11 eV (Drowart et al., 1972) and the theoretical result of 5.24 eV at the CCSD($T)/$(CBS + SO + DK + CV + ZPE) level (Bauschlicher, 1999). The dissociation energy of the P atom dissociated from PO₂(X²A₁) is predicted to be 6.23 eV, showing a good concordance of the predicted value of 6.15 eV at the CI/3-21G^* level (Kabbadj and Liévin, 1989).

Figure 2b and c shows the entrance channel of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π), i.e. the contour plots for the insertion of the P atom into the O₂ fragment with the insertion angle of 50^∘ and for P moving around one of the O atoms with R₁ fixed at 2.7 a₀, respectively. As illustrated in these two figures, the MEP of the entrance channel is connected by TS1, LM, and TS4 in turn, which is in accordance with the earlier prediction (Gomes et al., 2022; Douglas et al., 2019). The most important configuration for the reaction is the entrance barrier of TS1. The barrier height relative to the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) asymptote was theoretically predicted to be 0.032 eV at the B3LYP/AVQZ level (Douglas et al., 2019) and 0.158, 0.142, and 0.137 eV at the MRCI$\left(Q\right)/$AVXdZ (X=T, Q, 5) levels (Gomes et al., 2022), respectively. Our CHIPR PES is fitted by the energies at the MRCI$\left(Q\right)/$AV5Z level, and the barrier height of TS1 is 0.133 eV, which is particularly consistent with the result at the MRCI(Q)/AV5dZ level (Gomes et al., 2022). Also, the attributes of TS1, LM, and TS4 showed good agreement with previous theoretical results, which are presented in detail in Table 3.

Figure 2d shows the contour plot for bond stretching of R₁ and R₂ at colinear configuration, including the OOP (upper left corner) and OPO (lower right corner) configurations. The notable features here are the OPO colinear transition state TS2 and POO colinear one TS3. The TS3 is predicted to locate at R₁=2.475 a₀, R₂=2.917 a₀, and R₃=5.392 a₀ and is connected with two second-order saddle points (SP1 and SP3). The configuration of the colinear TS2 is $R_{\mathrm{2}}=R_{\mathrm{3}}=\mathrm{2.774}$ a₀ and very close to GM. The potential energy of TS2 is 0.816 eV higher than that of GM, which is close to previous theoretical results of 1.05 and 0.93 eV at the CI/3-21G^* level (Kabbadj and Liévin, 1989) and MP3/6-31G* level (Lohr, 1984), respectively. As shown in Table 3, the attributes of TS2 are similar to those of previous theoretical results (Kabbadj and Liévin, 1989; Lohr, 1984), except that the vibration frequencies are much less than those obtained at the SCF/3-21G^* level (Kabbadj and Liévin, 1989).

Figure 3a is the contour plot for the P atom moving around the centre of mass of O₂ with R₁ fixed at 2.285 a₀, which displays the smooth long-range behaviour of the CHIPR PES. There exist high barriers (SP1 and SP2) on the entrance channel of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) when the P atom inserts along the mid-perpendicular or molecular axis of O₂. When the Jacobi approaching angle is about 40–50^∘, it is much easier for the P atom to cross the barrier (TS1) and reach the LM, followed by the subsequent reaction process. At relatively low collision energies, therefore, the approaching P atom will bond first with one of the two O atoms rather than with them at the same time. Figure 3b shows the contour plot for O moving around the centre of mass of PO with R₂ fixed at 2.806 a₀, exhibiting the exit channel of the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) reaction. As shown in Fig. 3b, both TS3 and TS4 can be reached from the LM, but the dissociation of the system is relatively difficult to occur at TS3 due to the high-barrier SP3 on the dissociation path. The two wells LM and GM represent the POO and OPO isomers, which are separated by the isomerization barrier TS4. The minimum exit channel is connected by LM, TS4, GM and TS2 in turn, and then the system dissociates into the products of O(³P) + PO(X²Π).

Figure 4 depicts the relaxed triangular contour plot (Varandas, 1987) for the CHIPR PES of PO₂(X²A₁) using the scaled hyperspherical coordinates (${\mathit{\beta }}^{*}=\mathit{\beta }/Q$ and ${\mathit{\gamma }}^{*}=\mathit{\gamma }/Q$), given by

$$\begin{array}{}\text{(7)}& \left(\begin{array}{c}Q\\ \mathit{\beta }\\ \mathit{\gamma }\end{array}\right)=\left(\begin{array}{crr}\mathrm{1}& \mathrm{1}& \mathrm{1}\\ \mathrm{0}& \sqrt{\mathrm{3}}& -\sqrt{\mathrm{3}}\\ \mathrm{2}& -\mathrm{1}& -\mathrm{1}\end{array}\right)\left(\begin{array}{c}{R}_{\mathrm{1}}^{\mathrm{2}}\\ R_{\mathrm{2}}^{\mathrm{2}}\\ R_{\mathrm{3}}^{\mathrm{2}}\end{array}\right).\end{array}$$

The hidden coordinate Q (i.e. the sum of squares of the three bond distances) is allowed to relax to give the lowest potential energy, while β and γ define the shape of the triangle formed by the three atoms (see the work of Varandas, 1987, for details). It gives a better view of the major topographical features of the CHIPR PES, including the configurations and symmetries of all stationary points mentioned above. The MEP is connected by TS1, LM, TS4, GM and TS2 in turn. At high collision energies, the P atom is able to cross the C_2v barrier SP2 and reach the GM directly, which will be confirmed in the following dynamics calculations. In this condition, the P atom approaches along the mid-perpendicular of O₂ and the two PO bonds stretch symmetrically, accompanied by the progressively open OPO angle, as shown in Fig. 4. Thus, after the P atom crosses over the SP2, the system will reach the symmetric GM instead of the asymmetric TS1, LM, and TS4. Then, the system evolves through the linear transition state TS2 and finally dissociates into the products of O(³P) + PO(X²Π). There is also a collinear abstraction reaction path through TS3, which is difficult to occur due to the mutation of the molecular structure and the existence of collinear barriers (SP1 and SP3) on the entrance and exit channels. We should reiterate that the critical intermediates were carefully verified by the OPTG results, the relevant ab initio energy grids were constructed to be dense enough, and the CHIPR method reproduced them well. Hence, our CHIPR PES of PO₂(X²A₁) has reliable reaction channels and can be used to perform relevant dynamics calculations.

5 Dynamics calculations

Based on our CHIPR PES of PO₂(X²A₁), the reaction probability, ICS, and rate constants of the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) reaction were calculated using the QCT and TDWP methods. For each QCT calculation, the sampled trajectories, initial distance of the reactants, and integration time step of classical equations of motion were 100 000, 15 Å, and 0.2 fs, respectively. The state-specified QCT ICS of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v, j) → O(³P) + PO(X²Π) was calculated by (Li et al., 2014):

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\mathrm{r}}\left(E_{\mathrm{tr}};v,j\right)=\mathit{\pi }{b}_{\max }^{\mathrm{2}}{\displaystyle \frac{N_{\mathrm{r}}}{N}},\end{array}$$

where b_max is the maximum impact parameter, and N and N_r represent the total and reactive trajectories, respectively. The state-specified and thermal QCT rate constants, k(T), were obtained by (Li et al., 2014)

$$\begin{array}{}\text{(9)}& k\left(T\right)=g_{\mathrm{e}}\left(T\right){\left({\displaystyle \frac{\mathrm{8}{k}_{\mathrm{B}}T}{\mathit{\pi }{\mathit{\mu }}_{\mathrm{P}+{\mathrm{O}}_{\mathrm{2}}}}}\right)}^{\mathrm{1}/\mathrm{2}}\mathit{\pi }{b}_{\max }^{\mathrm{2}}{\displaystyle \frac{N_{\mathrm{r}}}{N}}\end{array}$$

and

$$\begin{array}{}\text{(10)}& \begin{array}{rl}k\left(T\right)& =g_{\mathrm{e}}\left(T\right){\left({\displaystyle \frac{\mathrm{2}}{k_{\mathrm{B}}T}}\right)}^{\mathrm{3}/\mathrm{2}}{\left({\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\right)}^{\mathrm{1}/\mathrm{2}}{Q}_{\mathrm{v}j}^{-\mathrm{1}}\left(T\right)\sum _{vj}(\mathrm{2}j+\mathrm{1})\\ & \times \exp \left(-{\displaystyle \frac{E_{\mathrm{v}j}}{k_{\mathrm{B}}T}}\right)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{E}_{\mathrm{tr}}{\mathit{\sigma }}^x\exp \left(-{\displaystyle \frac{E_{\mathrm{tr}}}{k_{\mathrm{B}}T}}\right)\mathrm{d}{E}_{\mathrm{tr}},\end{array}\end{array}$$

where ${\mathit{\mu }}_{\mathrm{P}+{\mathrm{O}}_{\mathrm{2}}}$ is the reduced mass of the reactants, k_B is the Boltzmann constant, Q_vj(T) is the rotational–vibrational partition function for all the rotational–vibrational states of O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$), E_vj is the energy of the (v, j) state, and E_tr is the translation energy. The rate constant was computed adiabatically for the ground-state PO₂(X²A₁) and the electronic degeneracy factor g_e(T) assumed the following form (Graff and Wagner, 1990):

$$\begin{array}{}\text{(11)}& g_{\mathrm{e}}\left(T\right)=g_{{\mathrm{PO}}_{\mathrm{2}}}\left(T\right){\left(q_{\mathrm{P}}\left(T\right)q_{{\mathrm{O}}_{\mathrm{2}}}\left(T\right)\right)}^{-\mathrm{1}},\end{array}$$

where gPO₂(T)=2 is the degeneracy of the ground-state PO₂(X²A₁), q_P(T)=4 is the electronic partition function accounting for the fine structure of P(⁴S), and qO₂(T) takes into account two spin-orbit states for O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$), given by

$$\begin{array}{}\text{(12)}& q_{{{\mathrm{O}}_{\mathrm{2}}}_{\left({}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}\right)}}\left(T\right)=g\left({}^{\mathrm{3}}{\mathrm{\Sigma }}_{{\mathrm{0}}^{+}}^{-}\right)+g\left({}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{1}}^{-}\right)\exp (-\mathrm{\Delta }/T),\end{array}$$

where $g{(}^{\mathrm{3}}{\mathrm{\Sigma }}_{{\mathrm{0}}^{+}}^{-})=\mathrm{1}$, $g{(}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{1}}^{-})=\mathrm{2}$ due to the doubly degenerate $\mathrm{\Omega }=\pm \mathrm{1}$, and Δ=2.88 K (2 cm⁻¹) is the spin-orbit splitting between ${}^{\mathrm{3}}{\mathrm{\Sigma }}_{{\mathrm{0}}^{+}}$ and ${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{1}}^{-}$ (Liu et al., 2014; Barrow and Yee, 1974).

Figure 4 The relaxed triangular contour plot for the ground-state PO₂ in hyperspherical coordinates (see the definition in the text). Contours are equally spaced by 0.008 E_h, starting from −0.41 E_h.

Table 4 Parameters used in the TDWP calculations. All parameters are given in atomic units, except for the numbers or quantities indicated otherwise.

During the TDWP calculations, the split-operator scheme was used to solve the Schrödinger equation. The Hamiltonian for the reactants (P and O₂) in the PO₂ system can be represented using the Jacobi coordinates:

$$\begin{array}{}\text{(13)}& H=-{\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{R}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial R^{\mathrm{2}}}}+{\displaystyle \frac{\left(\widehat{J}-\widehat{j}\right)}{\mathrm{2}{\mathit{\mu }}_{\mathrm{R}}{R}^{\mathrm{2}}}}+{\displaystyle \frac{{\widehat{j}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{r}}{r}^{\mathrm{2}}}}+V(R,r)+h\left(r\right),\end{array}$$

where R is the distance of the P atom relative to the centre of mass of O₂, r is the bond distance of O₂, μ_R is the reduced mass between P and O₂, μ_r is the reduced mass of O₂, $\widehat{J}$ is the total angular momentum and $\widehat{j}$ is the rotational angular momentum of O₂, V(R, r) is the Jacobi form of the CHIPR PES for PO₂, and h(r) is the diatomic reference Hamiltonian, given by

$$\begin{array}{}\text{(14)}& h\left(r\right)=-{\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{r}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial r^{\mathrm{2}}}}+V\left(r\right),\end{array}$$

where V(r) is the PEC of O₂. For the TDWP calculations on the adiabatic PES of PO₂(X²A₁), the time-dependent wave function was expanded by the body-fixed (BF) translational–vibrational–rotational basis:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathrm{\Psi }}_{v_{\mathrm{0}}{j}_{\mathrm{0}}{K}_{\mathrm{0}}}^{JM\mathit{\epsilon }}(R,r,t)=& \sum _{n,v,j,K}{F}_{nvjK,v_{\mathrm{0}}{j}_{\mathrm{0}}{K}_{\mathrm{0}}}^{JM\mathit{\epsilon }}\left(t\right)\\ & u_n^{\mathrm{v}}\left(R\right){\mathit{\phi }}_v\left(r\right)Y_{jK}^{JM\mathit{\epsilon }}(R,r),\end{array}\end{array}$$

where M and K are the projection of J onto the z axis of the space-fixed and BF coordinates, respectively; $u_n^{\mathrm{v}}\left(R\right)$, φ_v(r), and $Y_{jK}^{JK\mathit{\epsilon }}$ are the translational basis, reference vibration eigenfunction for O₂, and total angular momentum eigenfunction, respectively; ε is the parity of the system; n is the label of the translational basis; and v₀, j₀, and K₀ denote the initial rotational–vibrational state of O₂.

The dynamics information is extracted from the final wave packet after a long propagation time. The reaction probability $P_{v_{\mathrm{0}}{j}_{\mathrm{0}}{K}_{\mathrm{0}}}^J\left(E\right)$, total reaction cross section ${\mathit{\sigma }}_{v_{\mathrm{0}}{j}_{\mathrm{0}}}\left(E\right)$, and temperature-dependent rate constant $k_{v_{\mathrm{0}}{j}_{\mathrm{0}}}\left(T\right)$ can be calculated by

$$\begin{array}{}\text{(16)}& {\displaystyle }{P}_{v_{\mathrm{0}}{j}_{\mathrm{0}}{K}_{\mathrm{0}}}^J\left(E\right)={\displaystyle \frac{h}{{\mathit{\mu }}_{\mathrm{r}}}}\mathrm{Im}\left[\mathit{\psi }\left(E\right)\left|\mathit{\delta }\left(r-r_{\mathrm{0}}\right){\displaystyle \frac{\partial }{\partial r}}\right|\mathit{\psi }\left(E\right)\right],\text{(17)}& {\displaystyle }{\mathit{\sigma }}_{v_{\mathrm{0}}{j}_{\mathrm{0}}}\left(E\right)={\displaystyle \frac{\mathit{\pi }}{\left(\mathrm{2}{j}_{\mathrm{0}}+\mathrm{1}\right)k^{\mathrm{2}}}}\sum _{K_{\mathrm{0}}}\sum _J(\mathrm{2}J+\mathrm{1})P_{v_{\mathrm{0}}{j}_{\mathrm{0}}{K}_{\mathrm{0}}}^J,\text{(18)}& {\displaystyle }\begin{array}{rl}{k}_{v_{\mathrm{0}}{j}_{\mathrm{0}}}\left(T\right)=& g_{\mathrm{e}}\left(T\right)\sqrt{{\displaystyle \frac{\mathrm{8}{k}_{\mathrm{B}}T}{\mathit{\pi }{\mathit{\mu }}_{\mathrm{R}}}}}{\left(k_{\mathrm{B}}T\right)}^{-\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}E{\mathit{\sigma }}_{v_{\mathrm{0}}{j}_{\mathrm{0}}}\left(E\right)\\ & \exp \left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right)\mathrm{d}E,\end{array}\end{array}$$

where k is the wavenumber of the initial state with the fixed collision energy E, and k_B is the Boltzmann constant. Table 4 displays the parameters used in the TDWP calculations. The J-shifting method was used in the calculation, and the reaction probabilities of J=0 to 270 were obtained.

Figure 5 presents the CHIPR MEP for P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) along with the potential energies of the corresponding intermediates relative to the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) asymptote. This reaction is an exothermic reaction (−0.927 eV) with an intrinsic entrance barrier TS1 (0.133 eV) and two isomers (LM, GM) separated by the isomerization barrier TS4. Along the MEP, the P atom approaches O₂ at the Jacobi angle of about 40–50^∘ and bonds with one of the two O atoms to form an unstable transition state TS1, followed by a transitory POO isomer LM. The system then evolves through TS4, the OPO isomer GM, and the linear transition state TS2, accompanied by the progressively open OPO angle (i.e. 44.2^∘ for TS4, 134.5^∘ for GM, and 180^∘ for TS2). Finally, the linear PO₂ dissociates into the products of O(³P) + PO(X²Π). The structural diagrams of these intermediaries are also shown in Fig. 5, and the corresponding geometric parameters are given in Table 3.

Figure 5 The minimum energy path (MEP) and relative intermediates of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π). Energies of intermediates relative to the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) asymptote are given in brackets and in units of eV.

The TDWP reaction probability of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v=0, j=0) → O(³P) + PO(X²Π) at the total angular momentum values of J=0, 80, 160, and 240 are presented in Fig. 6. For J=0, the probability starts with a threshold and gradually rises to a peak. Subsequently, the probability decreases until the collision energy reaches 0.91 eV; then it climbs again and eventually stabilizes. The threshold is about 0.1 eV, resulting from the C_s barrier TS1 (0.133 eV) on the entrance channel. The secondary elevation after 0.91 eV is probably due to the opening of a new entrance channel, i.e. the P atom crosses the C_2v barrier SP2 (0.936 eV relative to the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) asymptote) along the mid-perpendicular of O₂ and then reaches GM directly, as discussed in Sect. 4. Both the threshold and the starting point of the secondary elevation are slightly smaller than the corresponding entrance barrier heights because of the quantum tunnelling effects. The trend of probability divides into three stages, as marked in Fig. 6. The first and third stages are the two ascending ones in the probability, showing the most common tendency of reactions controlled by barriers (TS1 and SP2). The second stage is the descending one after the peak, which is dominated by the exothermicity of the reaction.

Figure 6 The TDWP reaction probabilities for P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v=0, j=0) → O(³P) + PO(X²Π) for J=0, 80, 160, and 240 as a function of the collision energy.

For rotational cases (J>0), the centrifugal effect appears, and the rotational (or centrifugal) energy $J(J+\mathrm{1})h^{\mathrm{2}}/\left(\mathrm{2}{\mathit{\mu }}_{\mathrm{R}}{R}^{\mathrm{2}}\right)$ is added to the potential energy to obtain an effective potential (Waage and Rabinovitch, 1970). Since the centrifugal energy is positively correlated with J, the centrifugal barrier gradually increases with the increment of J, resulting in the need for more translational energy to push the reactants over the barrier. Hence, the threshold shifts to larger collision energies with the increment of J, as shown in Fig. 6.

The reaction probabilities obtained by the TDWP method are also oscillatory in all three stages, which is the typical characteristic of quantum resonances. There are numerous bound and quasi-bound states which exist in the LM and GM potential wells, so the temporary reaction complexes are formed there under the bondage of the potential well, leading to resonances. As shown in Fig. 5, the potential energy of the exit channel is 5.305 eV higher than GM leading to a deep potential well, although this reaction is exothermic. Hence, the bondage of the potential well is strong at low collision energy, resulting in the long lifetimes of the collision complexes and plenty of sharp and violent resonances before the peak (stage 1). Higher collision energy can help the complex get rid of the bondage of the potential well and form a short-lived complex, which weakens the quantum resonance and makes the curves of probability smoother in stages 2 and 3.

Figure 7 displays the state-specified QCT (v=0, 1, and 2) and TDWP (v=0) ICSs for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v, j=0) → O(³P) + PO(X²Π) reaction as a function of the collision energy. The TDWP ICSs of vibrational excitation were not calculated due to the extremely expensive cost for this reaction involving three heavy atoms and deep potential wells. As shown in Fig. 7, both QCT and TDWP ICSs rise rapidly from the threshold in stage 1 and gradually reach a plateau (stage 2); they then increase again and finally stabilize (stage 3). Also, the TDWP ICS at v=0 is larger than the QCT ICS at v=0 due to the quantum effects. The threshold is the minimum collision energy at which the reaction can occur (i.e. the point of intersection of ICS and the x axis); these thresholds are about 0.1 for TDWP ICS at v=0 and 0.133, 0.11, and 0.09 eV for QCT ICSs at v=0, v=1, and v=2, respectively. The threshold of QCT ICS at v=0 is consistent with the entrance barrier height, but the threshold of TDWP ICS at v=0 is less than it, because the TDWP method includes the tunnelling effect and the zero-point-energy (ZPE) correction.

Figure 7 The state-specified QCT (v=0, 1, and 2) and TDWP (v=0) integral cross sections for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v, j=0) → O(³P) + PO(X²Π) reaction as a function of the collision energy.

Furthermore, the vibrational excitation of O₂ has different effects on the reactivity for the three stages of this reaction. For stage 1, the reactivity increases with increasing vibrational excitation, because the increased vibrational energy facilitates the reaction through the path with a barrier, resulting in less collision energy required for high vibrational states. Also, the threshold tends to decrease for increasing vibrational excitations. For stage 2, where the exothermicity dominates, the vibrational excitations of the reactants suppress the reaction reactivity. For stage 3, both MEP and the new entrance channel are contributing to the reaction, in which the reaction through MEP is dominated by the exothermicity like stage 2, and the reaction through the new entrance channel is dominated by the entrance barrier (SP2) like stage 1. In other words, the reaction through the new entrance channel is promoted by the vibrational excitation, so the slopes of QCT ICSs at v=1 and v=2 are significantly greater than QCT ICS at v=0, as shown in Fig. 7. Moreover, the vibrational excitations of the reactants suppress the reactivity of the reaction through MEP, and the reaction mainly occurs in this channel. Hence, the combined effect of the vibrational excitation for the two channels in stage 3 is that the vibrational excitations of the reactants suppress the reaction reactivity.

Figure 8 presents the state-specified QCT and TDWP rate constants for P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v=0, j=0) → O(³P) + PO(X²Π) versus temperature ranging below 5000 K. This reaction is exothermic and dominated by the entrance barrier TS1, so the rate constants show a significant Arrhenius linear behaviour at relatively low temperatures, as shown in Fig. 8b. At high temperatures, the rate constants deviate from the Arrhenius behaviour, because the effect of the barrier is weakened and the reaction activation energy becomes temperature-dependent. As expected, the TDWP rate constant is higher than the QCT rate constant due to the quantum effects, especially at low temperatures.

Figure 8 (a) The state-specified rate constant k(T) for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$; v=0, j=0) → O(³P) + PO(X²Π) reaction as a function of temperature. (b) The corresponding Arrhenius plot. The green dashed-dotted line corresponds to the theoretical k(T) ($\mathrm{150}\le T/K\le \mathrm{1400}$) of P(⁴S) + O₂ based on the RRKM theory (Douglas et al., 2019). The dashed blue line corresponds to the theoretical rate constant of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) based on the TST and carried out on MEPs of the ²A^′ and ⁴A^′ electronic states of PO₂ (Gomes et al., 2022). Diamond symbols represent the experimental values using pulsed-laser photolysis and laser-induced fluorescence techniques (Douglas et al., 2019); circle symbols denote the experimental values using the discharge-flow method (Henshaw et al., 1987); triangle symbols represent the experimental values using resonance-fluorescence detection in a discharge-flow system (Clyne and Ono, 1982); quadrate symbols denote the experimental values using resonance fluorescence method (Husain and Slater, 1978); cross symbols represent the experimental values obtained from the attenuation of atomic resonance radiation in the vacuum ultraviolet (Husain and Norris, 1977).

Previous experimental rate constants are also given in Fig. 8 for comparison. Obviously, the rate constants at room temperature determined by various experiments differ by almost an order of magnitude. The earliest two experiments produced extremely large rate constants (Husain and Norris, 1977; Husain and Slater, 1978), and such large values were attributed to the large amounts of secondary dissociation products (Clyne and Ono, 1982). The subsequent experiments were carefully performed to minimize this possible effect and obtained relatively lower rate constants (Clyne and Ono, 1982; Henshaw et al., 1987; Douglas et al., 2019), in which the recent experiment presented the temperature-dependent rate constants (Douglas et al., 2019). As shown in the Arrhenius plot of Fig. 8b, however, the recent experimental rate constants do not exhibit a good linear behaviour, unless the values at maximum and minimum temperatures are excluded. Therefore, four experimental rate constants at about 300–600 K are expected to be more reliable. As shown in Fig. 8, the present QCT and TDWP rate constants are lower than the experimental values. It is partly because the experimental results may be suffering from secondary chemistry; that is, some O atoms and PO molecules could be produced from reactions of P(²P, ²D) + O₂, albeit at extremely low quantities.

The theoretical rate constants based on the RRKM theory (Douglas et al., 2019) and TST (Gomes et al., 2022) are also shown in Fig. 8. The RRKM rate constant available for temperatures of 150–1400 K was optimized by the experimental data, and the TST rate constant was fitted from the calculated values at the temperatures below 1000 K. The rate constant obtained by TST includes the contribution of the doublet and quartet electronic states of PO₂ corresponding to P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$), in which the quartet state works mainly at temperatures above about 600 K (Gomes et al., 2022). We did not consider the contribution from the quartet state due to the extremely high computational cost; therefore, our rate constants at very high temperatures may be slightly underestimated. However, the entrance barrier height including ZPE correction for the quartet state (0.3 eV) is about 3 times above that of the doublet state (0.105 eV) (Gomes et al., 2022), so the effect of quartet state on the rate constant may be small. At low temperatures, the RRKM rate constant is much higher than the TST rate constant and our QCT and TDWP rate constant, which can be explained by the fact that the low-temperature reactivity of such barrier-dominated reaction like P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) depends sensitively on the entrance barrier height. The entrance barrier height adopted when calculating the RRKM rate constant was obtained to be 0.032 eV at the B3LYP/AVQZ level (Douglas et al., 2019). This barrier height was also predicted to be 0.158, 0.142, and 0.137 eV at the MRCI$\left(Q\right)/$AVXdZ (X=T, Q, 5) levels (Gomes et al., 2022), in which the MRCI$\left(Q\right)/$AVTdZ result was modified to be 0.105 eV by including the ZPE correction, and then it was used in the TST study (Gomes et al., 2022). There is evidence that the B3LYP method is not good at predicting the barrier heights on the reaction path (Zhao and Truhlar, 2008; Peverati and Truhlar, 2012). Our CHIPR PES is fitted by the potential energies at the MRCI$\left(Q\right)/$AV5Z level, and the entrance barrier height is 0.133 eV, which is particularly consistent with the result at the MRCI$\left(Q\right)/$AV5dZ level (Gomes et al., 2022), so the CHIPR PES is accurate and reliable. The relatively small QCT rate constants at low temperatures are due to the ZPE problem in the QCT treatment. The TDWP rate constants consider the ZPE correction and the tunnelling through the entrance channel barrier, and the threshold of ICS is about 0.1 eV, agreeing with the modified barrier height (0.105 eV) used in the TST calculation. The difference in the barrier height (0.05 eV) may be due to the tunnelling effect.

At high temperatures, it is necessary to point out that the previous two theoretical Arrhenius formulas (Douglas et al., 2019; Gomes et al., 2022) were fitted without the backing of high-temperature rate constants. Since the rate constant will deviate from the Arrhenius behaviour at high temperatures, the RRKM and TST rate constants remains to be verified at high temperatures. Moreover, the rate constants obtained by RRKM theory are obviously lower than our QCT and TDWP rate constants at high temperatures, because the original Arrhenius formula ($\mathrm{4.2}\times {\mathrm{10}}^{-\mathrm{12}}\exp (-\mathrm{600}/T)$) applied for the RRKM rate constants assumes that the reaction activation energy is temperature-independent, which is not applicable at high temperatures. Note that the population of the rotational–vibrational excited O₂ increases at high temperatures, which could affect the reaction activity according to the analysis of ICS above. Figure 9 compares the state-specified (v=0, j=0) QCT and TDWP rate constants and thermal QCT rate constant for P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) versus the temperature ranging from 2000 to 5000 K. The thermal TDWP rate constant was not calculated due to the extremely expensive cost. As shown in Fig. 9, the rotational–vibrational ground-state (v=0, j=0) O₂ plays an absolutely dominant role at temperatures below 3000 K, while the rotational–vibrational excited O₂ appears at temperatures above 3000 K and promotes the reactive activity. Hence, our TDWP rate constant is reliable at temperatures below 3000 K, while it probably underestimates the rate constant at temperatures above 3000 K due to the neglect of the rotational–vibrational excitation of O₂.

Figure 9 The state-specified (v=0, j=0) QCT and TDWP rate constants and thermal QCT rate constant for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) reaction as a function of temperature.

Table 5 Parameters for the Kooij function obtained by fitting the computed rate constants.

The rate constants of TDWP (v=0, j=0), QCT (v=0, j=0), and QCT (thermal) for the P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) reaction can be approximated using the three-parameter Kooij function, given by (Laidler, 1996)

$$\begin{array}{}\text{(19)}& k\left(T\right)=A{\left({\displaystyle \frac{T}{\mathrm{300}}}\right)}^{\mathit{\alpha }}{e}^{-\mathit{\beta }/T},\end{array}$$

where A, α, and β are fitting parameters. The rate constant curves are divided into four temperature ranges, and the fitting parameters are summarized in Table 5. The fitted rate constants deviate less than 1 % from our calculated ones. These rate constants are calculated in the atom–diatom system and independent of pressure, which may be applicable to the ISM and upper planetary atmospheres.

6 Atmospheric implications

The ablation of IDPs in the upper atmosphere (mainly at heights between 70 and 110 km) of terrestrial planets delivers about 6200 kg yr⁻¹ of ablated phosphorus to the atmosphere (Carrillo-Sánchez et al., 2020), where the temperatures range from below 100 to over 2500 K (in ablation IDPs). Several reaction networks of meteor-ablated phosphorus in the Earth's upper atmosphere (Douglas et al., 2019, 2020; Plane et al., 2021) indicate that the initial oxidation of P will proceed through the successive oxidation by O₂ to produce OPO (i.e. Reactions R1 and R2). The oxidation by O₃ is not significant, because O₂ is 10⁵ times more abundant than O₃ at this altitude. Also, OPO is likely dissociated into PO and P as a result of hyperthermal collisions with air molecules (Carrillo-Sánchez et al., 2020). Therefore, the P + O₂ → O + PO reaction may occur throughout the upper mesosphere and thermosphere. Our rate constants are fitted by sufficient data below 5000 K, which is appropriate for most altitudes of the Earth's atmosphere and can be used to model its phosphorus chemistry.

A recent study developed a reaction network of phosphorus atmosphere chemistry (Plane et al., 2021), including the possible routes from P + O₂ → O + PO to the stable reservoirs (H₃PO₃ and H₃PO₄). Subsequently, they incorporated the rate constants of the associated reactions into the Whole Atmosphere Community Climate Model from the US National Center for Atmospheric Research (Gettelman et al., 2019) and then explored the vertical profiles of the P-containing species and the global mean P deposition flux. Also, they estimated that the fraction of the ablated phosphorus forming bioavailable metal phosphites was 11 %.

One of our concerns is that the theoretical predicted rate constants of P + O₂ → O + PO diverge from experimental result at 200 K, which is about the typical temperature of the upper mesosphere and lower thermosphere. For the reaction with entrance barrier, its rate constant usually has Arrhenius linear behaviour at relatively low temperatures. Our QCT and TDWP rate constants, as well as the TST rate constant (Gomes et al., 2022), exhibit highly consistent Arrhenius linear behaviour at relatively low temperatures, with slopes close to those of the experimental rate constants at about 300–600 K (Douglas et al., 2019), as shown in Fig. 8b. There is reason to believe that the experimental results of 300–600 K are very reliable, whereas the experimental rate constant at about 200 K diverges from this linear behaviour and seems to be slightly overestimated. Therefore, the lifetime of the P atoms in the atmosphere may be a little longer than previously thought (Carrillo-Sánchez et al., 2020). Further experiments for the rate constants at approximately 200 K of P + O₂ → O + PO are encouraged in the future.

7 Conclusions

In this work, we have constructed a global CHIPR PES of the ground state PO₂(X²A₁) based on a total of 6471 ab initio energy points computed at the MRCI$\left(Q\right)/$AV5Z level. The ab initio grids for the critical intermediates are constructed to be dense enough based on the OPTG results. The total RMSD of the CHIPR PES is 91.5 cm⁻¹, and the RMSDs for the critical intermediates are lower than 25 cm⁻¹. The PES presents complex topographical features with multiple potential wells and barriers. The long-range interactions, diatomic potentials, and dissociation energies of each asymptotic channel are reasonably reproduced. The attributes of the main intermediates agree well with available experimental and theoretical results as well as our OPTG results. The corresponding adiabatic MEP of P(⁴S) + O₂(X${}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) → O(³P) + PO(X²Π) features the barrier insertion of the P atom into the O₂ bond at the Jacobi angle of about 40–50^∘, and the reaction then evolves through LM, TS4, GM, and TS2 in turn. Based on the CHIPR PES, the state-specified reaction probability and ICSs for this reaction are calculated using QCT and TDWP methods. The results show three stages of reaction controlled by barrier TS1 (stage 1), exothermicity (stage 2), and barrier SP1 (stage 3) in turn. The vibrational excitation of O₂ promotes the reaction in stages 1 and 3 but suppresses the reactivity in stage 2. Meanwhile, the state-specified and thermal rate constants are predicted for the temperatures ranging from 200 to 5000 K and then compared with available experimental and theoretical results. The rate constants show a significant Arrhenius linear behaviour at relatively low temperatures, but they deviate from the Arrhenius behaviour at high temperatures. The rotational–vibrational ground-state O₂ plays an absolutely dominant role at temperatures below 3000 K, while the rotational–vibrational excitation of O₂ promotes the reactive activity at temperatures above 3000 K.

The presented CHIPR PES of PO₂(X²A₁) can be used for molecular simulations of reactive or non-reactive collisions and photodissociation of the PO₂ system in atmospheres and interstellar media. Moreover, this analytical PES of PO₂ can provide reliable two-body fragments and three-body fragment for the construction of PO₃ and HPO₂ PESs using MBE form so as to carry out a dynamics study of Reactions (R2) and (R3), which are also important for generating H₃PO₃ in the Earth's atmosphere. The computed reaction probability, ICSs, and rate constants may help to explain the relevant thermochemical reactions in related atmospheric and interstellar media.