General Comments

The manuscript presents a computational investigation into the reaction dynamics of the P(⁴S) + O₂ → PO + O system. Several thousand high level ab initio energy calculations, carried out using the multireference configuration interaction (MRCI(Q)) method with an aug-cc-pV5Z basis set, were carried out to map out the PES. These energy points were then fit using the combined-hyperbolic-inverse-power-representation (CHIPR) method in order to construct a global PES, from which features of the PES such as global and local minima and transitions states (TSs), as well as minimum energy pathway for the reaction, were identified. Using the CHIPR PES, rate coefficients for the reaction were calculated using the quasi-classical trajectory (QCT) and time-dependant wave-packet (TDWP) methods. Rate coefficients for the reaction of vibrationally excited O₂ with P(⁴S) were also calculated. 

The reaction P(⁴S) + O₂ → PO + O is important to both astrochemistry, where it is likely an important source of PO in the interstellar medium, and to atmospheric chemistry, being the first step in the oxidation of P atoms from the ablation of interstellar dust particles. As such, accurate rate coefficients over a wide range of temperatures are required to accurately model the P chemistry in these environments. There have been two previous computational studies of the reaction, by Douglas et. al. (2019) who calculate the surface at the B3LYP level and use transition state theory (TST) to fit the surface to experimental data, and by Gomes et. al. (2022), who calculate the surface at the MRCI(Q) level of theory, and also use TST to predict temperature dependant rate coefficients (I believe up to around 1000 K). The justification for the current work is twofold; firstly that the B3LYP method employed by Douglas et. al. (2019) is known to underestimate barrier heights, and secondly the TST theory employed by both studies may not provide accurate rate coefficients due to unincluded non-statistical effects, making the case for a dynamical study to be carried out.

The work presented in the manuscript is of a high quality, and the results and conclusions are generally well explained. I have several minor comments which I list below, as well several technical corrections I have spotted.

Specific Comments

Line 256 onwards and Figure 4 – I don’t fully understand what the x and y axis on Figure 4 physically represent. Some more detail on how to interpret the relaxed triangular contour plot and what physical dimensions the x and y axis relate to would be good.

Line 322 – ‘The system then evolves through TS4, the OPO isomer GM and the linear transition state TS2, accompanied by the progressively open θ.’ Do you mean to say that as you progress from the GM to TS2 the OPO angle increases?  I think this needs to be made clearer.

Line 331 – you describe the secondary elevation after 0.91 eV (stage 3) as probably due to the opening of a new entrance channel, with the P atom crossing over the second-order saddle point SP2 to reach the GM directly. I’m guessing this refers to the contour plot in Figure 2c, which I think you should refer to here. However, it appears that when the P atom approaches the O₂ at an angle of 90° (along the mid-perpendicular?) it reaches the LM rather than the GM. Is this the case? I think this requires slightly more explaining.

Line 338 – the last sentence of this paragraph I think need more explanation (i.e. how rotational excitation creates a barrier / increases the threshold energy).

Around line 343 – the resonances you are referring to are in stage 1 in Figure 6?  I think state this for clarity.

Around line 354 onwards – This sentence is a little confusing as you're initially referring to the absolute ICS values (for which the TDWP values are greater), and then later the threshold values (for which the QCT value is greater). Maybe split into two. When you state that the QCT ICS values are less than the TDWP ICS values, I think you need to state that this is for v=0. And when talking about the threshold values, I think maybe state what the threshold values are, as it is hard to read them off of Figure 7. Also you state that for stages 1 and 3, the reactivity increases with increasing vibrational excitation.  I can see this for stage 1, i.e. that at low collision energies the ICS is higher for v=1 and v=2 than for v=0. However, I’m insure where stage 3 begins, as in Figure 7 at higher collision energies the ICS for v=1 and v=2 is always lower than for v=0. You also state the threshold tends to decrease for increasing vibrational excitation, I can see this in Figure 7, but again it might be useful to give the threshold energies as they are difficult to read off the Figure.

Around line 380 – you state that ‘the experimental rate constants include other possible processes with excited states of O₂, O, and PO, that can be responsible for the depletion of P(⁴S).’ In the more recent study by Douglas et. al. (2019), 248 nm light is used, at which O₂ has as cross section of < 1e-24 cm² molecule^-1, suggesting interference from excited states of O₂ or O atoms are unlikely. However, there may be some O atoms and PO molecules produced from the reaction of P(²P, ²D) + O₂ present, although these would be present at very low concentrations compared to the O₂ co-reagent. I think is it safer to say the experimental results may be suffering from secondary chemistry, rather than stating that they definitely are.

Line 396 – This sentence suggests that both the Douglas et. al. (2019) results and the Gomes et. al. (2022) results are both fitted from calculated data below 1000 K (when the earlier work is fitted to experimental data).

Around line 405 – Here you discuss the barrier heights of TS1 from different levels of theory. You go on to mention that the barrier height used by Gomes et. al. (2022) is reduced from 0.137 eV down to 0.105 eV when including zero-point energy (ZPE), and it is this lower value that is used in the TST study to calculate rate coefficients. I’m not overly familiar with the QCT and TDWP methods for calculating rate coefficients, but I assume you do not need to include ZPE in your calculations as this is already accounted for? You do mention later on that the relatively small QCT rate constants at low temperature are due to the ZPE problem, whereas the TDWP rate constants are faster, with the threshold agreeing well with the 0.105 eV barrier height that includes ZPE. Does this suggest the TDWP rate coefficients take into account ZPE, while the QCT ones do not (and so a ZPE adjustment is required)?

Around line 425 – your results nicely show that above 3000 K the ro-vibrational excitation of the O₂ enhances the rate coefficient. However, your predicted rate coefficients at high temperatures are still lower than those of Gomes et. al. (2022) who also consider the ⁴A state in their calculations. Do you need to consider this state in your calculations too, and if not does this suggest your rate coefficients at very high temperatures are underestimated?

Line 433 – you give the parameterized rate coefficients in the supplemental, I think it would be good to have these in the main paper for ease of astro/atmospheric modellers requiring these rate coefficients.

I have another more general comment. The title of your paper includes ‘implications for atmospheric modelling’, however in the manuscript you don’t actually mention what the implications of the new rate coefficient are, say for Earth’s atmosphere or the interstellar medium. Thinking of implications for atmospheric modelling, you also don’t mention at what pressures your new rate coefficients are applicable. I’m assuming both the QCT and TDWP methods give low pressure limiting (or zero pressure) rate coefficients, which will be applicable to the ISM and upper planetary atmospheres where the pressure is low?

Technical Corrections

Line 33 – change planet to planets

Line 46 – change container to reservoir

Line 49 – add ‘the’ – ‘in the ISM’

Line 55 – strictly speaking, the rate coefficients for P(²D, ²P)  + O₂ are only loss rates, as the branching ratio between chemical reaction forming O + PO, and electronic relaxation forming P(⁴S) + O₂ is not known.

Line 78 / 79 – add ‘the’ – ‘that the potential energy surface’

Line 81 – add ‘of’ – ‘such as the reaction or non-reaction of collisions, and’

Line 169 – remove ‘the’ (to construct a global PES of PO₂)

Line 239 – change display to displays

Line 241 – rephrase sentence – e.g. change to ‘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’

Line 246 – maybe change evolved to reached – and add the – ‘both TS3 and TS4 can be reached from the LM’

Line 261 – change possible to able, change to (at end of sentence) to the – ‘At high collision energies, the P at is able to cross the C_2V barrier SP2 and reach the GM directly, …’

Line 280, equation 9 – change P+H₂ to P+O₂.

Line 298 – add ‘the’ – ‘is the Jacobi for of the CHIPR PES’

Line 337 – descending

Line 344 – add which – ‘bound and quasi-bound states which exist in the LM and GM potential wells’

Line 347 – remove ‘a’ – ‘resulting in plenty of sharp’

Line 406 – change ‘to’ to ‘in’ – ‘and then used in the TST study’

Line 417 – add the – ‘At high temperatures, the previous two …’

Line 444 – remove s from channels and change to reasonably – ‘The long-range interactions, diatomic potentials, and dissociation energies of each asymptotic channel are reasonably reproduced.’

Line 458 – remove ‘the’, and maybe change reaction to reactive? – ‘can be used for molecular simulations of reactive or non-reactive collisions and photodissociation of the PO₂ system in atmospheres and the interstellar medium.’

General comments: In this work, Chen, Qin et al performs a systematic and exhaustive quantum chemistry study on the reaction of the collision reaction between P and molecular oxygen giving rise to O and PO with relevance in the Earth’s atmospheric chemistry of phosphorus and its interstellar chemistry. While previous theoretical studies focused on providing “static information”, in this case, the authors performed both semi-classical and quantum dynamics simulations to provide with time-dependent properties. The quality of the study is high and the methods and techniques are the most state-of-the-art ones for this type of problems. The presentation of the manuscript is also very clear, systematic and precise.

Specific comments: Even though in the title is it stated “… implication for atmospheric modelling”, there is no much discussion on how the new data obtained from the dynamics simulations would impact the chemistry in the atmosphere and interstellar media of the studied system. More details and illustrative examples would help to enhance the scientific significance of the work.

Technical corrections:

• On page 4, equation (1), define R1, R2 and R3.
• On page 6, line 169, change “…to construct a global the PES…” by “…to construct a global PES…”
• In relation to Table 2, please define more clearly the meaning of “ascending ordered energies” and the meaning of the numbers.
• In Fig. 2, it would help the reader if Eh is defined.
• On page 8, line 210, change “Fig. 4(a)” by “Fig. 2(a)”.
• On page 8, line 211, rewrite this sentence since the O atom dissociation indeed requires energy and therefore it is not “barrierless”.
• On Fig. 4 caption, add at the end of the first sentence “(see definition in the text)” referring to the hyperspherical coordinates.
• In the last paragraph of Conclusions, re-formulate the grammar of the sentence “It can also be a reliable component for constructing other larger molecular systems containing PO2, such as PO3 and HPO2 correspond to 460 the reactions R2 and R3 for generating H3PO3 in the Earth's atmosphere.”
• It would also help the reader if Jacobian coordinated are briefly defined.