The absolute cross-section of the reactive collision of P⁺(³P) with D₂: a combined theoretical and experimental study

Abstract

In this work, we combine theoretical and experimental methods to study the P⁺(³P) + D₂ → PD⁺ + D reaction. As a result, the absolute cross-section as a function of collision energy is obtained. Experimentally, the cross-section is measured using the guided ion beam technique (GIB), where P⁺ is produced by VUV photons at the SOLEIL synchrotron using PCl₃ as a precursor. Theoretically, the cross-section is calculated from first principles. The potential energy surfaces of the three electronic states correlating with the P⁺(³P) are constructed by fitting MRCI points, and reaction dynamics are performed on each of them independently, hence neglecting couplings. The total cross-section is then obtained from the weighted contribution of each considered electronic state. Our findings show good agreement between the measured and calculated cross-sections, with a small discrepancy indicating that spin–orbit and non-adiabatic coupling, not considered in this work, may play a role in this reaction. The results hereafter presented demonstrate that the chemistry of the third-row atomic cations with molecular hydrogen is generally unfavoured, unlike their second-row homologues, thus manifesting the existence of boundaries to the applicability of the so-called chemical analogy (i.e., assuming the same chemical behaviour for elements belonging to the same group).

---

1. Introduction

Phosphorus is an element typically associated with refractory materials, such as phosphates, which are found in numerous minerals.¹ However, phosphorus may also be present in volatile molecules such as phosphine and its halogenated derivatives. Following the detection of phosphine in circumstellar envelopes in 2014,² the gas-phase chemistry of phosphorus-bearing compounds has attracted growing interest over the last decade.³

Indeed, the detection of PH₃ in such environments challenges astrochemical models, as no routes for its formation in the gas phase are currently known. The presence of phosphine in the atmospheres of Jupiter^4,5 and Saturn^6,7 is believed to be associated with photochemical formation routes, which also in this case are poorly understood. For this reason, it is sometimes assumed that the presence of phosphine (PH₃) may be associated with the existence of biological activity.^8,9 However, the great number of open questions about the chemistry of phosphorus in the gas phase,¹⁰ in particular the lack of knowledge on the chemistry of its hydrides, makes it very difficult to draw any conclusions on the role of phosphine as a biomarker.

To shed some light on the chemistry of phosphorus hydrides, in this work, we adopt a combined experimental and theoretical approach, investigating the possible formation of PD⁺ through the collision of the phosphorus cation in its ground state P⁺(³P) with molecular deuterium D₂, which is expected to behave in a similar way to H₂. To the best of our knowledge, limited data on this reaction are available in the literature. The first experimental study of this reactive system was carried out by Smith and coworkers through an ion flow tube setup,¹¹ where they measured the three-body association reaction P⁺ + H₂ + M → PH₂⁺. However, during the experiments they failed to observe PH⁺, suggesting a high endothermicity for such a channel. This was confirmed by Armentrout's preliminary results on the P⁺(³P) + H₂ → PH⁺ + H reaction obtained with a guided ion beam experiment,¹² although the definitive findings of this study were never published. Additionally, a recent theoretical study of the reaction was performed by Li and coworkers, by studying the reaction in the ground triplet state of the PH₂⁺ system.¹³

The manuscript is structured as follows. In Section 2, the experimental and theoretical methods are presented. Firstly, the experimental setup, the generation of the phosphorus cation and the procedure to measure the absolute reactive cross-section are described. Then, the theoretical approach to simulate the experimental results from the first principles will be presented. This includes the methodology to construct the potential energy surfaces (PESs) on which the reaction may occur and the details of how the reaction cross-section can be computed by probing the reaction dynamics. In Section 3, the theoretical and experimental results obtained are presented, compared and discussed, leading to the conclusions.

2. Methods

2.1. Experimental

The ion–molecule reaction P⁺(³P) + D₂ was studied using the Guided Ion Beam setup CERISES,^14,15 coupled to the DESIRS beamline¹⁶ at the French synchrotron SOLEIL. The setup is based on a sequence of four radio frequency (RF) devices, with a quadrupole–octopole–octopole–quadrupole (QOOQ) configuration. The system allows mass selection of the parent cations produced in the source (Quad 1), and the control of the collision energy between the ions colliding with the target gas (Oct1). Once the product ions (PD⁺) are obtained, they are mass-selected together with the parent ions (P⁺) before their detection (Quad 2). Absolute reaction cross-sections are extracted from the measured ion yields and the absolute pressure of the target gas, which is set so that the single-collision regime is ensured. The effective cell length is determined through the calibration reaction Ar⁺ + D₂ → ArD⁺ + D.¹⁷ Additional technical details on the experiment are reported in previous studies on the reactivity of the sulphur cation in its electronic ground (⁴S) or excited (²D) state with H₂,^18,19 and on the reactivity of the ethenium cations with ethyne.²⁰

The phosphorus cations (P⁺) are generated via dissociative photoionization of the PCl₃ precursor. Fig. 1 shows the ion yield as a function of the photon energy together with results from previous experiments by Au and Brion.²¹ Both datasets exhibit similar trends, but the higher energy resolution available at SOLEIL allows us to provide more accurate values for the appearance energies (AE) of the ion. The appearance energy of P⁺ is measured at 15.5 ± 0.1 eV photon energy and is associated with the generation of P⁺ with a low yield up to around 20.5 eV. Above this photon energy value, a significant increase of the ion yield is observed, reaching a maximum at 23.6 eV. Au and Brion assigned the first threshold to the dissociation leading to pathway (1), while the second threshold is ascribed to pathway (2).

PCl₃ + hν → P⁺(³P) + Cl₂ + Cl + e⁻ (AE = 18.02 eV) (1)

PCl₃ + hν → P⁺(³P) + 3Cl + e⁻ (AE = 20.53 eV) (2)

While the second threshold is in good agreement with its respective calculated AE of 20.53 eV (from thermodynamics tables), the first AE of P⁺ is considerably lower than the calculated value. To explain the discrepancy of the first AE, one should consider the possibility that the dissociative ionisation of PCl₃ might proceed via ion pair formation, according to the following pathways (indicated as vertical green bars in Fig. 1):

PCl₃ + hν → P⁺(³P) + Cl⁻ + Cl₂ (AE = 14.40 eV) (3)

PCl₃ + hν → P⁺(³P) + Cl₂⁻ + Cl (AE = 15.52 eV) (4)

PCl₃ + hν → P⁺(³P) + Cl⁻ + 2Cl (AE = 16.92 eV) (5)

Fig. 1 Photoion yields at m/z 31 (P⁺ products) from the dissociative photoionisation of PCl₃. Our results are represented by red triangles and the data from Au and Brion²¹ are represented by the blue circles. The vertical solid lines indicate the AE thresholds for the corresponding channels, as calculated from available thermochemical values (NIST and Active TC Tables). The vertical dashed lines represent the AE thresholds when considering P⁺(¹D) instead of P⁺(³P) for each channel. Green lines correspond to ion pair formation thresholds (pathways 3–5) while purple line correspond to dissociative ionization (pathways 1–2). The vertical dark blue arrow represents the photon energy at which data as a function of collision energy are collected.

The fragmentation pattern that matches our measured AE is described by (4), although the formation of Cl₂⁻ has not been reported in the electron ionization mass spectrometry experiment performed on PCl₃ by Halmann and Klein,²² while Cl⁻ has been observed. However, the P⁺ ion yield obtained from PCl₃ for photon energies below 20.5 eV was not high enough to carry out the experiments. In order to maximize the production of P⁺ and gain sensitivity for the study of its reactivity with a light target such as D₂, we therefore set the photon energy at 22.6 eV. Although at this photon energy the ion yield is high enough to perform the experiments, the generation of the phosphorus cation in its first electronic excited state is energetically open, and a small amount of P⁺(¹D), which cannot be discriminated from P⁺(³P), might contaminate the ion beam.

2.2. Theory

The reaction P⁺ + D₂ is studied theoretically from first principles within the Born–Oppenheimer approximation. P⁺(³P) is an atomic radical cation composed of three degenerate electronic states when neglecting the spin–orbit coupling. Since each of the three electronic states might interact in a different way with the D₂ molecule, it is necessary to consider the reaction dynamics on the three related PESs on which the reaction may take place. In this work, as a first approximation, no couplings between electronic states are considered, with the total reactivity given as the average of contributions of the three electronic states of P⁺.

2.2.1. Potential energy surfaces. High-level electronic structure calculations were performed over the whole configuration space of the PH₂⁺ system to make a global sampling of the three PESs correlating asymptotically to P⁺(³P) + H₂: 1³A′, 1³A″ and 2³A″. All calculations were performed with the internally contracted MRCI method²³ available within the MOLPRO suite of programs,²⁴ choosing aug-cc-pV5Z as a basis set.²⁵ The MRCI wavefunction is built over orbitals generated with the state-average CASSCF method²⁶ in which the three triplet states of interest to this work have been included, as well as the five singlet states correlating to P⁺(¹D) + H₂: 1¹A′, 2¹A′, 3¹A′, 1¹A″ and 2¹A″. By considering these eight electronic states in the state-average and an active space including the 3s and 3p valence orbitals of phosphorus together with the two 1s orbitals of both hydrogen atoms, the active space has been found to remain stable throughout the entire configuration space. Then, the three triplet states were computed at the MRCI level using the same active space, with the core orbitals of phosphorus (1s, 2s and 2p) being kept frozen. To minimize the size-consistency error inherent in the MRCI method, the Davidson correction²⁷ has been applied.

Jacobi coordinates of the reactants P⁺–H₂ and products H–PH⁺ were used to sample the geometries in the complex region as well as in both reactant and product channels. All the points are included within the PH internuclear distances between 2 and 24 bohr and H–H distances between 0.8 and 20 bohr. Then, a total of 20 000 geometries computed for the three triplet states have been considered for the fitting. The corresponding ab initio energies of the 1³A′, 1³A″ and 2³A″ electronic states have been fitted separately using the many-body expansion where the total energy is decomposed in a sum of one, two and three body terms:

where the one-body term V⁽¹⁾_A, represents the energies of the three separated atoms, the two-body term V⁽²⁾_D(r), provides the energies of the three diatomics constituting the system as a function of their internuclear distances, and the three-body term V⁽³⁾(r₁, r₂, r₃), depends on the three internuclear distances. With this procedure, the three-body term tends asymptotically to zero, providing a well-behaved extrapolation in the long-range region.

The two-body terms were fitted using the RKHS method,²⁸ in which the diatomic potential is described by the following expansion:

where q^n,m₁ is a distance-like kernel, N is the number of data points considered for the fitting and α_i are the optimized coefficients. To ensure an r⁻⁶ long-range smooth behaviour, the parameters n = 2 and m = 5 have been chosen for H₂, while n = 2 and m = 3 have been selected for PH⁺ to accurately reproduce the r⁻⁴ ion-induced dipole long-range interaction. More details on how the procedure is applied may be found elsewhere.^29–32

The three-body term, extracted from the MRCI energies by subtracting the one- and two-body contributions, is then fitted using the invariant polynomial procedure GFIT3C introduced in ref. 33–35, as it is well suited to fit triatomic systems with two equivalent atoms, as it was shown for SH₂⁺,^36–38 OH₂,^39,40 LiH₂⁺ (ref. 41) or AuH₂.^42,43 Briefly, within this procedure, the three-body term is expressed as a combination of modified Rydberg functions ρ_r:⁴⁴

where d_ijk are the optimized linear parameters and ρ_r = re^−βr with β > 0. To keep the invariability of the energy to the permutation of two equivalent atoms, the parameter β is set to be equal for the two PH⁺ diatomics. The indices ijk are constrained by i + j + k ≤ L, where L is the maximum order of the polynomial. In this work, we considered L = 10 for the three PESs.

2.2.2. Reaction dynamics. The reaction dynamics studies are performed on each of the three PESs using the quasi-classical trajectory (QCT) method as described by Karplus,⁴⁵ using the local code of the group, MDwQT.^37,46–48 Since the method is extensively described in other studies, we only report the parameters used to investigate the system of interest in this work. With the temperature of the chamber containing the target gas D₂ at about 300 K, the QCT simulations were performed considering the vibrational and rotational states for D₂ at v = 0 and j = 0–6. Having the fitted PESs provided by analytical derivatives, batches of 1 million trajectories were run for each rotational level and collision energy value, the latter scanned through a range from the threshold energy up to 2.5 eV, with a 0.1 eV step.

For each of the three PESs, the cross-sections as a function rotational level, j, and collision energy, E_kin, are computed as:

σ^S(E_kin, j) = πb_max²P_r(E_kin, j) (9)

where S refers to the electronic state, b_max is the maximum value of the impact parameter leading to a reactive collision and P_r = N_r/N, with N_r being the number of reactive collisions and N being the total number of trajectories with impact parameters ≤b_max. In this work, the standard histogram binning method is applied. The total cross-section as a function of energy is then recovered for each individual PES by a weighted average of the specific cross-sections and assuming a Boltzmann distribution of the rotational levels of D₂ at a temperature of 300 K. Additionally, since P⁺(³P) is produced by the dissociative photoionization of PCl₃, we assume that the three degenerate levels of P⁺(³P) are equally populated. Within this assumption, the integral cross-section is obtained as:

3. Results and discussion

3.1. PD₂⁺ triplet PESs

The three first triplet states of the PH₂⁺ triatomic system exhibit a relatively complex electronic structure. In particular, the two A″ states are strongly coupled via two conical intersections (CIs), one located in the P⁺ + D₂ entrance channel with a T-shape configuration (C_2v) and the other in the PD⁺ + D exit channel with a linear configuration (C_∞v), associated with the Σ–Π crossing. The positions of the two CIs obtained at the MRCI level are indicated in Table 1. As a consequence, the presence of these two CIs shapes the adiabatic PESs of the two ³A″ states, leading to relatively rich topographical features for a triatomic molecule in its triplet states. With respect to their asymptotic behaviour, both the 1³A′ and 1³A″ states correlate with the ground states of the reactants and products, P⁺(³P) + D₂(X¹Σ⁺_g) and PD⁺(X²Π) + D(²S), respectively. The 2³A″ state, degenerated with the other two in the reactant channel, correlates with the first excited state of the products, PD⁺(a⁴Σ⁻) + D(²S).

Table 1 Positions and energies of the two conical intersections present in the P⁺ + H₂ entrance channel (C_2v) and the PH⁺ + H exit channel (C_∞v)

	CIC2v	CIC∞v
E (eV)	0.007	1.381
R PH1 (bohr)	2.668	2.700
R PH2 (bohr)	2.668	7.015
R H1H2 (bohr)	3.020	4.315
	55.5	0

---

Except for the CIs, the fitted PESs reproduce very well the main topographical features of the ab initio points, as suggested by the small RMS error of the fits. The fit of the 1³A′ state present an RMS error of 17.5 meV, with a maximum deviation of 129.9 meV between the ab initio points and the fit. For the fit of the 1³A″ state, the RMS error is 23.4 meV, with a maximum deviation of 144.8 meV, while the RMS error of the fit of the 2³A″ state is 28.6 meV, with a maximum deviation of 277.0 meV. In all the cases, the regions affected by the highest uncertainty are located in the repulsive regions of the respective potentials and in the close proximity of the CIs. To assess the quality of the PESs and the ab initio calculations, the analysis of the asymptotic behaviour of the fitted PESs is insightful, also allowing us to probe the representations of the reactants and diatomic products. The derived spectroscopic constants of the asymptotic D₂ and PD⁺ are presented in Table 2. In order to make a comparison with previous data available in the literature, the constants associated with H₂ and PH⁺ are also reported. For both diatomics, the agreement obtained with previous experimental and theoretical data is good, showing that the asymptotic behaviour of our PESs is well described.

Table 2 Spectroscopic constants for PH⁺(X,a), H₂(X), PD⁺(X,a) and D₂(X) derived from the asymptotic diatomic molecules. Data are compared with values available in the literature. Values of D_e marked with (*) are obtained as a first approximation by summing half of the harmonic frequencies from the available D₀

State	R e (Å)	D e (cm^−1)	B e (cm^−1)	ω e (cm^−1)	ω e x e (cm^−1)
a This work. b Ref. 49. c Ref. 13. d Ref. 50. e Ref. 51.
PH^+
X^2Π^a	1.4260	3.558	8.4857	2380	50.32
theo.^b	1.4226	3.525	8.5369	2413	44.38
theo.^c	1.4270	3.415	8.4824	2386	48.97
exp.^d	—	3.496–3.517(*)	8.5080	2383	41.67
H&H^e	1.4352	3.503(*)	8.3851	2291	—
a^4Σ^− ^a	1.4799	1.881	7.8853	1892	60.73
theo.^b	1.4816	1.791	7.8819	1833	58.12
H2
X^1Σ^+g ^a	0.7426	4.696	60.651	4387	126.04
theo.^c	0.7415	4.723	60.834	4397	107.12
H&H^e	0.7414	4.751(*)	60.853	4401.21	121.34
PD^+
X^2Π^a	1.4265	3.558	4.3802	1710	25.98
H&H^e	1.4314	—	4.3505	1666	—
a^4Σ^− ^a	1.4799	1.881	4.0703	1360	31.38
D2
X^1Σ^+g ^a	0.7426	4.696	30.352	3103	63.06
H&H^e	0.7415	4.749(*)	30.444	3115	61.82

---

Since the fits well reproduce both the energies obtained from the electronic structure calculations, and the diatomic molecules in the asymptotic regions, the characteristics of the PESs of the different electronic states are taken from the features of the analytical PESs. The most relevant among such properties are summarized in Fig. 2. From the spectroscopic constants of the diatomic molecules, it can be determined that P⁺ + D₂ is an endothermic reaction, with a calculated threshold of 1.127 eV with respect to the zero-point energies (ZPE) of the reactants and products. The calculated second threshold, leading to the excited state of products is located at 2.762 eV.

Fig. 2 Adiabatic correlation diagram of the three first triplet states of PH₂⁺ according to the analytical PESs used in this work. The zero-point energy (ZPE) of D₂ and PD⁺ are also shown in parenthesis. The origin of energies is set to the P⁺(³P) + D₂ channel asymptote. The stationary points associated with HPH configurations are represented in blue (1³A′), red (1³A″) and brown (2³A″), while the stationary points associated with HHP configurations are not shown on this diagram and are only reported in Table 3.

The representative contour plots of the three PESs are shown in Fig. 3 to give an overview of their topography. All the corresponding stationary points are detailed in Table 3. The 1³A″ state is the most stable. Its global minimum, MIN1_1A″, lies at −1.882 eV below the entrance channel and correlates directly with the product channel without any additional barriers. In terms of the geometry of the system, MIN1_1A″ represents an insertion complex where the phosphorus center lies between the two deuterium atoms. The second well, MIN2_1A″, lies at −1.044 eV and is accessible from the entrance channel also through a barrierless path. MIN1_1A″ and MIN2_1A″ are separated by a submerged barrier, SP1_1A″, that lies below the entrance channel at −0.155 eV. It is worth noting that SP1_1A″ should correspond to the conical intersection CI_C2v that has been smoothed by the fitting procedure, thus appearing as a saddle-point in our PES. The other conical intersection, CI_C∞v, is also represented by a saddle-point, SP2_1A″, lying at 1.219 eV, slightly above the PD⁺ + D asymptotic energy (1.216 eV). Hence, SP2_1A″ acts as a tiny barrier in the exit channel, but since it lies only 3 meV above the asymptotic channel, it is not expected to play a significant role in the P⁺ + D₂ reaction. Additionally, the PES presents two more saddle-points, SP3_1A″ and SP4_1A″, not shown in Fig. 2 for clarity. SP3_1A″ corresponds to the bending inversion barrier associated with MIN1_1A″ and appears for a linear HPH configuration. This barrier lies at 1.520 eV above the entrance channel, and 3.40 eV above MIN1_1A″, and is not expected to play a significant role in the reaction dynamics. SP4_1A″ corresponds to the inversion barrier for MIN2_1A″, but unlike SP3_1A″, this saddle-point is not associated with any bending motions. Rather, it is related to the orbital motion of P⁺ around D₂, and while the well MIN2_1A″ is found for a T-shape, SP4_1A″ is linear with the HHP configuration and lies at −0.414 eV. It is noteworthy that like MIN2_1A″, SP4_1A″ correlates directly with the entrance channel without additional barriers. This aspect is quite relevant for the reaction dynamics, because starting from the entrance channel, P⁺ is subject to a relatively strong attraction to D₂, regardless its orientation. Considering such features, we might expect that reactive collisions would be possible on the 1³A″ state as soon as the collision energy is sufficient to overcome the tiny barrier in the exit channel, and that approximately corresponds to the reaction endothermicity.

Fig. 3 Contour plots of the three PESs as a function of the Jacobi distances r and R (where r denotes the D₂ internuclear distance and R denotes the distance between P⁺ and the center of mass of D₂) for the Jacobi angles of γ = 0°, 45° and 90°.

Table 3 Energies and geometries for the stationary points of the three fitted triplet electronic states of PH₂⁺

	E (eV)	R PH1 (bohr)	R PH2 (bohr)	R H1H2 (bohr)
1^3A′
SP11A′	1.258	2.825	5.652	2.827	0	180
SP21A′	1.520	1.572	1.572	3.144	180	0
1^3A″
MIN11A″	−1.882	2.671	2.671	4.696	123.0	28.5
MIN21A″	−1.044	3.258	3.258	1.595	28.3	75.8
SP11A″	−0.158	3.003	3.003	2.941	58.6	60.7
SP21A″	1.219	2.722	6.138	3.416	0	180
SP31A″	1.520	2.972	2.972	5.944	180	0
SP41A″	−0.414	3.110	4.609	1.499	0	180
2^3A″
MIN2A″	0.288	2.862	2.862	3.272	69.7	55.1
SP12A″	2.799	3.181	3.181	2.193	40.3	69.8
SP22A″	1.575	2.941	2.941	5.882	180	0

---

The 2³A″ state is the highest in energy. This state shows a metastable well, MIN_²A″, lying at 0.288 eV above the reactant channel, which can only be accessed by overcoming a remarkable barrier, SP1_²A″, of 2.799 eV. The well MIN_²A″ can be considered as the counterpart of CI_C2v, with a geometry resembling SP1_1A″. Due to the height of the barrier, it is nevertheless unlikely that this state may play an important role in the reaction dynamics, therefore even a non-rigorous treatment of CI_C2v is expected to be sufficient to get an accurate description of the reaction. Furthermore, the 2³A″ state correlates with the excited state of the product channel (PD⁺(⁴Σ⁻) + D(²S), lying at 2.875 eV), with no additional barriers, and it is also linked to a bending inversion saddle-point found in a linear HPH configuration. Finally, this state also includes a barrierless pathway directly connecting the reactants and products for PHH collinear configurations, although this state is expected to contribute minimally to the reaction due to the very large endothermicity involved.

An additional electronic state correlating reactants and products in their ground state is 1³A′. Unlike the two other states, it does not present any stable or metastable geometries outside the purely electrostatic wells arising from van der Waals interaction, and therefore their characterization is out of the scope of this work. The optimal reaction path on this PES is found for a direct PHH collinear mechanism, with only a small barrier in the exit channel, similar to the 1³A″ state. This is not surprising though, since for linear configurations this electronic state has a Π character, and becomes degenerate with 1³A″ in the exit channel and with 2³A″ in the entrance channel. The transition of degeneracies occurs at the conical intersection CI_C∞v, which, interestingly, coincides with the position of the barrier. In this region, strong couplings are expected as well as the Renner–Teller effect, which suggests that neglecting the correct description of CI_C∞v in our fully adiabatic approach may impact the reaction dynamics results. The acceptance cone of this pathway is, however, much lower on this state as compared to its homologue on the 1³A″ PES, as it suffers a sharper increase of the potential when linearity is broken. With only one pathway leading to the products, together with the reduction of the acceptance cone, this state is expected to contribute less to the reaction than the 1³A″ state, but at the same time, it could contribute more than the 2³A″ state because it is less endothermic.

To illustrate the different reaction paths of the three PESs, and see how the presence of CI_C∞v shapes the different PESs, the global minimum energy pathway and additional reaction routes determined at fixed angles are shown in Fig. 4. For the collinear reaction path, we clearly see how the degeneracy of the 1³A′ state with the two ³A″ states switches at the top of the barrier and how the avoided crossing increases when the angle starts to bend. It can also be appreciated how well the degeneracy is reproduced by the different fits, which were obtained separately. This can be interpreted as a further confirmation of the quality of the calculated PESs.

Fig. 4 Minimum energy path of the three electronic states (top left panel) and several reaction paths obtained at fixed θ values (with ). The reaction coordinate is here defined as the difference between the internuclear distances of H₂ and PH.

3.2. P⁺(³P) + D₂ reaction dynamics

The absolute cross-section σ^exp of the reaction P⁺ + D₂ → PD⁺ + D is measured with the CERISES setup for translational energies up to 20 eV in the laboratory frame, which corresponds approximately to 2.3 eV in the center-or-mass frame. The results are shown in the top panel of Fig. 5. It can be seen that at low collision energies, the measured cross-section is close to zero, and then increases monotonically as the collision energy increases. A non-negligible cross-section is obtained from around 1 eV to 2.3 eV. The fact that the cross-section is slightly above zero at low collision energies confirms the presence of P⁺(¹D) in the ion beam. However, the magnitude of the cross-section in this range, which remains apparently constant from 10 meV up to ≈0.95 eV, is ≈0.15 Ų, suggesting that the amount of excited state should be almost negligible. As a first approximation we assume that, according to the observed cross-section at low collision energies, the contribution of P⁺(¹D) to the reaction is constant over the whole energy range. An energy independent cross-section is rather unusual for a reaction involving an excited atom with molecular hydrogen, as generally, a Langevin behavior is expected, as in the case of S⁺ (²D) + H₂.¹⁹ However, the case of P⁺ (¹D) + D₂ → PD⁺ + H is slightly different. Unlike the sulfur cation, the reaction involving the excited P⁺ is slightly endothermic. A cross-section showing a threshold is thus expected, at least for the lower rotational levels. However, as mentioned previously, D₂ is at 300 K in the experiment, and rotational levels up to j = 6 may be populated. For higher rotational levels, the reaction becomes exothermic, and in this case, the cross-section is expected to follow Langevin behavior, assuming that there are no barriers in the entrance channel. The Boltzmann average of contributions showing threshold behavior and Langevin behavior may lead to a total cross-section nearly energy independent, as it is apparently observed. A detailed study of the reaction dynamics of P⁺(¹D) + D₂ would be necessary to confirm this hypothesis and provide more accurate results, but meanwhile, we have to settle for this approximation. The cross-section associated with P⁺(³P) is thus recovered shifting the entire cross-section and setting the base line to zero and it is shown in the bottom panel of Fig. 5, which includes a comparison with the calculated cross-sections on the three PESs obtained using the QCT method.

Fig. 5 (Top) Experimental absolute cross-section of P⁺ + D₂ → PD⁺ + D, where P⁺ is produced by the dissociative photoionization of PCl₃ at a photon energy of 22.6 eV. (Bottom) QCT absolute cross-section of P⁺(³P) + D₂ → PD⁺ + D and the cross-section calculated over each individual PES compared with the shifted experimental cross-section (see the text for details). The blue vertical line in the bottom panel represents the threshold of the reaction for D₂ (v = 0 and j = 6).

The experimental cross-section compares qualitatively well with the preliminary results published by Armentrout in Fig. 6 of ref. 12, with a cross-section slightly below 2 Ų at 2 eV collision energy. Due to the contamination of P⁺(¹D) and the approximation adopted to subtract it, the onset of the reaction is not perfectly defined in this experiment, but the reaction threshold is observed at a collision energy in the 0.8–0.9 eV range. Considering the dispersion in the energy of the ion beam and the Doppler broadening associated with the temperature of the chamber (300 K), this onset is compatible with the theoretical predictions, with an expected threshold of 0.974 eV when considering a rotational energy of 0.153 eV for D₂ (v = 0; j = 6), the highest level populated at 300 K, and the global energy profile provided by our analytical PESs.

It is insightful looking at the contributions of the different electronic states that are also plotted in the bottom panel of Fig. 5. It appears that up to a collision energy of 2 eV all the reactivity is provided by the 1³A″ state. The 2³A″ state remains unreactive up to 2.9 eV, consistent with its endothermic nature and its correlation with the ⁴Σ⁻ state of PH⁺, which lies higher in energy. In contrast, the unexpectedly low reactivity of the 1³A′ state, becoming reactive only slightly above 2 eV, is noteworthy. This behaviour is surprising, as the 1³A′ state correlates with the ground state of the products and, both SP1_1A′ and SP2_1A′ lie at considerably lower energies, 1.258 eV and 1.520 eV, respectively. This behaviour can be attributed to a dynamic effect arising from the inability to form a complex on this surface and from the large difference in the mass between P (31 amu) and D₂ (4 amu). When P⁺ approaches D₂ in a collinear geometry, the energetically favoured path for this state, the translational energy is directly transferred to the center of mass rather than into the D₂ bond. A minimum incident angle seems therefore necessary to transfer the energy into the bond, a necessary step to break it and be able to form the hydride bond. To get an idea of the incident angle required to produce the reaction, we represent in Fig. 6, the bending potential as a function of the Jacobi angle associated with SP1_1A′, as it provides a graphical view of the expected acceptance cone. By examining the intersection between the bending potential and the product appearance threshold energy given by the QCT simulation, it seems that the reaction can effectively take place only for incident angles larger than 48°. In contrast, on the 1³A″ surface insertion complex formation is possible, enabling the spontaneous breaking of the D₂ bond and a more efficient redistribution of energy, thereby allowing reactivity as soon as the collision energy becomes larger than the potential threshold to overcome.

Fig. 6 Potential of the 1A′ PES as a function of the Jacobi angle for r = 2.63 bohr and R = 4.14, corresponding to the position of SP1_1A′ in the Jacobi coordinates. The QCT threshold energy for which the reaction opens is marked by the red dashed line.

Since most of the reactivity is associated with the 1³A″ state, it is interesting to compare our results on this state with those of Li and coworkers, who recently studied theoretically the reaction dynamics of P⁺ + D₂ in this state.¹³ Globally, the results are very similar and both simulations agree on a calculated maximum cross-section value slightly above 4 Ų for a collision energy of 2.8 eV. A small fluctuation appears in the cross-section calculated by Li and coworkers between 2.8 eV and 4 eV, which is not observed in our case. This may be due to the fact that our cross-section is averaged over seven rotational states of D₂, while only D₂ (j = 0) was considered in their work. Another source of discrepancy may be related to the differences in the PESs used to simulate the reactive system between this work and that of Li and coworkers.

When comparing theory and experiment, it appears that the shapes of the experimental and the average QCT cross-sections (obtained from eqn (10)), presented in the bottom panel of Fig. 5, are very similar, although the experimental values lies clearly above the QCT results, which are underestimated by about 25%. This small discrepancy between theory and experiment may arise in part from the procedure applied to subtract the P⁺(¹D) contribution. Even so, as previously mentioned, the amount of P⁺(¹D) is expected to be small and the differences appear to be too large to be explained only by a problem of contamination by the excited state. It is thus more likely that the underestimation of the QCT cross-section is attributable to the adiabatic approximation applied in this work, suggesting that non-adiabatic and spin–orbit couplings may also play a non-negligible role in this reaction. In fact, as previously discussed, one of the conical intersection, CI_C∞v, may contribute to the reaction dynamics, as couplings can be very strong in its vicinity. Indeed, the results of this work show that, adiabatically, most of the reactivity is provided by the 1³A″ PES, while the initial flux is equally distributed among the three PESs (degenerate in the entrance channel). By coupling the 2³A″ PES to the 1³A″ PES (we remind that the 2³A″ state correlates with the excited state of the product, PD⁺ (a⁴Σ⁻) and it is therefore not reactive at lower energies), part of its flux may be transferred to the more reactive PES, enhancing the overall reactivity. In addition, the ³A′ state may also couple to both ³A″ states through Renner–Teller couplings, which may also be strong near colinear geometries where the Π character is well marked, and these couplings will therefore be efficient near CI_C∞v. Furthermore, nearby the threshold, around 1.1–1.2 eV, a tiny shoulder that resembles the one observed in the S⁺(⁴S) + H₂ reaction¹⁸ is observed, and may indicate a possible contribution of a spin–forbidden route. In the case of S⁺, the shoulder was attributed to a mechanism where part of the initial flux of the repulsive quartet state, of relatively low reactivity near the threshold, could be transferred to the doublet ground state of H₂S⁺ by intersystem crossing. This mechanism could be associated with a maximal contribution to the total cross-section of around 0.4 Ų. As the doublet ground state also correlates with the ground state products, but to an excited state of reactants, the probability to reach the products is greatly enhanced on this state. A larger effect is likely to occur in the P⁺ + D₂ reaction since two singlet states of D₂P⁺, one of them being the ground state, also correlates with the ground state of PD⁺. So several routes associated with the couplings between 1¹A′, 1¹A″, 1³A′ and 1³A″ states, which are all degenerate in the asymptotic region, may contribute to enhance the cross-section. Assuming that the sum of these contributions are likely to be larger than 0.4 Ų, and that the adiabatic contributions are estimated to be around 2 Ų, it is very plausible that spin–orbit and non adiabatic couplings may enhance the reactivity by 25% in this reaction. However, a much more complex model would be required to confirm this hypothesis. Nevertheless, since the mechanism of the dissociative photoionization of PCl₃ to produce P⁺(³P) is not known, we cannot be certain that the fine-structure levels of P⁺(³P) are equally populated as assumed in this work. So we cannot exclude that part of the discrepancy may also arise from the choice of the assumed partition function.

It is also interesting to compare the reactivity of P⁺ and molecular hydrogen with its neighbours in the periodic table. Phosphorus lies just below nitrogen in the 15th group, and since they share the same valence electronic configuration, at first glance one could expect a similar reactivity of P⁺ and N⁺. However, the much smaller reaction endothermicity in the N⁺ case compared to P⁺ (∼17 meV vs. 1.024 eV) leads to striking differences in the reactivity of nitrogen and phosphorus cations toward molecular hydrogen, as observed when comparing the cross-sections here reported with those of N⁺(³P) + H₂ and isotopologues.^52,53 While the maximum cross-section is of the order of 2 Ų for P⁺ + D₂ and appears at rather high collision energies, the cross-section of N⁺ + D₂ is much larger, reaching values higher than 50 Ų. Its trend with collision energies is also completely different, as it behaves as a capture process so that the highest values of the cross-section are found at very low collision energies. On the other hand, the reactivity of P⁺ with H₂ is found to resemble that of its neighbours along the third row of the periodic table, namely the sulphur and silicon cations. In each case, the reaction is markedly endothermic and characterized by relatively small cross-sections toward molecular hydrogen. With S⁺, the maximum cross-section is ≈1 Ų,^18,54 while in the case of Si⁺ (ref. 55), it hardly reaches ≈1.5 Ų. Since the bond type is, to some extent, similar for atoms along the same group, it is often assumed that their reactivity would be similar too (i.e., grouping the reactivity by the chemical analogy). Here, we clearly see that this is not the case. This is mainly ascribable to the difference in the involved bonding energetics, which may differ considerably between the second- and third-row atoms.

4. Conclusions

In this work, we present a combined experimental and theoretical study to explore the reactivity between the phosphorus cation in its ground state, P⁺(³P), with deuterium, D₂. The measured cross-section is simulated theoretically from first principles, by running quasi-classical trajectories over the three potential energy surfaces on which the reaction may take place. Despite a slight underestimation of the experimental reactive cross-section, attributed to a minor contamination from electronically excited P⁺ and to limitations of the Born–Oppenheimer approximation, the experimental and theoretical results are in good agreement. They confirm that the investigated reaction is highly endothermic, requiring large collision energies to proceed, and thus exhibits low overall reactivity. This low reactivity of P⁺ with hydrogen contrasts sharply with that of its second-row homologue, N⁺, which exhibits much higher reactivity. In contrast, the behaviour of P⁺ more closely resembles that of its third-row neighbours, Si⁺ and S⁺, which also react inefficiently with H₂.

From a mechanistic point of view, our results show that the reactivity arises predominantly from one of the three electronic states, the 1³A″ one, with only minor contributions from the other two. While the limited reactivity of the 2³A″ state is expected, as it correlates with an excited state of the products, the very small contribution from the 1³A′ state is more surprising. Indeed, the reaction threshold for this state is found to occur at a much higher collision energy than what would be anticipated from the topology of its adiabatic PES. This behaviour is attributed to the overall repulsive character of the PES, which hinders efficient access to the product channel. However, the 1³A′ and 2³A″ electronic states may contribute more to the overall reactivity than predicted in this work, and the rationale is twofold. First, the discrepancy between the experimental and theoretical cross-sections is larger than the estimated contribution from the contamination by P⁺(¹D) to the measured signal. Second, reaction dynamics calculations were performed under the assumption of no coupling between the electronic states. Although reactive collisions are not favoured on the 1³A′ and 2³A″ states under strictly adiabatic conditions, the inclusion of non-adiabatic and spin–orbit couplings may enable population transfer to other electronic states, where access to the product channel may occur more easily. The fact that 75% of the experimental cross-section is reproduced by the present approach, together with a crude estimate that a maximum of 10% of the cross-section may arise from P⁺(¹D) contamination that was not properly accounted for, suggests that over 15% of the cross-section may come from mechanisms involving non-adiabatic and spin–orbit effects. To achieve a comprehensive understanding of the reaction, additional model refinements are necessary and further studies will be required to assess the impact on the reaction of the various couplings neglected in the present theoretical treatment, and to refine the estimation of the contribution from singlet states, which was not considered in this work.

Author contributions

JAF: data curation – theory (lead) and experiment (equal); formal analysis – theory (lead); writing – review & editing (equal). MM: data curation – experiment (lead); formal analysis – experiment (lead); writing – review & editing (equal). CSS: data curation – experiment (equal); supervision – theory (equal); formal analysis – theory (supporting); writing – review & editing (equal). SGC: data curation – experiment (equal); supervision – theory (equal); formal analysis – theory (supporting); writing – review & editing (equal). CA: data curation – experiment (equal); resources (equal); writing – review & editing (equal); conceptualization (supporting); project administration (supporting). RT: data curation – experiment (equal); resources (equal); writing – review & editing (equal); conceptualization (supporting); project administration (supporting). NS: data curation – experiment (equal); resources (equal); writing – review & editing (equal). MP: data curation – experiment (equal); writing – review & editing (equal). DA: data curation – experiment (equal); supervision – experiment (lead); formal analysis – experiment (supporting); writing – original draft (supporting); writing – review & editing (equal); conceptualization (equal); project administration (lead). AZ: data curation – experiment (equal); supervision – theory (equal); formal analysis – theory (supporting); writing – original draft (lead) and review & editing (equal); conceptualization (equal); project administration (supporting).

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including the experimental raw data and the code to analyse them, as well as the data corresponding to the specific cross-sections for D₂ (j = 0–6) for each PES are available at the Digital CSIC repository at https://doi.org/10.20350/digitalCSIC/17915.

Acknowledgements

We acknowledge SOLEIL for providing synchrotron radiation necessary for the experiments and we would like to thank L. Nahon for assistance in using the DESIRS beamline under proposal 20231311. The research leading to these results has received funding from MICINN (Spain) under grant no. PID2021-122549NB-C21, PID2021-122549NB-C22 and PID2024-155352NB-C22. We also acknowledge the computing time on HPC DRAGO (CSIC). The COST Action CA21101 (COSY) is also acknowledged for its support and facilitating collaborations. MM acknowledges the Department of Physics at the University of Trento for PhD scholarship funding, and the COST Action CA22133 – PLANETS for the financial support while carrying out the experiments at the SOLEIL synchrotron. DA acknowledges financial support from the MUR PRIN 2020 project no. 2020AFB3FX “Astrochemistry beyond the second period elements”. MP acknowledges the COST Action CA21126 – Carbon Molecular Nanostructures in Space (NanoSpace) for financial support while conducting the experiments at the SOLEIL synchrotron.

Notes and references

1. E. Macia, Chem. Soc. Rev., 2005, 34, 691.
2. M. Agundez, J. Cernicharo, L. Decin and D. Teyssier, Astrophys. J. Lett., 2014, 790, L27.
3. F. Fontani, Front. Astron. Space Sci., 2024, 11, 1451127.
4. E. W. Weisstein and E. Serabyn, Icarus, 1994, 109, 367.
5. C. Gapp, M. Rengel, P. Hartogh, H. Sagawa, H. Feuchtgruber, E. Lellouch and G. L. Villanueva, Astron. Astrophys., 2024, 688, A10.
6. J. A. Kaye and D. F. Strobel, Icarus, 1984, 59, 314.
7. L. N. Fletcher, G. S. Orton, N. A. Teanby and P. G. J. Irwin, Icarus, 2009, 202, 543–564.
8. J. S. Greaves, A. M. S. Richards, W. Bains, P. B. Rimmer, H. Sagawa, D. L. Clements, S. Seager, J. J. Petlowski, C. Sousa-Silva, S. Ranjan, E. Drabek-Maunder, H. J. Fraser, A. Cartwright, I. Mueller-Wodarg, Z. Zhan, P. Friberg, I. Coulson, E. Lee and J. Hoge, Nat. Astron., 2021, 5, 655.
9. T. Encrenaz, T. K. Greathouse, R. Giles, T. Widemann, B. Bezard, F. Lefevre, M. Lefevre, W. Shao, H. Sagawa, E. Marcq and A. Arredondo, Astron. Astrophys., 2024, 690, A304.
10. J. Chantzos, V. M. Rivilla, A. Vasyunin, E. Redaelli, L. Bizzocchi, F. Fontani and P. Caselli, Astron. Astrophys., 2020, 633, A54.
11. D. Smith, B. J. Mcintosh and N. G. Adams, J. Chem. Phys., 1989, 90, 6213.
12. P. A. Armentrout, Int. Rev. Phys. Chem., 1990, 9, 115.
13. M. Li, X. Li, Z. Yang, Z. Zhu and W. Xing, Int. J. Quantum Chem., 2025, 125, e70004.
14. B. K. Cunha de Miranda, C. Romanzin, S. Chefdeville, V. Vuitton, J. Zabka, M. Polasek and C. Alcaraz, J. Phys. Chem. A, 2015, 119, 6082.
15. C. Alcaraz, C. Nicolas, R. Thissen, J. Zabka and O. Dutuit, J. Phys. Chem. A, 2004, 108, 9998.
16. L. Nahon, N. de Oliveira, G. Garcia, J. F. Gil, B. Pilette, O. Marcouille, B. Lagarde and F. Polack, J. Synchrotron Radiat., 2012, 19, 508.
17. K. M. Ervin and P. A. Armentrout, J. Chem. Phys., 1985, 83, 166.
18. A. Zanchet, O. Roncero, E. Karabulut, N. Solem, C. Romanzin, T. Thissen and C. Alcaraz, J. Chem. Phys., 2024, 161, 044302.
19. A. Zanchet, J. L. Chen-Qiu, P. Larregaray, L. Bonnet, C. Romanzin, N. Solem, T. Thissen and C. Alcaraz, ACS Earth Space Chem., 2025, 9, 738.
20. V. Richardson, M. Polášek, C. Romanzin, P. Tosi, T. Thissen, C. Alcaraz, J. Zabka and D. Ascenzi, Molecules, 2024, 29, 810.
21. J. W. Au and C. E. Brion, Chem. Phys., 1997, 218, 87.
22. M. Halmann and Y. Klein, J. Chem. Soc., 1964, 4324.
23. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803.
24. H.-J. Werner, P. J. Knowles, R. Lindh, M. Schütz, P. Celani, T. Korona, F. R. Manby, G. Rauhut, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, A. W. Lloyd, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R. Pitzer, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and M. Wang, MOLPRO, version 2022, a package of ab initio programs, 2022, see https://www.molpro.net.
25. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007.
26. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053.
27. E. R. Davidson, J. Comput. Phys., 1975, 17, 87.
28. T.-S. Ho and H. Rabitz, J. Chem. Phys., 1996, 104, 2584.
29. A. Zanchet, B. Bussery-Honvault and P. Honvault, J. Phys. Chem. A, 2006, 110, 12017.
30. A. Zanchet, B. Bussery-Honvault, M. Jorfi and P. Honvault, Phys. Chem. Chem. Phys., 2009, 11, 6182.
31. P. Gamallo, A. Zanchet, F. J. Aoiz and C. Petrongolo, Phys. Chem. Chem. Phys., 2020, 22, 17091.
32. O. Denis-Alpizar, A. Zanchet and T. Stoecklin, Phys. Chem. Chem. Phys., 2024, 26, 13432.
33. A. Aguado and M. Paniagua, J. Chem. Phys., 1992, 96, 1265–1275.
34. A. Aguado, C. Suarez and M. Paniagua, J. Chem. Phys., 1993, 98, 308–315.
35. A. Aguado, C. Tablero and M. Paniagua, Comput. Phys. Commun., 1998, 108, 259–266.
36. A. Zanchet, M. Agúndez, V. J. Herrero, A. Aguado and O. Roncero, Astronom. J., 2013, 146, 125.
37. A. Zanchet, O. Roncero and N. Bulut, Phys. Chem. Chem. Phys., 2016, 18, 11391.
38. A. Zanchet, F. Lique, O. Roncero, J. R. Goicoechea and N. Bulut, Astron. Astrophys., 2019, 626, A103.
39. P. G. Jambrina, A. Zanchet, J. Aldegunde, M. Brouard and J. Aoiz, Nat. Commun., 2016, 7, 13439.
40. A. Zanchet, M. Menendez, P. G. Jambrina and J. Aoiz, J. Chem. Phys., 2019, 151, 094307.
41. J. Hernández-Rodríguez, C. Sanz-Sanz, P. A. Enríquez, M. González and M. Paniagua, Phys. Chem. Chem. Phys., 2023, 25, 28052.
42. A. Dorta-Urra, A. Zanchet, O. Roncero, A. Aguado and P. B. Armentrout, J. Chem. Phys., 2011, 135, 169901.
43. A. Dorta-Urra, A. Zanchet, O. Roncero and A. Aguado, J. Chem. Phys., 2015, 142, 154301.
44. R. Rydberg, Z. Phys., 1931, 73, 25.
45. M. Karplus, R. N. Porter and R. D. Sharma, J. Chem. Phys., 1965, 43, 3259.
46. C. Sanz-Sanz, A. Aguado, O. Roncero and F. Naumkin, J. Chem. Phys., 2015, 143, 234303.
47. A. J. Ocaña, E. Jimenez, B. Ballesteros, A. Canosa, M. Antiñolo, J. Albaladejoa, M. Agundez, J. Cernicharo, A. Zanchet, P. del Mazo, O. Roncero and A. Aguado, Astrophys. J., 2017, 850, 28.
48. P. del Mazo, A. Aguado, A. Zanchet and O. Roncero, ChemPhysChem, 2023, 24, e202300291.
49. X. Li, X. Zhang and B. Yan, Spectrochim. Acta, Part A, 2015, 142, 1–7.
50. C. P. Edwards, P. A. Jackson, P. J. Sarre and D. J. Milton, Mol. Phys., 1986, 57, 595–604.
51. K. P. Huber and G. Herzberg, Molecular spectra and molecular structure: IV. Constants of diatomic molecules, Springer Science & Business Media, 2013.
52. L. S. Sunderlin and P. A. Armentrout, J. Chem. Phys., 1994, 100, 5639.
53. S. Gomez-Carrasco, D. Felix-Gonzalez, A. Aguado and O. Roncero, J. Chem. Phys., 2022, 157, 084301.
54. G. F. Stowe, R. H. Shultz, C. A. Wight and P. B. Armentrout, Int. J. Mass Spectrom. Ion Processes, 1990, 100, 177.
55. J. L. Elkind and P. A. Armentrout, J. Phys. Chem., 1984, 88, 5454.

---