© J. Chantzos et al. 2020

1 Introduction

Phosphorus is essential for biotic chemistry, since it is a fundamental component of many important biological molecules, such as nucleic acids and phospholipids. Phosphorus is therefore vital to life on Earth and can consequently play an important role in exoplanets (Schaefer & Fegley 2011). Despite its importance, the chemistry of P-bearing molecules is in its infancy and remains poorly understood. The aim of this work is to add an important missing piece to the puzzle: unveiling the first steps of P chemistry via observations and chemical simulations of simple P-bearing molecules in diffuse clouds.

The ion P⁺ was detected in several diffuse clouds by Jura & York (1978), where an elemental abundance of ~2 × 10⁻⁷ with a low P depletion factor of between approximately two and three was derived. However, a more recent study by Lebouteiller et al. (2006) showed that phosphorus remains mostly undepleted towards diffuse clouds. In addition, P has been identified towards dwarf and giant stars (Maas et al. 2017; Caffau et al. 2016), while detections of simple P-bearing molecules (PN, PO, HCP, CP, CCP, NCCP, PH₃) have been done towards the circumstellar material of carbon- and oxygen-rich stars (Agúndez et al. 2007, 2014a,b; Tenenbaum et al. 2007; De Beck et al. 2013; Ziurys et al. 2018). The species PN and only very recently PO are the only P-bearing molecules to have been discovered towards dense star-forming regions (Turner & Bally 1987; Fontani et al. 2016, 2019; Rivilla et al. 2016; Lefloch et al. 2016; Mininni et al. 2018) and molecular clouds in the Galactic Center (Rivilla et al. 2018). The limited number of available observations hinders our understanding of the chemical pathways involved in P chemistry. The main uncertainty in P chemistry is the unknown depletion factor of P in molecular clouds. In general, chemical models of dark clouds start with the so-called “low-metal abundances”, where the elemental abundances of heavy elements (such as P, S, Fe, Mg) are reduced by orders ofmagnitude to reproduce molecular observations (e.g., Agúndez & Wakelam 2013), but with poor understanding of the chemical processes at the base of such depletions. In the case of P, the level of depletion is still very uncertain. While Turner et al. (1990) and Wakelam et al. (2015) used high depletion factors of 600–10⁴ with respect to cosmic P abundance, recent works have shown that it could be as low as ~ 100 (Rivilla et al. 2016; Lefloch et al. 2016).

As only a very limited number of P-bearing molecules have been detected in star-forming regions, it is very hard to put constraints on the elemental abundance of P in the gas phase and on the major chemical pathways. In order to elucidate the interstellar P chemistry we focus on diffuse clouds, which represent the first steps of molecular-cloud evolution. Diffuse clouds can provide important constraints on P chemistry, since P in these objects is not strongly affected by depletion, meaning that the initial P abundance that can be used for chemical simulations is well constrained (Lebouteiller et al. 2006). With this approach we are able to remove an important uncertainty in our model and use a reliable starting point for our chemical simulations.

The existing chemical and physical models focus solely on diffuse clouds (e.g., Dalgarno 1988; Le Petit et al. 2004; Cecchi-Pestellini et al. 2012; Godard et al. 2014). For example, in Godard et al. (2014), a model including dissipation of turbulence was applied to reproduce the observed molecular abundances in the diffuse interstellar medium (ISM). The main results showed that chemical complexity is strongly linked to turbulent dissipation, which was able to reproduce the high abundances of CO and other species (such as C⁺ and HCO⁺) observed towards Galactic diffuse clouds. Le Petit et al. (2004) describe the development of a chemical model of the diffuse cloud towards ζ Persei that was able to reproduce the abundance of ${\text{H}}_{\text{3}}^{\text{+}}$ and other species, like CN and CO. This was achieved by modeling two phases, namely a small dense phase (~ 100 au) with a density of n(H) = 2 × 10⁴ cm⁻³ and a larger diffuse region (4 pc) with n(H) = 100 cm⁻³. In addition, the reproduction of the CH⁺ abundance and that of the rotationally excited H₂ required the inclusion of shocks into the model. Similar results were achieved by Cecchi-Pestellini et al. (2012) when including the injection of hot H₂ into the model.

Previous observations (Corby et al. 2018; Liszt et al. 2018; Thiel et al. 2019) prove the chemical complexity and the wide range of densities, temperatures, and visual extinctions of diffuse and translucent clouds, making them promising targets for observations of P-bearing molecules. Diffuse clouds are characterized by low densities with n(H) = 100−500 cm⁻³ and are therefore more exposed to interstellar radiation, which can destroy molecules. Translucent clouds on the other hand are an intermediate state between diffuse and dense molecular clouds, being more protected from UV radiation (1 mag < A_V < 5 mag). They are denser with typical densities of n(H) = 500−5000 cm⁻³ and are consequently cooler (T_gas = 15−50 K), showing higher chemical complexity (Snow & McCall 2006; Thiel et al. 2019). One prominent candidate that has been widely studied in previous works (e.g., Liszt et al. 2018, and references therein) is the gas that lies along the line of sight to the compact extragalactic continuum source B0355+508. This strong blazar is located at a very low latitude in the outer Galaxy (b = −1.6037°), meaning that the way through the Galactic disk is long and therefore gathers a significant amount of distributed Galactic diffuse gas (Pety et al. 2008). Indeed, the line of sight towards B0355+508 shows a complex kinematic structure which incorporates several diffuse and translucent clouds. The detections of numerous molecules like S- and CN-bearing species as well as small hydrocarbons towards B0355+508 also indicate the rich chemistry present in this diffuse and translucent gas (e.g., Liszt et al. 2018, and references therein). The substantial velocity structure coupled with a high chemical complexity of this line of sight enables us to adjust our chemical and physical model to every cloud component and find which physical conditions most favor the abundances of P-bearing molecules. Other background sources that have previously been studied are either lacking the chemical (like B0224+671) or the velocity (such as B0415+479) features which are essential for the present work.

In this paper, we present single-pointing observations of the (2–1) transitions of HCP, CP, PN, and PO and chemical simulations of their molecular abundances towards the line of sight to B0355+508 in order to investigate P-bearing chemistry within diffuse and translucent clouds, the precursors of molecular clouds. In Sect. 2 we describe the observational details. Section 3 summarizes the results of the observations. In Sect. 4 we describe our updated phosphorus chemical network as well as the grid of models that we apply in order to reproduce the observations of HNC, CN, CS, and CO towards every cloud component along the line of sight. Furthermore, in Sect. 5 we focus on the P-bearing chemistry based on our best-fit model (which was determined in Sect. 4). In particular, we report the predicted molecular abundances of HCP, CP, PN, PO, and PH₃ and we study their dependence on visual extinction, cosmic-ray ionization rate, and diffusion-to-desorption energy ratio on dust grains. A future outlook and conclusions are summarized in Sects. 6 and 7.

2 Observations

The observations of the HCP (2–1), CP (2–1), PN (2–1), and PO (2–1) transitions in the 3 mm range were carried out at the IRAM 30 m telescope located at Pico Veleta (Spain) towards the line of sight to the compact extragalactic quasar B0355+508. Table 1 lists the observed transitions, the spectroscopic constants, and the telescope settings at the targeted frequencies: the upper state energy is described by E_up, the upper state degeneracy is given by g_u, while A_ul stands for the Einstein coefficient of the transition u → l. The main beam efficiency and the main beam size of the telescope at a given frequency are denoted by the parameters B_eff and θ_MB, respectively.For our observations we used the EMIR receiver with the E090 configuration (3 mm atmospheric window). We applied three observational setups, in which every setup covered a total spectral coverage of 7.2 GHz (each sub-band covered 1.8 GHz). As a backend we used the Fast Fourier Transform Spectrometer with a frequency resolution of 50 kHz (0.15 km s⁻¹ at 100 GHz). In addition, we applied the wobbler switching mode with an amplitude offset of ± 90″. Pointing and focus of the telescope was performed every 2 h on the background source B0355+508 itself and was found to be accurate to within 2″.

The intensity of the obtained spectra was converted from antenna ($T_{\text{A}}^{*}$) to main beam temperature (T_mb) units, using the following relation: $T_{\text{mb}}=\frac{F_{\text{eff}}}{B_{\text{eff}}}\times T_{\text{A}}^{*}$, where F_eff is the forward efficiency. F_eff is equal to 95% in the targeted frequency range.

3 Results

The compact extragalactic source B0355+508 is located at α = 3h 59m 29.73s, δ = 50°57′50.2″ with a low galactic latitude of b = −1.6037°, incorporating a large amount of Galactic gas along the line of sight that harbors up to five diffuse and/or translucent clouds at velocities of −4, − 8, − 10, − 14 and −17 km s⁻¹ (e.g., Liszt et al. 2018, and references therein). The flux of the blazar B0355+508 is variable over time and has been measured at ~ 3 mm to be on average equal to (4.62 ± 1.02) Jy, after averaging the flux of 76 different observations (Agudo 2017). This corresponds to a temperature T_c of (0.96 ± 0.21) K at a beam size of 27″ by taking into account the Rayleigh-Jeans-Approximation. The obtained spectra were reduced and analyzed by the GILDAS software (Pety 2005). Every detected line was fitted via the standard CLASS Gaussian fitting method. For the derivation of the peak opacity we use the radiative transfer equation $$\begin{array}{lll}{T}_{\text{mb}}\hfill & =\hfill & (J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})-J_{\nu }(T_{\text{c}}))\times (1-\exp (-\tau ))\Rightarrow \hfill \\ \tau \hfill & =\hfill & -\ln \left(1-\frac{T_{\text{mb}}}{J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})-J_{\nu }(T_{\text{c}})}\right),\hfill \end{array}$$(1)

where T_ex is the excitation temperature, T_bg is the cosmic background temperature, and $J(T)=(\frac{h\nu }{k_B}){(e^{\frac{\text{h}\nu }{k_{B}T}}-1)}^{-1}$ describes the Rayleigh-Jeans temperature in Kelvin¹. After obtaining the peak opacity τ, the column density N is then estimated by following the relation: $$N=\tau \sqrt{\frac{16{\pi }^3}{\ln 2}}\frac{{\nu }^3{Q}_{\text{rot}}(T_{\text{ex}})\Delta v{e}^{E_u/k_{\text{B}}{T}_{\text{ex}}}}{c^3{A}_{ul}{g}_u(e^{h\nu /k_{\text{B}}{T}_{\text{ex}}}-1)},$$(2)

with k_B being the Boltzmann constant, Δv the line width (FWHM), ν the transition frequency, c the speed of light, and h the Planck constant. Q_rot(T_ex) gives the partition function of a molecule at a given excitation temperature T_ex.

The (2–1) transitions of HCP, CP, PN, and PO were not detected within our observations (see Fig. 1). We derive 3σ upper limits for the opacities and column densities of the P-bearing species using Eqs. (1) and (2). Due to the low densities in diffuse clouds, molecules are expected to show no collisional excitation. The column densities were calculated assuming T_ex = T_bg = 2.7 K, which simplifies Eq. (1) to: $$\begin{array}{lll}\tau \hfill & =\hfill & -\ln \left(1+\frac{T_{\text{mb}}}{J_{\nu }(T_{\text{c}})}\right).\hfill \end{array}$$(3)

The results are summarized in Table 2.

We detected the HNC (1–0), CN (1–0), and C³⁴S (2–1) transitions in absorption as well as the ¹³CO (1–0) in emission at the 3 mm range with a high signal-to-noise ratio (S/N), ranging from 6 to 80². Figure 2 shows all the detected spectra towards the line of sight to B0355+508. In the case of CN we were able to detect and resolve four hyperfine components from 113.123 to 113.191 GHz (see Fig. 3). Every hyperfine component was detected in the three velocity components at − 8, − 10, and − 17 km s⁻¹ except for the one weak transition at 113.123 GHz, which was identified only in two clouds (at − 10, − 17 km s⁻¹). The molecule HNC was identified in all five cloud components, while C³⁴ S (in absorption) and ¹³ CO (in emission) were detected solely towards the densest features, at − 10 and − 17 km s⁻¹. Table 3 lists the identified species and the corresponding spectroscopic parameters.

For estimating the CN column density we use the hyperfine component at 113.170 GHz. Our derived opacities and column densities of CN agree within a factor of two to three with previous results (Liszt & Lucas 2001), while the HNC results are well reproduced within a factor of 1.5. Table 4 summarizes the derived opacities and column densities of the detected species, as well as the obtained line intensities and rms levels.

Lucas & Liszt (2002) reported the detection of the main isotopolog C³²S with the IRAM Plateau de Bure interferometer (PdBI) and estimated a column density of (4.27 ± 0.16) × 10¹² cm⁻² at the − 10 km s⁻¹ component and (3.06 ± 0.32) × 10¹² cm⁻² at − 17 km s⁻¹. With the above values and the column densities of C³⁴S calculated in this work, we derive a sulfur isotopic ratio ³²S∕³⁴S of 12.8 ± 4.8 and 18.7 ± 9.5 for the components at − 10 and − 17 km s⁻¹, respectively. The latter value is in good agreement with the ³²S∕³⁴S ratio for the local ISM of 24 ± 5 (Chin et al. 1996). However, the isotopic ratio determined for v = −10 km s⁻¹ is significantly lower than the local interstellar value, which could be the result of opacity effects of the C³²S line. In addition, the determination of the sulphur isotopic ratio was based on just one spectral line of the main species C³²S and its isotopolog, which also yields a high uncertainty. To our knowledge this is the first detection of C³⁴S towards this line of sight, owing to the high spectral resolution of ~50 kHz and high sensitivity (rms of ~4 mK) achieved with our observations.

In Liszt & Lucas (1998), detections of the main species ¹²CO and its isotopolog ¹³CO are reported, which were obtained with the PdBI as well as the NRAO 12 m telescope. The single-dish observations covered a large beam of 60′′, thus seeing CO and its isotopolog in emission, while the interferometric observations were sensitive only to the very narrow column of gas towards the strong background blazar, giving rise to absorption lines. For deriving the excitation temperatures and the column densities, both emission (single-dish data) and absorption lines (interferometric data) were considered.

In Fig. 2 it is clearly visible that the strong ¹³ CO emission line at − 10 km s⁻¹ overlaps with an absorption feature at around − 8 km s⁻¹. This is probably due to the fact that absorption is present close to the background source, meaning that emission and absorption lines are merged together in our observations with the IRAM 30 m telescope. This contamination effect is influencing the line profile at − 10 km s⁻¹ which subsequently results in an unreliable fit. This could possibly explain why the ¹³ CO column density derived at − 10 km s⁻¹ deviates by a factor of about four from previous results (Liszt & Lucas 1998), while towards − 17 km s⁻¹, N(¹³CO) is well reproduced within 10% (see Table 4)³. Our derived isotopic ratio ¹²CO∕¹³CO at − 17 km s⁻¹ is equal to 16.7 ± 1.4. For this calculation, we used the column density of ¹²CO derived in Liszt & Lucas (1998) with N(¹²CO) = (6.64 ± 0.47) × 10¹⁵ cm⁻². The resulting CO isotopic ratio is almost a factor of four lower than the local interstellar ratio ¹² C∕¹³C = 60 (Lucas & Liszt 1998). This was already confirmed by previous studies (Liszt 2007, 2017) that show an increased insertion of ¹³ C into CO towards clouds in the translucent regime with elevated densities and/or smaller radiation fields, which lead to an enhanced abundance of ¹³CO by a factor of two to four. Under these conditions isotope exchange fractionation (¹³ C⁺ + ¹²CO →¹²C⁺ + ¹³CO + 35 K) is more dominant than selective photodissociation.

Fig. 1 Spectra of the nondetected (2–1) transitions of PO, PN, HCP, and CP. The upper x-axis shows the rest frequency (in MHz) and the lower one is a velocity axis (in km s−1). The red dashed line indicates the 3σ level and the blue dashed line shows the transition frequency of the corresponding molecule. In the case of PO, we show as an example one of the observed transitions at 108.998 GHz.

Fig. 2 Spectra of the detected species HNC, CN, C34S, and 13 CO in the 3 mmrange towards the line of sight to the extragalactitc source B0355+508. The red line represents the CLASS Gaussian fit.

Fig. 3 Detected hyperfine components of the CN (1–0) transition between 113.12 and 113.20 GHz. The three strongest hyperfine components were detected in the three clouds with vLSR =−8, − 10, − 17 km s−1 except for the one weak transition (N, F) = (1, 1∕2) - (0, 1∕2), which was identified only in the two densest clouds (at − 10, − 17 km s−1).

4 Chemical modeling

The goal of the present study is to constrain and improve our model of diffuse and translucent clouds to make reliable predictions regarding the abundances of P-bearing species (and also others). For this reason, we used the observations of HNC, CN, CS, and CO in order to constrain the physical parameters in our model. The chemical code that we applied was developed by Vasyunin & Herbst (2013) with an updated grain-surface chemistry (Vasyunin et al., in prep.). The model includes a gas-grain chemical network with 6000 gas-phase reactions, 200 surface reactions, and 660 species. Accretion and desorption processes regulate and connect the gas-phase and grain surface chemistry. The code numerically solves coupled differential equations (chemical rate equations) and computes a set of time-dependent molecular abundances. Since the observations were carried out towards diffuse and translucent clouds, we considered as initial elemental abundances the standard Solar elemental composition (see Table 5). We note that our initial elemental abundances are significantly different compared to the low metal abundances used in Wakelam et al. (2015) for dark clouds (200 times more abundant S and up to 10⁴ more abundant Fe, Cl, P, and F). In particular, the initial abundance of P is 2.6 × 10⁻⁷ and is therefore well constrained unlike in dense molecular clouds. This approach will help us better elucidate the chemistry ofP since a key parameter for the chemical model is well determined. In addition, we begin our chemical simulations with the entirety of hydrogen being in its atomic form in order to start with pure atomic diffuse cloud conditions.

4.1 The chemical network of phosphorus

The phosphorus chemical network that has been used in previous studies (Fontani et al. 2016; Rivilla et al. 2016) has been extended with new available information in the literature (new reactions, updated reaction rates, desorption energies, etc.). In particular, chemical reactions of several P-bearing species, such as PN, PO, HCP, CP, and PH₃, were included and/or updated in our chemical network. The reaction rates were taken from the online chemical databases KInetic Database for Astrochemistry (Wakelam et al. 2015, KIDA)⁴ and the UMIST Database for Astrochemistry (McElroy et al. 2013, UDfA)⁵, as well as from numerous previous papers (Thorne et al. 1984; Adams et al. 1990; Millar 1991; Anicich 1993; Charnley & Millar 1994; Jiménez-Serra et al. 2018). In particular we included several reactions involving the formation and destruction of PH_n (n = 1, 2, 3) and their cationic species from Charnley & Millar (1994) and Anicich (1993), along with the chemical network proposed by Thorne et al. (1984) that contains production and loss routes for P, PO, P⁺, PO⁺, PH⁺, HPO⁺ and H₂ PO⁺. In addition, weextended the PN chemical network based on the work by Millar et al. (1987), and we took into account the gas-phase reaction P + OH → PO + H proposed byJiménez-Serra et al. (2018), as well as two formation routes of PN in the gas phase, N +CP → PN + C and P + CN → PN + C, by Agúndez et al. (2007). Finally, we included the photodissociation reactions of PN, PO, HCP, and PH_n (n = 1, 2, 3) based on the reaction rates given in KIDA and UDfA. The reaction rates of the photodissociation of PH_n were assumed to be equal to the analogous reactions for NH_n.

Concerning the chemistry taking place on grain surfaces, we took into account the hydrogenation reactions of P-bearing species (where the letter “g” denotes a grain surface species) as well as their corresponding desorption reactions:

• gH + gP → gPH,

• gH + gPH → gPH₂,

• gH + gPH₂ → gPH₃.

The desorption energy of PH₃ was calculated based on that of NH₃ and amounts to ~ 5800 K. This corresponds to an evaporation temperature of ~100 K, which is in good agreement with the value of ~90 K proposed by Turner et al. (1990)⁶. The reactive desorption efficiency in our chemical model is set equal to 1%. An increased reactive desorption of 10% changes the predicted abundances of the aforementioned P-bearing molecules by less than a factor of two. Another nonthermal desorption mechanism included in our model is the cosmic-ray desorption, which is fully described in Hasegawa & Herbst (1993). Based on this study, dust grains are heated upon impact with cosmic rays reaching a peak temperature T_dust of 70 K, which subsequently leads to preferential desorption of molecules from grain surfaces. This type of desorption is however negligible in diffuse clouds, where photodesorption dominates. In our model we adopt for all species a photodesorption rate of 3 × 10⁻³ molecules per incident UV photon, as was determined in Öberg et al. (2007) based on laboratory measurements of pure CO ice. This stands in good agreement with the photodesorption yield of ~ 10⁻³ molecules per UV photon found for other species, such as H₂O, O₂ and CH₄ (Öberg et al. 2009; Fayolle et al. 2013; Dupuy et al. 2017).

4.2 Comparison to observations

In order to reproduce the observed abundances of HNC, CN, CS, and CO in every cloud towards the line of sight to B0355+508 we produce a grid of models applying typical physical conditions for diffuse or translucent gas (Snow & McCall 2006; Thiel et al. 2019). We note here that, since our chemical model is not treating isotopic species, we are using as a reference for our comparison the main species ¹²CO and C³² S instead of ¹³CO and C³⁴S. For the fractional abundances of ¹²CO and C³²S, we are adopting the column densities determined in Liszt & Lucas (1998) and Lucas & Liszt (2002). In addition, for the clouds at − 14 and − 4 km s⁻¹ we use for CN and CS the upper limits derived in this work (N(CN) < 10¹² cm⁻²) and in Lucas & Liszt (2002). The parameter space that we investigate is listed below:

• n(H) = 100−1000 cm⁻³, spacing of 100 cm⁻³,

• A_V = 1−5 mag, spacing of 1 mag,

• T_gas = 20−100 K, spacing of 10 K.

The chemical evolution in each model is simulated over 10⁷ yr (100 time steps) assuming static physical conditions. For the cosmic-ray ionization rate ζ(CR) we use a value of 1.7 × 10⁻¹⁶ s⁻¹ (Indriolo & McCall 2012, see Sect. 5 for further explanation). This also corresponds to the values applied in Godard et al. (2014) and Le Petit et al. (2004), where the best-fit models provided a ζ(CR) of 10⁻¹⁶ s⁻¹ and 2.5 × 10⁻¹⁶ s⁻¹, respectively.

Given the above parameter space, we calculate the level of disagreement D(t, r) between modeled and observed abundances (for the species HNC, CN, CS and CO), which, following Wakelam et al. (2010) and Vasyunin et al. (2017), we define as $$D(t,r)={\sum }_{j=1}^{N_{\text{species}}}\log (x_{\text{mod}}^j(t,r))-\log (x_{\text{obs}}^j),$$(4)

with r = (n(H), A_V, T_gas) and $x_{\text{obs},\text{mod}}^j$ is the observed or modeled abundance of species j, respectively. We then determine the minimal value of D(t, r), (noted as D_min(t, r)), that corresponds to the best-fit model; the parameters (t, r) provided by the best-fit model are the ones giving the smallest deviation between observations and predictions. The smaller the D_min(t, r), the better the agreement.

According to Pety et al. (2008) and references therein, the clouds with v_LSR = −10 and − 17 km s⁻¹ towards B0355+508 show strong ¹² CO emission lines that originate from dense regions (with n(H) = 300−500 cm⁻³ and N(H₂) > 10²¹ cm⁻²) which are just outside the synthesized beam (combination of the 30 m and PdBI telescopes), but still within the IRAM 30 m beam (Pety et al. 2008). Based on the CO (2–1) maps shown in Pety et al. (2008) with 22′′ and 5.8′′ resolutions, respectively, the cloud at v_LSR = −10 km s⁻¹ shows the most pronounced, dense substructure. This is confirmed by the fact that this particular cloud component produces the most detectable amounts of observed species in the data presented in this work and previous works by Liszt & Lucas (2001) and Lucas & Liszt (2000), and is therefore chemically the most complex one. The component at v_LSR = −8 km s⁻¹ shows a similar structure to the one at v_LSR = −10 km s⁻¹, incorporating a dense region as well. Pety et al. (2008) suggest that the two components are part of the same cloud, even though they are distinguishable in absorption and show different levels of chemical complexity. With a higher spatial resolution of 5.8′′, the CO emission at v_LSR = −8 km s⁻¹ is separated from the one at v_LSR = −10 km s⁻¹ and is clearly visible. The diffuse gas seen at velocities − 14 and − 4 km s⁻¹ shows barely any ¹³ CO (Liszt & Lucas 1998, and this work) or ¹²CO (Pety et al. 2008; Liszt & Lucas 1998) in emission, which suggests that the density of these clouds is too low to sufficiently excite CO. Pety et al. (2008) estimated a low to moderate density of these clouds to be ~ 64−256 cm⁻³ with A_V < 2 mag.

Since we are performing single-dish observations with a beam of ~ 22′′, we also cover the high-density regions that produce significant ¹² CO emission (Pety et al. 2008). For this reason we constrain our grid of models to high densities of ≥ 300 cm⁻³ for the denser clouds (− 8, − 10, − 17 km s⁻¹), while for the low-density objects at − 14 km s⁻¹ and − 4 km s⁻¹ we restrict our input parameters to n(H) ≤ 200 cm⁻³ and A_V < 2 mag. The calculation of the molecular abundances was done with respect to the H₂ column densities (N(H₂) ~ 4−5 × 10²⁰ cm⁻²) that were derived towards every cloud by Liszt et al. (2018). However, our model provides the fractional abundance of a species X with respect to the total number of hydrogen nuclei, as in n(X)∕n(H), with the total volume density of hydrogen defined as n(H) = n(H I + 2 × H₂). The surface mobility parameters that we set as default values in our model (see Sect. 5.3) enable fast and effective formation of H₂ on the surfaces of grains. At the end of our simulations (at t = 10⁷ yr), the H₂ abundance reaches a value of 40−50% (depending on the set of parameters). This means that almost the entire hydrogen is predicted to be in its molecular form at the late phases of the chemical evolution. Following this consideration, we divide all the observed abundances and their upper limits mentioned in this paper by a factor of two, since n(X)∕n(H I + 2 × H₂) ≃ n(X)∕2 × n(H₂) = 0.5 × n(X)∕n(H₂). We note that this expression applies to the high-density parts of the clouds and does not account for the low-density (and H I rich) gas along the line of sight.

Figure 4 shows the results of the grid of models that was applied to reproduce the observations towards v_LSR = −17 km s⁻¹. In particular, we plot D(t_best, r), which describes the deviation between observed and modeled abundances at the time of best agreement t_best, versus the density, temperature, and visual extinction. Between an A_V of 1 and 3 mag the smallest level of disagreement D(t_best, r) reduces by 13%. The main discrepancy between observed and modeled abundances at low A_V comes from the fact that a high visual extinction results in higher molecular abundances and is therefore able to reproduce the chemical complexity seen towards the translucent clouds. For models with A_V > 3 mag the minimal D(t_best, r) barely changes (less than 1% of increase). The smallest D(t_best, r) increases with respect to the density and temperature up to 3 and 2%, respectively. This is a clear indication that the most influential physical parameter in our analysis is the visual extinction.

For the cloud component at v_LSR = −17 km s⁻¹ the best-fit model with D_min(t, r) is reached ata time t_best = 6.2 × 10⁶ yr and has the parameters: r_best = (n(H), A_V, T_gas) = (300 cm⁻³, 3 mag, 40 K). At t_best = 6.2 × 10⁶ yr we also fulfill the assumption of having most of the hydrogen in molecular form, as the H₂ abundance reaches a value of 0.45. Based on this model, we show in Fig. 5 the time dependent abundances of CO, CN, CS and HNC over 10⁷ yr as well as the corresponding observed abundances towards v_LSR = −17 km s⁻¹. Our chemical model reproduces the observed species CO, CN, and CS very well within a factor of ~ 1−1.4 at t_best = 6.2 × 10⁶ yr. As can be seen in Fig. 5, the predicted abundances follow the order of the observed quantities. The most significant discrepancy is found in case of HNC, where the chemical model underestimates the observed abundance by a factor of four at the time of best agreement (see Table 6). According to the model, one of the main destruction mechanisms of HNC is: C⁺ + HNC → C₂N⁺ + H; based on the online chemical databases, its reaction rate remains uncertain. In UDfA the given reaction rate was determined theoretically (Leung et al. 1984) and therefore entails a high uncertainty. Experimental studies of the above chemical route are still needed to make reliable predictions of the HNC abundance.

The smallest deviation with the observations towards the cloud with v_LSR = −10 km s⁻¹ was produced by the same set of parameters: r_best = (n(H), A_V, T_gas) = (300 cm⁻³, 3 mag, 40 K). However, in this case the best-fit model gives a D_min(t, r) that is slightly larger by a factor of approximately 1.4. The molecular abundances observed towards v_LSR = −8 km s⁻¹ are best reproduced with an A_V of 5 mag, a density of 400 cm⁻³ and a gas temperature of 40 K. The smallest level of disagreement between an A_V of 3 and 5 mag differs by less than 1%.

For the two remaining clouds with v_LSR = −4 km s⁻¹ and v_LSR = −14 km s⁻¹ the best-fit model in both cases is given by the parameters: (n(H), A_V, T_gas) = (200 cm⁻³, 1 mag, 30 K) at t_best = 10⁷ yr. Here, the discrepancy in both clouds arises mostly from the fact that the model underestimates the CS abundance by a factor of between approximately six and nine. In our model, CS is being effectively destroyed via photodissociation due to the low A_V. Table 7 lists the best-fit parameters that were determined towards every cloud component.

We note that towards the same line of sight there have been detections of several other molecules, as reported in Liszt et al. (2008). The best-fit model determined towards v_LSR = −17 km s⁻¹ and v_LSR = −10 km s⁻¹ is able to reproduce within one order of magnitude the species OH, C₂H, H₂ CO, NH₃ and CH, while other species such as HCN, SO, H₂S and C₃H₂ are strongly underestimated by up to two orders of magnitude. This is a clear indication that the chemical network of certain molecules (other than P-bearing ones) still needs to be extended and updated. This however will be addressed in future work, as the present paper focuses mainly on P chemistry.

Fig. 4 Results of the grid of models applying typical conditions for diffuse or translucent clouds in order to reproduce the observations towards the cloud at vLSR =−17 km s−1. The deviation between observations and model at the time of best agreement tbest is given by D(tbest, r), which is plotted versus density, temperature, and visual extinction. The best-fit model is given at a time tbest = 6.2 × 106 yr and has the following parameters: (n(H), AV, Tgas) = (300 cm−3, 3 mag, 40 K).

Fig. 5 Chemical evolution of the abundances of CO, CN, CS and HNC over 107 yr predicted by our best-fit model with the parameters (n(H), AV, Tgas) = (300 cm−3, 3 mag, 40 K). The colored horizontal bands correspond to the observed abundances towards the cloud with vLSR =−17 km s−1, including the inferred uncertainties. The vertical dashed line indicates the time of best agreement (t = 6.2 × 106 yr) between observations and model results.

5 Discussion: the chemistry of phosphorus

Based on the above results we can conclude that the molecular abundances observed at v_LSR = −8, − 10, − 17 km s⁻¹ can be best reproduced by a more “shielded” (A_V > 1 mag) interstellar medium that allows the build-up of molecules to occur more efficiently. The resulting visual extinction A_V of 3 mag should be viewed as an average value over the region covered by our beam (~ 22′′). Within this region the denser clumps are most likely translucent in nature. Hence, the observed cloud components are probably heterogeneous clouds, incorporating diffuse and translucent material, filled with relatively abundant molecules. This result stands in good agreement with a study by Liszt (2017), which involved modeling the CO formation and fractionation towards diffuse clouds. One of the main results of this latter study was that strong ¹³ CO absorption lines observed in the millimeter- and UV-range can be explained by higher densities (≥ 256 cm⁻³) and weaker radiation (and thus higher visual extinction), as already mentioned in Sect. 3. Our conclusions also agree well with the work done by Thiel et al. (2019) in which the physical and chemical structure of the gas along the line of sight to SgrB2(N) was studied; here, complex organic molecules such as NH₂CHO and CH₃CHO were detected in the majority of the clouds, which at the same time proved to have relatively high visual extinctions (A_V = 2.5−5 mag with N(H₂) > 10²¹ cm⁻²), thus consisting mainly of translucent gas. According to Thiel et al. (2019) the column density of H₂ that corresponds to an A_V of 3 mag is ~3 × 10²¹ cm⁻². This is also consistent with the study by Pety et al. (2008), which states that the bright ¹² CO emission originates from dense regions with N(H₂) > 10²¹ cm⁻². The gas observed at velocities of − 14 and − 4 km s⁻¹ on the other hand, corresponds mainly to a “classical” diffuse cloud with a visual extinction of ~ 1 mag according to the above analysis. These clouds also yielded the smallest amounts of the detected molecular abundances. Since chemical complexity seems to be favored towards translucent rather than diffuse gas, for the following discussion we use the model that provided the best fit towards the dense clouds with v_LSR = −17 km s⁻¹ and v_LSR = −10 km s⁻¹ as a reference.

According to our best-fit model, P⁺ has a gas-phase abundance of 1.8 × 10⁻⁷ at the end of our simulations, being a factor of approximately 1.4 lower than its cosmic value, which indicates that little depletion takes place. The main reservoir of phosphorus other than P⁺ is atomic P, having an abundance of 7.4 × 10⁻⁸ at 10⁷ yr. Atomic P is formed mainly through the electronic recombination of P⁺. For our models with elevated densities (10³ cm⁻³) we reach high elemental depletion (such as for C⁺, S⁺ and P⁺) after running the code for 10⁷ yr. This is consistent with the results presented by Fuente et al. (2019), which show significant depletion of C, O, and S happening already towards translucent material at the edge of molecular clouds (3–10 mag) with 1− 5 × 10³ cm⁻³. In Appendix A we investigate further the expected depletion of phosphorus when transitioning from diffuse- to dense-cloud conditions. We find that there is a significant depletion of atomic P on dust grains after the final density of 10⁵ cm⁻³ is reached. This in turn leads to a strong increase of gPH₃, that becomes the main carrier of phosphorus in the dense phase. We also find a considerable decrease of the molecules HCP, CP, PN, PO and PH₃ due to freeze-out on grains and their destruction route with ${\text{H}}_{\text{3}}^{\text{+}}$ after the final density is attained at t ~ 10⁶−10⁷ yr.

The most abundant P-bearing molecules in the gas phase are HCP and CP, with maximal abundances of 3.4 × 10⁻¹⁰ and 2.1 × 10⁻¹⁰, respectively⁷. The formation and destruction pathways of both HCP and CP are strongly related to the electron fraction, as they are mainly produced (throughout the entire chemical evolution) by dissociative recombination of the protonated species ${\text{PCH}}_{\text{2}}^{\text{+}}$ and destroyed by reacting with C⁺, the main carrier of positive charge in diffuse clouds. Two additional P-bearing species that are predicted by our model to have “observable” abundances in the gas phase are PN and PO, with respective maximal abundances of 4.8 × 10⁻¹¹ and 1.4 × 10⁻¹¹. The most productive formation pathways for PN start with P + CN → PN + C, N + PH → PN + H and end with N + CP → PN + C. In the late stage of evolution (~ 0.5 × 10⁶ − 10⁷ yr) PN is primarily being destroyed by He⁺ + PN → P⁺ + N + He. The species PO is mainly produced over the entire chemical evolution of 10⁷ yr by the dissociative recombination of HPO⁺: HPO⁺ + e⁻→ PO + H, and is mostly destroyed by reactions with C⁺ and H⁺. On the other hand, HPO⁺ is efficiently formed via P⁺ + H₂O → HPO⁺ + H⁸. An additional reaction that becomes relevant at progressive times (~ 10⁶−10⁷ yr) is O + PH → PO + H with a ~ 10% reaction significance.

Another relatively abundant P-bearing species in the gas-phase based on our best-fit model is phosphine, PH₃, with a maximal abundance of ~1.6 × 10⁻¹¹ at a late time of 10⁷ yr. We note here that the species PH is also predicted to be detectable with a maximal abundance of ~ 3.6 × 10⁻¹¹. Unlike PH however, PH₃ has already been detected in circumstellar envelopes of evolved stars (Agúndez et al. 2014a), indicating that it could be an important P-bearing species in interstellar environments such as diffuse and translucent clouds. We therefore focus in the following sections on the PH₃ rather than the PH chemistry. Based on our chemical model, PH₃ is formed most efficiently on dust grains in the early phase, being released to the gas-phase via reactive desorption: gH + gPH₂ → PH₃. Its formation proceeds after 1.4 × 10³ yr with the photodesorption process gPH₃ → PH₃ being the most effective reaction. Since the evaporation temperature of PH₃ lies at ~ 100 K, the main mechanism driving the desorption of PH₃ at low temperatures is photodesorption (instead of thermal desorption). Switching off the photodesorption in our model leads to a decrease of the PH₃ gas-phase abundance of two orders of magnitude. Once in the gas phase, PH₃ is mostly destroyed by reactions with C⁺ and H⁺ as well as through the photodissociation reaction: PH₃ + hν → PH₂ + H. The most abundant species on grains is gPH₃, with a maximal abundance of 7.2 × 10⁻¹⁰. Almost all the atomic P that depletes onto the dust grains reacts with gH and forms gPH (gP + gH → gPH), which subsequently forms gPH₃ through further hydrogenation. Table 8 summarizes all the main formation and destruction pathways for the molecules PN, PO, HCP, CP and PH₃ at three different times (t = 10³, 10⁵, 10⁷ yr). The last column shows the significance of the given reaction in the total formation or destruction rate of the species of interest.

Figure 6 depicts the time-dependent abundances of PN, PO, HCP, CP and PH₃ over 10⁷ yr predicted by the best-fit model along with the computed 3σ upper limits. The predicted abundance for PO lies a factor of about 40 below the observational upper limit at t = 10⁷ yr, while the current upper limits of HCP and CP are about one order of magnitude higher than the model predictions. Finally, for PN the modeled abundance almost reaches the observed value at the end of our simulations. This means that in all cases the predicted abundances of P-bearing species are lower than the derived upper limits. Future observations of the ground-energy transitions (1–0) will help us to constrain these upper limits even more (see Sect. 6 for further justification). Table 9 lists the predicted abundances of the above species given by our chemical model at t = 10⁷ yr along with the corresponding upper limits. In the case of PO we show only the lowest value of the four upper limits that were derived for each transition.

In the following discussion we focus on how deviations from our best-fit model can affect the chemistry of P-bearing species. In particular, we examine the dependence of the abundances of HCP, CP, PN, PO and PH₃ on increasing visual extinction A_V, increasing cosmic-ray ionization rate ζ(CR), and alternating the surface mobility constants (diffusion-to-desorption ratio E_b ∕E_D and possibility of quantum tunneling for light species).

5.1 Effects of visual extinction on the P-bearing chemistry

In this section, we analyze how an increase of A_V is affecting the predicted abundances of P-bearing species. For this purpose we consider the parameters of the best-fit model with n(H) = 300 cm⁻³ and T_gas = 40 K, while varying the A_V from 1 to 10 mag. By keeping the density constant, we avoid high levels of elemental depletion. The increase in visual extinction can then be explained by a figurative increase of the size of the source. Figure 7 (left panel) shows the predicted abundances of P-bearing species at the end of our simulations (t = 10⁷ yr) under the effect of varying the visual extinction. All species reach a maximal abundance at an of 4 mag. The abundances of HCP, CP, and PN barely change for A_V > 4 mag, while for the rest of the molecules the abundances drop; especially in the case of PH₃ where we see a substantial decrease of almost two orders of magnitude. As already mentioned in Sect. 5, the most effective formation process of PH₃ is the photodesorption gPH₃ → PH₃. Thus, a high visual extinction attenuates the incoming UV-field and therefore the desorption of gPH₃.

In order to better understand the A_V dependence of the remaining molecular abundances, we plotted in Fig. 7 the predicted abundances of the species that mainly form and destroy HCP, CP, PN, and PO (see Table 8) as a function of the visual extinction. In particular, we simulated the abundances of ${\text{PCH}}_{\text{2}}^{\text{+}}$, C⁺, P⁺, He⁺, and H⁺ as well as HPO⁺. In the case of HCP (and also CP) its abundance increases up to an A_V of 4 mag and then remains constant above that value. This behavior is correlated with the increase of the ${\text{PCH}}_{\text{2}}^{\text{+}}$ abundance up to an A_V of 3 mag as well as the decreaseof C⁺ up to a visual extinction of 4 mag. The species PO seems to be more strongly affected by the increasing A_V. Its abundance will also increase for A_V ≤ 4 mag which again stands in correlation with the decrease of the C⁺ abundance (the main “destroyer” of PO), followed by a drop in abundance up to 7 mag. This on the other hand results from the decrease in HPO⁺ abundance (the main precursor of PO) in the same A_V range. An increase in A_V will decrease the P⁺ abundance (due to the decrease of the total ionization rate), as can be seen in Fig. 7. In addition, an enhanced A_V slightly decreases the H₂O abundance (by a factor of two), since the most effective formation for H₂ O at late times is the photodesorption gH₂O → H₂O (see footnote 7). Therefore, for higher A_V, both P⁺ and H₂ O decrease, meaning that HPO⁺ and subsequently PO reduce in abundance as well.

5.2 Effects of the cosmic-ray ionization rate on the chemistry of P-bearing species

As already mentioned in Sect. 4.2, for all the applied models we use a value of 1.7 × 10⁻¹⁶ s⁻¹ for the cosmic-ray ionization rate ζ(CR), as was derived by Indriolo & McCall (2012). This is also consistent with previous work in which diffuse and translucent clouds were studied as well (Fuente et al. 2019; Godard et al. 2014; Le Petit et al. 2004). However, we should note here that in Indriolo & McCall (2012) several cosmic-ray ionization rates were derived towards 50 diffuse lines of sight, ranging from 1.7 × 10⁻¹⁶ to 10.6 × 10⁻¹⁶ s⁻¹ with a mean value of 3.5 × 10⁻¹⁶ s⁻¹. Due to the complex and not yet fully understood nature of our observed clouds, we test our chemical model by also applying the elevated values of ζ(CR) = 3.5 × 10⁻¹⁶ s⁻¹ and 10.6 × 10⁻¹⁶ s⁻¹ in order to examine the influence of the cosmic-ray ionization rate on P-bearing chemistry. As for the remaining parameters of the code (such as A_V and T_gas) we use the values given by our best-fit model (see Sect. 4.2). Figure 8 shows the chemical evolution of the species PN, PO, HCP, CP and PH₃ over 10⁷ yr for ζ(CR) = 1.7 × 10⁻¹⁶ s⁻¹ and 10.6 × 10⁻¹⁶ s⁻¹ in the left and right panels, respectively.

Table 10 summarizes the predicted abundances of P-bearing species for the three different cosmic-ray ionization rates given in Indriolo & McCall (2012). As one can see that PN shows the most substantial decrease in abundance with increasing ζ(CR). From the lowest (ζ(CR) = 1.7 × 10⁻¹⁶ s⁻¹) to the highest (ζ(CR) = 10.6 × 10⁻¹⁶ s⁻¹) cosmic-ray ionization rate, the PN abundance decreases by a factor of approximately 730, while for HCP, CP, and PH₃ we see a drop by a factor of about 85 and 50, respectively. As already mentioned, PN is heavily destroyed by He⁺ with a ~40% reaction significance. An increase of ζ(CR) up to a value of 10.6 × 10⁻¹⁶ s⁻¹ significantly enhances the ionization of He and H by a factor of about 20 and 30, respectively, via cosmic-ray-induced secondary UV photons: He + CRP → He⁺ + e⁻ and H + CRP → H⁺ + e⁻ (the abundance of C⁺ increases by a factor of ~ 6). Therefore, the destruction path with H⁺ becomes relevant for all P-bearing species showing a 10–40% loss efficiency. The effect is the strongest in the case of PN, because PN is mainly formed through CP which drastically decreases and is also efficiently destroyed by both He⁺ and H⁺. The PO abundance is only reduced by a factor of about 20 after increasing ζ(CR) up to 10.6 × 10⁻¹⁶ s⁻¹, despite being heavily destroyed by H⁺. On the other hand, the significance of the dissociative recombination of HPO⁺ increases up to 50% which in turn counterbalances the loss through H⁺. An increased ζ(CR) of 10.6 × 10⁻¹⁶ s⁻¹ enhances the abundance of P⁺ up to ~2.5 × 10⁻⁷, nearly reaching its cosmic value of ~2.6 × 10⁻⁷ (Asplund et al. 2006), while the abundance of atomic P decreases down to ~ 9.5 × 10⁻⁹ via the enhanced reaction with C⁺ and H⁺.

Fig. 6 Variation of the predicted abundances of PN, PO, HCP, CP and PH3 over 107 yr in our best-fit model. The dashed lines represent the 3σ upper limits derived from the observations at vLSR = −17 km s−1. In the case of PO we use 5 × 10−10 as an upper limit (see Table 9 and text for explanation).

5.3 Effects of the diffusion-to-desorption ratio on the chemistry of P-bearing species

The chemistry in the ISM is heavily influenced by the presence of dust grains (Caselli & Ceccarelli 2012). The mobility of the depleted species on the surface of dust grains depends on two mechanisms: thermal hopping and quantum tunneling for the lightest species H and H₂ through potential barriers between surface sites (Hasegawa et al. 1992). Without the possibility of tunneling, the species are not able to scan the grain surface quickly at low temperatures and the total mobility decreases. The parameters that strongly determine the surface chemistry are the diffusion-to-desorption energy ratio E_b ∕E_D as well as the thickness of the potential barrier between adjacent sites. Based on previous studies (Hasegawa et al. 1992; Ruffle & Herbst 2000; Garrod & Herbst 2006), Vasyunin & Herbst (2013) proposed three different values for the E_b ∕E_D ratio: 0.3, 0.5, and 0.77. In the case of low ratios (E_b∕E_D = 0.3) we activate in our model the possibility of quantum tunneling for light species, while for the other two cases, surface mobility is only controlled by thermal hopping (and quantum tunneling is deactivated). The potential barriers are assumed to have rectangular shape and a thickness of 1 Å (Vasyunin & Herbst 2013). In our model we utilize the first set of parameters (E_b ∕E_D = 0.3, with tunneling), nevertheless, since the chemistry of P-bearing species is still highly uncertain, we examine how the remaining two sets of parameters (E_b∕E_D = 0.5, 0.77, no tunneling) influence the predicted abundances. Table 11 lists the predictions for PN, PO, HCP, CP, and PH₃ as well as H₂ at t = 10⁷ yr for the three different sets of surface mobility parameters proposed in Vasyunin & Herbst (2013). As Table 11 shows, the H₂ abundance decreases by a factor of four by switching from setup 1 (E_b∕E_D = 0.3 with tunneling) to setup 2 (E_b∕E_D = 0.5 no tunneling), and finally experiences a dramatic drop of a factor 50 when increasing the E_b ∕E_D up to 0.77 (overall change of a factor 200 between setups 1 and 3).

The reduction of the H₂ abundance has a significant impact on the formation of ${\text{PCH}}_{\text{2}}^{\text{+}}$ and PH, which affects the PN, PO, HCP and CP abundances through the following reactions:

PN

• ${\text{P}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{PH}}_{\text{2}}^{\text{+}}$

• ${\text{PH}}_{\text{2}}^{\text{+}}{\text{+e}}^{\text{-}}\to \text{PH+H}$

• N + PH → PN + H

PO

• ${\text{P}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{PH}}_{\text{2}}^{\text{+}}$

• ${\text{PH}}_{\text{2}}^{\text{+}}{\text{+e}}^{\text{-}}\to \text{PH+H}$

• O + PH → PO + H

HCP

• ${\text{HCP}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{PCH}}_{\text{2}}^{\text{+}}$

• ${\text{PCH}}_{\text{2}}^{\text{+}}{\text{+e}}^{\text{-}}\to \text{HCP+H}$

CP

• ${\text{HCP}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{PCH}}_{\text{2}}^{\text{+}}$

• ${\text{PCH}}_{\text{2}}^{\text{+}}{\text{+e}}^{\text{-}}\to {\text{CP+H}}_{\text{2}}$

Both PH and ${\text{PCH}}_{\text{2}}^{\text{+}}$ decrease by a factor of approximately 20 when increasing the E_b∕E_D up to 0.77. In addition, the abundance of H⁺ is increased by a factor of about 25, since the reduction of H₂ formation leads to more atomic hydrogen and

subsequently more H⁺. The enhanced H⁺ abundance results in a stronger destruction of all P-bearing species through their reaction with H⁺. The species HCP and CP are also strongly affected by changing the surface mobility parameters, with an overall decrease by a factor of about 70 and 65 in abundance, respectively. In both cases the dissociative recombination of ${\text{PCH}}_{\text{2}}^{\text{+}}$ is essential during the whole chemical evolution for the formation of HCP and CP showing a reaction significance of 30–99%. A decrease of ${\text{PCH}}_{\text{2}}^{\text{+}}$ due to lower H₂ abundance therefore results in reduced HCP and CP formation. The largest effect is seen for PN, where a diffusion-to-desorption ratio of 0.77and no quantum tunneling of light species reduces the PN abundance by a factor of 300 (see Fig. 9).

Besides the effective loss through H⁺, the substantial decrease in PN is also related to the reduction of CP, which is the main precursor of PN at late times. In addition, the reaction N + PH → PN + H is important for PN formation over the entire chemical evolution of 10⁷ yr with a 10–50% formation efficiency (for E_b∕E_D = 0.77 and no tunneling). This means that the reduction of the H₂ abundance decreases PH, which in turn produces less PN. In the case of PO however, the change in abundance between the two extreme cases is only a factor of about 20. Here, the route O +PH → PO + H increases in significance only up to 3% at late times (4 × 10⁶−10⁷ yr), indicating that the decrease of PH will not considerably affect PO production. Furthermore, the reduction of PO due to H⁺ is compensated through its effective formation via the dissociative recombination of HPO⁺. Finally, the abundance of PH₃ decreases only by a factor of 13 in total when changing the surface chemistry constants. Despite being heavily destroyed by H⁺, PH₃ is still sufficiently formed through the photodesorption of gPH₃.

Fig. 7 Predicted abundances of P-bearing molecules as a function of visual extinction AV. The molecular abundances shown here are computed at t = 107 yr. Right panel: predicted abundances of ${\text{PCH}}_{\text{2}}^{\text{+}}$, C+, P+, He+, H+, and HPO+ as they are contributing the most to the formation and destruction of HCP, CP, PN, PO, and PH3 (left panel).

Fig. 8 Chemical evolution of P-bearing molecules as a function of time under the effects of cosmic-ray ionization rates of ζ(CR) = 1.7 × 10−16 s−1 (left panel) and 10.6 × 10−16 s−1 (right panel).

6 Future observations

Thanks to the sensitive observations (rms of ~6 mK) of the (2–1) transitions of HCP, CP, PN and PO we were able to obtain good upper limits for the column densities and abundances of the above species (see Tables 2 and 9) and thus constrain P chemistry. The observations of HNC, CN, CS, and CO helped us to put important constraints on the main physical parameters of the targeted diffuse and translucent clouds, that is, the visual extinction, the density, and the gas temperature. For the prospect of future observations we want to estimate the expected line intensities of the (1–0) transitions of HCP, CP, PN and PO (at ~ 40−65 GHz) based on our new and improved diffuse-cloud model. Since the densities present in diffuse and translucent clouds are too low to show any collisional excitation (T_ex = T_bg = 2.7 K), the (1–0) transitions are expected to be more strongly populated than the higher-energy transition levels. For these calculations, we take into account the nonthermal nature of the blazar emission, meaning that the flux increases with decreasing frequency. In particular, we apply a power law to the emission of the blazar with $\frac{F}{F_0}={(\frac{\nu }{{\nu }_0})}^{-\alpha }$, where F is the flux, ν is the corresponding frequency, and α is the spectral index. By using the fluxes determined in Agudo (2017) at 3 and 1.3 mm we infer a spectral index of α ~ 1.06. Following this, we determine the flux at 7 mm to be ~11 Jy, which in turn corresponds to a temperature of ~26 K with a beam size of 17′′ (at 7 mm with the Green Bank Telescope). As Table 12 shows, the derived peak intensities of the species PN, PO, HCP and CP vary from 10 to 200 mK, making these lines “detectable” with radio telescopes, such as the Green Bank Telescope (GBT) and the Effelsberg Telescope. The capabilities of these instruments will allow us to reach rms levels down to 4 mK and enable possible detections up to the 50σ level. The only exception is PH₃ with a (1–0) transition at 266.944 GHz. The flux of the background source at that frequency based on the above power law is equal to 1.91 Jy. This corresponds to a background temperature T_c of 0.4 K with a beam size of9′′ (with the IRAM telescope), which in the end results in a very weak, nondetectable absorption line.

7 Conclusions

The aim of this work is to understand through observations and chemical simulations which physical conditions favor the production of P-bearing molecules in the diffuse ISM and to what degree. Observing diffuse clouds offers us the opportunity to constrain animportant parameter in our chemical simulations, namely the depletion level of phosphorus (and in general the initial elemental abundances).

We performed single-pointing observations (IRAM 30 m telescope) of the (2–1) transitions of the species PN, PO, HCP and CP at 3 mm towards the line of sight to the bright continuum source B0355+508. None of the above transitions were detected. Nevertheless, the sensitive observations yielding an rms level of ~ 6 mK allowed us to obtain reliable upper limits (see Tables 2 and 9).

We have obtained high S/N detections of the (1–0) lines of HNC, CN, and ¹³ CO between 80 and 110 GHz. We also show a first detection of C³⁴S (2–1) at 96 GHz towards the two densest cloud components at − 10 and − 17 km s⁻¹. Following this, we were able to derive a sulfur isotopic ratio ³²S∕³⁴S of 12.8 ± 4.8 and 18.7 ± 9.5 towards the − 10 and − 17 km s⁻¹ features, with the latter being close to the local interstellar value of 24 ± 5 (Chin et al. 1996). The detected molecular species show the highest abundances towards the two components at − 10 and − 17 km s⁻¹, as already shown in previous studies (e.g., Liszt et al. 2018, and references therein).

Based on the detected molecular abundances, we updated our chemical model in order to provide reliable predictions of abundances and line intensities of P-containing molecules that will serve as a guide for future observations. For this purpose we ran a grid of chemical models, with typical physical conditions of diffuse or translucent clouds, trying to reproduce the observed abundances and upper limits of HNC, CN, CO, and CS in every cloud component along the line of sight (at − 4, −8, −10, −14 and −17 km s⁻¹). For the clouds with v_LSR = −10 km s⁻¹ and −17 km s⁻¹, the best agreement between observed and modeled abundances is reached at a time t_best = 6.2 × 10⁶ yr and at r_best = (n(H), A_V, T_gas) = (300 cm⁻³, 3 mag, 40 K). We chose this set of parameters as a reference for modeling the phosphorus chemistry.

According to our best-fit model mentioned above, the most abundant P-bearing species are HCP and CP (~ 10⁻¹⁰) at a time of t = 10⁷ yr. The species PN, PO, and PH₃ also show relatively high predicted abundances of 4.8 × 10⁻¹¹ to 1.4 × 10⁻¹¹ at the end of our simulations. All species are effectively destroyed through reactions with C⁺, H⁺, and He⁺. The molecules HCP, CP, and PO are efficiently formed throughout the entire chemical evolution via the dissociative electron recombination of the protonated species ${\text{PCH}}_{\text{2}}^{\text{+}}$ and HPO⁺, respectively. In addition, the species PH₃ is mainly formed on dust grains through subsequent hydrogenation reactions of P, PH, and PH₂ and then released to the gas-phase via photodesorption. Finally, PN is formed at late times (10⁵ −10⁷ yr) mainly through the reaction N + CP → PN + C.

We also examined how the visual extinction A_V, the cosmic-ray ionization rate ζ(CR), and the surface mobility on dust grains affect the chemistry of P-bearing species. We found that all P-bearing species are strongly sensitive to the visual extinction: low A_V values of 1 and 2 mag lead to very low P-bearing molecular abundances of ~ 10⁻¹⁴−10⁻¹², indicating that a translucent region rather than a diffuse one is needed to produce observable amounts of P-containing species. All examined species in our study are influenced by the cosmic-ray ionization rate as well. An increasing ζ(CR) enhances the abundance of He⁺, H⁺ and C⁺, which in turn are effectively destroying all P-bearing species. A similar conclusion was found when changing the diffusion-to-desorption ratio to E_b ∕E_D = 0.77 and deactivating the possibility of quantum tunneling of light species on grain surfaces. This setup increases the H⁺ abundance, which in turn efficiently reacts with and destroys PN, PO, HCP, CP, and PH₃. Finally, we performed a study of the P-depletion level by tracing the phosphorus chemistry from a diffuse to a dense cloud with the application of a dynamical model that varies the density, the gas and dust temperature, the cosmic-ray ionization rate, and the visual extinction with time (see Appendix A). We came to the main conclusion that at high densities of ~ 10⁵ cm⁻³, atomic P is strongly depleted through freeze-out on dust grains, resulting in a significant increase of the gPH₃ abundance. The molecules PN, PO, HCP, CP, and PH₃ are also affected by freeze-out on grains and are destroyed by their reaction with ${\text{H}}_{\text{3}}^{\text{+}}$ when reaching the dense phase at timescales of ~10⁶−10⁷ yr.

Based on the predictions of our improved diffuse-cloud model, the (1–0) transitions of HCP, CP, PN, and PO are expected to be detectable with estimated intensities ranging from 10 to 200 mK. A possible detection of the above species will help us to further constrain the physical and chemical properties of our model and help us to better understand interstellar phosphorus chemistry.

Fig. 9 Chemical evolution of P-bearing molecules as a function of time for a diffusion-to-desorption ratio Eb ∕ED of 0.3 (with quantum tunneling) shown in the left panel and for a Eb∕ED of 0.77 (without quantum tunneling) in the right panel.

Acknowledgements

We thank the anonymous referee for his/her comments that significantly improved the present manuscript. The authors also wish to thank the IRAM Granada staff for their help during the observations. V.M.R. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 66493. Work by A.V. is supported by the Latvian Science Council via the project lzp-2018/1-0170. J.C. acknowledges Dr. J. C. Laas for his support with the Python programming.

Appendix A The depletion of phosphorus

The main advantage of studying the early phases of star formation is to avoid high levels of elemental depletion and thus to constrain the initial abundances used in our model to their cosmic values. This is crucial especially for phosphorus, as the small number of detections of P-bearing species in the ISM makes the determination of the P-depletion level quite difficult. In order to obtain an approximate estimation of the expected depletion level, we apply a dynamical model with time-dependent physical conditions that allows us to follow the chemical evolution of P-bearing species from a diffuse to a dense cloud. In particular we simulate a “cold” stage in which a free-fall collapse takes place within 10⁶ yr (Vasyunin & Herbst 2013; Garrod & Herbst 2006). During that time the density increases from n(H) = 300 to 10⁵ cm⁻³ and the visual extinction rises from 1 to 40 mag. The gas temperature decreases from 40 to 10 K, while the dust temperature drops slightly, from 20 to 10 K. Finally, the cosmic-ray ionization rate also changes from 1.7 × 10⁻¹⁶ s⁻¹ to 1.3 × 10⁻¹⁷ s⁻¹. We note here that the changes in the above-mentioned physical constants happen within 10⁶ yr, while the total chemical evolution is over 10⁷ yr. This means that between 10⁶ and 10⁷ yr the model becomes static with the above parameters retaining the values they reached at 10⁶ yr. That way, we simulate a long-lived collapse that provides enough time for chemical processes such as depletion to evolve.

As a first step, we plot the chemical evolution of the sum of abundances of gas-phase and solid-phase P-bearing species separately (see lower left panel of Fig. A.1). It is clearly visible how at late times, the gas-phase species decrease, and in return the grain species increase in abundance due to depletion. In particular, the sum of the gas-phase abundances of P-bearing species reduces by a factor of ~ 3000 at t = 10⁷ yr. This does not correspond to the elemental depletion, but it indicates the redistribution of phosphorus between the gas phase and the dust grains. The right-hand panel of Fig. A.1 shows the time-dependent abundances of the main carriers of phosphorus in the gas phase and on grains. The species that experience the largest change during the transition from diffuse to dense cloud are P⁺ and gPH₃. The P⁺ abundance strongly decreases down to ~10⁻¹⁶, mainly through its destruction reactions with OH, CH₄, S, and H₂. Atomic P decreases significantly because of freeze-out on dust grains, which is also evident through the increase in gP. According to the model almost all P that freezes out, quickly reacts with hydrogen on grains, and finally forms gPH₃ (after successive hydrogenation), which reaches a high abundance of ~2.5 × 10⁻⁷ at the end of our simulations.

Finally, Fig. A.2 shows the time-dependent abundances of PN, PO, HCP, CP, PH₃ in the gas phase (left panel) and the corresponding grain species (right panel). All species reach their peak abundances at around 10⁶ yr, followed by a strong decrease due to freeze-out on dust grains as well as through their reaction with ${\text{H}}_{\text{3}}^{\text{+}}$ (at t = 10⁶−10⁷ yr). The species PN, PH₃, and HCP show a more significant freeze-out than CP and PO, as they are the most abundant molecules in the gas-phase at t = 10⁶ yr. The freeze-out process is also clearly evident from the substantial increase of the corresponding grain species once high densities of ~ 10⁴−10⁵ cm⁻³ are reached (see right panel of Fig. A.2.)

Fig. A.1 Results of our dynamical model that simulates the transition from a diffuse to a dense cloud. Left panel: sum of abundances of all P-bearing species in the gas phase (red line) and the solid phase (blue line) as a function of time. Right panel: chemical evolution of the main carriers of phosphorus in the gas and solid phase: P+, P, gP and gPH3. In both figures the density profile of the free-fall collapse is depicted as a black dashed line.

Fig. A.2 Chemical evolution of PN, PO, HCP, CP and PH3 (left panel) and the corresponding grain species (right panel) as a function of time based on our dynamical model (diffuse to dense cloud). The black dashed line illustrates the density profile of the free-fall collapse. The gPH3 abundance is shown in Fig. A.1.

References

All Tables

All Figures

	Fig. 1 Spectra of the nondetected (2–1) transitions of PO, PN, HCP, and CP. The upper x-axis shows the rest frequency (in MHz) and the lower one is a velocity axis (in km s−1). The red dashed line indicates the 3σ level and the blue dashed line shows the transition frequency of the corresponding molecule. In the case of PO, we show as an example one of the observed transitions at 108.998 GHz.
In the text

	Fig. 2 Spectra of the detected species HNC, CN, C34S, and 13 CO in the 3 mmrange towards the line of sight to the extragalactitc source B0355+508. The red line represents the CLASS Gaussian fit.
In the text

	Fig. 3 Detected hyperfine components of the CN (1–0) transition between 113.12 and 113.20 GHz. The three strongest hyperfine components were detected in the three clouds with vLSR =−8, − 10, − 17 km s−1 except for the one weak transition (N, F) = (1, 1∕2) - (0, 1∕2), which was identified only in the two densest clouds (at − 10, − 17 km s−1).
In the text

Fig. 4 Results of the grid of models applying typical conditions for diffuse or translucent clouds in order to reproduce the observations towards the cloud at vLSR =−17 km s−1. The deviation between observations and model at the time of best agreement tbest is given by D(tbest, r), which is plotted versus density, temperature, and visual extinction. The best-fit model is given at a time tbest = 6.2 × 106 yr and has the following parameters: (n(H), AV, Tgas) = (300 cm−3, 3 mag, 40 K).

In the text

Fig. 5 Chemical evolution of the abundances of CO, CN, CS and HNC over 107 yr predicted by our best-fit model with the parameters (n(H), AV, Tgas) = (300 cm−3, 3 mag, 40 K). The colored horizontal bands correspond to the observed abundances towards the cloud with vLSR =−17 km s−1, including the inferred uncertainties. The vertical dashed line indicates the time of best agreement (t = 6.2 × 106 yr) between observations and model results.

In the text

	Fig. 6 Variation of the predicted abundances of PN, PO, HCP, CP and PH3 over 107 yr in our best-fit model. The dashed lines represent the 3σ upper limits derived from the observations at vLSR = −17 km s−1. In the case of PO we use 5 × 10−10 as an upper limit (see Table 9 and text for explanation).
In the text

	Fig. 7 Predicted abundances of P-bearing molecules as a function of visual extinction AV. The molecular abundances shown here are computed at t = 107 yr. Right panel: predicted abundances of ${\text{PCH}}_{\text{2}}^{\text{+}}$, C+, P+, He+, H+, and HPO+ as they are contributing the most to the formation and destruction of HCP, CP, PN, PO, and PH3 (left panel).
In the text

	Fig. 8 Chemical evolution of P-bearing molecules as a function of time under the effects of cosmic-ray ionization rates of ζ(CR) = 1.7 × 10−16 s−1 (left panel) and 10.6 × 10−16 s−1 (right panel).
In the text

	Fig. 9 Chemical evolution of P-bearing molecules as a function of time for a diffusion-to-desorption ratio Eb ∕ED of 0.3 (with quantum tunneling) shown in the left panel and for a Eb∕ED of 0.77 (without quantum tunneling) in the right panel.
In the text

Fig. A.1 Results of our dynamical model that simulates the transition from a diffuse to a dense cloud. Left panel: sum of abundances of all P-bearing species in the gas phase (red line) and the solid phase (blue line) as a function of time. Right panel: chemical evolution of the main carriers of phosphorus in the gas and solid phase: P+, P, gP and gPH3. In both figures the density profile of the free-fall collapse is depicted as a black dashed line.

In the text

	Fig. A.2 Chemical evolution of PN, PO, HCP, CP and PH3 (left panel) and the corresponding grain species (right panel) as a function of time based on our dynamical model (diffuse to dense cloud). The black dashed line illustrates the density profile of the free-fall collapse. The gPH3 abundance is shown in Fig. A.1.
In the text