NCCP collisions with para -H 2 : Accurate potential energy surface and quantum dynamics at interstellar temperatures

The effect of para -hydrogen ( j p = 0) collisions on the rotational de-excitation transitions of molecule NCCP is investigated in this study. The scattering information is obtained by spherically averaging a four-dimensional potential energy surface (4DPES) o v er various H 2 molecule orientations. The calculations used the CCSD(T)-F12a method and aug-cc-pVTZ basis set to generate a 4DPES for the NCCP–H 2 van der Waals system. Within the NCCP– para -H 2 4DPES, a minimum energy point of 191.82 cm − 1 is attained at a distance of 3.6 Å from the centre of mass of H 2 and NCCP. To compute cross-sectional data for NCCP interacting with para -H 2 ( j p = 0), close coupling calculations are employed, encompassing total energies up to 600 cm − 1 . The resulting rate coefﬁcients [ k j → j (cid:3) ( T )] are calculated across a temperature range spanning from 5 K to 200 K. In accordance with propensity, even (cid:2) j = − 2 transitions are highly preferred. Comparatively, the derived k j → j (cid:3) ( T ) for NCCP–H 2 are determined to be 1.5– 4.5 times of NCCP–He. This observation implies that relying on a scaling factor of 1.38 to extrapolate rate coefﬁcients for H 2 collisions from those of NCCP–He is not a reliable approach.

CCP species, one can adopt a rotational temperature of 10-50 K for NCCP (Halfen, Clouthier & Ziurys 2008 ).This temperature is a common parameter for species found in the outer layers.The column density values of NCCP and CCP indicated a slightly lower abundance of NCCP than CCP in the IRC + 10216 star.Given the sufficient abundance of reactants in the ISM, it could be formed through the reaction between HNC and CP (Pham-Tran et al. 2001 ).Among the C 2 NP isomers, NCCP is the most stable one (Kwon & McKee 2001 ) and its internal state distribution can be comprehended by facilitating collisions with the most pre v alent gases in the ISM.
Despite the crucial role of phosphorus in biological systems, our understanding of phosphorus gas-phase chemistry is limited (Pasek & Lauretta 2005 ), but the researchers are nowadays focused on interstellar phosphorus chemistry.Ho we v er, the y make conflicting predictions regarding the abundance of the primary phosphoruscontaining molecules (Charnley & Millar 1994 ;Maci á, Hern ández & Or ó 1997 ).These models can only be advanced by determining the molecule abundance based on the observations.To model molecular emission from interstellar clouds, it is mandatory to calculate the rate coefficients resulting from collisions with hydrogen, which is the most pre v alent species in ISM (Ritika, Kumar & Dhilip Kumar 2022 ).Hydrogen exists in two nuclear spin states, i.e. para -hydrogen ( j p = 0) and ortho -hydrogen ( j o = 1) with a ratio of 1:3 at room temperature (Lee et al. 2006 ;Kushwaha et al. 2023 ).Ho we ver, at extremely low temperatures, hydrogen predominantly exists as para -H 2 ( j p = 0).This paper focuses on determining rotational deexcitation rate coefficients of NCCP in the presence of para -H 2 at temperatures between 5 and 200 K.Throughout the collision process, it is assumed that para -H 2 remains in its lowest rotational state ( j p = 0).The assumption is valid because within this temperature range, the probability of exciting of H 2 molecule to the rotational level j p ' = 2 is negligible (Roueff & Lique 2013 ;Ritika & Dhilip Kumar 2023 ).
In this work, we have computed a new truncated ab initio fourdimensional potential energy surface (PES) for the NCCP-para -H 2 , considering three different orientations of H 2 .To achieve this, we employ the explicitly correlated coupled-cluster singles and doubles with perturbative triples with F12 approximation, i.e. the CCSD(T)-F12a method.
The calculations are performed using an augmented correlation consistent polarized valence triple zeta (aVTZ) basis set, which has been observed to be the most accurate approximation to CCSD(T)/CBS (complete basis set) method.The close-coupling approach (CC) is further used to investigate the transitions in the collision of NCCP with para -H 2 .The results are then compared to the previous findings for the NCCP-He complex.To e v aluate the reliability of approximating k j → j ( T ) for para -H 2 based on rate coefficients obtained from helium collisions, we use a square-root formula that takes into account the ratio of their reduced masses.The paper's subsequent sections elaborate on the content presented.Section 2 provides details about the computational aspects of the averaged 2D PES calculations.In Section 3 and 4 , we delve into the quantum dynamical calculations and results of the inelastic rotational cross-sections and rate coefficients.The paper is finally concluded in Section 5 .

AB INITIO C A L C U L AT I O N S O F N C C P -H 2
In this study, we focus on rotational collisions between NCCP and para -H 2 at interstellar temperatures.Both molecules, NCCP and H 2 are treated as rigid rotors during the PES calculations.According to a recent investigation by Faure et al. ( 2016 ), using state-averaged geometries yield scattering data for the CO-H 2 system that is similar to data obtained through full-dimensional calculations.Therefore, we fix the H 2 bond length at the ground vibrational level's average distance, r HH = 0.766 Å.
For NCCP, we use the experimental equilibrium distance (Pham-Tran et al. 2001 ) with its bond distances r NC = 1.159Å, r CC = 1.378Å, r CP = 1.554Å.The geometric parameters of the NCCP-H 2 system are described using certain angles ( θ a , θ b , and δ) and a distance parameter ( R ), which represents the separation between the centre of mass of both NCCP and H 2 as shown in Fig. 1 .The θ a and θ b denote the rotations of NCCP and H 2 , respectively, relative to the z -axis and the symbol δ indicates the dihedral angle that describes the relative orientation between NCCP and H 2 .The adiabatic 4DPES, which relies on the parameters θ a , θ b , δ, and R is
In the F12 calculations, F12 integrals are computed using density fitting approximations.We have used the 'ri basis = aVTZ/MP2FIT' in PES calculations.We present a comparison in Table 1 , where we hav e e xamined different CCSD(T)-F12 approximations and basis sets.The errors with respect to CCSD(T)/CBS are e v aluated and presented in a table.It is observed that CCSD(T)-F12a/aVTZ proved to be sufficiently accurate in describing van der Waals intermolecular interactions (Tew & Klopper 2006 ;Roueff & Lique 2013 ).All calculations are conducted using MOLPRO (Werner et al. 2015 ) package.The values of R ranges from 3.0 to 26 Å and there are 19 discrete values of θ a in the grid, ranging from 0 • to 180 • in increments of 10 • .For each combination of R and θ a , we consider three distinct orientations of H 2 molecules, which are given as: ( Consequently, 2964 geometries are computed for the NCCP-H 2 system in the C s symmetry group, represented by ( R , θ a , θ b , and δ).In all surface computations, we apply the standard Boys & Bernardi approach to rectify the basis set superposition error: (Boys & Bernardi 1970 ) The collinear configuration of NCCP •••H-H occurs when θ a is equal to 0 • and θ b is also equal to 0 • and H-H ••• NCCP corresponds to the configuration when θ a = 180 • and θ b = 0 • .Fig. 2 displays the onedimensional potential energy curves of NCCP-H 2 concerning the variable R .It illustrates these curves for three different orientations of H 2 at θ a = 0 • , 90 • , and 180 • .As evident from this figure, the potential exhibits strong anisotropy due to the orientation of H 2 .
The most fa v ourable interaction arises from a collinear (H and N head-to-head ), which is influenced by the H 2 quadrupole and NCCP dipole orientations.In the linear arrangement with angles θ a = 0 • and θ b = 0 • , the interaction is weaker due to a repulsive long-range electrostatic force between the unfa v ourable NCCP dipole and H 2 quadrupole orientations.
When θ a = 90 • , the interaction potential shows minor dependency on the rotation of the H 2 .
Since our objective involves solving the CC equations to analyse scattering phenomena, expanding the PES function, V , at a given R value using angular functions is necessary.In the context of scattering between two linear rigid rotors, the potential can be represented as follows (Launay 1977 ;Flower et al. 1979 ).
(2)    where with μ ≤ min ( l 1 , l 2 ).Here, inde x es l 1 and l 2 are connected, respectively, with the rotational motion of NCCP and H 2 .Specifically, index l 2 is even, reflecting the homonuclear geometry of the molecule H 2 .
In specific scenario of collisions involving para -H 2 molecules (where j p = 0), with l 2 = 0, and reduced mass μ = 0, the corresponding wavefunction simplifies to Y 00 ( θ b , δ) = 4 × π −1/2 .The expectation value of the potential is then deduced by averaging the total PES o v er the angular motion ( θ b , δ) of the H 2 entity: where P l 1 denotes the Legendre polynomials.Equation ( 4) has the appropriate form to expand the interaction potential between the fixed rotor and spherical projectile.As mentioned MNRAS 527, 9826-9832 (2024)   earlier, in the context of low-energy rotational excitation, the impact of the j p ≥ 2 channels of H 2 are ignored.Additionally, the energy difference between j p = 0 and j p = 2 levels (510 K) is significantly larger than the thermal kinetic energy ( T ≤ 150 K), suggesting that the influence of closed channels ( j p = 0) on cross-sections is expected to be minor.Therefore, we have chosen to work with the minimal rotational basis of j p = 0 for para hydrogen.The PES for NCCP-H 2 ( j p = 0) is simplified into a 2DPES, V av , computed by averaging the total PES o v er angles θ b and δ, where The spherically averaged 2DPES has a global minimum of 125.40 cm −1 at θ a = 110 • and R = 3.6 Å.It is interesting to contrast the differences and similarities between the PES for NCCP-He and our computed reduced PES for NCCP-H 2 (Ritika & Dhilip Kumar 2022 ).The o v erall behaviour of both PESs is similar, with global minima at 110 • , but reduced NCCP-H 2 has a deeper well measuring 46.40 cm −1 at R = 3.5 Å. Fig. 3 presents the spherically averaged 2D contour plot for NCCP-H 2 .For fitting the averaged surface, we have included 29 different radial ( v l 1 ) coefficients.The multipole expansion coefficients are presented in Fig. 4 as a function of radial coordinates.Due to the asymmetric nature of NCCP, both even and odd v l 1 values persist.These calculations offer insightful information regarding the rotational de-excitation rates in NCCP-H 2 collisions.

Computational details
The CC method formulated by Arthurs, Dalgarno & Bates ( 1960 ) is employed for the determination of the integral inelastic crosssection between different quantum states for NCCP-para -H 2 system using the MOLSCAT quantum package (Hutson & Green 1994 ).The calculations are done using the AIRY propagator, an approach advanced by Alexander & Manolopoulos ( 1987 ) in 1987.The NCCP possesses spectroscopic attributes denoted as B 0 (0.0902 cm −1 ) and MNRAS 527, 9826-9832 (2024) D 0 (6.708 × 10 −9 cm −1 ), and it has a linear configuration in its electronic state (X 1 + ).The collisional problem is limited to a rigid rotor by disregarding the H 2 rotational motion and the NCCP vibrational motion.The range of total energy considered for the dynamic computations spans from 0.5 to 600 cm −1 , encompassing all resonances.The parameters of the propagator are adjusted to strike a balance between accuracy and computational efficiency .Notably , R min is set at 3.0 Å while R max is set at 26 Å.The calculations are conducted using different step sizes depending on the total energy range: E t = 0.1 cm −1 for E t ≤ 60 cm −1 , E t = 0.2 cm −1 for 60-120 cm −1 , E t = 0.5 cm −1 for 120.5-160 cm −1 , E t = 1 cm −1 for 161-210 cm −1 , E t = 5 cm −1 for 215-250 cm −1 and E t = 50 cm −1 for 250-600 cm −1 .This range ensures the inclusion of all resonances with lower energy levels.The STEPS parameter, dictating the grid magnitude for propagation of short range, is set sufficiently high to enable a thorough analysis of low-energy resonances.Typically, STEPS = 10 is set for E t levels greater than 100 cm −1 , and STEPS = 20 is used for E t levels less than 100 cm −1 .The j max , rotational basis, is determined at a level that allows all coupling coefficients to be taken into account in dynamical calculations.This consideration guarantees the inclusion of both open and closed channels in the analysis, yielding accurate scattering data.For instance, the parameter j max is initially configured at 20 for energy levels ranging up to 60 cm −1 , with a gradual increment to 60 for energy levels extending up to 600 cm −1 .The determination of j max values is achieved through optimization, where the thresholds for inelastic and elastic transitions are defined as OTOL = 0.005 Å 2 and DTOL = 0.05 Å 2 , respectively.

De-excitation cross-sections
Figs 5 (a) and 5 (b) indicate the de-excitation cross-section variation with respect to the total energy when it interacts with the perturber, H 2 .The analysis discloses that rotational shifts are feasible for transitions of e ven-to-e ven, odd-to-odd, odd-to-e ven, and e ven-toodd nature.The molecular asymmetry of NCCP results in transitions that produce both even and odd j values and its energy levels allow both odd and even j values.Both illustrations reveal the presence of shape and Feshbach resonances (FR) within energy limits up to E t = 300 cm −1 (Christoffel & Bowman 1983 ).The FR are pronounced at E t ≤ 130 cm −1 but wane as E t rises, disappearing past 300 cm −1 .
Quasi-bound states are detected as a consequence of the H 2 projectile being confined within a potential well.This confinement gives rise to FR, which persists until the H 2 possesses sufficient energy to get out of the trap.In contrast, shape resonances emerge as a result of tunnelling effects that enable the H 2 projectile to pass through the centrifugal barrier.Moreo v er, as the j value increases, the cross-section trend generally ascends, consistent with typical rotational de-excitation behaviour.
Fig. 5 (a) and 5 (b) depict the rotational de-excitation crosssections for both j = −1 and j = −2 transitions among the eight rotational levels.As j rises (from j = 1 to 4 for j = −1 and 3 to 5 for j = −2), the cross-sections in both cases show an increase.The information on variation of final state ( j ') concerning crosssections are provided in Fig. 6 .The cross-sections decrease as the collision energy is increased from 100 cm -1 to 400 cm -1 .These results demonstrate a notable inclination towards fa v ouring e ven v alues of j compared to odd ones.The same propensity is observed in other phosphorus-containing systems such as PN-H 2 by Najar et al. ( 2017 ) and HPO-He by El Hanini et al. ( 2021), where j = −2 fa v ours.

R AT E C O E F F I C I E N T S O F N C C P -PARA -H 2 A N D C O M PA R I S O N W I T H N C C P -H E
Here, the k j → j ( T ) for NCCP-para -H 2 interaction is calculated by thermally averaging the cross-sections obtained from CC calculations at various temperatures.The temperature range used for the measurements is 5-200 K.The rate coefficients between state ( j ) and ( j ') are calculated using the following Boltzmann average formula: The kinetic energy is given as E k = E t − E j , with E j showing the rotational energy and k B denoting the Boltzmann constant.
The reduced mass and temperature are symbolized by μ and T , respectively.Fig. 7 illustrates the rate coefficients of interest as a function of temperature.The observed trends confirm the same propensity that fa v ours even j = −2 transitions over j = −1 transitions.The obtained values of k j → j ( T ) for NCCP-para -H 2 are compared with previously determined rate coefficients of NCCP-He.It is generally assumed that k H 2 can be deduced from k He by using a scaling factor (SF) that takes into account the difference in their reduced masses.The SF for NCCP can be determined using the following equation: The abo v e approximation suggests that the cross-sections involving para -H 2 ( j p = 0) and He are the same and this SF of 1.38 is only attributed to the reduced mass factor incorporated in equation ( 6), which determines the associated rate coefficients.Since a limited number of systems have been examined with both He and H 2 , it's necessary to verify the accuracy of this approximation.Fig. 8 provides a comparison of k j → j ( T ) data, comparing the results for NCCP-He with NCCP-H 2 .This indicates that accurate collisional rate coefficients using helium collision cannot be achieved for para -H 2 interactions.The observed differences in the results are expected due to the difference in the well-depth of both NCCP-He and NCCP-H 2 .NCCP-H 2 has three times deeper well than NCCP-He resulting in a different range of resonances presented in Fig. 5 .The ratios of  the k j → j ( T ), i.e. [NCCP-para -H 2 ( j p = 0)/NCCP-He], are tabulated in T able 2 .T able 2 shows that the ratio is inconsistent and exhibits fluctuations at T = 10 K and T = 50 K.At low-temperature values, the differences between the two sets of systems are more pronounced, and as the temperature rises further, the correspondence between the two sets of coefficients w ould lik ely become quite fa v ourable.The outcomes align with the earlier research by Lique et al. ( 2007 ) concerning collisions involving SO and H 2 .In an attempt to examine these differences, we have also graphed the rate coefficients for various transitions at T = 10 K and T = 50 K as illustrated in Fig. 9 .The estimated abundance of NCCP in molecular clouds could have been revised in light of these recently calculated rate coefficients.The bending frequency ( v 2;  2014 ), where minor differences between rigid and non-rigid models were found, almost identical results were obtained without and with considering the vibrations.For the NCCP molecule, this should be the case.

S U M M A RY A D C O N C L U S I O S
The collisional cross-sections involving NCCP and para -H 2 are computed using the exact CC approach.These computations span a wide spectrum of E t ranging from 0.5-600 cm −1 .Using Boltzmann velocity distribution, the cross-sections are averaged to calculate j → j ( T ) for T up to 200 K. accuracy of dynamical data highly depends on the interaction potentials of colliding entities developed for the collision NCCP and H 2 .To model this collision, a new NCCP-H 2 4DPES is computed the F12a/aug-cc-pVTZ level of theory.This PES is subsequently averaged to 2DPES using three distinct angular orientations of H 2 .Within the averaged NCCP-H 2 PES, a global energy minimum of 125.40 cm −1 is identified, which is notably deeper than the 46.40 cm −1 minimum energy observed in the NCCP-He collision.The point of maximum depth was found at an angle of 110 • for both collision systems.Upon expansion of PES, both systems demonstrate analogous trends and shapes.The isotropic term v 0 increases by about ∼20 cm −1 due to the lower minimum energy of the NCCP-H 2 collision.Furthermore, the spectrum of resonances in the NCCP-para -H 2 collision spans a wide range, extending to around 200 cm −1 .In contrast, the NCCP-He collision's range only encompassed approximately 50 cm −1 .Comparing the newly calculated k j → j ( T ) for NCCP-para -H 2 collisions with the existing MNRAS 527, 9826-9832 (2024) NCCP-He rate coefficients reveals notable distinctions, particularly at lower temperatures.The rate coefficients for collisions with He appear to provide a plausible model for interactions with para -H 2 at moderately higher temperatures.It's crucial to highlight that the SF cannot be deemed a dependable approach for replicating para -H 2 collisions.Analysing the data, it becomes evident that the rate coefficients for NCCP-H 2 collisions are 1.5-4.5 times greater than those for NCCP-He collisions.Modelling the NCCP abundance in the ISM will be simpler owing to new insights into collision rates and other dynamic properties between H 2 and NCCP.

AC K N OW L E D G E M E N T S
The access to HPC resources is acknowledged by the authors with gratitude to IIT Ropar.Ritika is grateful for the support received from the MoE for the research funding through the PMRF fellowship (PMRF-192002-244).Additionally, the authors would like to acknowledge the SERB, New Delhi funding under grant no.SERB-CRG/2020/003006, which supported the core research activities.

Figure 3 .
Figure 3. Spherically averaged 2D contour plot for NCCP-H 2 as function of θ a and R .

Figure 7 .
Figure 7. Rate coefficients of rotational de-excitation of NCCP-para -H 2 for selected transitions at 5-200 K.

Figure 9 .
Figure 9.Comparison of rate coefficients of NCCP with para -H 2 and He along with the SF (1.38) at T = 10 K and 50 K.