Deep Search for Phosphine in a Prestellar Core

Understanding in which chemical forms phosphorus exists in star- and planet-forming regions and how phosphorus is delivered to planets are of great interest from the viewpoint of the origin of life on Earth. Phosphine (PH3) is thought to be a key species to understanding phosphorus chemistry, but never has been detected in star- and planet-forming regions. We performed sensitive observations of the ortho-PH3 10 − 00 transition (266.944 GHz) toward the low-mass prestellar core L1544 with the Atacama Compact Array stand-alone mode of the Atacama Large Millimeter/submillimeter Array. The line was not detected down to 3σ levels in 0.07 km s−1 channels of 18 mK. The nondetection provides the upper limit to the gas-phase PH3 abundance of 5 × 10−12 with respect to H2 in the central part of the core. Based on the gas-ice astrochemical modeling, we find the scaling relationship between the gas-phase PH3 abundance and the volatile (gas and ice with larger volatility than water) P elemental abundance for given physical conditions. This characteristic and well-constrained physical properties of L1544 allow us to constrain the upper limit to the volatile P elemental abundance of 5 × 10−9, which is a factor of 60 lower than the overall P abundance in the interstellar medium. Then the majority of P should exist in refractory forms. The volatile P elemental abundance of L1544 is smaller than that in the coma of comet 67P/C-G, implying that the conversion of refractory phosphorus to volatile phosphorus could have occurred along the trail from the presolar core to the protosolar disk through, e.g., sputtering by accretion/outflow shocks.


INTRODUCTION
Phosphorus-bearing molecules, such as nucleic acids, are essential for life on Earth.Then, understanding in which chemical forms phosphorus exist in star-and planet-forming regions and how phosphorus is delivered to planets are of great interest from the viewpoint of the origin of life on Earth.The recent detection of phosphorus-bearing molecules, PO, in the coma of comet 67P/Churyumov-Gerasimenko (Altwegg et al. 2016;Rivilla et al. 2020) indicates that at least some fraction of phosphorus was present as volatiles (i.e., gas and ice) in the protosolar disk, whereas in meteorites, phosphorus is identified as inorganic minerals (Pasek & Lauretta 2005) and alkyl phosphonic acids (Cooper et al. 1992).
In low-mass and high-mass star-forming regions, PO and PN in the gas phase are the only neutral P-bearing molecules detected so far (e.g., Turner & Bally 1987;Ziurys 1987;Fontani et al. 2016;Lefloch et al. 2016;Rivilla et al. 2016;Mininni et al. 2018;Rivilla et al. 2020;Bergner et al. 2022;Koelemay et al. 2023).In addition, the detection of PO + was reported in a molecular cloud located in the Galactic center (Rivilla et al. 2022).These gas-phase molecules are detected in shocked regions, but not in hot core/corino regions, implying that P in the solid phase is sputtered by shocks, and PO and PN are formed by gas-phase reactions, assisted by photochemistry (Jiménez-Serra et al. 2018;Rivilla et al. 2020).
However, to date, no infrared detection of solidstate P-bearing species are reported in any interstellar sources; hence, there is no observational constraint on solid-state P. Astrochemical models have predicted that phosphine (PH 3 ) ice is the dominant carrier of volatile P in molecular clouds (e.g., Charnley & Millar 1994;Chantzos et al. 2020).Then there are several theoretical studies to study the phosphorus chemistry in the shock regions, assuming that P is initially locked up in PH 3 (e.g., Aota & Aikawa 2012;Jiménez-Serra et al. 2018).
Despite the possible importance of PH 3 to understand the phosphorus chemistry, very little is known about PH 3 in star-forming regions.Because PH 3 is relatively volatile species (the sublimation temperature is ∼40 K (Molpeceres & Kästner 2021) and PH 3 ice can be nonthermally desorbed even at ∼10 K, see below), PH 3 in the gas-phase may be detectable by radio observations.However, the PH 3 gas has never been detected in molecular clouds and star-forming regions (e.g., Turner et al. 1990;Lefloch et al. 2016).To the best of our knowledge, the only detection of PH 3 has been reported in the circumstellar envelope of a carbon star (Agúndez et al. 2014).
In contrast to the poor observational understanding of PH 3 , its chemistry has been studied in detail by quantum chemical calculations and laboratory experiments.The formation of PH 3 via gas-phase reactions is inefficient at low temperatures (∼10 K), because of the endothermicity of relevant reactions (Thorne et al. 1984).Thus, PH 3 is thought to be formed on grain surfaces by the sequential hydrogenation of atomic P through barrierless reactions on grain surfaces (Molpeceres & Kästner 2021).Once PH 3 is formed on grain surfaces as ice, PH 3 ice can be sublimated into the gas-phase at ≳40 K (Molpeceres & Kästner 2021) or non-thermally desorbed by chemical desorption (desorption induced by the energy released by chemical reactions; Nguyen et al. 2020Nguyen et al. , 2021;;Furuya et al. 2022), by photodesorption or by shock-induced desorption (e.g., Lefloch et al. 2016;Jiménez-Serra et al. 2018).Thermal desorption should be the dominant desorption mechanism at ≳40 K, while at lower temperatures, non-thermal desorption processes would be more efficient than thermal desorption.
Then there would be two possible approaches to detect PH 3 : observing thermally desorbed PH 3 in the warm (≳40 K) inner envelopes of protostellar sources (hot corinos/cores) or observing non-thermally desorbed PH 3 in prestellar cores.Considering the non-detection of PH 3 , but the detection of PO and PN in the outflow shock region L1157-B1 (Lefloch et al. 2016), PH 3 might be short-lived in the warm gas and quickly converted into PO and PN as suggested by astrochemical models (Aota & Aikawa 2012;Jiménez-Serra et al. 2018).Therefore, we focus on the latter approach in this work.
We report ACA stand-alone observations of the ortho-PH 3 1 0 − 0 0 transition (266.944GHz) toward the dust continuum peak of low-mass prestellar core L1544 in Taurus at a distance of 135 pc (Schlafly et al. 2014).L1544 is one of the best-studied prestellar cores from both the physical and chemical points of view and is a good target for the search of PH 3 (e.g., Crapsi et al. 2007;Spezzano et al. 2017;Furuya et al. 2018;Caselli et al. 2019;Redaelli et al. 2021).The detection of hydrides in the gas phase toward the dust continuum peak of the core has been reported: H 2 O (Caselli et al. 2012), NH 3 (Caselli et al. 2017), and H 2 S (Vastel et al. 2018).These hydrides are likely formed on the grain surface via a sequence of hydrogenation reactions of atoms and released into the gas phase by non-thermal desorption mechanisms.Judging from the large Einstein A coefficient and the small upper state energy (12.8 K; Müller et al. 2005) of the ortho-PH 3 1 0 − 0 0 , the transition should be the strongest one in the cold gas (<10 K).For example, the upper state energy of ortho-PH 3 2 0 − 1 0 transition at 533.795 GHz is 38.4K (Müller et al. 2005), which is much higher than the gas temperature of the prestellar core.As the critical density of ortho-PH 3 1 0 − 0 0 is high (>10 5 cm −3 , Appendix A), the central region of L1544 inside a radius of ∼1000 au, where the density is larger than 10 6 cm −3 (Caselli et al. 2019), is the suitable place to observe the transition.In addition, high spatial resolution observations with ACA allow us to selectively probe PH 3 in the central region of the core, where PH 3 is expected to be non-thermally desorbed by chemical desorption.Gas-phase PH 3 desorbed by the chemical desorption is a good probe of the elemental abundance of volatile phosphorus (see Section 4.2).
The source had been observed with the ACA standalone mode to map the dust continuum emission (Caselli et al. 2019) and NH 2 D lines (Caselli et al. 2022).There has been no detection of P-bearing species in L1544, while Turner et al. (1990) and Rivilla et al. (2018) reported the non-detection of PN with the NRAO 12 m telescope and IRAM 30 m telescope, respectively.

OBSERVATIONS
Observations were carried out with the ACA standalone mode of ALMA in July and August 2023 as a Cycle 9 program (2022.1.01490.S, PI: K. Furuya).The total number of 7 m antennas is 9-11, where the minimum-maximum baseline lengths are 8.9-48.9m.The total on-source integration time is 7. ing center of antennas is RA = 05 h 04 m 17. s 21 and Dec = +25 • 10 ′ 42.′′ 8 (International Celestial Reference System, ICRS), which corresponds to the dust continuum peak of L1544 (Spezzano et al. 2016).The sky frequencies of 250.41-250.47, 250.79-250.85, 265.86-265.92, 266.89-267.01, 267.53-267.59GHz, are covered by five spectral windows with a velocity resolution of 0.07 km s −1 , where the forth window is used for the analysis of the PH 3 transition.The sky frequency of 252.06-253.94GHz is also observed with a continuum mode.
Raw data is processed with the Common Astronomy Software Applications (CASA) package (version 6.4.1).The CASA task tclean is used for imaging, and the masking is done using the auto-multithresh algorithm.The continuum emission is subtracted from the spectral data using the uvcontsub task in CASA.The synthesized beam size at the frequency of the PH 3 transition is 4. ′′ 8 × 5. ′′ 8 with the Briggs weighting and a robustness parameter of 0.5.The beam size corresponds to 650 au × 780 au at the distance of L1544.The maximum recoverable scale is about 25 ′′ .The synthesized image is corrected for the primary beam pattern using the impbcor task in CASA.The aggregate continuum image is constructed by using the continuum window in addition to the line-free channels selected from the other five spectral windows.The spectrum and continuum flux are extracted from the 5. ′′ 0 × 9. ′′ 0 (700 au × 1260 au) elliptical region centered at RA = 05 h 04 m 17. s 01 and Dec = +25 • 10 ′ 42.′′ 60, which corresponds to the continuum peak in our observations.The elliptical region roughly corresponds to the brightest regions in the dust continuum emission (Fig. 1).

OBSERVATIONAL RESULTS
The ortho-PH 3 1 0 − 0 0 transition was not detected as shown in Figure 2. We obtain the upper limit to the PH 3 column density from the 3σ upper limit to the integrated intensity 3σ √ ∆vδv, where σ is the root-meansquare noise level of the spectra, ∆v is the FWHM line width of the spectra, and δv is the velocity resolution.∆v is assumed to be 0.3 km s −1 , similar to those of the high-density tracers, H 2 D + and D 2 H + , obtained with IRAM 30 m single-dish telescope toward L1544 (Vastel et al. 2004).Assuming local thermal equilibrium (LTE), the 3σ intensity upper limit is converted to the upper limit to the PH 3 column density (N (PH 3 )) following Mangum & Shirley (2015).See Appendix A for the discussion on the validity of the LTE assumption.The parameters of the observed transition are taken from the Cologne Database for Molecular Spectroscopy (Müller et al. 2005) and Leiden Atomic and Molecular Database (Schöier et al. 2005).The excitation temperature (T ex ) is assumed to be the same as the gas temperature of the core center, 6 K (Keto et al. 2015).As the production of PH 3 by the gas-phase chemistry is negligible (Thorne et al. 1984), gas-phase PH 3 is most likely originated from ices.Then we assume that the ortho-to-para ratio of PH 3 is unity (i.e., the statistical ratio; see Hama et al. 2016).As a result, the upper limit to N (PH 3 ) is estimated to be 1.2 × 10 11 cm −2 .Assuming the dust temperature of 6 K, We estimate the H 2 column density of the examined region to be 2.6 × 10 22 cm −2 based on the dust continuum data (see Appendix B for details).The optical depth of the dust emission is estimated to be ∼0.001.Accordingly, the upper limit to the gas-phase PH 3 abundance is 5 × 10 −12 with respect to H 2 .If T ex is assumed to be 9 K (see Section 4), the upper limits to N (PH 3 ) and the PH 3 abundance are 7.6 × 10 10 cm −2 and 6.7 × 10 −12 , respectively.
As many previous observations toward L1544 were done with single-dish telescopes, IRAM 30m telescope in particular, it may be useful to provide the upper limit to N (PH 3 ) with the corresponding beam size.We perform a similar analysis to the spectrum extracted from the 9. ′′ 2 × 9. ′′ 2 region, which corresponds to IRAM 30m's beam at 267 GHz, centered to the continuum peak.When T ex = 6 K and 9 K, the upper limit to N (PH 3 ) is constrained to be 9.6 × 10 10 cm −2 and 5.9 × 10 10 cm −2 , respectively.The upper limit to the PH 3 abundance is 4 × 10 −12 and 6 × 10 −12 at T ex = 6 K and 9 K, respectively.

VOLATILE PHOSPHORUS ABUNDANCE IN L1544
To constrain the volatile P abundance from the upper limit to the PH 3 abundance, we conduct gas-ice . The spectrum centered on 266.944GHz.The vertical line marks the VLSR of L1544 (7.2 km s −1 ).The ortho-PH3 10 − 00 transition was not detected down to 3σ levels in 0.07 km s −1 channels of 18 mK (horizontal dotted line).
astrochemical simulations under the physical conditions appropriate for L1544.

Model description
We utilize the gas-ice astrochemical code based on the rate equation approach (Rokko code; Furuya et al. 2015).We adopt a three-phase model, where the gas phase, a surface of ice, and the chemically inert bulk ice mantle are considered (Hasegawa & Herbst 1993a).The top four ice layers are assumed to be the surface (Vasyunin & Herbst 2013), and the rest is considered as the bulk ice mantle.Our gas-ice chemical network is based on that of Garrod (2013) with some modification to the P chemical network as explained below.The gas phase chemical network relevant to PH 3 is updated referring to previous studies (Charnley & Millar 1994;García de la Concepción et al. 2021;Fernández-Ruz et al. 2023;Gomes et al. 2023;Millar et al. 2023).As the surface chemical network relevant to PH 3 , we consider a sequence of hydrogenation addition reactions of atomic P to form PH 3 and the hydrogen abstraction reaction from PH 3 by atomic H (Nguyen et al. 2020).The hydrogen addition reactions are barrierless, while the hydrogen abstraction reaction has a barrier with the tunneling corrected rate constant of 8.7 × 10 7 s −1 (Molpeceres & Kästner 2021).Laboratory experiments by Turner et al. (2016) showed that electron irradiation to mixed ices of PH 3 and CH 4 leads to the formation of phosphanes and methylphosphanes as large as P 8 H 10 and CH 3 P 8 H 9 , whose sublimation temperature is higher than that of water.Although it is unclear whether such a mechanism can work in the water-ice-dominated ISM ice, the current chemical network is limited to only small P-bearing molecules, such as PN, PO, CP, and PH n , where n is 1-3.The volatility of these species is higher than that of water (Sil et al. 2021;Molpeceres & Kästner 2021).From this point of view, our model only considers the volatile P chemistry.
We consider a binding energy distribution for all surface species using the method developed by Furuya (2023), assuming the binding energy distribution follows a Gaussian distribution.In particular, the mean binding energy of atomic H is set to be 350 K, while that of atomic P and PH n , where n is 1-3, is taken from Molpeceres & Kästner (2021).Following Cazaux et al. (2017), the activation energy for surface hopping from a site with the binding energy of E to another site with the binding energy of E ′ is given by where χ is set to 0.65 for atomic H and H 2 and to 0.4 for the other species.E hop for atomic H in our model is consistent with that constrained by laboratory experiments; there are deep potential sites with E hop ≳ 350 K, while the majority of sites are shallow sites with E hop ≲ 250 K on amorphous solid water (Hama et al. 2012).
Due to the low temperature in the core (<10 K), thermal desorption is negligible for PH n (Molpeceres & Kästner 2021).As non-thermal desorption, we consider desorption by the stochastic heating by CRs (Hasegawa & Herbst 1993b), sputtering by CRs (Dartois et al. 2018), photodesorption, and chemical desorption.To the best of our knowledge, there is no study on the CR sputtering nor the photodesorption of P-bearing species.We calculate the CR sputtering yield for all species in a way proposed by Dartois et al. (2021), using the parameters appropriate for water ice.We assume the yield of 10 −3 per incident UV photon for all P-bearing species.The chemical desorption probability of PH 3 upon reaction between PH 2 and atomic H is set to 4 % as constrained by laboratory experiments and their modeling (Nguyen et al. 2021;Furuya et al. 2022).Note that those studies on chemical desorption assumed that the grain surfaces are covered by water ice.The detection of water vapor at the center of L1544 likely indicates that water ice is present on the ice surface and the water vapor is produced via photodesorption by cosmic-ray induced UV photons (Caselli et al. 2012).
There are various studies on the physical structure of L1544 based on dust continuum and molecular line observations.As discussed below, the two of the most important physical parameters for determining the gasphase PH 3 abundance in our models are the dust temperature and the ratio of the CR ionization rate (ξ) to the H 2 gas density (n H2 ).According to the physical model of L1544 proposed by Keto et al. (2015), the  Keto et al. (2015).Because the abundances of the major molecular ions depend on ξ/n H2 (Dalgarno 2006) like the gas-phase PH 3 abundnance, we assume ξ = 3 × 10 −17 s −1 throughout this work.
We consider four different models, varying in temperature, gas density, and initial abundances, which are summarized in Table 1.In each model, we vary the volatile P elemental abundance from 2.57 × 10 −7 (i.e., the solar value; Asplund et al. 2021) to 2.57 × 10 −9 , while the abundances of other elements are taken from Aikawa & Herbst (1999).In Models A, C, and D, initially, hydrogen is assumed to be present as H 2 , while the other elements are assumed to be present as atoms or atomic ions.In Model B, we first evolve the gas-ice chemistry for 1 Myr under the dense cloud condition (n H2 = 1 × 10 4 cm −3 , the temperature of 10 K, and A V =10 mag), and the final abundances are used as the initial abundances of the L1544 model.

Model results
The left and right panels of Figure 3 show the temporal evolution of P-bearing species in Model A and Model B. While the abundances of ice species depend on the choice of initial abundances, the gas-phase PH 3 abundance is similar between the two models.In our models, the chemical desorption of PH 3 is the dominant route for the supply of gas-phase PH 3 , and the other non-thermal desorption mechanisms are negligible.PH 3 ice is formed by a sequence of hydrogenation of atomic P on grain surfaces.Through the hydrogen abstraction reaction on the grain surface, PH 3 ice is converted to PH 2 ice, and the hydrogenation of PH 2 ice produces PH 3 ice again, PH 2 + H → PH 3 . (3) When PH 3 ice is produced by Reaction 3, 4 % of the produced PH 3 ice desorbs to the gas phase (the rest remains as ice) by the chemical desorption as assumed in our models.The PH 3 -PH 2 loop with PH 3 chemical desorption has been observed in laboratory experiments (Nguyen et al. 2020(Nguyen et al. , 2021)).In our models, the PH 3 -PH 2 loop keeps a non-negligible fraction of PH 3 in the gas phase, because PH 3 can desorb in each cycle of the PH 3 -PH 2 loop, and because the loop continues as long as atomic H is available on grain surfaces.The abundance of atomic H in the gas phase depends on the ξ/n H2 ratio (e.g., Tielens 2005), and thus the gas-phase PH 3 abundance depends on the ξ/n H2 ratio.We also confirmed that the gas phase PH 3 abundance linearly depends on the chemical desorption probability, by additionally running the modes with the desorption probability of 2 % and 8 %.In our models, gas-phase PH 3 is mainly destroyed by PH 3 + H 3 + → PH 4 + + H 2 , followed by the electron dissociative recombination to form PH 2 + H 2 .The dissociation rate through this pathway is similar to the adsorption rate of PH 3 onto dust grains.Photodissociation by CR-induced UV photons and the neutralneutral reaction PH 3 + H → PH 2 + H 2 are included in our models, but these are negligible compared to the pathways mentioned above.
In our models, PH 3 is the most abundant gas-phase P-bearing molecule over PO and PN.PN is mainly produced by PH + N → PN + H, while PO is mainly produced by P + OH → PO + H as in previous models (Charnley & Millar 1994;Aota & Aikawa 2012;Jiménez-Serra et al. 2018).However, as the abundance of N atoms and OH in the gas phase is very low, as oxygen is locked up in H 2 O ice and nitrogen is locked up in N 2 and NH 3 ices, a negligible fraction of atomic P and PH are converted to PO and PN.Rivilla et al. (2018) reported that the PN/C 34 S column density ratio toward L1544 is lower than ∼ 8 × 10 −3 using the data obtained with IRAM 30m telescope (Jiménez-Serra et al. 2016).According to Vastel et al. (2018) who observed L1544 with IRAM 30m telescope, the C 34 S column density toward the dust continuum peak is (7.8 − 8.3) × 10 11 cm −2 .Taken together, the upper limit to the PN column density is ∼ 6×10 9 cm −2 , which is around one order of magnitude lower than that of PH 3 and corresponds to the upper limit abundance of 10 −12 -10 −13 .Such a low abundance of PN in L1544 is consistent with our model predictions, where the gas-phase PN abundance is ∼10 −15 after around 10 4 yr. Figure 4 shows the gas-phase PH 3 abundance at 10 5 yr as a function of assumed volatile P elemental abundance.10 5 yr was chosen, because at that timescale the observationally derived CO-to-H 2 column density ratio of L1544 is reproduced by a gas-ice astrochemical model, considering 1D physical structure of the core (Vasyunin et al. 2017).The gas-phase PH 3 abundance is almost proportional to the assumed total volatile P elemental abundance in all four models.This characteristics allows us to constrain the volatile P abundance from the upper limit to the gas-phase PH 3 abundance.Then we conclude that the volatile P elemental abundance in L1544 is less than 5 × 10 −9 ., i.e., lower than the solar value by a factor of at least ∼60.We confirmed that this conclusion holds, even if a different time other than 10 5 yr (3 × 10 4 yr or 10 6 yr) is chosen.

DISCUSSION AND CONCLUSIONS
In L1544, some hydrides have been detected.According to previous observations with single-dish telescopes and modeling studies, the gas-phase abundances of NH 3 and H 2 O with respect to H 2 in the central 1000 au region of L1544 are ∼8×10 −9 (Crapsi et al. 2007;Caselli et al. 2017) and ∼3×10 −10 (Caselli et al. 2012;Keto et al. 2015), respectively.H 2 S was detected with IRAM 30m, but due to the complicated spectral shape, abundance in the central part was not constrained well (See Appendix B of Vastel et al. 2018, for more details).The abundances of NH 3 and H 2 O are more than an order of magnitude higher than the upper limit to the abundance of gas-phase PH 3 constrained in this study.Assuming the solar abundance represents the elemental abundance in the local ISM (see Przybilla et al. 2008), the proportion of oxygen, nitrogen, and phosphorus locked up in the gas-phase hydrides in L1544 are 6 × 10 −7 , 1 × 10 −4 , and < 2×10 −5 , respectively.Despite nitrogen and phosphorus being in the same family, the proportion of nitrogen locked up in NH 3 is much higher than the proportion of phosphorus locked up in PH 3 .This difference may partly stem from the fact that in addition to the non-thermal desorption, gas-phase NH 3 can be formed via the gas-phase chemistry in the CO-depleted regions (e.g., Aikawa et al. 2005), while the gas phase formation of PH 3 is inefficient at low temperatures (≲10 K) (Thorne et al. 1984).
The solar abundance of phosphorus is 2.57×10 −7 with respect to hydrogen (Asplund et al. 2021).According to the observations of UV absorption lines proving diffuse clouds, phosphorus is depleted in the gas phase, indicating that some part of phosphorus is already incorporated in refractory dust grains in diffuse clouds (Jenkins 2009;Ritchey et al. 2023).The gas phase phosphorus abun-dance is smaller than the solar value by a factor of up to ten, depending on sight lines (Ritchey et al. 2023).
1 In this work, we showed that volatile (gas + ice) P abundance in prestellar core L1544 is less than 5 × 10 −9 and is smaller than the solar value by a factor of at least 60.Then, almost all (∼99 % or even higher) phosphorus may be locked up in refractory dust grains during the evolution from diffuse clouds to prestellar cores.
In low-mass Class 0/I protostars, PO and PN have been detected in the outflow-shocked regions, while they are not detected in hot corinos (e.g., Yamaguchi et al. 2011;Lefloch et al. 2016;Bergner et al. 2022).This is also true for high-mass star-forming regions (Rivilla et al. 2020;Fontani et al. 2024).These suggest that shock-induced sputtering of solid P carriers followed by the gas-phase chemistry is important for the formation of those species.In the outflow shock region L1157-B, volatile P abundance is still a factor of ∼100 lower than the solar value, based on the observations of PN and PO and shock-chemical models (Aota & Aikawa 2012;Lefloch et al. 2016), implying only a fraction of solid P is released into the gas phase.The observational studies of PN and PO toward high-mass star-forming regions also indicate a factor of ∼100 lower volatile P abundance than the solar value (Fontani et al. 2016;Rivilla et al. 2016).More observations of P-bearing molecules in prestellar cores and protostellar sources would help to understand how much P is released into the gas phase at which conditions.
In the coma of comet 67P/C-G, PO is the main reservoir of phosphorus and the P/O elemental ratio is (0.5 − 2.7) × 10 −4 (Rubin et al. 2019;Rivilla et al. 2020), close to the Solar value of 5.2 × 10 −4 (Asplund et al. 2021), implying that volatile P is not significantly depleted compared to L1544.Note that our astrochemical model predicts that PO is not a major reservoir of volatile P in prestellar cores (Fig. 3).Taken together, the substantial conversion of solid P in dust grains to volatile P could have occurred along the trail from the presolar core to the protosolar disk through e.g., sputtering by accretion/outflow shocks as suggested by Bergner et al. (2022) based on the PN and PO observations in the protostellar source.

ACKNOWLEDGMENTS
We are grateful to the referee for providing useful comments that helped to improve the manuscript considerably.This work is supported in part by JSPS KAKENHI Grant numbers 20H05845, 20H05847, 21H01145, and 21K13967.KF was supported by the ALMA Japan Research Grant of NAOJ ALMA Project, NAOJ-ALMA-324.This paper makes use of the following ALMA data: ADS/JAO.ALMA#2022.1.01490.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile.

A. ESTIMATION OF THE CRITICAL DENSITY AND NON-LTE ANALYSIS
To check the validity of the LTE assumption, the critical density of the ortho-PH 3 1 0 − 0 0 transition is estimated and the non-LTE analysis is performed.We assume that the main collisional partner of PH 3 for excitation and de-excitation is H 2 , and H 2 is in the para form, because the ortho-to-para ratio of H 2 is typically much less than unity in prestellar cores (e.g., Pagani et al. 2009).Hereafter, let us denote ortho-PH 3 as o-PH 3 and para-H 2 as p-H 2 .Because the rate coefficients for collisions between o-PH 3 and p-H 2 are unavailable 1 According to Ritchey et al. (2023), the transition oscillator strengths used in earlier studies (e.g., Lebouteiller et al. 2005) to derive PII column density are in need for revision and those earlier studies overestimate the gas-phase P abundance (i.e., underestimate the degree of P depletion).More details can be found in Ritchey et al. (2023).
in the literature, we estimate them in the following two ways.Let us denote the rate coefficients for collisions between species X and species Y as γ(X, Y).
The critical density of p-H 2 for the transition is estimated to be A 10 /γ, where A 10 is the Einstein A coeffi-cient for the transition.The critical density is ∼4 × 10 5 cm −3 and ∼3 × 10 5 cm −3 , when γ 1 and γ 2 are adopted, respectively.As the gas density of the central regions of L1544 (r ≲ 1000 au) is ∼10 6 cm −3 or even higher (Keto et al. 2015;Caselli et al. 2019), PH 3 1 0 − 0 0 is expected to be thermalized.
To confirm the validity of the LTE assumption, we also perform non-LTE analysis with RADEX (van der Tak et al. 2007).As in our LTE analysis presented in Section 3, we assume T ex of 6 K, the line width of 0.3 km/s, ortho-to-para ratio for PH 3 of unity, and assume that H 2 is solely present as the para-form.When γ 1 (γ 2 ) is adopted, the upper limit to N (PH 3 ) is estimated to be 1.7 × 10 11 cm −2 (1.4 × 10 11 cm −2 ) and 1.4 × 10 11 cm −2 (1.3 × 10 11 cm −2 ) for the H 2 gas density of 1 × 10 6 cm −3 and 3 × 10 6 cm −3 , respectively.These results are similar to that obtained with the LTE approximation (< 1.2 × 10 11 cm −2 ).
While we estimated γ(o-PH 3 , p-H 2 ) in the ways described above, LAMBDA database provides an approximation for γ(o-PH 3 , H 2 ) (not γ(o-PH 3 , p-H 2 )) scaling γ(o-PH 3 , He) considering the mass difference between H 2 and He.If we use γ(o-PH 3 , H 2 ) in the LAMBDA database as γ(o-PH 3 , p-H 2 ), the critical density is ∼ 5 × 10 6 cm −3 , which is an order of magnitude higher than the value estimated with γ 1 or γ 2 .The difference comes from the fact that γ(o-NH 3 , p-H 2 ) is ∼10 times higher than γ(o-NH 3 , He) (Bouhafs et al. 2017).This indicates that He is not a good proxy for H 2 at least for o-NH 3 .If the collisional rate coefficient in the LAMBDA database is used, the upper limit to N (PH 3 ) is 5×10 11 cm −2 for the gas density of 10 6 cm −3 and 2 × 10 11 cm −2 for the gas density of 3 × 10 6 cm −3 .These values are a factor of a few larger than the value obtained with the LTE approximation.The collisional rate coefficient for the PH 3 -H 2 system is highly desired.

B. DERIVATION OF THE H 2 COLUMN DENSITY
Based on the standard treatment of optically thin dust emission, we apply the following equation to calculate the H 2 column density (N (H 2 )): where F ν /Ω is the continuum flux density per beam solid angle as estimated from the observations, κ ν is the mass absorption coefficient of dust grains coated by thick ice mantles taken from Ossenkopf & Henning (1994) and we use 1.21 cm 2 g −1 at 1.16 mm, T d is the dust temperature and we assume the gas and dust temperatures are the same in this work, B ν (T d ) is the Planck function, Z is the dust-to-gas mass ratio and we use a canonical value of 0.008 for the solar neighborhood, µ is the mean atomic mass per hydrogen (1.41, Cox 2000), and m H is the hydrogen mass.
Our N (H 2 ) is a factor of ∼4 lower than the value obtained by Crapsi et al. (2005) based on IRAM 30m observations at 1.2 mm (2.6 × 10 22 cm −2 versus ∼ 1 × 10 23 cm −2 ).The dust opacity adopted in Crapsi et al. (2005) is a factor of ∼2 larger than our adopted value.On the other hand, Crapsi et al. (2005) assumed the dust temperature of 10 K, which is higher than our assumed value of 6 K.As B ν (10 K)/B ν (6 K) ∼ 1.7 at 1.2 mm, the difference due to the dust opacity and the dust temperature is almost canceled out.The contribution from the largely-extended (>25") continuum emission that was resolved out in the ACA observations may account for the above discrepancy.
Figure 1.ALMA 1.2 mm continuum emission toward L1544.The black solid ellipse represents the region where we extract the spectrum and continuum flux discussed in the main text.The synthesized beam size is shown by the gray-filled circle.North is up, and east is to the left.

Figure 3 .Figure 4 .
Figure 3. Temporal evolution of the abundances of P-bearing species in the models with the elemental P abundance of 1×10 −8 .Solid lines indicate gas-phase species, while dashed lines indicate icy species.Purple lines represent the sum of PH and PH2 abundances to increase the visibility.The left panel shows the results from Model A, while the right panel shows the results from Model B. The horizontal dotted line indicates the upper limit to the gas-phase PH3 abundance in L1544.

Table 1 .
Redaelli et al. (2021)arameters.Atomic H 2 gas density in inner ∼600 au of the core, where our ALMA observations prove, is ∼3×10 6 cm −3 , while the temperature is ∼6-9 K.Redaelli et al. (2021)constrained the CR ionization rate to be 3 × 10 −17 s −1 based on the observations of N 2 H + , N 2 D + , HC 18 O + , and DCO + and astrochemical and radiative transfer models, adopting the physical model by