The chemistry of phosphorus-bearing molecules under energetic phenomena

For decades, the detection of phosphorus-bearing molecules in the interstellar medium was restricted to high-mass star-forming regions (as e.g. SgrB2 and Orion KL) and the circumstellar envelopes of evolved stars. However, recent higher-sensitivity observations have revealed that molecules such as PN and PO are present not only toward cold massive cores and low-mass star-forming regions with PO/PN ratios>1, but also toward the Giant Molecular Clouds in the Galactic Center known to be exposed to highly energetic phenomena such as intense UV radiation fields, shock waves and cosmic rays. In this paper, we carry out a comprehensive study of the chemistry of phosphorus-bearing molecules across different astrophysical environments which cover a range of physical conditions (cold molecular dark clouds, warm clouds, hot cores/hot corinos) and are exposed to different physical processes and energetic phenomena (proto-stellar heating, shock waves, intense UV radiation and cosmic rays). We show how the measured PO/PN ratio (either>1 as in e.g. hot molecular cores, or<1 as in UV strongly illuminated environments) can provide constraints on the physical conditions and energetic processing of the source. We propose that the reaction P + OH -->PO + H, not included in previous works, could be an efficient gas-phase PO formation route in shocks. Our modelling provides a template with which to study the detectability of P-bearing species not only in regions in our own Galaxy but also in extragalactic sources.


INTRODUCTION
Phosphorus is considered as one of the main biogenic elements due to its key role in biochemical processes in living organisms (Maciá et al. 2005). Phosphorus is believed to be synthesised in massive stars and subsequently ejected into the interstellar medium (ISM) by supernovae explosions (Koo et al. 2013). Although its cosmic abundance is relatively high (fractional abundance of ∼3×10 −7 with respect to H; Asplund et al. 2009), it is difficult to detect in the ISM since it is likely to be heavily depleted onto dust grains (see Turner et al. 1990).
Indeed, phosphorus-bearing species (hereafter Pbearing) have mainly been detected toward diffuse clouds (in the form of P + ; Jura et al. 1978) and circumstellar envelopes around evolved stars (in the form of HCP, CP, CCP, NCCP, PO, PN and PH 3 ; Agúndez et al. 2007Agúndez et al. , 2008Agúndez et al. , 2014aTenenbaum et al. 2007;Halfen et al. 2008;Milam et al. 2008;De Beck et al. 2013). More recently, higher-sensitivity observations have revealed the presence of P-bearing molecules such as PN and PO in regions with active low-mass and high-mass star formation (Yamaguchi et al. 2011;Caux et al. 2011;Rivilla et al. 2016;Lefloch et al. 2016). PN has also been reported toward massive starless cores , which suggests that this molecule could form in relatively quiescent and cold gas (Mininni et al. 2018). In all cases, PO seems to be slightly more abundant than PN (by factors 2-3; see e.g. Rivilla et al. 2016) 3 , and the PO/PN ratio can potentially be used as a length indicator of the pre-stellar collapse phase (see Aota & Aikawa 2012;Lefloch et al. 2016). The analysis of the P-bearing content of the Giant Molecular Clouds (GMCs) in the Galactic Center shows that PN and PO can also be detected in molecular clouds subject to extreme conditions such as those found in galactic nuclei (Rivilla et al. 2018).
The lack of many detections of P-bearing molecules has indeed hampered the study of the chemistry of these species in the ISM. Initial theoretical works focused on the chemistry of cold and warm molecular clouds (Millar 1991) and massive hot molecular cores (Charnley & Millar 1994). With the detection of P-bearing molecules in the circumstellar envelopes of evolved stars, some more modelling was carried out by MacKay & Charnley (2001) and Agúndez et al. (2007), exploring the non-equilibrium chemistry of phosphorus in the outer envelope of these objects. The recent detections of PN and PO in star-forming regions has triggered a new interest for the chemistry of phosphorus in the ISM. While Yamaguchi et al. (2011) and Lefloch et al. (2016) have proposed that P-bearing species are produced in shocked gas after the sputtering of dust grains, Rivilla et al. (2016) suggest that P-bearing molecules are formed during the cold collapse phase and subsequently frozen onto dust grains (Caux et al. 2011;Fontani et al. 2016;Rivilla et al. 2016). The detection of PN and PO in the Galactic Center poses even more unknowns to the chemistry of phosphorus (Rivilla et al. 2018). These authors have indeed found a positive trend between the abundance of PN and the optically-thin 29 SiO isotopologue, which suggests that PN is generated in shocks. However, it is unclear whether P-bearing molecules can survive additional energetic processing due to the intense become available for PO with He and/or H 2 (see Lique et al. 2018). UV radiation and/or enhanced cosmic rays, known to be present in this region of the Galaxy.
In this paper, we present a comprehensive study of the chemistry of phosphorus in the ISM covering a wide range of physical conditions (from molecular dark clouds to cold/warm/hot cores) and of physical processes (protostellar heating, shock waves, UV radiation and cosmic rays). Our goal is to establish the most likely conditions under which P-bearing molecules form, and to what extent. Special focus is put on the analysis of the PO/PN abundance ratio to determine why this ratio is found to be constant (within a factor of 2-3) across different sources. In Section 2, we describe the details of our chemical modelling, while in Section 3 we present the results. In Section 4, we discuss the implications of our modelling for Galactic and extragalactic studies, and we summarised our conclusions in Section 5.

CHEMICAL MODELLING OF P-BEARING MOLECULES
To model the chemistry of P-bearing species, we have used the gas-grain chemical code uclchem (Holdship et al. 2017) 4 , which has been recently updated with a new treatment for the grain surface chemistry (Quénard et al. 2018). The code includes thermal and non-thermal desorption processes as described in Viti et al. (2004) and Roberts et al. (2007) respectively, and considers the grain surface processes of radical diffusion, chemical reactive desorption and reaction-diffusion competition as described in Quénard et al. (2018).
The chemical network is based on the one developed by Quénard et al. (2018), which has been expanded to include P-bearing species. The gas-phase reactions were taken from the UMIST database (McElroy et al. 2013), which includes the original network for phosphorus of Millar (1991). This network is based on the early experimental work of Smith et al. (1989) and Thorne et al. (1983Thorne et al. ( , 1984, and only ∼30% of the reactions have rate coefficients measured in the laboratory. We have also added the chemical reactions involving PH 3 from the chemical network of Charnley & Millar (1994), who took the rate coefficients of the neutralneutral reactions H + PH n (n=1-3) from Kaye & Strobel (1983) following the experimental data of Lee et al. (1976). Other ion-neutral reactions involving PH n products and species such as PS or HPS were also included in our network following Anicich (1993) and according to the experiments of Smith et al. (1989). The two gas-phase reactions N+CP→PN+C and P+CN→PN+C proposed by Agúndez et al. (2007), with reaction rates of 3×10 −10 cm 3 s −1 , were also included in our network. Note that these reaction rates are highly uncertain since no experimental values are available. However, several models were run without these reactions and they yielded similar results.
For the grain surface chemistry, we have considered the hydrogenation reactions of P-bearing species into their non-saturated and saturated forms, as well as their associated diffusion and chemical reaction desorption reactions. Binding energies were taken from the KIDA database (Wakelam et al. 2017) 5 . In total, the network 1.58×10 −9 3 P 2.57×10 −9 4 F 6.7×10 −9 5 a Extracted from Asplund et al. (2009). b As in the low-metal abundance case from Graedel et al. (1982) and Morton et al. (1974).
contains 3695 reactions involving 401 species, 265 out of which are gas phase species while the remaining 136 are grain surface species. The initial elemental abundances considered in our models were taken from Asplund et al. (2009, see the Solar values in their Table 1) except for Mg, Si, Cl, and F, which have been depleted by factors between 5.4 and 5700 (see Table 1 and references therein). For P, its elemental abundance has been depleted by a factor of 100 with respect to its Solar value (P/H=2.57×10 −7 ), to simulate the fact that phosphorus is not detected in quiescent molecular clouds (Turner et al. 1990). Note that in the models of Aota & Aikawa (2012), a higher initial abundance of P/PH 3 of ∼10 −8 is assumed in the ices, which corresponds to a depletion factor of ∼10. However, a better agreement with the observations toward the B1 shocked region in the L1157 molecular outflow is obtained when this initial abundance is reduced to ∼3×10 −9 (see their Section 3.2). This is also supported by recent findings by Lefloch et al. (2016) and Rivilla et al. (2016) in, respectively, L1157 B1 and a sample of massive star-forming regions. Therefore, it is reasonable to assume that P is likely depleted by a factor of 100 in molecular dark clouds.
uclchem is run in three phases. Phase 0 simulates the chemistry in a diffuse cloud with density n(H)=100 cm −3 and temperature 20 K for 10 6 yrs. We have also tested higher temperatures for the Phase 0 stage (T= 100 K; see Turner et al. 1999) and the abundances of the P-bearing species at the end of Phase 1 (the collapse phase, see below) differ by less than 10%. In Phase 1, the cloud undergoes free-fall collapse keeping the temperature constant at 10 K until the final density (n(H)=2×10 4 , 2×10 5 and 2×10 6 cm −3 ) is reached. As explained in Section 3.1, we have investigated the effects of a short-lived/long-lived collapse phase since it has been found to have important implications for the chemistry of PN and PO (see e.g. Aota & Aikawa 2012;Lefloch et al. 2016). Finally, Phase 2 simulates the physical processes associated with star formation such as warm-up of the protostellar envelope (up to 50, 100 and 300 K, as described in Viti et al. 2004), interaction of C-type shocks (for shock velocities of v s =20 and 40 km s −1 and using the parametrization of the physical structure of C-shock waves formulated by Jiménez-Serra et al. 2008), intense UV radiation (with χ=1, 10 2 and 10 4 Habing, as explored in Viti & Williams 1999) and cosmic rays (with ζ=1.3×10 −17 , 1.3×10 −15 and 1.3×10 −13 s −1 ; see Harada et al. 2015).

RESULTS
3.1. Effects of collapse on the chemistry of P-bearing species In Figure 1 we show the results obtained for the collapse phase (Phase 1 in uclchem ) for final gas densities n(H)=2×10 4 , 2×10 5 , and 2×10 6 cm −3 considering a short-lived collapse (the code stops once the final density is reached, i.e. at time-scales of 5.2-5.4×10 6 yrs de-pending on the final density) and a long-lived collapse (chemistry runs until a final time-scale of 6×10 6 yrs, i.e. between 6-8×10 5 years longer after the final density is attained; see also Table 2). In earlier works (see Aota & Aikawa 2012;Lefloch et al. 2016), the length of the collapse was found to have a strong impact on the subsequent evolution of the chemistry of PO and PN, and on their predicted abundance ratio. Indeed, in the short-lived scenario the gas-phase abundance of atomic nitrogen by the end of the collapse is large, which enhances the formation of PN and yields PO/PN ratios ≤1. On the contrary, if the collapse is long-lived, atomic nitrogen has enough time to deplete on grain surfaces and get converted mainly into NH 3 , giving as a result PO/PN ratios ≥1 (Lefloch et al. 2016). We therefore explore the consequences of the length of the collapse in our models by making the core become static for an extra 6-8×10 5 years after the end of the collapse. This extra time is consistent with the estimated lower-limit dynamical age of pre-stellar cores (≥1-2×10 5 years; see Section 4 in Pagani et al. 2012) and provides a good match to observations toward this type of cores (Quénard et al. 2018).
In our models, we assume a physical radius for the core of 0.3 pc, which implies visual extinctions of 13 mag, 120 mag and 1200 mag for hydrogen densities of n(H)=2×10 4 , 2×10 5 and 2×10 6 cm −3 respectively, according to the following expression (Bohlin et al. 1978;Frerking et al. 1982): where r out is the radius of the core, n(H) is the number density of hydrogen nuclei, A v0 is the background extinction assumed to be 2 mag, and A v is the total extinction within the core. As shown in Figure 1, for all models the most abundant P-bearing species during the collapse are atomic phosphorus (P) and PN, with maximum abundances of ∼5×10 −10 -10 −9 at time-scales of ∼5×10 6 yrs (i.e. when the final density is reached). The formation of PN in detectable abundances starts with the gas-phase reactions N + CP → PN + C and P + CN → PN + C, and proceeds through the reaction N + PO → PN + O in the gas phase by the end of the collapse. The formation of PO and PS is delayed with respect to PN and it occurs in the gas-phase via reactions O + PH → PO + H and HPS + + e-→ PS + H respectively. PH 3 dramatically increases its gas-phase abundance by the end of the collapse, because it is efficiently formed on grain surfaces via hydrogenation and subsequently non-thermally desorbed. Note also the dramatic drop of P + with density as a result of the enhanced freeze out for ions.
At the highest densities in the collapse, the gas-phase abundances of all P-bearing species start to decrease as a result of freezing out onto dust grains. This depletion is more extreme in the long-lived case since species can freeze out for longer reaching gas-phase abundances even lower than 10 −13 . Except for models with n(H)=2×10 6 cm −3 where all P-bearing species are fully frozen out in both the long-lived and short-lived cases, these differences in depletion by the end of the collapse have an impact on the PO/PN abundance ratio predicted during , and 2×10 6 cm −3 (right). The length of the collapse is varied so that the chemistry stops when the final density is reached (short-lived collapse; upper panels) or the chemistry is let to evolve for a few 10 5 yrs more after the final density is reached (long-lived collapse; lower panels).

Effects of proto-stellar heating
In this Section, we analyze the effects of the increase of dust and gas temperature due to the protostar on the chemistry of phosphorus. We consider a range of hydrogen gas densities (n(H)=2×10 4 , 2×10 5 , and 2×10 6 cm −3 ) and final temperatures (T = 50, 100 and 300 K) that are typically found in warm and hot cores (see Table 2). We simulate the effect of the presence of a protostar at the center of the core by subjecting it to an increase in the gas and dust temperature. The temperature reaches its maximum value at the contraction time of the core (occurring at time-scales ≤10 5 yrs), following a power law derived by Viti et al. (2004) from the observational luminosity function of Molinari et al. (2000). For simplicity, in our models we only consider the case of proto-stellar heating produced by a 15 M ⊙ proto-star.
The results of the Phase 2 of our models are shown in Figure 2. For T = 50 K, the only P-bearing that is released into the gas phase from the mantles is atomic P (at ∼3×10 4 yrs). PH 3 is not thermally desorbed at such low temperatures and its solid abundance increases steadily during the warming-up phase thanks to hydrogenation (in particular, via #PH 2 + #H → #PH 3 , where # denotes a grain mantle species). However, once the temperature reaches its maximum and is kept constant to T = 50 K, hydrogenation ceases because atomic P is thermally desorbed at this temperature (see above) and PH 3 in the gas phase is slowly destroyed via proton transfer reactions with H + 3 , HCO + and H 3 O + . PO is formed in the gas phase via the reaction PH + O → PO + H while the formation of PN depends on the abundance of PO through N + PO → PN + O. In fact, these two species reach an equilibrium at time-scales ∼3-4×10 5 yrs for all densities showing similar abundances of a few 10 −11 . PS remains very low in all models with abundances ≤10 −12 .
For T = 100 and 300 K, the evolution of the P-bearing species are similar presenting a first jump in their abundances at ∼5×10 4 yrs once the temperature of dust/gas reaches 100 K and all ices are thermally desorbed into the gas phase. Since PH 3 is the main reservoir of P in the mantles, this species is the most abundant for the first few 10 5 yrs. After that, it gets destroyed via proton transfer reactions with mainly H 3 O + . The decrease in the abundance of PH 3 is faster in models with T=300 K because the destruction reaction PH 3 + H → PH 2 + H 2 becomes efficient at higher temperatures due to its activation barrier of 735 K. Also, note that PH 3 is destroyed in the gas phase faster than other species such as PN and PO because, while PN/PO can re-form in the gas phase via reactions between P, PH and PH 2 and O, O 2 and N, the only gas-phase formation routes for PH 3 in our network are PH + 4 + NH 3 → PH 3 + NH + 4 and PH + 4 + e − → PH 3 + H (see also Thorne et al. 1984). For PO, this molecule tends to increase its abundance with time in most models and, except for densities n(H)=2×10 6 cm −3 , PO always stays above PN. Once PS is desorbed from the mantles of dust grains, its abundance remains practically constant at ∼5-6×10 −12 . In all these models the abundance of P + remains below 10 −13 .
The effects of the length of the collapse on the subsequent chemistry of P-bearing species during Phase 2 is illustrated in Figure 3. For warm temperatures (T=50 K), the abundances in Phase 2 derived after considering a short-collapse phase, are enhanced by at least a factor of 10 with respect to the long-collapse phase (see left panel in Figure 3). However, for higher temperatures, -Abundances of P-bearing species as a function of time simulated during the warming-up of the molecular envelope by a central protostar up to temperatures of 50 K (lower panels), 100 K (middle panels) and 300 K (upper panels). The final densities of the envelope are n(H)=2×10 4 (left), 2×10 5 (middle), and 2×10 6 cm −3 (right). These models have been run considering a long-lived collapse phase (Section 3.1). Fig. 3.-Examples of the evolution of the abundances of P-bearing species as a function of time simulated during the warming-up of the molecular envelope by a central protostar considering a short-lived collapse phase. The abundances for models with n(H)=2×10 4 cm −3 and T=50 K are at least a factor of 10 higher than those for the long-lived collapse case. For higher temperatures, however, the behaviour of the P-bearing species is similar for both the short-lived and long-lived collapse scenarios.
although the main P reservoir is PN instead of PH 3 at the beginning of Phase 2, the overall behaviour of all P-bearing species is similar to that of the long-collapse scenario since the gas phase chemistry clearly dominates.

Effects of UV-photon radiation
The effects of intense UV illumination on the chemistry of phosphorus are evaluated by varying the interstellar radiation field within ranges typical of photondominated regions (PDRs). We consider UV-photon radiation fields of χ= 100 and 10 4 Habing, visual extinctions of A v = 3, 7.5 and 13 mag, and temperatures T = Fig. 4.-Evolution of P-bearing species under the effects of a radiation field of χ= 1, 100 and 10 4 Habing and for different visual extinctions (Av=3 and 7.5 mag) after a long-lived collapse. Solid lines denote models with T= 50 K, while the results from models with T=100 K are shown in dashed lines. The results from models with Av= 13 mag are identical to those obtained with Av= 7.5 mag, and they are thus not shown. 50 and 100 K. These values are similar to those found in PDRs such as the Horsehead Nebula and the Orion Bar (Goicoechea et al. 2009;Schilke et al. 2001). In the models, the temperature is kept constant and the ices are instantaneously evaporated at the beginning of Phase 2. For completeness, we also present the results for models with χ= 1 Habing (see Table 2). The results are reported in Figures 4 and 5.
From Figure 4, we find that at low extinctions (A v = 3 mag) PN and P are the most abundant P-bearing species (with abundances ≥5×10 −10 -10 −9 ) for low and intermediate UV radiation fields (χ≤100 Habing), while PO is efficiently destroyed. This is partly explained by the larger photo-destruction rate of PO (with α=3×10 −10 s −1 and γ=2.0, where the rate is calculated as k=αe −γAv ) than for PN (with α=5×10 −12 s −1 and γ=3.0). Note that these photo-destruction rates have been estimated from the ones for NO and N 2 (see Harada et al. 2010) and, therefore, may be inaccurate. However, even in the case that PN were destroyed at higher rates via UV photodissociation, it would re-form again in the gas phase thanks to the large amounts of atomic N and P available in the gas phase, which react via N + PO → PN + O and P + CN → PN + C yielding PN. Like PO, PH 3 and PS are also increasingly photo-dissociated with increasing UV radiation field. Note that for even higher values of the UV field (χ= 10 4 Habing), all molecular material is photo-dissociated including PN, which shows abundances ≤10 −12 .
For higher extinctions (A v =7.5 mag), the behaviour of all P-bearing species is similar to those obtained for the warming-up case with T= 50 and 100 K (see Figure  2) since UV photo-dissociation becomes practically ineffective at A v ≥5 mag. As expected, the results of models with A v =13 (not shown in Figure 4) are identical to those with A v =7.5 mag. The differences between models with T=50 K and T=100 K are minimal except for PS, whose chemistry is more sensitive to higher temperatures (see dashed lines in Figure 4). In models with a short-lived collapse phase, PN is the most abundant P-bearing molecule with an abundance of ∼10 −9 , following P closely. As for the long-lived case, only small differences are found between models with T= 50 K and T= 100 K. The results with A v = 7.5 mag are also very similar to the short-lived warm-up models from Figure 3.
We finally note that in all UV-photon illuminated models at low extinctions (A v =3 mag) or under high UV radiation fields (χ= 10 4 Habing), the abundance of PN always remains above that of PO.

Effects of cosmic-rays
We now evaluate the effects of cosmic-rays on the chemistry of P-bearing species. Cosmic-rays are present in a variety of environments from shocks in molecular outflows (Podio et al. 2014), to star-forming regions Fontani et al. 2017), to the GMCs in the CMZ (Goto et al. 2014). In Figures 6  and 7, we present the evolution of the P-bearing species within a molecular cloud with hydrogen volume densities  Figure 4, but assuming a short-lived collapse phase. Solid lines indicate models with T= 50 K, dashed lines denote models with T=100 K. The abundance variations due to the change in temperature, are minimal. The results from models with Av= 13 mag are identical to those obtained with Av= 7.5 mag, and they are thus not shown. n(H)=2×10 4 and 2×10 5 cm −3 and gas temperatures T = 10 and 100 K, under the influence of cosmic-rays whose ionization rates have been increased by factors 100, 1000 and 10 4 (i.e. ζ=1.3×10 −15 , 1.3×10 −14 and 1.3×10 −13 s −1 , respectively). For the T=100 K case we consider that the ices are instantaneously evaporated at the beginning of Phase 2, as in Section 3.3. The physical conditions assumed for these models are presented in Table  2.
For the models with n(H)=2×10 4 cm −3 and T=10 K, Figure 6 shows that a higher ζ (by factors 100 and 1000) enhances the abundances of species such as PN and PO with respect to those predicted by the collapse models at the end of the collapse (see lowermost left panel of Figure 1). This is due to i) the non-thermal desorption of these species by cosmic-ray induced secondary UV photons, and to ii) the rapid photo-dissociation of PH 3 by this secondary UV field. As already noted in Section 3.2, once PH 3 is destroyed, it cannot re-form efficiently in the gas phase (Thorne et al. 1984), which unlocks high abundances of gas-phase P that can then be used to produce PN and PO.
For ζ=1.3×10 −13 s −1 (i.e. an enhancement in ζ by a factor of 10 4 ), molecules are rapidly photo-dissociated by the strong secondary UV-photon radiation field, yielding large abundances of P + and P. A similar trend is seen for higher temperatures (T=100 K), where PH 3 is rapidly destroyed in the gas phase in favour of PN and PO. In all models, PS is either destroyed or its abundance remains below ≤10 −12 . PN lies above PO with PO/PN ratios ∼0.5-1 except for the case with T=100 K and ζ=1.3×10 −13 s −1 where PO is a factor of ∼10 more abundant than PN.
For higher densities (n(H)=2×10 5 cm −3 ; see Figure 7), PN and PO are largely enhanced for ζ=1.3×10 −13 s −1 with respect to the low-density case thanks to the higher abundance of atomic N and O in these models as they are more protected by the higher extinction. Atomic N and O yield higher abundances of PN and PO via reactions N + PO → PN + O and O + PH → PO + H, respectively. PN also tends to be more abundant than PO in the highdensity models.
We note that all models above consider a long-lived collapse phase (see Section 3.1). If the collapse phase is short-lived, the largest discrepancies are generally found for models with ζ≥1.3×10 −14 s −1 , where the PN and PO abundances decrease by factors 10-100 with respect to the long-lived collapse case (see Figures A12 and A13). This is explained by the fact that in the short-lived collapse, the ice reservoir of P-bearing molecules is not as large as in the long-lived case, and the gas-phase Pbearing species are rapidly destroyed by the cosmic-ray induced secondary UV photons.
3.5. Effects of C-type shock waves General behaviour of P-bearing species in C-shocks P-bearing species such as PN and PO have recently been observed in shocked regions in molecular outflows (Yamaguchi et al. 2011;Lefloch et al. 2016). In addition, these species have also been detected in the turbu-  a This refers to the length of the collapse during Phase 1 in our modelling (see also the caption in Table 2). lent GMCs in the Central Molecular Zone (Rivilla et al. 2018), believed to be affected by widespread low-velocity shocks (Requena-Torres et al. 2006) and likely under the influence of enhanced cosmic-ray ionization rates (Goto et al. 2014;Harada et al. 2015).
In Table 3 we present the physical parameters assumed for the C-type shock models used in our calculations. The pre-shock density is denoted by n(H), v s is the shock speed, T n,max refers to the maximum temperature of the neutral fluid attained within the shock and t sat is the saturation time at which most of the material within the icy mantles of dust grains is sputtered into the gas phase (for more details, see Jiménez-Serra et al. 2008;Holdship et al. 2017). Different values of the cosmic-ray ionization rate (ζ=10 −17 , 10 −15 , and 10 −13 s −1 ) are also considered. Figures 8 and 9 show the evolution of P-bearing molecules across the C-type shock structure for hydrogen gas densities of n(H)=2×10 4 cm −3 and n(H)=2×10 5 cm −3 , respectively. All species are enhanced early in the shock (at ∼60 yrs for n(H)=2×10 4 cm −3 and at ∼6 yrs for n(H)=2×10 5 cm −3 ) due to sputtering. In particular, PH 3 and PN present enhancements larger than fac- tors 100 and 10, respectively. For models with shock speeds v s =20 km s −1 , PH 3 stays relatively abundant throughout the shock after the sputtering of the ices (abundance between 10 −10 -10 −9 ), except in the presence of a high cosmic-ray ionization rate when PH 3 is more efficiently destroyed thanks to reactions with He + and H + 3 (see model with ζ=1.3×10 −13 s −1 in Figure 8). For higher shock velocities (v s = 40 km s −1 ), PH 3 is destroyed due to the endothermic reaction PH 3 + H → PH 2 + H 2 , whose activation energy is 735 K (see Table  1 in Charnley & Millar 1994). The destruction of PH 3 is even faster at higher temperatures (as in models with n(H)=2×10 5 cm −3 and v s =40 km s −1 ) or under the effects of cosmic rays.
From Figures 8 and 9, we find that PN remains relatively constant throughout the shock with abundances typically falling between 1-5×10 −10 . In contrast to PN, PO shows a different behaviour depending on the model, becoming more abundant than PN in shocks with n(H)=2×10 4 cm −3 and ζ≥10 −15 s −1 , or in shocks with n(H)=2×10 5 cm −3 and ζ≥10 −13 s −1 . PS experiences little changes in the gas phase after its injection into the gas phase by sputtering, except in models with ζ=1.3×10 −13 s −1 where it is destroyed by H + 3 and He + . For the short-lived collapse models, Figures A14 and  A15 show that the most abundant P-bearing species within the shock are atomic P and PN, since they have not had enough time to deplete onto dust grains and to subsequently hydrogenate during the collapse phase. PO becomes comparable in abundance to PN only for high cosmic-ray ionization rates (ζ∼10 −13 s −1 ). Lefloch et al. (2016) also presented models of the chemistry of P-bearing species in C-type shock waves using uclchem and the C-shock parametric approximation of Jiménez-Serra et al. (2008). The chemical network used in these models, however, did not include either PH 3 or the endothermic reactions PH 3 + H → PH 2 + H 2 , PH 2 + H → PH + H 2 , and PH + H → P + H 2 (with energy barriers ranging from 318 to 735 K; Charnley & Millar 1994). This has important consequences in the evolution of the molecular abundances of PO and PN and its ratio in the shock since, as opposed to the results of Lefloch et al. (2016), PN is always more abundant than PO in our models with standard values of the cosmic-ray ionization rate (i.e. ζ=10 −17 s −1 ; see Figures 8 and 9). If the exothermic reactions of PH 3 , PH 2 and PH with H are switched-off in our models, we recover the results of Lefloch et al. (2016).

Comparison with the Lefloch et al. (2016)'s models
Millimeter observations of PN and PO in outflows (Lefloch et al. 2016) and in Galactic Center Giant Molecular Clouds (Rivilla et al. 2018), however, indicate that the PO/PN abundance ratio is always ≥1, which can be reconciled with our models only if i) cosmic rays are significantly enhanced; or ii) the chemical network of PO is incomplete. In the following section (Section 4), we explore the possibility that the chemical network of PO lacks a few key formation reactions. Codella et al. (2018) have recently carried out the modelling of the chemistry of NO (a species with a molecular structure similar to that to PO) in hot corinos and shocked gas in molecular outflows. These authors have noted that the enhancement of NO in the L1157-B1 shocked region, is mainly produced by the reaction of atomic nitrogen, N, with OH. The equivalent reaction for P (P + OH → PO + H), however, does not exist in Fig. 8.-Evolution of the abundances of P-bearing species as a function of time across two C-type shocks with a pre-shock density of n(H)=2×10 4 cm −3 and shock speeds of vs=20 km s −1 (left panels) and vs=40 km s −1 (right panels). Dashed lines indicate the evolution of the temperature of the neutral fluid within the shock (Tn). To simulate the extreme conditions in the Galactic Center, we also consider that the shocked gas is affected by enhanced cosmic ray ionization rates of ζ= 10 −15 and 10 −13 s −1 . These models have been obtained considering a long-lived collapse. Table 4. Reaction rate assumed for the models of Section 4, as extracted from UMIST.

Reaction)
α (cm 3 s −1 ) β γ (K) P + OH → PO + H 6.1×10 −11 -0.23 14.9 any database and, given the similarity between NO and PO, the reaction between P and OH may also be efficient (note that the activation barrier is only of ∼15 K; see Table 4).
We have therefore re-run all our models after including the reaction P + OH → PO + H in our chemical network. As reaction rate, we have assumed the one from the analogous reaction with NO (see Table 4). Our results show that for most models the inclusion of this reaction does not affect appreciably the final abundance of PO, since it is enhanced only by factors ∼2-3. However, for models with C-type shock waves, the abundance of PO can be enhanced by several orders of magnitude (see Figure  10). From this Figure, it is clear that thanks to the new reaction P + OH → PO + H, PO even becomes more abundant than PN within the shock for standard cosmic-ray ionization rates, in contrast to our results from Section 3.5. This occurs thanks to the high temperatures (of thousands of K) attained within the shock. In fact, if this reaction were as efficient as proposed here, the PO/PN ratio in shocks would be dominated by this reaction and not by enhanced cosmic rays any longer (see lower panels of Figure 10). Theoretical calculations and/or experiments are needed to clearly establish the efficiency of the PO gas-phase formation route P + OH → PO + H.

PN and PO
In this section we compare our modelling results to the abundances of PN and PO measured toward low-mass/high-mass starless/pre-stellar cores (Turner et al. 1990;Mininni et al. 2018), massive hot cores (Turner et al. 1990;Rivilla et al. 2016), GMCs in the Galactic Center (Rivilla et al. 2018), and molecular outflows (Yamaguchi et al. 2011;Lefloch et al. 2016): • Cold Cloud and Starless/Pre-stellar Cores: Early millimeter observations of PN toward a sample of low-mass cold cloud cores (with T∼10 K) provided upper limits to the abundance of PN of ≤0.2-4×10 −11 in cold cloud cores (Turner et al. 1990).  -Comparison of models obtained for n(H)=2×10 5 cm −3 and vs=40 km s −1 using the original chemical network (as presented in Section 2; left panels) and a new network where the proposed PO formation route P + OH → PO + H is included. Models with comic-ray ionization rates of ζ∼10 −17 s −1 (upper panels) and ζ∼10 −15 s −1 (lower panels) are presented. All models shown here consider a long-lived collapse phase.
Recent higher-sensitivity observations performed toward the L1544 pre-stellar core reveal that PN remains undetected with upper limits of ≤4.6×10 −13 (as derived from the dataset of Jiménez-Serra et al. 2016). As shown in Figure 1 for the long-lived col-lapse case (see bottom panels), these upper limits are consistent with the low gas-phase PN abundances predicted at a few 10 5 yrs after the end of the collapse in cold high-density cores (with hydrogen gas densities ≥2×10 5 cm −3 ). For massive starless cores, which show slightly warmer temperatures (with T∼25-30 K; Fontani et al. 2016), the derived abundances of PN are ∼10 −11 and 5×10 −12 (Mininni et al. 2018), in agreement with those predicted by our models under warm temperatures for n(H)≥2×10 5 cm −3 and time-scales of a few 10 5 yrs (see models with T=50 K in Figure 2). No upper limits are available for PO toward low-mass/highmass cold cores.
• Galactic Center GMCs: Rivilla et al. (2018) have recently carried out observations of PN and PO toward a sample of quiescent GMCs in the Galactic Center that present either shock-dominated or radiation-dominated chemistries. Rivilla et al. (2018) found a positive trend between the abundances of PN and 29 SiO, which suggests that PN may be formed in shocks as SiO. The derived PN abundances toward shock-dominated, quiescent GMCs are ∼0.5-5×10 −11 (Rivilla et al. 2018), i.e. within a factor of 10 those predicted by our models with v s =20 km s −1 , n(H)=2×10 4 cm −3 and high cosmic ray ionization rates (see Figure  8, left panels). The measured PO/PN ratio of 1.5 in the quiescent GMC G+0.693-0.03 also supports the models with enhanced cosmic rays, although, as discussed in Section 4, a PO/PN ratio ≥1 could also be obtained in shocks if the reaction P + OH → PO + H were efficient. This reaction is currently not included in any chemical database. For the quiescent GMCs affected by radiation, only upper limits to the abundance of PN are available. This is consistent with our results of the chemistry of P-bearing species in PDRs, for which UV radiation fields χ∼10 4 Habing destroy PN and PO down to abundances ≤10 −12 (for A v =3 mag; see Figures 4 and 5). The lack of detection of PN in the Orion Bar (Cuadrado et al. 2015;Rivilla et al. 2018), confirms this hypothesis (the Orion Bar is affected by a UV radiation field of χ=5×10 4 Habing; see e.g. Schilke et al. 2001).
• Molecular outflows: In Lefloch et al. (2016), the preferred model to explain the PO/PN abundance ratio of ∼2 measured in the L1157-B1 shock cavity, involves a shock with gas densities of 10 5 cm −3 , shock speeds of 40 km s −1 , and a standard cosmicray ionization rate (1.5×10 −17 s −1 ). Unless the chemical network of PO is incomplete (see Section 4), our shock model results suggest that the chemistry of P-bearing species in L1157-B1 is likely affected by cosmic rays. This is consistent with the detection of molecular ions such as HOCO + and SO + in L1157-B1, which requires a cosmicray ionization rate of ∼3×10 −16 s −1 (Podio et al. 2014). Theoretical studies and/or new experiments are needed to establish whether the reaction P + OH → PO + H is indeed as efficient as proposed in Table 4.
5.2. Phosphine (PH 3 ) Although phosphine (PH 3 ) has been detected in the atmospheres of Jupiter and Saturn (Larson et al. 1977;Weisstein & Serabyn 1994), there is no evidence of its presence in star-forming regions. According to our models in Section 3, PH 3 tends to be destroyed rapidly in the gas phase (for time-scales shorter than a few 10 4 years) after its ejection from dust grains, as also noted by Charnley & Millar (1994). Nevertheless, for some cases PH 3 could remain in the gas phase for longer time-scales in hot molecular cores and in shocked regions, reaching gas-phase abundances as high as ∼10 −9 (see Figures 2,8 and 9). PH 3 could then be observed via its J=1→0 and J=2→1 rotational transitions at 266.9 and 533.8 GHz respectively, with instruments such as the IRAM 30m telescope, ALMA and SOFIA. To our knowledge, no observations (or reported upper limits) do exist for this molecule toward hot molecular cores. However, upper limits to the abundance of PH 3 have been reported toward the B1 shocked region in the L1157 molecular outflow (of ≤10 −9 ; Lefloch et al. 2016). These upper limits, which were inferred from high-sensitivity IRAM 30m telescope data 6 , are consistent with the abundances predicted by our models considering the interaction of shock waves (Figures 8 and 9).

THE PO/PN RATIO
For the sources where PO has been detected in the ISM (e.g. hot cores, molecular outflows and the quiescent Galactic Center GMC G+0.693 in the Galactic Center; see Rivilla et al. 2016;Lefloch et al. 2016;Rivilla et al. 2018), observations have revealed that the PO/PN ratio is always ≥1. Our models of Section 3 have however shown that the PO/PN ratio may become ≤1 for certain physical conditions and energetic phenomena. Figure 11 presents the PO/PN ratios obtained by our models for a time-scale of 10 5 yrs (typical of hot molecular cores and Galactic Center GMCs; Requena-Torres et al. 2006;Rivilla et al. 2016) for the models of proto-stellar heating, cosmic rays and UV radiation, and of 10 3 yrs (typical of outflows; see e.g. Podio et al. 2014) for the models with shocks and shocks with cosmic rays. All models consider a long-lived collapse phase. As seen from Figure 11, the different energetic processes tend to cluster the PO/PN ratios toward different parts of the diagram. While proto-stellar heating and cosmic rays typically give PO/PN≥1, this ratio Fig. 11.-Compilation of the PO/PN abundance ratios obtained within our models as a function of temperature. The ratios corresponding to two different time-scales are shown for every model: 10 5 yrs for models considering protostellar heating (red squares), cosmic rays (blue pentagons) and UV radiation (green and black triangles); and 10 3 yrs for models with shocks and shocks+cosmic rays (magenta and black empty stars). Horizontal dashed magenta line denotes the PO/PN ratio of 1.
is expected to be <<1 in PDRs (with A v ∼3 mag). This explains why the observed PO/PN ratios in hot molecular cores are ≥1 . For PDR regions with higher extinction (with A v ∼7 mag), the chemistry of P-bearing species behaves in a similar way to that affected by photo-stellar heating, since the efficiency of molecular photo-dissociation drastically decreases at extinctions A v ≥5 mag.
For models with shocks, the PO/PN ratio is ≤1 in the absence of cosmic rays while this ratio is boosted above 1 under their effects. However, note that this trend may be a result of the incompleteness of the PO chemical network (see Section 4). If the PO formation reaction proposed in Table 4 is experimentally or theoretically confirmed, the PO/PN ratios in shocks are also expected to be ≥1.
Finally, for models with a short-lived collapse phase, the derived PO/PN ratios typically fall below 1, except for a few exceptions which involve enhancements of the cosmic-ray ionization rate by factors of 10 2 and 10 4 .

CONCLUSIONS
P-bearing species such as PN and PO have started to be routinely detected in regions of the ISM affected by different energetic phenomena such as proto-stellar heating, cosmic rays, UV radiation and shock waves. In this paper, we have re-visited the chemistry of phosphorus in the ISM under energetic processing, putting special emphasis on the predicted ratio between PN and PO.
Our models show that a long-lived collapse is required to reproduce the PO/PN ratios ≥1 observed in hot molecular cores (see Section 6). The models considering UV illumination confirm that P-bearing species are expected to be destroyed under the effects of strong UV radiation fields, as observed in the Orion Bar and Galactic Center GMCs (Section 3.3). Moderate UV radiation fields (χ∼100 Habing), however, should yield detectable abundances of P-bearing species such as PN, PO and PS. Models with enhanced Cosmic ray ionization rates provide a wide range of PO/PN ratios, although PN tends to be more abundant than PO. For models with C-type shock waves, an observed ratio PO/PN≥1 is only obtained under the effects of enhanced rates of cosmic rays, unless the chemical network of PO is incomplete (Section 3.5). We propose that the reaction P + OH → PO + H (currently missing in all chemical databases) could be as efficient as its analogue reaction involving NO. Theoretical/experimental investigations are needed to establish the actual efficiency of this gas-phase formation route of PO.