Laboratory and Astronomical Detection of the SiP Radical (X²Π_i): More Circumstellar Phosphorus

The initial search for SiP was based on the previous laser induced fluorescence (LIF) study of Jakubek et al. (2002) and focused on the lower-lying spin-orbit ladder, Ω = 3/2, which was predicted to consist of closely spaced lambda doublets. After over 20 GHz of scanning (∼2.5 B), a repeating set of doublets was found within ∼1 GHz of the predictions, separated by a few MHz. As scanning continued to lower frequency, each doublet then further split into two lines, creating a closely spaced quartet. As the frequency further decreased, the features evolved into a broad doublet.

This behavior can be explained by a combination of lambda doubling and hyperfine interactions, the latter arising from phosphorus nuclear spin of I = 1/2. The lambda doubling for the Ω = 3/2 state increases with J, from ∼0.2 MHz at 151 GHz to ∼6 MHz at 485 GHz. The hyperfine splitting does the opposite. As a consequence, at higher frequency, only a lambda doublet is observed, as the hyperfine splitting collapses. As the frequency decreases, the hyperfine components from the doublets cross and then blend together, creating another "doublet," but with much broader features. Along with chemical tests, which showed the lines depended on both phosphorus and silane, this pattern securely identified the spectra as arising from the Ω = 3/2 ladder of SiP. Lambda doublets arising from the Ω = 1/2 spin-orbit ladder were subsequently found; their splitting was considerably larger, as predicted, but the spectra were obviously perturbed. These data are the subject of another paper (Burton et al. 2022).

The pattern observed for the Ω = 3/2 spin-orbit component of SiP is illustrated in Figure 1. Here, measured spectra of the J = 30.5 ← 29.5, J = 14.5 ← 13.5, and J = 9.5 ← 8.5 rotational transitions of SiP in the Ω = 3/2 ladder near 485 GHz, 231 GHz, and 151 GHz are shown. The lambda doublets are labeled by e and f , and the sets of hyperfine components by F_e or F_f , depending on their lambda doublet, where F = J + I. The spectrum in the upper panel (J = 30.5 ← 29.5) only shows the lambda doublets, as the hyperfine splitting is collapsed. The quartet pattern is clearly visible in the middle panel (J = 14.5 ← 13.5), as the hyperfine interaction increases to lower frequency. The splitting continues to increase until two hyperfine components of opposite lambda doublets cross and then blend together, creating two broad features, as shown in the lower panel (J = 9.5 ← 8.5).

Zoom In Zoom Out Reset image size

Overall, 15 rotational transitions (J + 1 ← J) were measured for SiP in the Ω = 3/2 ladder in the frequency range 151–533 GHz, as listed in Table 1. The data were analyzed using a Hund's case (c) Hamiltonian consisting of rotational, effective lambda doubling (q_eff ) and magnetic hyperfine (h) terms. Note that h = aΛ + (b + c)Σ = a + 1/2 (b + c) in this case. The results of this analysis are given in Table 2; the rms of the fit was 55 kHz. For comparison, the rotational constant for the Ω = 3/2 ladder was estimated from the LIF study of Jakubek et al. (2002). The derived constant, B_3/2 ∼ 7917 MHz, is within 50 MHz of the measured mm-wave value, in good agreement.

Following the laboratory measurements, a search was conducted for the 2 and 1.2 mm transitions of SiP toward the envelope of AGB star IRC+10216. Unfortunately, four of the eight transitions were completely contaminated by other strong lines. The J = 14.5 → 13.5 transition near 230 GHz was masked by CO, while the J = 9.5 → 8.5, J = 10.5 → 9.5, and J = 15.5 →14.5 transitions were obscured by C₄H v₇ = 1, C¹³CCN, and C₄H v₇ = 1, respectively. Only the J = 8.5 → 7.5, J = 13.5 → 12.5, J = 16.5 → 15.5, and J = 17.5 → 16.5 were not completely hidden. As shown in Figure 2, the J = 13.5 → 12.5 and J = 16.5 → 15.5 were clean of contamination and appear as individual lines in IRC+10216 (middle panels), although slightly broadened by lambda doubling and phosphorus hyperfine structure. The J = 8.5 → 7.5 line is a "shoulder" on the J = 48.5 → 47.5 transition of C₆H (top panel), while the J = 17.5 → 16.5 transition appears as a "trough" between features of Si¹³CC and SiC₂ v₃ = 1 (bottom panel); emission is clearly present at the respective frequencies. There are no other molecules that can account for the SiP features. Integration times are 113, 112, 103, and 61 hr at 278, 262, 215, and 135 GHz, respectively.

Zoom In Zoom Out Reset image size

The line parameters T_A*(K), V_LSR and ΔV_1/2 for SiP were derived from the line profiles and are listed in Table 3. As shown in the table, the lines are quite weak with T_A*(K) ≤ 1.2 mK, and appear at the expected velocity, V_LSR = −26 km s⁻¹. The nominal linewidths (∼36–39 km s⁻¹) are slightly broader than those typically found for this object (ΔV_1/2 ∼ 30 km s⁻¹) because of hyperfine splitting.

The abundances for SiP were derived using the non-LTE radiative transfer code developed in the Ziurys group, ESCAPADE (Adande et al. 2013). The code, based on the Sobolev approximation, employs the escape probability formalism. Levels are populated by collisions and radiation from dust, and spectra are predicted as a function of r, the distance from the star. A density profile was assumed based on a mass-loss rate of ∼10⁻⁵ M_⊙ yr⁻¹ (Anderson & Ziurys 2014), and the temperature law T_K (r) = 200 K (10¹⁵ cm/r)^{0 .73} was employed (Crosas & Menten 1997). Here 10¹⁵ cm is the distance from the star where the gas kinetic temperature is T_K = 200 K. The dust was modeled as a blackbody at a temperature of 400 K (Crosas & Menten 1997). A stellar radius of R_* ∼ 3.5 × 10¹³ cm and an average expansion velocity of V_exp = 14.5 km s⁻¹ were assumed. The calculation was started at r ∼ 7.0 × 10¹⁴ cm ("r_inner "), the edge of the dust formation zone (∼20 R_*) where the envelope reaches its terminal velocity (Men'shchikov et al. 2001).

The line profiles for SiP were modeled with a shell distribution following the equation $f\left(r\right)=\ {f}_{0\ }\exp \left(-{\left(\tfrac{r-{r}_{{shell}}}{{r}_{{outer}}}\right)}^{2}\right)\ \$(Adande et al. 2013; Anderson & Ziurys 2014), with variables peak abundance f₀ , r_shell , and r_outer . The radius r_shell is the distance from the star where the abundance reaches a maximum, and the radius r_outer is when the abundance drops by a factor of 1/e from its peak value, defining the width of the shell with maximum radius r_max = r_shell + r_outer . The fitting parameters f₀ , r_shell , and r_outer were varied over a grid of values. The ranges were 10⁻¹⁰ –10⁻⁸, 200–800 R_*, and 20 R_*–500 R_*, for f₀ , r_shell , and r_outer , respectively. The best fit was determined by a reduced chi-squared analysis, along with visual inspection. Initially, the unblended J = 13.5 → 12.5 and J = 16.5 → 15.5 lines, along with the shoulder feature of the J = 8.5 → 7.5 transition, were fit. Then the blended J = 17.5 → 16.5 transition was modeled with the same parameters, and the predicted line profile was consistent with the observed spectra. Neighboring lines were not fit in the analysis. A purely spherical distribution was also attempted (here r_shell = 0), but with poorer results, based on the above criteria.

The collision cross-sections were scaled from those of SiS (Schöier et al. 2005); infrared excitation was considered for the v = 1 level of SiP, which lies ∼615 cm⁻¹ above the v = 0 state (Jakubek et al. 2002). The dipole moment derivative for SiP is not known and was estimated from that of SiS. Only a theoretical dipole exists for SiP: 0.9 D (Ornellas et al. 2000). Uncertainty in the modeling results is estimated at ∼30%. It should be noted that the A²Σ⁺ excited state of SiP lies ∼427 cm⁻¹ in energy above the ground state, lower than the v = 1 level. However, Frank-Condon factors favor the v = 1 level.

For abundance comparisons, CP and PN spectra were also modeled, based on three, relatively contaminant-free transitions at 1, 2 and 3 mm, measured at ARO (see Table 3). PN is a well-characterized molecule, with a known transition rates and cross-sections (see EXOMOL, LAMBDA, BASECOL). These data do not exist for CP, except for a theoretically derived dipole moment (0.86 D; Rohlfing & Almöf 1988). These quantities were therefore scaled from those of CS (Schöier et al. 2005). The same model parameters were used as for SiP, and, again, a shell distribution gave a better fit. Modeled abundances are listed in Table 3. Uncertainties are estimated at ∼20%.

The modeling suggests that the abundance of SiP is f₀ ∼ 2 × 10⁻⁹ at its peak value, relative to H₂, located in a shell centered at ∼300 R_*, extending out to ∼550 R_* (1/e point). Moving the shell maximum further from the star resulted in less accurate fits. However, the SiP lines are very weak and the modeling only provides a very rough estimate, given the assumptions and the lack of collisional and radiative parameters for this molecule.

The PN and CP spectra, in comparison, had better signal-to-noise ratios as they are ∼3 times stronger. Two of the four lines analyzed for both molecules were uncontaminated, while the third feature was partially blended. The peak abundance determined for PN in IRC+10216 was f₀ ∼ 9 × 10⁻¹⁰, and the molecule was found in a shell at ∼550 R_*, extending out to ∼700 R_*. The line profiles could not be reproduced with a spherical model, which had previously been used by Milam et al. (2008). These authors, however, derived a similar abundance of f₀ ∼ 3 × 10⁻¹⁰, with an outer radius of ∼600 R_*, based on the same line measured here but with less sensitivity. For CP, the modeling resulted in f₀ ∼ 7 × 10⁻⁹ for a shell at ∼650 R_*, extending out to ∼900 R_*. Milam et al. determined an abundance of f₀ ∼ 10⁻⁸ for this radical, in close agreement, determined from a single transition. In this case, a shell distribution was assumed, with the shell maximum at ∼460 R_*, extending out to ∼770 R_*.

Using the derived abundances, ESCAPADE was used to predict the intensities of the 3 mm transitions (J = 4.5 → 3.5, 5.5 → 4.5 and 6.5 → 5.5 at 71.7, 87.6 and 103.5 GHz, respectively) in IRC+10216 at the ARO 12 m telescope. Values obtained were T_A* < 0.3 mK—consistent with the upper limits of 1–2 mK derived from 12 m observations of this source at 3 mm.

SiP is generated in a shell that is closer into the star than either CP or PN at r ∼ 300 R_*, SiP formation may be linked to the liberation of both phosphorus and silicon from dust grains. Once formed, SiP is then channeled into other P-bearing molecules, such as PN and CP, which have shells near 550–650 R_*.

Thermodynamic calculations predict that phosphorus should be condensed into the dust grains as the mineral schreibersite, (Fe,Ni)₃P, in circumstellar material (e.g., Lodders & Fegley 1999; Lodders 2003). The condensation temperature of schreibersite is T ∼ 1200 K at 10⁻⁴ atm Lodders 2003). Schreibersite is commonly found in meteorites, including primitive carbonaceous chondrites (Pasek 2019). It is therefore possible that at least some phosphorus condenses into schreibersite in the circumstellar dust formation zone. However, some phosphorus may be liberated from the solid state by shocks in the inner envelope. The phosphorus may also remain in the gas phase and avoid condensation, which may not be an equilibrium process (Höfner & Olofsson 2018).

Unlike the other observed P-bearing species in IRC+10216, SiP also contains the refractory element silicon. There are a wide range of Si-containing molecules in IRC+10216, so gas-phase silicon is readily available, perhaps arising from the shock processing of SiC grains (e.g., Bernal et al. 2019). If such shocks also disrupt schreibersite, then SiP could form in the heated material. This scenario could explain its relatively confined spatial distribution.

Recently, Golubev et al. (2013) conducted a theoretical study of the interstellar formation of SiP. Based on a previous work by Andreazza et al. (2006), these authors examined the formation of the radical by radiative association of silicon and phosphorus atoms. They suggest that the reaction Si(³P) + P(⁴S) creates SiP in its excited state A²Σ⁺, which, for stabilization, decays to the nearby ground state, releasing a photon:

$$\begin{eqnarray}&&{\rm{S}}{\rm{i}}{(}^{3}{\rm{P}})+{\rm{P}}{(}^{4}{\rm{S}})\,\to \,{\rm{S}}{\rm{i}}{\rm{P}}({{\rm{A}}}^{2}{{\rm{\Sigma }}}^{+})\,\to \,{\rm{S}}{\rm{i}}{\rm{P}}({{\rm{X}}}^{2}{\rm{\Pi }})+h\nu .\end{eqnarray} \tag{ 1 }$$

This reaction is feasible because both the X and A states share the same dissociation limit and are separated in energy by only 427 cm⁻¹ (Jakubek et al. 2002). Golubev et al. (2013) calculated the rate of this reaction to be ∼10⁻²⁰ cm³ s⁻¹ for 10–1000 K.

There are no reactions for the formation or destruction of SiP listed in the KIDA and UMIST databases (Wakelam et al. 2012; McElroy et al. 2013). However, in analog to SiN, there are two possible alternative gas-phase synthesis processes:

$$\begin{eqnarray}&&{\rm{P}}+\mathrm{SiC}\,\to \,\mathrm{SiP}+{\rm{C}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{P}}+{\rm{S}}{\rm{i}}{\rm{H}}\,\to \,{\rm{S}}{\rm{i}}{\rm{P}}+{\rm{H}}.\end{eqnarray} \tag{ 3 }$$

Based on the nitrogen analogs, these reactions are likely to be exothermic. The first reaction can also lead to CP + Si, with an unknown branching ratio. Assuming an equal product distribution, this reaction has a projected rate of k₂ ∼ 2 ×10⁻¹¹ cm³ s⁻¹, based on SiN. The second pathway is estimated to have k₃ ∼ 2 × 10⁻¹⁰ cm³ s⁻¹, also extrapolated from SiN. Both processes may be sufficiently fast to be the major gas-phase pathways to SiP. Measurements and/or theoretical estimates of these reaction rates are certainly needed.