New Detections of Phosphorus Molecules towards Solar-type Protostars

Phosphorus is a necessary element for life on Earth, but at present we have limited constraints on its chemistry in star- and planet-forming regions: to date, phosphorus carriers have only been detected towards a few low-mass protostars. Motivated by an apparent association between phosphorus molecule emission and outflow shocking, we used the IRAM 30m telescope to target PN and PO lines towards seven Solar-type protostars with well-characterized outflows, and firmly detected phosphorus molecules in three new sources. This sample, combined with archival observations of three additional sources, enables the first exploration of the demographics of phosphorus chemistry in low-mass protostars. The sources with PN detections show evidence for strong outflow shocks based on their H$_2$O 1$_{10}$-1$_{01}$ fluxes. On the other hand, no protostellar properties or bulk outflow mechanical properties are found to correlate with the detection of PN. This implies that gas-phase phosphorus is specifically linked to shocked gas within the outflows. Still, the PN and PO line kinematics suggest an emission origin in post-shocked gas rather than directly shocked material. Despite sampling a wide range of protostellar properties and outflow characteristics, we find a fairly narrow range of source-averaged PO/PN ratios (0.6-2.2) and volatile P abundances as traced by (PN+PO)/CH$_3$OH ($\sim$1-3%). Spatially resolved observations are needed to further constrain the emission origins and environmental drivers of the phosphorus chemistry in these sources.


INTRODUCTION
Along with H, C, O, N, and S, phosphorus is considered to be a key biogenic element and is found in all known forms of life on Earth. It is a critical component in DNA and RNA, ATP, and a host of other biological molecules (Macia 2005). While the biological importance of this element is well established, there remain many questions about how phosphorus was distributed during the formation of our Solar System. Given its importance for terrestrial life, phosphorus may also play a central role in the development of life on other planets. Ultimately, understanding the speciation and distribution of phosphorus in other nascent planetary systems can help us to contextualize Solar System constraints on the phosphorus chemistry, and to explore prospects for prebiotic chemistry in systems other than our own. * NASA Sagan Fellow The first phosphorus molecule detected in a star forming region was phosphorus nitride (PN). Turner & Bally (1987) presented PN detections in Ori (KL), W51M, and Sgr B2-three high-mass star-forming regions (SFRs)-and proposed that the higher-thanpredicted gas-phase abundances were the result of grain disruption. While phosphorus molecules were subsequently detected around evolved stars (e.g. Guelin et al. 1990;Agúndez et al. 2007;Tenenbaum et al. 2007), it was not until 2011 when Yamaguchi et al. presented the first detection of PN in a low-mass (i.e. solar-type) SFR, L1157. This was soon followed by the first detection of PO in a low-mass SFR, towards the same source L1157 (Lefloch et al. 2016). In recent years, PN and PO have been detected in a host of other massive SFRs and molecular clouds (e.g. Rivilla et al. 2016;Mininni et al. 2018;Rivilla et al. 2018;Rivilla et al. 2020;Bernal et al. 2021). In low-mass SFRs, PN detections were reported towards three sources as part of the TIMASS and ASAI surveys (Caux et al. 2011;Lefloch et al. 2018), and PN and PO were both detected towards the B1-a protostar (Bergner et al. 2019). Characterizing the phosphorus molecule emission in low-mass star forming regions is of particular importance because these systems are analogs to the proto-Solar Nebula.
While the chemistry of PO and PN in SFRs remains poorly understood, several common themes have emerged. In all low-and high-mass SFRs where PO and PN have both been detected, the source-averaged PO/PN ratio is within the range of 1-3 (Lefloch et al. 2016;Rivilla et al. 2016Rivilla et al. , 2018Bergner et al. 2019;Bernal et al. 2021). Correlations between the line profiles and column densities of PN and PO with the shock tracer SiO support that gas-phase P molecules are typically a product of grain sputtering in shocked regions (Lefloch et al. 2016;Mininni et al. 2018;Rivilla et al. 2018;Fontani et al. 2019;Bergner et al. 2019;Bernal et al. 2021). While phosphorus chemical models predict that PO and PN both form in the gas phase following the desorption of PH 3 from grains, no other P molecules including PH 3 have been detected in SFRs with PO and PN detections (Yamaguchi et al. 2011;Lefloch et al. 2016;Jiménez-Serra et al. 2018).
In the Solar System, the phosphorus abundance in primitive meteorites is close to Solar, indicating a large reservoir of refractory phosphorus in the Solar Nebula (most commonly in the form of phosphate minerals; Lodders 2003). Additionally, the volatile phosphorus carrier PO was recently detected towards comet 67P/Churyumov-Gerasimenko (hereafter 67P) by the Rosetta mission (Altwegg et al. 2016;Rivilla et al. 2020) with a high elemental P/O ratio ∼0.1-0.5× the Solar value. This is especially noteworthy given that comets may have played a role in the delivery of prebiotic material to young terrestrial planets (e.g. Rubin et al. 2019a). The PO/PN ratio in comet 67P was found to be at least 10. Whether the Solar System is typical in its volatile P abundance and composition remains an open question.
The small sample size of low-mass SFRs with P molecule detections is a major obstacle to constraining the phosphorus astrochemistry relevant to planet formation. Informed by growing evidence that outflow shocks play a key role in releasing phosphorus in the gas in low-mass star-forming regions (Yamaguchi et al. 2011;Lefloch et al. 2016;Bergner et al. 2019Bergner et al. , 2022, we targeted PO and PN lines towards seven low-mass protostars with well-characterized outflows using the IRAM 30m telescope. In this paper, we present detections of phosphorus carriers towards three new low-mass protostars (NGC1333-IRAS 3, Ser-SMM1, Ser-SMM4), with one additional tentative detection (L723). This sam-ple is combined with archival observations of B1-b and NGC1333-IRAS 4a from the ASAI program (Lefloch et al. 2018), for which PN detections were previously reported but not characterized. Section 2 describes our observations, including the source and line targets. In Section 3 we outline our analysis methods, and Section 4 presents the resulting PO and PN column densities and column density ratios with respect to CH 3 OH and H 2 . With our expanded sample of phosphorus molecules in Solar-type star forming regions, in Section 5 we discuss the trends identified across our sample, and comparisons between the P inventories in low-mass protostars and primitive Solar System materials.

Source targets
Our source sample (Table 1) is selected from the embedded low-mass protostars observed by the Herschel WISH program. For these sources, the outflow properties have been previously characterized based on the H 2 O ground-state 1 10 -1 01 transition and auxiliary CO observations (Kristensen et al. 2012;Yıldız et al. 2015). Of the sources observable by the IRAM 30m telescope, our targets were chosen to span a range of physical properties (bolometric luminosities from ∼2-40 L , envelope masses from ∼1-50 M ) as well as evolutionary stages (bolometric temperatures from ∼20-150 K).
We note that the Herschel H 2 O 1 10 -1 01 observations correspond to a beam size of 39 , which is larger than the IRAM 30m beam FWHM (16-27 at 150-90 GHz). Higher-frequency water lines (corresponding to smaller beam sizes of ∼13-28 ) were also detected by Herschel towards the WISH sources (Mottram et al. 2014), confirming that shocks are present on the smaller spatial scales accessed by the IRAM 30m beam.  Mottram et al. (2014). Properties for B1-b and B1-a are taken from Pezzuto et al. (2012) and Hatchell et al. (2007). Because B1-b is a binary system, we list the average T bol and total L bol of the two components. Sources with new coverage of phosphorus molecule lines are listed above the horizontal line, and those with archival coverage are below the line.

Description of observations
Observations of our source targets were taken with the IRAM 30m telescope using the EMIR 90 GHz (3mm) and 150 GHz (2mm) receivers with the Fourier Transform Spectrometer backend. The 3mm setting was observed with a spectral resolution of 200 kHz, and the 2mm setting was observed with a resolution of 50 kHz. Observations were taken between 2020 March 03 and 2020 March 06. Additional details on the observing parameters and quality can be found in Appendix A. Initial data reduction steps were performed using CLASS 1 , assuming a main-beam efficiency described by the Ruze equation with a surface RMS of 66µm and B eff,0 =0.863, as appropriate for the IRAM 30m telescope. Spectral line fitting and subsequent analysis were performed with Python.

Line Targets
Every source was observed in the 2mm setup, the 3mm setup, or both the 2mm and 3mm setups, as a result of variable weather conditions (Table 6). Each setup covers one PN line and four PO hyperfine components. In addition, C 17 O 1-0 was covered in the 3mm setup, and a range of CH 3 OH and SO lines were covered in both the 2mm and 3mm setups. The line parameters used in our analysis are summarized in Table 2. Additional 1 https://www.iram.fr/IRAMFR/GILDAS/ details on the spectral setups can be found in Appendix A.

Additional sources
In addition to the new observations described above, we also made use of archival IRAM 30m observations of NGC1333-IRAS 4a and B1-b from the ASAI program (Lefloch et al. 2018). PN was previously reported as detected towards these sources, however there has not yet been an analysis of the PN (or PO) column densities or line kinematics. ASAI provided coverage of the 3mm and 2mm PN and PO lines targeted by our survey, as well as higher-frequency (1.3mm) lines. The observational details can be found in Lefloch et al. (2018).
We also include the previous constraints on PN and PO emission in B1-a from Bergner et al. (2019) when comparing across the source sample (Section 4.3). Note that B1-a and B1-b were not part of the WISH program.

Line Fitting and Detections
For each line target, we subtracted a local baseline and fit a Gaussian to the spectrum. In the case of PO, we stacked the four hyperfine components of each rotational level to increase the SNR. This was done by producing spectra centered on each hyperfine line and rebinned to a common velocity grid. Each hyperfine spectrum was scaled based on its intrinsic line strength relative to the strongest hyperfine component (marked with * in Table 2), then all spectra were averaged. The PN and Line parameters are from CDMS (Müller et al. 2001(Müller et al. , 2005 based on data from Cazzoli et al. (2006), Bailleux et al. (2002), Kawaguchi et al. (1983), Xu et al. (2008), Clark & De Lucia (1976), and Klapper et al. (2003). + indicates lines that are covered only for the ASAI sources. Note that we include only E-CH3OH lines in our analysis. PO hyperfine lines of a given J-level were stacked after normalizing to those marked with a *. f indicates an f -type PO transition, while all other PO lines are e.
PO Gaussian fit results, including velocity-integrated intensities, are listed in Table 7, and those for CH 3 OH are in Table 8. Uncertainties on the integrated intensities consist of a 10% calibration uncertainty added in quadrature with the Gaussian fit uncertainty, and are propagated through subsequent analysis steps. We consider a line to be firmly or tentatively detected if the velocity-integrated intensity is >3σ or 2-3σ, respectively. Figure 1 shows the observed 3mm and 2mm PO and PN spectra along with the Gaussian fits for sources with firm or tentative detections. Note that several of the 2mm PN lines appear to be composed of multiple velocity components, however due to the low SNR we fit only a single Gaussian to each feature. Spectra for the sources with no PO or PN detections, as well as 1.3mm spectra for the ASAI sources, are shown in Appendix B.
We firmly detected at least one PN line towards five sources: NGC1333-IRAS 3, Ser-SMM4, Ser-SMM1, NGC1333-IRAS 4a, and B1-b. PN is also tentatively detected towards L723 (2.7σ), though follow-up observations are needed to confirm if this feature is real. PO is firmly detected at 2mm in Ser-SMM4, and tentatively in Ser-SMM1 and NGC1333-IRAS 3. These tentative PO lines exhibit similar velocity offsets and line-widths as narrow features that are seen in the PN 2mm line profiles, and we therefore expect that they are not spurious. This can be more clearly seen in Figure 4, in which the PN and PO spectra are overlaid. The 1.3mm PO line is also tentatively detected in B1-b ( Figure 5).
Note that for NGC1333-IRAS 3 and L723, PN is only detected in one of the two observed setups. While a nondetection of the PN 2mm feature in L723 is readily explained by a low excitation temperature, it is surprising to detect the 2mm feature but not the 3mm feature in NGC1333-IRAS 3. However, there are no other plausible line carriers nearby the 2mm PN line frequency, and we have ruled out that ghost lines could be contributing. Moreover, the narrow spike in the PN 2mm line near 2.5 km s −1 overlaps precisely with the PO 2mm feature, as is also seen in Ser-SMM1. We therefore treat the 2mm PN line as a detection, and discuss the excitation conditions implied by the non-detected 3mm PN line in Section 3.3. The PO 3mm lines in Ser-SMM1 and NGC1333-IRAS 3 are also not detected despite tentative detections in the 2mm setup. This is likely due to the very narrow PO linewidths and low SNRs, combined with the coarse (200 kHz) spectral resolution of the 3mm setting. Similarly, PO is tentatively detected at 1mm in B1-b, but not detected in the 2mm or 3mm settings. Again, this could be due to a non-LTE excitation and/or high noise levels in the lower-frequency settings.  (Table 1).

Column Density Calculations
To calculate the column densities of our detected molecules, we followed the population diagram method laid out by Goldsmith & Langer (1999). This method assumes that the observed transitions are optically thin and well described by local thermodynamic equilibrium (LTE). We discuss the effects of non-LTE conditions in Section 3.3.
For a given transition, we find the upper-state population (N u ) divided by the upper-state degeneracy (g u ) using the velocity-integrated main-beam temperature T mb dV : where ν is the line frequency and A ul is the transition Einstein coefficient. To go from the upper-level populations to the total column density (N T ), we then use: where E u is the upper state energy (K), T r is the rotational temperature, and Q(T r ) is the partition function.
If the emitting region does not fill the beam, beam dilution (η bf ) must be accounted for. We calculate the beam dilution factor by: where θ is the beam FWHM in arcsec and R src is the source emitting radius. The IRAM 30m beam FWHM is 27 at 90 GHz and 16 at 150 GHz. We divide each upper-level population by the appropriate beam dilution factor before calculating the column density. When multiple transitions of a given molecule are detected, we fit Equation 2 to the observed upper-level populations to derive the total column density and the rotational temperature. When only a single transition is available, we assume a rotational temperature to calculate the column density.

PN and PO
In the case of Ser-SMM1, B1-b, and NGC1333-IRAS 4a, we detected two PN lines and could fit for both the rotational temperature and column density. For PN and PO in all other sources, we detected only a single line, and solved for the column density assuming a rotational temperature of 10 K based on the values found in Ser-SMM1, B1-b, and NGC1333-IRAS 4a, and previously in the low-mass sources B1-a and L1157-B1 (Bergner et al. 2019;Lefloch et al. 2016). Because the source size is unknown, and cannot be determined from only one or two transitions, we adopt an emission radius corresponding to a linear scale of 1000 au for each source. This is based on imaging of PN and PO in the low-mass protostar B1-a, which shows the phosphorus molecules emit from distinct clumps with radii of this scale (Bergner et al. 2022). Note that the single-dish line fluxes for B1-a were fully recovered in the interferometric observations, ruling out a more diffuse emission component. The resulting column densities are shown in Table 4.
When we detected no line of PN or PO in a source, we solved for the column density upper limit using the 3σ upper limit on the integrated intensity, where σ = RMS × FWHM / √ n ch , and n ch is the number of channels spanned by the FWHM. We found the RMS in a 30 km s −1 spectral window around the line position, and assumed the average FWHM of CH 3 OH lines in the same source.
If the emitting radius is in fact larger than our assumed value, then our reported column densities are over-estimates. While this introduces a large uncertainty in our derived column densities, it is not currently possible to make a better estimate since the origin of phosphorus molecules within the protostellar environment remains poorly understood. Additionally, the beam dilution correction (Equation 3) assumes a Gaussian emission distribution, although the distributions may be more complicated in reality. Ultimately, spatially resolved observations of PN and PO are needed to obtain more robust column densities.

CH3OH
Our observations also covered multiple CH 3 OH lines in each of the 2mm and 3mm setups. To derive column densities, we fitted rotational diagrams as described in Section 3.2. We used exclusively transitions of the E isomer of CH 3 OH in calculating column densities, since our observations did not cover enough bright A-CH 3 OH transitions to independently fit rotational diagrams for this isomer. We then solved for the total CH 3 OH column density assuming a CH 3 OH E/A ratio of 0.5 based on previous surveys of CH 3 OH towards protostellar outflow sources (Holdship et al. 2019).
Given the larger number of CH 3 OH transitions detected compared to the phosphorus molecules, our fitting approach differed in several ways. First, we attempted to fit for the source size in addition to the column density and rotational temperature, but found there was no preference for a small source size. Given this, we assumed that the CH 3 OH emission fills the beam for subsequent fitting. Note that JCMT observations of similar sources (including the Serpens sources in our sample) reveal large-scale CH 3 OH emission (tens of au) surrounding each protostellar core (Kristensen et al. 2010). Therefore, it is reasonable that CH 3 OH emission would fill the IRAM 30m beam for our observations (16-27"). Second, we found that some sources required both a hot and cold component to fit the observed upper level populations, as described in detail in Appendix C. Table  4 shows the column densities and excitation temperatures of both components when relevant. Finally, given the larger number of data points, we performed the rotational diagram fitting using the Markov Chain Monte Carlo code emcee (Foreman-Mackey et al. 2013) as this provides more robust estimates of the fit uncertainties. CH 3 OH rotational diagrams are shown in Appendix C.

C 17 O and H2
For sources observed in the 3mm setup, our observations covered the C 17 O 1-0 hyperfine complex at 112 GHz. To constrain the C 17 O column density in these sources, we fit the C 17 O hyperfine spectrum following the procedure described in Bergner et al. (2019). The C 17 O spectra are shown in Appendix C. We assume the C 17 O emission fills the beam since CO should be present throughout the protostellar envelope. The derived column densities are therefore beam-averaged. The C 17 O column density was then used to estimate the H 2 column density by assuming a CO/C 17 O ratio of 2005 (Wilson 1999), and a H 2 /CO ratio of 10 4 . The resulting H 2 column densities are listed in Table 4.

Non-LTE Effects
As noted in Section 3.1, the non-detection of the 3mm PN line towards NGC1333-IRAS 3 is not consistent with LTE emission: for any positive excitation temperature, the 3mm line should be detected above the observed noise levels. To explore whether non-LTE PN emission could explain this, we ran a grid of RADEX models (van der Tak et al. 2007), using the PN collisional rates from the LAMDA database (Schöier et al. 2005) with data taken from Tobo la et al. (2007). Figure 2 shows the excitation landscape for the PN 2-1 (3mm) line, assuming a line width of 1 km s −1 and a column density of 10 12 cm −2 . For low gas temperatures and high gas densities, the line is thermalized. At high temperatures and high densities, the line emits superthermally, while at low densities and all temperatures it emits subthermally. Importantly, for gas densities around 10 6 cm −3 and temperatures >40 K, the 2-1 line is characterized by a negative excitation temperature, i.e. a population inversion. The PN 3-2 line shows a similar overall excitation map as the 2-1 line, but the population inversion region is shifted to higher temperatures (>70 K).
Given that phosphorus molecule emission is associated with outflow shocks, it is certainly plausible that the emitting region of NGC1333-IRAS 3 would be warm and dense. Thus, we expect that the non-detection of the PN 3mm line can be explained by a level population inversion. That this behavior is not seen for other sources in our survey, or for previously observed low-mass protostars (Lefloch et al. 2016;Bergner et al. 2019) may reflect that the emitting gas is warmer and/or denser in NGC1333-IRAS 3. Indeed, with a non-LTE analysis of four PN lines, Lefloch et al. (2016) infer a gas density of ∼10 5 cm −2 and a kinetic temperature around 40-80 K, i.e. the 'subthermal' regime of Figure 2. The Figure 2. Difference between the RADEX-calculated excitation temperature of the PN 2-1 line and the gas kinetic temperature, shown for a range of gas densities and temperatures. Different excitation regimes are indicated: the line is thermalized in the white region, superthermally excited in the red region, and subthermally excited in the blue region. Interior to the white contour, Tex is <0, reflecting a population inversion (nu/gu > n l /g l ).
low (but positive) PN excitation temperatures found in Ser-SMM1, B1-b, and NGC1333-IRAS 4a also point to emission either within the subthermal or thermalized regimes, corresponding to lower densities or lower temperatures compared to the population inversion regime.
The complex excitation behavior of PN illustrated in Figure 2 underscores the need for multi-line coverage to obtain good constraints on the column densities and emitting conditions. Indeed, with only two lines it is not possible to perform a full non-LTE fit with our current observations, and coverage of higher-J PN lines is needed to fully understand the excitation of the emitting molecules. Still, we note that for the cool gas temperatures (<40 K) tested in Bergner et al. (2019), the derived PN column densities are not significantly affected by non-LTE effects.

Column Densities
The column density findings for each source are presented in Table 4. For sources with detections, the column densities of both PN and PO range are around 10 12 cm −2 . Non-detection upper limits are on the order <10 11 -10 12 cm −2 . The phosphorus molecule column densities in our sample are about one order of magnitude lower than was reported in Ori (KL), W51M, and Sgr B2 (Turner & Bally 1987), several high mass protostars. The PO and PN column densities in our sample are in good agreement with the low-mass protostars L1157 B1 (Yamaguchi et al. 2011) and B1-a (Bergner et al. 2019).
Because the degree of beam dilution remains uncertain, it is important to note that our PN and PO column densities depend on the actual source size (R src ) of the emitting region. If R src is smaller than estimated (1000 au), the true column densities will be higher than listed, and the opposite for a larger R src . If T rot is lower (e.g. 5 K rather than 10 K), the PN and PO column densities could be underestimated by a factor of 2-3. A higher T rot (20 K) does not strongly affect the PN or PO column density.
CH 3 OH column densities span ∼10 13 -10 14 cm −2 in our sample, and fall within the range found in a survey of embedded low-mass protostars by Graninger et al. (2016). H 2 column densities in our sample are ∼10 22 cm −2 , which is an order of magnitude higher than the H 2 column density found in B1-a (Bergner et al. 2019).

Column density ratios
In Table 5, we present the column density ratios of PO and PN relative to other molecules. The PO/PN ratio is of particular interest for comparing the phosphorus chemistry in our source sample with other star-forming regions. Intriguingly, to date all source-averaged measurements of PO/PN in dense interstellar regions are between 1 and 3, despite probing a wide range of physical environments (Lefloch et al. 2016;Rivilla et al. 2016Rivilla et al. , 2018Bergner et al. 2019;Rivilla et al. 2020;Bernal et al. 2021). Across our sample, we find PO/PN ratios within this range or somewhat lower (0.6-2.2). Note that we exclude NGC1333-IRAS 3 from this comparison, since it likely suffers from non-LTE effects for at least the PN emission (Section 3.3). Coverage of additional PN and PO lines is needed to perform a full non-LTE fit and confirm the PO/PN ratio in NGC1333-IRAS 3.
In some sources (e.g. Ser-SMM1, NGC1333-IRAS 3), the kinematics of the PN and PO lines are different, suggesting that they may emit from different gas. This is consistent with previous spatially resolved observations of phosphorus molecules: in both AFGL 5142 and B1-a (Rivilla et al. 2020;Bergner et al. 2022), PN and PO emit from several distinct emission clumps which have different velocity offsets and also different PO/PN ratios. This is likely true for these sources as well. We emphasize that the derived PO/PN ratios are averages across the source, and that there is likely small-scale variation in the PO/PN ratio across these sources.
We also present PN+PO H abundances as an indicator of the total phosphorus abundance in the gas. For sources with H 2 column density estimates, we find PN+PO H abundances of the order 10 −11 . This is several orders of magnitude below the Solar P abundance of 3×10 −7 (Asplund et al. 2009), and on the low end but consistent with the range of P/H ratios previously inferred in the dense ISM (Lefloch et al. 2016;Rivilla et al. 2016;Bergner et al. 2019;Bernal et al. 2021). Note that we use beamaveraged H 2 column densities to estimate this ratio. Spatially resolved observations are needed to constrain the H 2 column density co-spatial with the phosphorus molecules.
Lastly, we list the PN+PO CH3OH ratio to enable a comparison with cometary ice compositions (Section 5.3). Because CH 3 OH does not form efficiently in the gas, the gasphase CH 3 OH content in these sources should originate from ice desorption. Assuming that phosphorus is released into the gas from the same sputtering process, the observed gas-phase P CH3OH ratio should reflect the underlying ice-phase P CH3OH ratio. For sources with both a hot and cold CH 3 OH component (Table 4), we use only the cold component to find this ratio given that the excitation conditions are closer to the phosphorus molecules. We find PN+PO CH3OH ratios of ∼0.7-2.7% in our sample, with a low upper limit (<0.3%) found for NGC1333-IRAS 4a.

Trends across the source sample
With five firm detections of phosphorus molecules towards low-mass star forming regions, we can now search for relationships between the source properties and the gas-phase phosphorus inventory. Our source sample was selected to span a range of physical properties (e.g. envelope mass) and evolutionary stages (traced by bolometric temperature). Additionally, most sources in our sample were part of the Herschel WISH program and have auxiliary constraints on the outflow characteristics, traced by H 2 O and CO emission (Kristensen et al. 2012;Yıldız et al. 2015). In general, we expect the H 2 O emission to correlate with shocked gas, whereas CO emission traces the bulk mechanical properties of the outflows. Figure 3 shows the the PN column density as a function of the protostar and outflow properties. We focus on PN since it is firmly detected towards more sources than PO. Comparisons between the PN column density and the protostar properties (top row) reveal no trends. PN is detected towards both Class 0 (T bol < 70K) and Class I (T bol >70K) sources, and towards sources with envelope masses ranging from 3-52 M and bolometric luminosities from 1-100 L . Thus, the release of phosphorus into the gas does not appear correlated with any intrinsic physical property of the protostar.
The middle row of Figure 3 shows comparisons between the PN column density and various features of the H 2 O 1 10 -1 01 emission: the line full-width at zero-intensity (∆v max ), velocity-integrated line intensity (T mb dv), and kinematic components (broad, medium, or both, as defined in Section 2.1), all taken from Kristensen et al. (2012). While we do not see any monotonic relationships, we do find similarities among the sources where PN is detected. The sources where PN is firmly detected have the highest velocity-integrated H 2 O intensities of the sample (>5 K km s −1 ). These same sources also have broad H 2 O linewidths (∆v max >50 km s −1 ). Together, this indicates that phosphorus is released into the gas when outflow shocking is relatively strong, producing broad H 2 O lines and a high H 2 O content in the gas. Figure 3 (bottom) shows comparisons with the outflow mass (M outf low ), outflow force (F CO ), and dynamical age (T dyn ) measured from CO 3-2 emission in Yıldız et al. (2015). Note that we adopt the average values between the red and blue outflow lobes. In all cases, we do not see any correlations: PN is both detected and not detected over a range of outflow masses, forces, and ages. This lack of correlation with CO-derived mechanical properties, taken together with the association between PN and H 2 O emission, implies that gas-phase phosphorus is specifically related to shocked gas, rather than the bulk outflowing material.
We also attempted to look for trends between the source or outflow properties and the PO/PN ratio. However, given the narrow range in PO/PN ratios across the sample (Table 5) Table 1), H2O 110-101 emission properties (middle row; Kristensen et al. 2012), and CO-derived outflow properties (bottom row; Yıldız et al. 2015), described in detail in the text. Purple represents a firm detection, pink represents a tentative (2-3σ) detection, and yellow represents an upper limit non-detection. B1-a is shown as a grey star. Note that the NGC1333-IRAS 3 column density, plotted as a diamond, is particularly uncertain due to non-LTE effects (Section 3.3).
In Section 4.3, we showed that the detection of phosphorus molecules appears to trend with the presence of shocked gas as traced by the H 2 O 1 10 -1 01 line. However, the PO and PN lines in our sample are significantly narrower than the H 2 O lines presented in Kristensen et al. (2012), a few km s −1 compared to tens of km s −1 , indicating that P molecules emit from a much narrower range of physical environments within the outflow compared to H 2 O. Moreover, the PO and PN lines are offset by just a few km s −1 from the source rest velocities, reflecting an origin in either weakly shocked or postshocked gas. To further explore the origin of phosphorus molecule emission within the protostellar environment, Figure 4 shows the line profiles of PN and PO with one another and with the weak-shock tracers CH 3 OH and SO.
In both Ser-SMM1 and NGC1333-IRAS 3, the PO line is a narrow feature red-shifted by ∼3 km s −1 relative to the source rest velocity. The PN lines in both sources exhibit a broader component red-shifted by ∼1-2 km s −1 , along with a narrow component with a similar width and offset to the PO feature. This may reflect that PO preferentially emits from higher-velocity material within the shock, while PN emits from multiple environments.
In Ser-SMM4, a similar pattern is seen, in which the PO line is narrower but overlaps with the PN feature. In this case, the PO line is much broader than those seen in Ser-SMM1 and NGC1333-IRAS 3, and the PO line is blue-shifted rather than red-shifted relative to the source velocity.
We also compare the PN line profiles to the weakshock tracers CH 3 OH and SO. In the low-mass protostar B1-a, the kinematics of CH 3 OH and the sulfur carrier SO 2 were found to closely resemble the phosphorus molecules (Bergner et al. 2019), and similarities between PN and SO abundances have been previously identified in numerous massive star-forming regions (Mininni et al. 2018;Rivilla et al. 2020). Among our sources, the CH 3 OH and SO lines generally peak close to the source rest velocity, with SO characterized by a narrower line width compared to CH 3 OH. Ser-SMM1, Ser-SMM4, and B1-b exhibit similarities in the line shapes and velocity offsets of PN compared to CH 3 OH and SO. This is particularly true when comparing SO and PN in Ser-SMM1 and B1-b. The PN line profiles do not agree well with CH 3 OH or SO in NGC1333-IRAS 3 and NGC1333-IRAS 4a. In both cases, there is a strong PN emission component red-shifted with respect to the CH 3 OH and SO lines. Thus, while there is some overlap, neither CH 3 OH nor SO appear to uniquely trace the same material as the phosphorus molecules in these sources.

DISCUSSION
Our survey, combined with archival observations (Lefloch et al. 2018), has enabled us to characterize the emission of phosphorus molecules towards five new Solar-type star forming regions. Here, we discuss the implications of our findings with respect to the interstellar phosphorus chemistry and the volatile phosphorus reservoir in the young Solar System.

Phosphorus molecules and outflow properties
Prior detections of PO and PN in the low-mass star forming regions B1-a (Bergner et al. 2019) and L1157-B1 (Yamaguchi et al. 2011;Lefloch et al. 2016) suggest that the phosphorus emission in those sources originates from outflow shocks. However, the exact relationship between the phosphorus chemistry and outflow physics currently remains ambiguous.
Of the sources targeted in this work, we firmly detect phosphorus molecule emission only towards sources with the highest H 2 O 1 10 -1 01 fluxes (Figure 3), reflecting robust shock-induced sputtering of icy material from the grains. We do not identify any correlations between phosphorus molecules and either the protostellar properties or the bulk mechanical properties of the outflowing gas. This implies a particular association between gas-phase phosphorus and shocked gas. This is consistent with the emission patterns seen in the low-mass protostar B1-a, in which PN and PO emit not from the SiO outflow, but from distinct shock positions associated with ice sublimation (Bergner et al. 2022).
The PO and PN lines are characterized by narrow linewidths (<2 km s −1 ) and small offsets from the source systemic velocities (a few km s −1 ). This indicates that the phosphorus molecules do not originate in strongly shocked material itself, but rather in a post-shocked stage. In this scenario, shocking is required to release the parent phosphorus carriers from the surface of grains, followed by gas-phase chemistry in the postshocked gas to form PO and PN. Indeed, Rivilla et al. (2020) and Bergner et al. (2022) infer a similar origin in post-shocked gas based on spatially resolved observations of PO and PN towards a high-mass and low-mass protostar, respectively. Notably, however, the PN line in L1157-B1 extends to velocities as high as 20 km s −1 (Lefloch et al. 2016), suggestive that the emission in this source may be more closely connected to the shocked material.
In our sources, we see interesting velocity structure in the phosphorus molecule line profiles (Figure 4). In Ser-SMM 4, NGC1333-IRAS 3, and NGC1333-IRAS 4a, the PN lines are fairly broad and appear to be composed of several velocity components. In Ser-SMM1 and NGC1333-IRAS 3, PO emits as a very narrow feature a few km s −1 offset from the main PN line, coincident with a very narrow PN feature. This velocity structure likely corresponds to emission from different shock positions along the outflow. Indeed, spatially resolved observations towards both AFGL 5142 and B1-a show that PO and PN emit from several distinct clumps adjacent to the outflow, which are characterized by different velocity offsets (Rivilla et al. 2020;Bergner et al. 2022). That PO generally exhibits fewer velocity components compared to PN in our sources suggests that it emits from a smaller number of clumps along the outflow. This is analogous to the emission pattern seen in B1a, where PN emits brightly from two emission clumps, while PO emits brightly from only one clump (Bergner et al. 2022). Resolved observations are needed to further explore the the progression of post-shock chemistry in these sources, and the origins of chemical differentiation at different shock positions.

Volatile P in low-mass protostars
In the two sources for which we could measure the H 2 column density, we estimate source-averaged PN+PO H ratios of 1.5-2.7 ×10 −11 (Table 5). This corresponds to <0.01% of the Solar P abundance (Asplund et al. 2009).  (Table 1).
For reference, the combined PO and PN abundances in L1157-B1 and B1-a are 3.6×10 −9 and 10 −10 -10 −9 , respectively. Thus, our measured abundances support previous findings that volatile P is only a trace component of the total P inventory in planetary system progenitors.
It is important to note that our PO and PN column densities are fairly uncertain due to unknown beam dilution factors. Still, even adopting a very small emitting area of 1 in radius would only increase the derived abundances by an order of magnitude, leading to an abundance that is still significantly lower than the Solar P/H ratio. The column density ratios of PO/PN in our sample range from 0.6 to 2.2 (Section 4.2). These values are generally consistent with the source-averaged PO/PN ratios reported in the other low-mass protostars B1-a and L1157-B1, which range from 1 to 3 (Lefloch et al. 2016;Bergner et al. 2019). Shock chemistry models predict that PO/PN should be sensitive to factors like the preshock duration, shock density, shock velocity, time postshock, and gas-phase O/N ratio (Aota & Aikawa 2012;Lefloch et al. 2016;Jiménez-Serra et al. 2018). However, we find a small range in source-averaged PO/PN ratios across our sample: a factor of 3.5, as compared to 1-2 orders of magnitude of variations in the source and outflow properties (Section 2.1). Thus, the sourceaveraged PO/PN ratios appear only weakly sensitive to environmental factors.
Spatially resolved observations will likely provide a clearer path to determining environmental drivers of the phosphorus chemistry: existing studies have revealed a wider range of local PO/PN ratios from ∼1-8 (Rivilla et al. 2020;Bergner et al. 2022). Follow-up mapping towards the sources with PN and PO detections would provide valuable new constraints on how and why the phosphorus chemistry varies within a protostellar environment.

Comparison to Comet 67P
Comets represent the most minimally processed relics from the planet formation epoch in our Solar system. Volatile P was recently detected in the Solar System comet 67P/Churyumov-Gerasimenko (Altwegg et al. 2016), with a reported PO/PN ratio of at least 10 (Rivilla et al. 2020). The source-averaged PO/PN ratios previously measured towards the proto-Solar analogs L1157 and B1-a are low (<3) compared to comet 67P (Lefloch et al. 2016;Bergner et al. 2019). With our expanded sample, we have shown that B1-a and L1157 are not unusual in their low PO/PN ratios, and that sourceaveraged PO/PN ratios in low-mass star forming regions appear in general to be lower than the cometary measurement. This difference between the protostellar and cometary PO/PN ratios may be explained if the solar system exhibits an unusual volatile P chemistry compared to typical low-mass protostars, or if processing between the protostellar and disk stages alters the volatile P composition in the ice.
The (PO+PN)/CH 3 OH ratio can be used to compare the volatile P content between low-mass protostars and comet 67P, assuming that the observed CH 3 OH, PO, and PN emission all trace material sputtered from the grains by shocks. In our sources, we find sourceaveraged (PO+PN)/CH 3 OH ratios of 0.7-2.7%, as well as a low upper limit (<0.3%) for NGC1333-IRAS 4a ( Table 5). The ratio in B1-a is typically 1-3% (Bergner et al. 2022), and comet 67P was measured to have a (PO+PN)/CH 3 OH ratio of ∼5% (Rubin et al. 2019b). Thus, most measurements of (PO+PN)/CH 3 OH along the low-mass star formation sequence are on the order of a few percent. This suggests that the volatile P content of the young Solar system is fairly typical. Still, there is currently only one comet for which P/CH 3 OH ratios have been measured, and better constraints on the volatile P inventory of primitive icy material in the Solar System would improve this comparison.
To date, it is unknown whether prebiotically accessible P on the early Earth was sourced from a soluble mineral like schreibersite (e.g. Pasek & Lauretta 2005), or from the delivery of volatile icy material through impact (e.g. Rubin et al. 2019a;Rivilla et al. 2020). The similarity in volatile P content between a Solar System comet and protosolar analogs implies that, if impact delivery was indeed an important source of P for origins of life chemistry on Earth, then planets in other systems could also be seeded with prebiotic P in the same way.

CONCLUSIONS
We have targeted PN and PO lines towards seven lowmass protostars with the IRAM 30m telescope in an effort to expand our understanding of P astrochemistry in proto-Solar analogs. We report new detections of phosphorus molecules towards NGC1333-IRAS 3, Ser-SMM4, Ser-SMM1. Combined with archival observa-tions of NGC1333-IRAS 4a and B1-b (Lefloch et al. 2018), we have characterized the PN and PO column densities and line kinematics towards five low-mass protostars, enabling explorations of the demographics of phosphorus chemistry in proto-Solar analogs.
The sources with PN detections show evidence for powerful outflow shocks based on their H 2 O 1 10 -1 01 emission, supporting the central role of shocks in releasing P into the gas in dense SFRs. Correlations are not found with the protostar properties or the bulk mechanical properties of the outflow, implying a specific association between gas-phase phosphorus and shocked gas. Based on line kinematics, the PO and PN emission seems to arise from quiescent post-shocked gas rather than directly from the strongly shocked gas. Ultimately, spatially resolved observations are required in order to unambiguously determine the origin of the P-bearing species in these protostellar regions.
While our source sample spans a wide range of protostellar and outflow characteristics, we find a fairly narrow range of source-averaged PO/PN ratios (0.6-2.2). This suggests that the phosphorus chemistry is not particularly sensitive to environmental factors. However, we expect the PO/PN ratios to vary locally across these sources, and spatially resolved observations are needed to further explore the drivers of this phosphorus chemistry.
With this analysis, we can begin to compare how the 'typical' protostellar phosphorus chemistry compares to the young Solar System. We find that the volatile P content in our sources is consistent with that inferred for comet 67P, with (PN+PO)/CH 3 OH ratios around 1-3%. However, the source-averaged PO/PN ratios in low-mass protostars are universally lower than the ratio >10 found in the comet. Spatially resolved observations of additional Solar-type star forming regions, and further constraints on volatile P chemistry within the Solar System, are needed to further constrain how this biocritical element is inherited to the planet-formation epoch.  Table 6 provides an overview of the observations taken at the IRAM 30m telescope, including for each source the spectral setups, time on-source, system temperature, and resulting spectral RMS. The listed image frequencies correspond the following spectral coverages: 110  GHz. All observations were performed in position switching mode, with the 'off' positions listed in Table 6. Pointing corrections were performed every 2h, and focus corrections every 4h. The primary line targets are listed in Table 2. B. PN AND PO LINE FITS AND ADDITIONAL SPECTRA Figure 5 shows the ASAI 1.3mm spectra of PN and PO, analogous to the 2mm and 3mm spectra shown in Figure  1. The PO line is tentatively detected towards B1-b, and the other lines are not detected.
The line fitting results for PN and PO towards all sources are given in Table 7. We did not detect phosphorous emission towards three of the seven sources with new observations (L1527, L1551, Ser-SMM3). For these sources, Note-Parameters derived from Gaussian fitting of each PO and PN line. 1σ uncertainties on the velocity-integrated intensities are listed in brackets. Tentative detections (2-3σ) are listed with an *. Source rest velocities have not been subtracted from the listed velocity offsets. For non-detections, the integrated intensity upper limit is found from the RMS in a 30 km s −1 window centered on the line position (Section 3.2.1). The line widths used to calculate upper limits are taken to be the average of detected CH 3 OH lines towards the same source, shown in italics. Figure 6 shows each spectrum centered around the expected transitions of PO and PN. For each transition, our fitting routine did not identify any emission line at a significance >2σ. For comparison, we also show SO lines observed towards these sources. . Non-detections of PN and PO in the sources where no lines were detected. We also include SO lines for comparison: 99.30 GHz for L1527 and L1551, and 138.18 GHz for Ser-SMM3. Ser-SMM3 and L1551 were observed in only one setup each. Velocities are centered on the rest velocity of each source (Table 1).

C. CH 3 OH AND C 17 O FITTING
The results from CH 3 OH line fitting are listed in Table 8. To find the CH 3 OH column density in each source, we used the Markov Chain Monte Carlo (MCMC) package emcee (Foreman-Mackey et al. 2013) to fit for both the total column density (N T ) and rotational temperature (T R ), as described in Section 3.2. In three sources (NGC1333-IRAS 10 10 10 11 T R -22.6 K N T -3.4e+14 Ser-SMM4 10 11 T R,cold -5.3 K T R,hot -67.0 K N T,cold -3.9e+14 N T,hot -5.1e+14 Ser-SMM1 10 10 10 11 T R,cold -6.8 K T R,hot -75.4 K N T,cold -6.5e+13 N T,hot -4.5e+14 NGC1333-IRAS 3 10 9 10 10 10 11 T R,cold -6.8 K T R,hot -61.8 K N T,cold -8.1e+13 N T,hot -7.2e+13 T R,cold -5.7 K T R,hot -57.0 K N T,cold -2.1e+14 N T,hot -9.8e+13 B1-b Figure 7. CH3OH rotational diagrams for the six sources with firm or tentative phosphorus molecule emission. Ser-SMM4 is well fit with a single temperature component, while the other five sources required a two-component fit. The resulting column densities (cm −2 ) and rotational temperatures (K) are listed in the top right corner of each panel, as well as in Table 4.
3, L723, Ser-SMM1) the CH 3 OH upper-level populations are not well described by a single rotational temperature. For these sources we fitted the observed upper-level populations as the sum of a hot and cold temperature component with independent column densities. The resulting rotational diagrams for the sources with firm or tentative P molecule detections are shown in Figure 7.
C 17 O 1-0 hyperfine lines were covered for five sources, and were fitted as described in Bergner et al. (2019) to derive the C 17 O column density and rotational temperature. The observed spectra and best-fit models are shown in Figure  8.   Figure 8. C 17 O 1-0 observed spectra (grey) and best-fit spectral models (blue). The resulting column density (cm −2 ) and rotational temperature (K) for each source is listed.