Relevance of the P+O2 Reaction for PO Formation in Astrochemical Environments: Electronic Structure Calculations and Kinetic Simulations

Phosphorus monoxide (PO) is a key brick of prebiotic chemistry since it is a potential precursor of phosphates, which are present in all living systems. Prompted by the lack of information on the different processes involved in the formation of PO, we have revisited and analyzed in detail the P(4 S) + O2(3Σ−) and P(4 S) + O2(1Δ) reactions leading to PO. The former process has been widely studied from both experimental and theoretical points of view, however, with contradictory results. We have used high-level quantum-chemical calculations to accurately describe the reaction mechanisms. Next, rate constants have been computed using a master equation approach based on ab initio transition state theory. By incorporating the P(4 S) + O2(3Σ−) reaction in an astrochemical model, we have found that this reaction cannot be overlooked when aiming at a complete understanding of the PO abundance in regions dominated by shocks with speeds below 40 km s−1.

Phosphorus is one of the most important elements for life because it is present in all living organisms in the form of phosphates (PO - 4 2 ) and their derivatives.PO is thus considered one of the most important prebiotic molecules detected in the ISM so far, because it is the only one showing the P-O bond and has been suggested as the starting compound in a reaction scheme leading to phosphates (Douglas et al. 2019).Among the phosphorus-bearing molecules detected in the ISM, PO is also the most abundant one, followed by PN.Therefore, these two species play a key role in explaining the chemistry of phosphorus in the Universe.
PN is slightly less abundant than PO; the PO/PN ratio ranges in the ∼1-3 interval, and its origin has been highly debated (Lefloch et al. 2016;Rivilla et al. 2016Rivilla et al. , 2018Rivilla et al. , 2020;;Bergner et al. 2019;Bernal et al. 2021).Interestingly, the measured PO/PN ratio is very similar across a number of starforming regions.Prior to the recent work by García de la Concepción et al. (2021), there was not a clear explanation for the observed PO/PN ratios, especially in the presence of energetic phenomena, such as shocks (Jiménez-Serra et al. 2018).García de la Concepción et al. (2021) computed the kinetic parameters for the reaction between atomic phosphorus (both in the ground ( 4 S) and in the first excited ( 2 P) electronic state) with the hydroxyl radical (OH( 2 Π)).In particular, the P( 4 S) + OH( 2 Π) reaction was demonstrated to be one of the most important sources of PO in the ISM.The incorporation of these new kinetic parameters in astrochemical models allowed the prediction of PO/PN ratio values consistent with those observed in shocked regions (García de la Concepción et al. 2021).In that work, the reaction between P( 4 S/ 2 P) and H 2 O was also investigated, but it was found to affect the chemistry of phosphorus in shocked regions in a negligible manner.It is notable that Douglas et al. (2022) performed an experimental revision of reaction rates for processes with phosphorus in an excited state and obtained smaller rate coefficients for the reaction between P( 2 P) and H 2 O with respect to the estimate by García de la Concepción et al. (2021).
The reaction between atomic P( 4 S) and O 2 ( 3 Σ − ) is another potential source of interstellar PO, with important implications also in atmospheric chemistry (see Douglas et al. 2019;Chen et al. 2023, and references therein).This reaction has been widely studied experimentally under conditions mimicking those typical of the Earth's atmosphere (Husain & Norris 1977;Husain & Slater 1978;Aleksandrov et al. 1982;Clyne & Ono 1982;Sausa et al. 1986;Douglas et al. 2019), obtaining, however, contradictory results.The combined experimental/theoretical study performed by Douglas et al. (2019) pointed out an emerged barrier in the entrance channel, with an estimated value of 3.1 kJ mol −1 at the (U)B3LYP/aug-cc-pVQZ level.In Douglas et al. (2019), the pressure-and temperature-dependent global rate constants were evaluated with the MESMER software (Glowacki et al. 2012) by fitting the RRKM parameters to the experimental data.Since a quantitative description of the strong static correlation present in this system can be hardly achieved by hybrid density functionals, Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.Gomes et al. (2022) reconsidered the P + O 2 reaction, employing multireference ab initio methods (CASSCF and MRCI(Q) for geometry optimizations and energy refinements, respectively).These computations pointed out an emerged barrier of 10.1 kJ mol −1 in the entrance channel, an outcome far from the theoretical estimate by Douglas et al. (2019) and their interpretation of the experimental RRKM parameters.Very recently, an analytical global potential energy surface (PES) for the electronic ground state of PO 2 , characterizing in detail the P( 4 S) + O 2 ( 3 Σ − ) reaction, has been published (Chen et al. 2023).In Chen et al. (2023), an even higher emerged barrier is reported: 12.8 kJ mol −1 at the MRCI(Q)/aug-cc-pV5Z level.In our opinion, these persisting disagreements call for a reexamination of the P( 4 S) + O 2 ( 3 Σ − ) reaction with the aim of obtaining a quantitative description of the PO production under the very different conditions characterizing the Earth's atmosphere and the ISM.
Based on the premises above and despite the modeling based on a large phosphorus chemical network reported in Douglas et al. (2022), in our opinion, the reaction with molecular oxygen still requires a thorough analysis.For this reason, in the present work, we perform a detailed study of the PES governing the P + O 2 reaction both in the doublet and quartet spin states.These can both be accessed following O 2 ( 3 Σ − ) recombination with P( 4 S) (Gomes et al. 2022).In addition, we also investigated the reaction between P( 4 S) and O 2 ( 1 Δ), which is expected to take place on the excited quartet PES.In our study, we put much effort into elucidating the reaction mechanism, thus focusing on a detailed sampling of the reactive PES together with the accurate structural determination of all intermediates and transition states.This led to the characterization of new stationary points and provided new kinetic results, which were found to be in good agreement with the available experimental data (Douglas et al. 2019).We then used the kinetic data estimated in this way to study the role played by the P( 4 S) + O 2 ( 3 Σ − ) reaction in the formation of PO in shocked regions.Our astrochemical model showed that there are two regimes for the chemistry of PO: one characterized by moderate shock speeds (20-30 km s −1 ), where P( 4 S) + O 2 ( 3 Σ − ) dominates, and another one with higher shock speeds (40 km s −1 ), where the reaction P( 4 S) + OH( 2 Π) takes over.The competition between these two reactions directly affects the production of PO, thus making them essential processes to be incorporated in astrochemical models in order to correctly describe the chemistry of phosphorus in the ISM.

Computational Details
The first step of this computational study is the characterization of all the stationary points of the reactive PESs to build the energy profiles.Kinetic calculations are then carried out in the same temperature range used in the experiment (Douglas et al. 2019) (between 187 and 732 K) and in the low-pressure limit.
In the description of the reactive PESs, we denote the doublet and quartet states as Dx and Qx, respectively, where x is an integer number indicating the energy order.Accordingly, D0 is the electronic ground state and Q1 is the lowest quartet state.The reaction between P( 4 S) and O 2 ( 3 Σ − ) can take place on the ground doublet D0 ( 2 A') PES, on the quartet Q1 ( 4 A'), or on the sextet ( 6 A') PES.However, since the sextet state cannot lead to the ground-state products, PO( 2 Π) + O( 3 P), through adiabatic processes, it is not considered in the following.At the same time, the O 2 ( 1 Δ) excited electronic state of the oxygen molecule can react adiabatically only on the excited quartet Q2 ( 4 A") PES.These PESs are sketched in Figure 1.

Electronic Structure Calculations
Since the system under study has a high multiconfigurational character (Gomes et al. 2022), the structures of all the stationary points on the ground doublet (D0) as well as quartet (Q1 and Q2) PESs have been optimized without any constraint using the fully internally contracted CASPT2 (FIC-CASPT2) model (Andersson et al. 1990) in conjunction with the Dunning aug-cc-pVQZ basis set (Dunning 1989;Kendall et al. 1992;Woon & Dunning 1993).In all cases, the CAS wave functions have been obtained without state averaging because all test calculations performed indicated a significant energy gap between the ground and excited states.Since the inclusion of additional tight d functions (+d basis set; Dunning et al. 2001) on the phosphorus atom, which is expected to be important whenever third-row elements are considered, was found to provide no significant improvement, to save computational cost, the standard aug-cc-pVQZ basis set has been systematically employed.The selected active space includes 5 electrons in 4 atomic orbitals for the phosphorus atom (3s 2 3p 3 ) and 12 electrons in 8 molecular orbitals for molecular oxygen , leading to a full valence active space containing 17 electrons in 12 orbitals.All FIC-CASPT2 calculations have been carried out with the ORCA program (Neese et al. 2020) (v.5.0.2).
The nature of all the stationary points has been validated by Hessian evaluations, which have also been used to compute the harmonic zero-point energy (ZPE) corrections.The connectivity between wells and transition states was determined by inspection of the transition vector (i.e., the normal mode corresponding to the imaginary frequency at the TS) and, when needed, by following the intrinsic reaction coordinate (Fukui 1981).
To validate the FIC-CASPT2 approach, selected structures and energies of key stationary points of the investigated PESs have been determined using the MOLPRO CASPT2 and MRCI implementations (Werner et al. 2020).The same approach has also been used to determine minimum energy crossing points (MECPs) between PESs with different spin multiplicities or different symmetry.As detailed in the following, this was required because the reactive PESs were found to be more complicated than expected.MECP structures have been determined at the CASPT2/aug-cc-pVTZ level, while their energies have been refined at the CASPT2/CBS level.This is obtained by extrapolating to the complete basis set (CBS) the CASPT2/aug-cc-pVnZ energies (with n = 3, 4), thereby exploiting the n −3 extrapolation formula (Helgaker et al. 1997).The same approach has been applied to the structural determination of the reactants, the entrance well, and TS1 as well as to investigating the reactivity on the quartet Q2 PES, which is accessed from the O 2 ( 1 Δ) interaction with P( 4 S).

Kinetics Calculations
Pressure-and temperature-dependent phenomenological rate coefficients have been calculated by solving the one-dimensional master equation using the MESS software (Georgievskii et al. 2013).The global rate coefficients have been evaluated in the 187-732 K temperature range and using pressure conditions that make collisional stabilization negligible (P = 1×10 −7 atm).The temperature dependence of the global rate constant has been modeled using the Arrhenius-Kooij formula (Kooij 1893): The α, β, and γ parameters have been determined through nonlinear least-squares fits.Finally, the fitted rate constants have been incorporated in the phosphorus chemical model built by Jiménez-Serra et al. (2018), which has been recently updated in García de la Concepción et al. (2021).

Results
In the following, the reactive PESs for both P( 4 S) + O 2 ( 3 Σ − ) and P( 4 S) + O 2 ( 1 Δ) are described in detail and compared to the available literature.In the discussion, we always refer to ZPEcorrected relative energy values with respect to the reactants in their electronic ground state.Subsequently, the kinetics is addressed, with the emphasis being on the temperature and pressure dependence of the global rate constants.A comparison with previous results is also performed.

Mechanism of the P(
As shown in Figure 1(a), the P( 4 S) + O 2 ( 3 Σ − ) reaction in the ground D0 state (red profile) starts with a barrierless approach of the reactants, which leads to the weak van der Waals prereactive complex, denoted as 2 C_P-O 2 (A').This lies only 2.4 kJ mol −1 below the reactants.The reaction then proceeds toward the 2 TS1_P-O 2 transition state, which lies 8.5 kJ mol −1 above the reactants and corresponds to the formation of the P-O bond.Consequently, the conversion of 2 C_P-O 2 (A') into 2 POO(A') through 2 TS1_P-O 2 (A') governs the rate-determining step of the entire process.As mentioned in the "Computational Details" section, its energy barrier has also been determined at the CASPT2/CBS level, which provides a value of 8.4 kJ mol −1 , in very good agreement with the result above.The 2 POO(A') intermediate, lying 61.0 kJ mol −1 below the reactants, is connected to the deep potential well formed by 2 OPO(A') (−570.7 kJ mol −1 ) through 2 TS2_P-O 2 (A'), which leads to a submerged barrier of 31.9 kJ mol −1 with respect to 2 POO(A').The 2 POO(A') intermediate and the 2 OPO(A') well are also connected through an allowed crossing (AC) between the A' and A" doublet states (D0/D1 AC).The corresponding MECP is denoted as 2 AC_P-O 2 and lies 16.8 kJ mol −1 below 2 TS2_P-O 2 (A').The D0/D1 AC leads to 2 OPO(A") (D1 PES), which is separated from 2 OPO(A') on the D0 ground state by 360.6 kJ mol −1 .Two pathways are then accessible from 2 OPO(A"): internal conversion to 2 OPO(A') or direct decomposition to 2 PO and O( 3 P).As the main reaction pathway of 2 OPO(A') is also its decomposition to 2 PO and O( 3 P), the 2 OPO(A") reactivity was not further analyzed.The last step of the reaction mechanism on the D0 PES is the dissociation of the 2 OPO(A') intermediate to give the bimolecular products, that is, 2 PO and atomic oxygen in the ground electronic state (O( 3 P)), through a barrierless process.The overall process is exothermic, showing a reaction energy of −83.2 kJ mol −1 at the CASPT2/CBS level.Note that we carefully checked that the reaction on the D0 PES proceeds from 2 TS2_P-O 2 (A') to 2 OPO(A') and not directly to the bimolecular products, as previously suggested (Gomes et al. 2022).
The reaction on the A' Q1 quartet electronic PES (blue profile in Figure 1(a)) is also characterized by a barrierless entrance channel, which leads to the formation of the 4 C_P-O 2 (A') prereactive complex lying only 0.8 kJ mol −1 above 2 C_P-O 2 (A').Therefore, in this region, the two PESs are nearly degenerate.Furthermore, the geometries of the two prereactive complexes are very similar, with the P-O distance showing the largest difference: 3.53 and 3.29 Å for 4 C_P-O 2 (A') and 2 C_P-O 2 (A'), respectively.The subsequent reaction step is also very similar to that of the doublet PES in terms of structural parameters.However, the transition state and the intermediate have energies significantly higher than their counterparts on the electronic ground state.4 TS1_P-O 2 (A') lies 20.0 kJ mol −1 above 2 TS1_P-O 2 (A') and leads to a reaction barrier of 30.1 kJ mol −1 . 4POO(A'), with an energy of −16.4 kJ mol −1 , lies 44.6 kJ mol −1 above 2 POO(A').From 4 POO(A'), the reaction proceeds differently on the quartet Q1 PES with respect to the doublet D0 PES.In fact, the last step of the mechanism is the dissociation of 4 POO(A') into the bimolecular products through the 4 TS2_P-O 2 (A') transition state, whose energy barrier is 37.8 kJ mol −1 with respect to 4 POO(A').By analogy with the D0 ( 2 A') PES, the rate-determining step of the reaction on the A' Q1 PES is the first transition state ( 4 TS1_P-O 2 (A')).However, this process should be considerably slower than the former one on the D0 PES because the emerged transition state, which governs the rate-determining step, lies 28.5 kJ mol −1 above the reactants, with this making the process very unlikely under the conditions of interest in the present work.
The potential energy profile of the P( 4 S) ), which takes place on the excited A" Q2 PES, has been investigated at the CASPT2/CBS level.It was found that the entrance channel is barrierless and leads to the formation of the 4 C_P-O 2 (A") prereactive complex, lying 1.5 kJ mol −1 below the energy of the corresponding reactants.The successive step, leading to the formation of 4 POO(A"), requires overcoming a significant energy barrier of 37.7 kJ mol −1 .This indicates that the reaction of P( 4 S) with O 2 ( 1 Δ) is significantly slower than that with O 2 ( 3 Σ − ).
Gomes et al. (2022) investigated the reactive PES at the MRCI (Q)/aug-cc-pV(T+d)Z level on top of stationary points optimized at the CASSCF level (in conjunction with a triple-zeta basis set) and incorporated the ZPE corrections at the same level of theory.Interestingly, despite the high level of theory, their PES does not show any prereactive complex for both the doublet and quartet PESs.Furthermore, they found that the second transition state of the doublet PES ( 2 TS2_P-O 2 (A') in our PES) leads directly to the products through the dissociation of the O-O bond without passing through the second intermediate (denoted as 2 OPO(A') in Figure 1(a)).The impossibility of optimizing a smooth saddle point connecting 2 POO to 2 OPO pointed out in Gomes et al. (2022) might be related to their use of geometries optimized at the CASSCF level in place of the more physically sound CASPT2 structures employed in the present work.This difference in the doublet PES is reflected in the relative energy of 2 TS2_P-O 2 (A'), which is −29.1 kJ mol −1 in our study, to be compared to a value of −40.3 kJ mol −1 in Gomes et al. (2022).On the other hand, the relative energy of the emerged barrier of the doublet PES obtained by Gomes et al. (2022) (10.1 kJ mol −1 ) is close to that computed in the present work (8.5 kJ mol −1 , which becomes 10.9 kJ mol −1 when referring to the prereactive complex), and the same applies to the other stationary points characterized in Gomes et al. (2022).In Chen et al. (2023), despite the high level of theory employed (MRCI(Q)/aug-cc-pV5Z), no prereactive complex was identified as well; this is probably related to the ab initio grids being not sufficiently dense in the description of the P + O 2 approach.Interestingly, Chen et al. (2023) found the reaction to proceed from the intermediate denoted as 2 OPO(A') in Figure 1(a) to the products via a linear OPO transition state.However, we were not able to find any linear OPO structure connecting 2 OPO(A') to PO( 2 Π) + O( 3 P).
In view of the sensitivity of the system to the level of theory employed for the characterization of structures and energetics, additional calculations have been carried out for the rate-determining barrier (thus implying calculations for 2 TS1_P-O 2 (A') and reactants).These are summarized in the Appendix (Table A1).In particular, the geometries were also optimized at the CASPT2(17e,11o)/aug-cc-pVTZ and MRCI (17e,11o)/aug-cc-pV(T+d)Z levels.For 2 TS1_P-O 2 (A'), consistent structural results have been found at both levels of theory: the P-O, O-O, and P-O-O distances and angles are 2.26 Å, 1.24 Å, and 118.3°at the CASPT2 level and 2.29 Å, 1.23 Å, and 118.4°at the MRCI level, respectively.Energies were then determined at the MRCI(17e,11o)/CBS level, with extrapolation from the aug-cc-pVTZ and aug-cc-pVQZ levels and using the Davidson correction.The ZPE-corrected relative energy of 2 TS1_P-O 2 determined at the MRCI/CBS level (on top of CASPT2/aug-cc-pVTZ geometries) is 11.2 kJ mol −1 , thus about 3 kJ mol −1 higher than its CASPT2/CBS counterpart (using CASPT2/aug-cc-pVTZ geometries as well).Such an energy difference is within the uncertainties of the adopted computational protocols, which is relatively large considering the high multireference character of the system.For what concerns the comparison with previous calculations, it is interesting to note that we were not able to reproduce the 10.1 kJ mol −1 barrier obtained in Gomes et al. (2022) at the MRCI/aug-cc-pV(T+d)Z level of theory.Indeed, employing the same level and active space used by Gomes et al. (2022), we determined a relative energy of 14.4 kJ mol −1 when using the MRCI/aug-cc-pV(T+d)Z geometry, and 7.6 kJ mol −1 when resorting to the CASSCF/aug-cc-pV(T+d)Z geometry (i.e., the same level used in Gomes et al. 2022).These results indicate a quite relevant impact of the structure on the 2 TS1_P-O 2 (A') energy barrier.Moving to the quartet Q1 PES, it should be noted that the relative energy of the second transition state, that is, 4 TS2_P-O 2 (A'), obtained in this work differs noticeably from that derived by Gomes et al. (2022): 21.4 kJ mol −1 versus 8.4 kJ mol −1 .

Kinetics of the P(
Channel-specific rate constants for the investigated system have been determined by integrating the one-dimensional master equation employing the PES reported in Figure 1(a) (red profile), but without including the D0/D1 crossing.The reason behind this choice is twofold: the MESS code does not include a reliable model to describe nonadiabatic transitions and, above all, the contribution of this pathway to the system reactivity is minor.We have in fact performed a preliminary investigation of the competition between the two reaction channels leading to 2 OPO(A'), that is, the pathways through 2 TS2_P-O 2 (A') and 2 AC_P-O 2 , using the MC-RRKM stochastic code, which implements nonadiabatic transition state theory (Barbato et al. 2009;Vanuzzo et al. 2021).It was found that, despite the higher energy, the reactive flux through 2 TS2_P-O 2 (A') is larger than that through 2 AC_P-O 2 by a factor of at least 2 in all the investigated temperature and pressure conditions.As the calculated rate constants are insensitive to a change in the formation rate of 2 OPO(A') up to a factor of 2 and the reactive flux is higher through 2 TS2_P-O 2 (A'), we computed the phenomenological rate constants for the P( 4 S) + O 2 ( 3 Σ − ) reaction considering only the D0 (A') PES.
As mentioned in Section 3.1, the rate-determining step of the P( 4 S) + O 2 ( 3 Σ − ) reaction is the formation of 2 POO(A') through 2 TS1_P-O 2 (A'), which is the only emerged stationary point of this reaction path and determines a barrier of 10.9 kJ mol −1 (with respect to the prereactive complex at the ZPE-corrected FIC-CASPT2/aug-cc-pVQZ level).For this reason and because of the low-pressure conditions enforced in the simulations, the rate constant decreases when decreasing the temperature, thus following an Arrhenius behavior in the whole range of temperatures studied.The temperature dependence is well reproduced by the Arrhenius-Kooij model, with the fitted parameters being α = 4.0 × 10 −12 cm 3 molecule −1 s −1 , β = 0.9, and γ = 814.0K (see Equation (1)).The resulting fitting function is shown in Figure 2 (dashed red line).
The rate constant calculated for this reaction varies from 3.5 × 10 −14 cm 3 molecule −1 s −1 at 187 K to 2.9 × 10 −12 cm 3 molecule −1 s −1 at 732 K.The computed rate constants are in good agreement with those experimentally determined by Douglas et al. (2019)   computed rate constant is a little more than a factor of three lower than the experimental value (1.4 × 10 −13 versus 3.5 ×10 −14 cm 3 molecule −1 s −1 ; see Douglas et al. 2019), while at 732 K the experimental rate is overestimated by about a factor of two.However, it should be noted that the hightemperature experimental rate constant is smaller than that measured by the same authors at 583 K, thus leading to an inconsistent Arrhenius behavior.The rate constant we computed at 298 K, 2.6 × 10 −13 cm 3 molecule −1 s −1 , is not only in quite good agreement with the experimental value from Douglas et al. (2019) (2.7 × 10 −13 cm 3 molecule −1 s −1 at 291 K), but it is also consistent with the experimental determinations (at 298 K) by Henshaw et al. (1987) (1.1 × 10 −13 cm 3 molecule −1 s −1 ), Clyne & Ono (1982) (0.99 × 10 −13 cm 3 molecule −1 s −1 ), Husain & Slater (1978) (2.1 × 10 −13 cm 3 molecule −1 s −1 ), and Husain & Norris (1977) (2.0×10 −13 cm 3 molecule −1 s −1 ).In view of the agreement with the available experimental data (see Table 1), it is reasonable to conclude that the rate constant computed in the present work has a confidence interval of about 2. Unfortunately, the single experimental rate constant available at temperatures lower than 291 K (Douglas et al. 2019) shows the largest disagreement with our estimate.Therefore, new experimental determinations of rate constants at low temperatures would be particularly valuable.Should the experimentaltheoretical discrepancy be confirmed, then additional theoretical calculations, focused on the determination of anharmonic or multidimensional quantum tunneling contributions to the rate constant, might be needed.This is, however, beyond the scope of the present work, as we do not expect this latter process to play a major role in the chemistry of PO in cold molecular dark clouds at T ∼ 10-20 K, since both theory and experiment indicate that the rate constant monotonically decreases when decreasing the temperature.
In Figure 2, in addition to the experimental rate constants by Douglas et al. (2019) (black dots) and our results (red dashed line), the Arrhenius-Kooij fit of the rate constants from Gomes et al. (2022) (green solid line) is also shown.These were calculated by the authors by integrating the master equation (using the same MESS code employed in the present work) over an ab initio PES determined at the MRCI(Q) level, as explained in Section 3.1 (Gomes et al. 2022).In this respect, we note a potential typo in the β Arrhenius-Kooij parameter reported by Gomes et al. (2022) (−1.66) because we could reproduce the rate constants reported in their paper only when employing a positive value (1.66).From an inspection of Figure 2, it is noted that the results from Gomes et al. (2022) underestimate the temperaturedependent experimental rate constants measured by Douglas et al. (2019), while they are in good agreement with the room temperature measurements performed by Henshaw et al. (1987) and Clyne & Ono (1982).A nearly constant underestimation with respect to our ZPE-corrected FIC-CASPT2/aug-cc-pVQZ plot is also observed.This can be traced back to the different relative energy of 2 TS1_P-O 2 : 10.1 kJ mol −1 (Gomes et al. 2022) versus 8.5 kJ mol −1 (present work).The situation changes when moving to our ZPE-corrected MRCI/CBS results (using CASPT2/augcc-pVTZ geometries): 2 TS1_P-O 2 lying at 11.2 kJ mol −1 above the reactants leads to a further underestimation of the computed rate constants with respect to the experimental data.The resulting temperature dependence is reported in Figure 2 as a blue dashed line, and the corresponding Arrhenius-Kooij parameters are α = 4.2 × 10 −12 cm 3 molecule −1 s −1 , β = 0.9, and γ = 1144.5K.In our opinion, also based on the thorough investigation of the 2 TS1_P-O 2 relative energy summarized in Table A1, the ratedetermining barrier cannot be estimated with an accuracy better than 3 kJ mol −1 .Therefore, the FIC-CASPT2/aug-cc-pVQZ and MRCI/CBS results of Figure 2 can be considered as providing the confidence interval for the rate constant and its temperature behavior.The kinetic results shown in Figure 2 as a yellow solid line are from Chen et al. (2023).They were obtained with timedependent wave packet methods and exhibit even larger deviations from experiment than our MRCI/CBS results and those by Gomes et al. (2022) (green solid line), despite the improved level of theory, at least in the structural determination.Finally, in Figure 2, the Arrhenius behavior obtained in Douglas et al. (2019) by fitting the RRKM parameters to the experimental data is reported as an orange line.This is not able to reproduce the experiment within the given uncertainties, thus reinforcing the need of additional experiments, in particular at low temperatures.
As already mentioned in Section 3.1, the D0 and Q1 reactive PESs sketched in Figure 1(a) are almost degenerate in the region close to the prereactive complexes ( 2 C_P-O2(A') and 4 C_P-O2(A')).An MECP between the two PESs has therefore been searched using the automatic procedure implemented in EStokTP (Cavallotti et al. 2019) and described in Cavallotti et al. (2020).As mentioned in the "Computational Details" section, the calculations have been performed at the CASPT2/aug-cc-pVTZ level using a full valence active space.The MECP structure has been found at O-O and O-P bond lengths similar to those of the quartet entrance complex (i.e., 1.22 and 3.68 Å, respectively), but for an O-O-P angle of 158°.The spin-orbit coupling (SOC) between the two PESs has been computed using a CASSCF wave function and a full valence space and exploiting the Breit-Pauli Hamiltonian, as implemented in the MOLPRO program.The computed SOC is quite small (0.1 cm −1 ), which is a clear diagnostic of a negligible rate for the intersystem crossing between the two PESs.Furthermore, the master equation simulation clearly shows that for low pressures there is no stabilization of 2 C_P-O2(A') at any temperature.Indeed, at a temperature of 187 K, the 2 C_P-O2(A') complex is stabilized by molecular collisions only at pressures above 1000 atm.At such high pressures, reachable only in situations such as collisions of meteorites with atmospheres or the surface of rocky exoplanets, the kinetics of the P( 4 S) + O 2 ( 3 Σ − ) reaction is expected to change noticeably.Under those conditions, since both the energy and the geometry of the 2 C_P-O2(A') and 4 C_P-O2(A') prereactive complexes are very similar, the Boltzmann population of the quartet complex ( 4 C_P-O2(A')) would increase, thus decreasing the total rate constants due to the participation of the quartet Q1 PES.Thus, different physical conditions would lead to different reaction rates.Given that the 2 C_P-O 2 (A') prereactive complex is not stabilized when the physical conditions considered in our simulation apply and that our kinetic calculations for the doublet spin state are in very good agreement with the most recent experimental results (Douglas et al. 2019; see Table 1), there is no evidence of the participation of the Q1 PES in the system reactivity at the temperatures considered herein and in previous works.However, due to the high temperatures reached in shocked regions, we also computed the rate constants for the P( 4 S) + O 2 ( 3 Σ − ) reaction in the Q1 state in order to evaluate its impact on our astrochemical model (see below).The Arrhenius-Kooij parameters obtained for this reaction in the 187-732 K temperature range are α = 1.1 × 10 −11 cm 3 molecule −1 s −1 , β = 1.3, and γ = 2988.0K.The high-energy emerged barrier for the Q1 PES ( 4 TS1_P-O 2 (A") lying at 28.5 kJ mol −1 ) is reflected in the large γ parameter, which leads to a different kinetic behavior with respect to that found for the D0 state.
To sum up, at low temperatures, the fastest reaction rate is obtained for the reaction occurring on the D0 PES, but when raising the temperature, the contribution of the reaction on the Q1 PES rapidly increases until becoming similar to the D0 one at 1400 K and dominant above 2000 K.This is consistent with the results by Gomes et al. (2022).Since the D0 and Q1 PESs are degenerate in the reactant asymptote and their crossings are ineffective, the production of PO depends on the sum of the rate constants in the D0 and Q1 states.The fitting of the new rate constants to the Arrhenius-Kooij expression leads to α, β, and γ parameters of 1.5 × 10 −12 cm 3 molecule −1 s −1 , 1.8, and 516.5 K, respectively.In the temperature range considered in this work, the new rate constants are very similar to those calculated for the D0 PES, displaying the most significant deviation at 732 K with a difference of 0.6 × 10 −12 cm 3 molecule −1 s −1 .On the other hand, the rate constant for the P( 4 S) + O 2 ( 1 Δ) reaction (Q2, Figure 1(b)) was not determined in this work, as we expect it to have an entirely negligible impact on the P( 4 S) reactivity under the typical ISM conditions.

Discussion: Astrophysical Implications
Molecular oxygen has been unambiguously detected toward the ρ Oph A and Orion star-forming regions (see Goldsmith et al. 2002;Liseau et al. 2012;Chen et al. 2014).Its abundance ranges from some 10 −8 to a few 10 −7 , and thus it should be considered as an important oxidizing reactant of other chemical species.Due to the relevance of PO as a prebiotic molecule, the reaction between P( 4 S) and O 2 ( 3 Σ − ) has been widely studied (see Douglas et al. 2019;Gomes et al. 2022;Chen et al. 2023, and references cited therein).However, the P( 4 S) + O 2 ( 3 Σ − ) reaction presents an emerged barrier and this process is not viable under the typical low temperatures found in quiescent and cold molecular dark clouds.This process is, thus, not expected to contribute to the formation of PO in these cold objects.
In contrast, this reaction might contribute to the production of PO in regions of the ISM affected by energetic phenomena such as shocks in molecular outflows.It has been found that the abundance of PO with respect to H 2 in shocked regions is relatively high (∼1×10 −9 ).In addition, the measured PO/PN abundance ratio lies above 1 (Aota & Aikawa 2012;Lefloch et al. 2016;Rivilla et al. 2016Rivilla et al. , 2018Rivilla et al. , 2020;;Bergner et al. 2019;Bernal et al. 2021).By performing astrochemical modeling using the newly calculated rate constant of the P( 4 S) + OH( 2 Π) reaction, García de la Concepción et al. (2021) have recently shown that the PO/PN ratio 1 measured toward the L1157-B1 outflow region could be reproduced thanks to the inclusion of this reaction rate in the chemical network of phosphorus.Here we study the effects of the newly calculated rate constant of the P( 4 S) + O 2 ( 3 Σ − ) reaction on the chemistry of shocked gas with temperatures considerably higher than those found in cold and quiescent molecular dark clouds (see Section 4.1).
We note that the reaction between P( 4 S) and O 2 ( 1 Δ) is also an activated process.Thus, we do not expect it to contribute to the formation of PO in molecular clouds.In addition, it should be considered that, as already mentioned, in the ISM the majority of molecular oxygen is expected in its electronic ground state.Even in shocked regions where the temperature can reach 4000 K, the population of O 2 ( 1 Δ) is negligible with respect to that of O 2 ( 3 Σ − ).Furthermore, although in molecular outflows there are regions of shocked gas, the density is still too low to allow an increase in the population of O 2 ( 1 Δ) by molecular collisions, so its reaction will not contribute to the formation of PO.

Chemical Modeling in a Shocked Region
Based on what was introduced above, in this section, we evaluate the effects of the P( 4 S) + O 2 ( 3 Σ − ) reaction in the phosphorus chemistry of shocked gas, while we do not consider the P( 4 S) + O 2 ( 1 Δ) reaction.In particular, we carry out the modeling of the phosphorus chemistry of an outflow shocked region with and without considering the PO formation through the P( 4 S) + O 2 ( 3 Σ − ) reaction.For this purpose, we consider the C-type shock models that have been explored to explain the molecular abundances measured toward the B1 shocked region in the L1157 molecular outflow (see Viti et al. 2011;Holdship et al. 2016;Lefloch et al. 2016).
For the modeling of the phosphorus chemistry, we use the UCLCHEM chemical code (Holdship et al. 2017) and the phosphorus chemical network built by Jiménez-Serra et al. (2018) and recently improved by García de la Concepción et al. (2021).This chemical network has been updated by incorporating the rate constants of the P( 4 S) + O 2 ( 3 Σ − ) reaction (see Table 1; FIC-CASPT2/aug-cc-pVQZ results have been employed because they give the best agreement with experimental data).UCLCHEM is run in three phases: Phase 0 simulates the chemistry of a diffuse molecular cloud for 10 6 yr assuming an n (H) density of 10 3 cm −3 and a temperature of 20 K.In phase 1, the cloud undergoes freefall collapse at a constant temperature of 10 K until the final density of 2 × 10 5 cm −3 is reached.Finally, phase 2 simulates the physical processes associated with the passage of a magnetohydrodynamic (C-type) shock, which are commonly found in molecular outflows.
For the physical structure of the C-type shock, we employ the parametric approximation developed by Jiménez-Serra et al. (2008) and implemented in UCLCHEM (Holdship et al. 2017).Table 2 reports the input parameters used to model the physical structure of the L1157-B1 shock, where n(H) refers to the initial volume density of the gas (2 × 10 5 cm −3 ), v s is the shock speed, set to 20, 30, and 40 km s −1 , T n,max is the maximum temperature reached by the neutral gas in the postshock region (800, 2000, 4000 K, respectively), B 0 is the magnetic field strength (450 μG), and t sat is the time at which the majority of the ice mantles of dust grains are released into the gas phase due to sputtering (5.7, 4.4, and 4.6 yr, respectively; see Tables 4 and 5 in Jiménez-Serra et al. 2008).
We note that the model that best reproduces the molecular abundances measured toward the L1157-B1 shock region is the one with a gas density of n(H) = 2 × 10 5 cm −3 and shock speed v s = 40 km s −1 (Viti et al. 2011;Lefloch et al. 2016;García de la Concepción et al. 2021).However, in this work, we also explore the other two shock speeds, v s = 20 and 30 km s −1 , to illustrate the effects of the inclusion of the P( 4 S) + O 2 ( 3 Σ − ) reaction on the chemistry of PO at different T n,max temperatures of the shocked gas.
Figure 3 shows the results of the models of the phosphorus chemistry for the C-shock parameters reported in Table 2.The  2. These shock parameters are those that have been explored in the past to reproduce the molecular abundances measured in the L1157-B1 shocked region (Viti et al. 2011;Holdship et al. 2016;Lefloch et al. 2016).Dashed red lines show the results for the model without the new rate constants for the formation of PO.The abundance of O 2 is the same in both models.Solid red lines refer to the results from the model with the new PO rate constants calculated in this work for the P( 4 S) considering the sum of the rate constants in the D0 and Q1 states.
abundances of PO obtained with the newly calculated rate constants of the P( 4 S) + O 2 ( 3 Σ − ) reaction are shown with red solid lines, while the results obtained without the new data are represented with red dashed lines.Note that, in Figure 3, we show the amount of PO produced by the title reaction considering the sum of the rate constants in the D0 and Q1 states (see the end of Section 3.2).The model in the top panel is calculated for a shocked speed of 40 km s −1 , which is the same used for the simulation reported in García de la Concepción et al. (2021).The middle panel reports the model for a shock speed of v s = 30 km s −1 , and the lower panel for the model with vs = 20 km s −1 .
The model with v s = 40 km s −1 does not present any noticeable difference with respect to the previous simulation reported in García de la Concepción et al. (2021).In this model, the abundances of H 2 O, P, PN, and O 2 increase at the beginning of the shock (at about time ∼5 yr) due to the full injection of the icy mantles from dust grains.At ∼10 yr, the abundance of OH drops, whereas the abundance of O 2 remains constant.In this period, the production of PO is ruled by the P( 4 S) + O 2 ( 3 Σ − ) reaction.In contrast, the P + OH process does not significantly contribute to the formation of PO (see dashed red line in Figure 3, top panel).At longer times (between 20 and 100 yr), and due to the high temperatures reached in the postshocked gas (∼4000 K), O 2 and H 2 O are destroyed in the gas phase, while OH is formed in large amounts.This is the shock regime where the P( 4 S) + O 2 ( 3 Σ − ) reaction does not influence the chemistry of PO because the latter is formed almost exclusively through the P + OH reaction, as reported in previous simulations (García de la Concepción et al. 2021).Overall, the effects of the P( 4 S) + O 2 ( 3 Σ − ) reaction on the total production of PO in a shock with V s = 40 km s −1 are minimal.
The effects of this reaction, however, are more important for the models with shock velocities v s = 30 and 20 km s −1 (see Figure 3, middle and lower panels, respectively).For lower shock speeds, neither O 2 nor H 2 O get destroyed in the shock.This situation prevents the formation of large amounts of OH in the postshocked gas, and, thus, the P( 4 S) + O 2 ( 3 Σ − ) reaction governs the formation of PO at all timescales in the postshocked gas.The inclusion of the new kinetic data in our model with v s = 20 km s −1 for the P( 4 S) + O 2 ( 3 Σ − ) reaction increases the abundance of PO in the shock by more than one order of magnitude, impacting strongly on the predicted abundance of PO at low shock velocities (see red solid line).Note that the higher the shock velocity, the higher the T n,max and the more abundant the PO in the postshocked gas (see Figure 3), as a result of the increase in the rate of the P + O 2 reaction at the higher gas temperatures of the model with v s = 30 km s −1 .
To summarize, the formation of PO is governed in shocks by the P + OH and P + O 2 reactions.However, their relative efficiency depends on the velocity of the shock.While at low shock speeds (v s 30 km s −1 ) the PO formation is ruled by the P + O 2 reaction, at higher shock velocities (v s 40 km s −1 ) its formation occurs mainly through the P + OH reaction.These results show that both processes are essential in our understanding of the production of the prebiotic PO molecule in the ISM.

Conclusions
In this work, we have revisited the reactions of P( 4 S) with O 2 ( 3 Σ − ) and O 2 ( 1 Δ) from a theoretical point of view.We have confirmed that multireference methods are needed for a correct description of the geometry of the stationary points, energetic, and reaction mechanism, and we have pointed out that additional experiments at low temperatures would be valuable.
The key step of the P( 4 S) + O 2 ( 3 Σ − ) reaction occurring in the D0 ground state is governed by the first transition state, which corresponds to an emerged barrier (with respect to the reactants) of 8.5 kJ mol −1 and is thus the rate-determining step of the whole process.The quartet Q1 PES presents a similar behavior but with a higher barrier to be overcome.D0 and Q1 prereactive complexes are nearly degenerate and an MECP has been found.However, the intersystem crossing rate is expected to be small because of a negligible SOC.Simulations demonstrated that the prereactive complexes are never stabilized in low and medium pressure regimes, and thus they do not have any influence on the kinetics of the P( 4 S) + O 2 ( 3 Σ − ) reaction, at least in the conditions of interest to this work.In the reaction with O 2 ( 1 Δ), we also note the formation of the prereactive complex through a barrierless association followed by the rate-determining transition state, which lies 37.7 kJ mol −1 above the complex.However, because of the high energy involved, we did not investigate the reactive PES beyond the formation of 4 POO(A").
The P( 4 S) + O 2 ( 3 Σ − ) reaction in the D0 and Q1 states can occur only in high-temperature interstellar environments such as shocks in molecular outflows, because the emerged barriers make the reactions too slow at low temperatures, with this being especially true for the Q1 state, which presents a barrier ∼3 times higher than that of the D0 state.Our chemical model demonstrated that in shocked regions where the gas velocity is high (40 km s −1 ), the production of PO is not influenced by the P( 4 S) + O 2 ( 3 Σ − ) reaction because the high temperatures (∼4000 K) enable the formation of OH from the destruction of H 2 O.As a consequence, the abundance of PO in high-velocity shocks is ruled by the reaction between P( 4 S) and OH.In contrast, in regions where shock speeds are slower (20-30 km s −1 ), the gas temperatures attained in the shock are not high enough to allow the destruction of H 2 O and O 2 .Consequently, without a large production of OH, it is the P(

Figure 3 .
Figure 3. Abundances of P, PO, PN, OH, O 2 , and H 2 O predicted for phase 2 of our model and employing the shock parameters collected in Table2.These shock parameters are those that have been explored in the past to reproduce the molecular abundances measured in the L1157-B1 shocked region(Viti et al. 2011;Holdship et al. 2016;Lefloch et al. 2016).Dashed red lines show the results for the model without the new rate constants for the formation of PO.The abundance of O 2 is the same in both models.Solid red lines refer to the results from the model with the new PO rate constants calculated in this work for the P( 4 S)+ O 2 ( 3 Σ − ) → 2 PO + O( 3 P)considering the sum of the rate constants in the D0 and Q1 states.
(see Figure 2, black dots, and Table 1), with the exception of the highest and lowest temperature measurements.At the lowest temperature, the