Protonation of apolar species: From Cl₂H⁺ to (E)-NCCHCHCNH⁺ through computational investigations

1 Introduction

The overwhelming majority of chemical species in the interstellar medium (ISM) have been discovered via radioastronomy, which exploits their rotational signatures. However, the observation of new astronomical species comes to a setback when the species of interest are weakly polar or apolar (Yamamoto, 2017). In the former case, a possible way out is offered by the line stacking technique (Yen et al., 2016), while for apolar molecules, one has to rely on the observation of the so-called molecular proxies. These are polar systems whose presence in the astronomical object under consideration can be directly related to non-observable species. Molecular proxies can be obtained in several ways, with some of them being well established in the literature.

A first possibility is to claim the presence of an apolar species based on the observation of its cyano-derivatives, which usually have a significant electric dipole moment. Indeed, the high abundance of the CN radical together with its highly reactive nature leads to species containing the CN moiety through formation pathways that are usually accessible even in the conditions of the ISM. A significant example is provided by the observation of benzonitrile (C₆H₅CN) (McGuire et al., 2018), whose detection in the ISM has been undoubtedly related to the presence of benzene in the same environment. Indeed, different experiments have demonstrated the effectiveness of the reaction between the CN radical and benzene C₆H₆ (Lee et al., 2019; Cooke et al., 2020). Similarly, the observation of cyanonaphthalene towards TMC-1 (McGuire et al., 2021) suggests the presence of naphthalene in this cold core. Other than the CN radical, proxies can also be formed by the CCH radical. However, the observation of its derivatives is more challenging because the ethynyl group leads to molecular proxies with a smaller electric dipole moment when compared to the CN-containing counterpart (Cernicharo, J. et al., 2021; McCarthy et al., 2021).

Another important way to observe apolar species is to consider their protonated form. Protonation in the ISM occurs owing to the presence of ${{\text{H}}_3}^{+}$, which in turn is obtained from the reaction between H₂ and ${{\text{H}}_2}^{+}$, with the latter being formed via cosmic ray ionization (Yamamoto, 2017). The ${{\text{H}}_3}^{+}$ initiates the protonation via this general mechanism:

which proceeds only if the proton affinity (PA) of the generic species A is greater than that of H₂. Based on the PA values reported by Yamamoto (2017), ${{\text{H}}_3}^{+}$ is the primary source of N₂H⁺ and HCO⁺/HOC⁺, but it does not lead to the formation of O₂H⁺ because O₂ has a smaller PA than H₂.

In the present paper, we propose the use of protonated species as suitable proxy for the radio-observation of apolar species of astrochemical interest. To this end, their computed rotational spectra are reported here for the first time. This study will involve simple diatomic molecules, such as P₂, Si₂, and Cl₂, a linear centrosymmetric species, i.e., CNNC, and a C_2h species, the (E) isomer of NC − CH = CH − CN. The first three are simple diatomic molecules that can be considered as reservoirs of atomic P, Si, and Cl, respectively. This is particularly interesting because the chemistry of third-row atoms in astronomical environments is still in early stages, with only a few species containing P, Si, or Cl being detected so far (McGuire, 2022). Since both vibrational (infrared) and rotational spectra of homonuclear diatomics cannot be observed, the only way-out for their detection (if avoiding electronic spectroscopy) is offered by their protonated forms. As an example, the N₂H⁺ molecule–observed for the first time in 1974 by Turner (1974)—is considered a good tracer of N₂, which is otherwise not observable at radio frequencies. For this reason, the study of the simple ${{\text{HP}}_2}^{+}$, ${{\text{HSi}}_2}^{+}$, and Cl₂H⁺ ions can potentially provide useful data to guide laboratory and/or astronomical observations. While the P₂ and Si₂ precursors might be considered somewhat exotic, Cl₂ is a simple species stable under Earth conditions. To assess the accuracy of our work, the N₂H⁺ and HCO⁺ ions will be used as references to test our computational strategy.

Another interesting molecule is the linear species diisocyanogen (CNNC) because two isomers of its family, namely CNCN and NCCN, have already been detected in the ISM. While the former was observed directly via its rotational transitions (Agúndez et al., 2018), the presence of NCCN was inferred thanks to the observation of its protonated counterpart, i.e NCCNH⁺ (Agúndez et al., 2015). Thus, the computational spectroscopic characterization of CNNCH⁺ surely represents a preliminary step toward its astronomical detection. In this case, the well-characterized NCCNH⁺ species will be used as reference to understand the uncertainties associated with our methodology.

The last molecule considered is the (E) isomer of 1,2-dicyanoethene, which is the main product of the reaction between the CN radical and vinyl cyanide (CH₂ = CH − CN, Marchione et al. (2022)), also known as acrylonitrile. While the (Z) isomer of 1,2-dicyanoethene, being a polar molecule, might be detected directly, without the use of any molecular proxy, the (E) form lacks of a permanent electric dipole moment. Thus, we aim at providing good estimates for the rotational parameters of its counterpart, i.e. (E)-NC − CH = CH − CNH⁺.

This work will also report the fundamental vibrational frequencies and intensities of all the species mentioned above. These data could provide reliable starting points for experimental measurements of their vibrational spectra as well as will be benchmark references for future computations.

The manuscript is organized as follows. In the next Section 2, we detail the computational strategy adopted to describe the protonated species, with specific focus on the computation of proton affinity and formation routes as well as spectroscopic parameters. Then, the results are reported in Section 3 where they are discussed together with those of chemically related species. After briefly addressing the PA values, the formation routes for the species considered are detailed, and then the computational spectroscopic characterization is discussed. The last Section 4 will draw the main conclusions of the present study.

2 Computational details

For most of the species considered in this work, experimental and computational data are missing, thus a preliminary study using density functional theory (DFT) was carried out to understand their stable form. For this purpose, the double-hybrid rev-DSD-PBEP86 functional (Santra et al., 2019), combined with the GD3BJ empirical dispersion (Grimme et al., 2011; 2010), was used in conjunction with the jun-cc-pVTZ basis set (Dunning Jr, 1989; Papajak and Truhlar, 2011). Overall, this level of theory is shortly denoted as revDSD/junTZ. Noted is that, for third-row elements, the d-augmented jun-cc-pV(T+d)Z set has been employed (Dunning et al., 2001). The revDSD/junTZ level has been employed to derive preliminary values for equilibrium geometries, electronic energies and (harmonic) zero-point energy corrections (ZPEs) of the protonated species here considered.

2.1 Proton affinity and formation routes

To understand if a protonated molecule can be formed in the ISM starting from ${{\text{H}}_3}^{+}$, its PA value has to be considered in comparison with that of the H₂ molecule. The PA at T = 0 K is computed as (McNaught and Wilkinson, 1997; Marino et al., 2000):

where, in this study, both ΔH₀ values are obtained by combining the electronic energy obtained using the jun-Cheap composite scheme (hereafter junChS) (Alessandrini et al., 2019) with anharmonic ZPE. The latter correction was computed at the ae-CCSD(T)/cc-pwCVQZ for protonated diatomics, while the cc-pwCVTZ basis set was used for NCCNH⁺ and CNNCH⁺. For protonated 1,2-dicyanoethene, the ZPE was evaluated at the revDSD/junTZ level. Here, ae-CCSD(T) stands for the coupled-cluster (CC) method including a full treatment of single and double excitations and a perturbative estimate of triples (Raghavachari et al., 1989), with all electrons (ae) correlated. Anharmonic computations were carried out at the corresponding reference geometry, on top of which junChS energy evaluations were also performed. Within the junChS approach, the electronic energy, E(junChS), is obtained as follows (Alessandrini et al., 2019):

Here, the first term on the right-hand side is the fc-CCSD(T)/jun-cc-pVTZ energy (with jun-cc-pV(T+d)Z being used for third-row elements), where fc indicates that the frozen-core approximation is employed. The second term accounts for the extrapolation to the complete basis set (CBS) limit at the MP2 level, with MP2 denoting Møller-Plesset perturbation theory to the second order (Møller and Plesset, 1934). This extrapolation exploits the n⁻³ formula by Helgaker et al. (1997), and requires energy computations with two consecutive members of a hierarchical family of basis sets. In this case, the jun-cc-pV(T+d)Z and jun-cc-pV(Q+d)Z sets were used. The third term in Eq. 3 is the core-valence (CV) correction and incorporates the effect of the correlation of inner-shell electrons, which is not included in the previous computations. This term is computed as the difference between all-electron and frozen-core calculations, both with the cc-pwCVTZ basis set (Peterson and Dunning Jr, 2002).

For three species, namely Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and CNNCH⁺, the formation route was also investigated. To this end, the revDSD/junTZ level of theory was employed to locate minima (MINs) and transition states (TSs) on reactive potential energy surface (PES). The correct assignment of two MINs to the connecting TS was assured by exploiting the intrinsic reaction coordinate (IRC) analysis (Fukui, 1981), performed at the same level of theory. To refine the electronic energy of each stationary point, the junChS approach was employed, while the harmonic ZPE was evaluated at the revDSD/junTZ level.

2.2 Spectroscopic characterization

Focusing on rotational spectroscopy, the main parameters influencing the spectrum are the rotational constants (Barone et al., 2021). For a given inertia axis (γ = a, b, or c), in the framework of vibrational perturbation theory to second order (VPT2) (Mills, 1972), the rotational constant of the vibrational ground state ($B_0^{\gamma }$, corresponding to A₀ for γ = a, B₀ for γ = b, and C₀ for γ = c), can be written as:

The first term on the right-hand side of Eq. 4 is the equilibrium rotational constant, which entirely depends on the isotopic composition of the molecule and equilibrium structure. Thus, $B_e^{\gamma }$ is straightforwardly obtained from a geometry optimization and accounts for about 99% of the value of $B_0^{\gamma }$ (Puzzarini et al., 2008; 2010; Alessandrini et al., 2018). The second term on the right-hand side of Eq. 4 is the vibrational contribution to the rotational constant which is obtained as the sum over all r vibrational modes of the vibration-rotation interaction constants $({\alpha }_r^{\gamma })$. In Eq. 4, d_r denotes the degeneracy of the rth vibrational mode. The vibrational correction represents a minor contribution to $B_0^{\gamma }$, but its consideration is crucial to reach accurate computational results, comparable with experimental ones (Puzzarini et al., 2008; 2010; Alessandrini et al., 2018).

Since $B_e^{\gamma }$ is the major contribution to $B_0^{\gamma }$, it should be computed at the highest possible level of theory when aiming at high accuracy. For this reason, for Cl₂H⁺, ${{\text{HP}}_2}^{+}$, ${{\text{HSi}}_2}^{+}$, and CNNCH⁺, the equilibrium geometry was obtained by exploiting a composite scheme entirely based on coupled-cluster (CC) theory (Heckert et al., 2005; 2006), which requires the minimization of the following energy gradient:

This composite scheme is denoted as “CBS+CV+fT+fQ” and it consists of 5 different terms.

• The extrapolation to the CBS limit, which is obtained by the sum of the HF-SCF and fc-CCSD(T) extrapolated gradients, i.e. $\frac{d{E}^{\infty }(\mathrm{H}\mathrm{F}-\mathrm{S}\mathrm{C}\mathrm{F})}{dx}$ and $\frac{d\mathrm{\Delta }{E}^{\infty }(\mathrm{C}\mathrm{C}\mathrm{S}\mathrm{D}(\mathrm{T}))}{dx}$, respectively. The former is evaluated using the three-point formula by Feller (1993) and the cc-pVTnZ basis sets, with n = T,Q,5. For the CCSD(T) correlation energy, ΔE (CCSD(T)), gradient the n⁻³ formula (Helgaker et al., 1997) is exploited in conjunction with the cc-pVTZ and cc-pVQZ basis sets.

• The CV contribution, which is obtained as difference between ae- and fc-CCSD(T)/cc-pwCVTZ computations.

• The fT term, which accounts for the full treatment of triples. This is obtained as difference between two computations carried out with the cc-pVTZ basis set, the first employing the CCSDT (Noga and Bartlett, 1987; Scuseria and Schaefer, 1988) method (the acronym standing for CC singles, doubles and triples) and the second one with CCSD(T), both within the fc approximation.

• The fQ term, which is the contribution due to quadruple excitations. This is computed as difference between CCSDTQ (CC singles, doubles, triples and quadruples; Kállay and Gauss (2008); Kucharski and Bartlett (1991)) and CCSDT calculations, both with the cc-pVDZ basis set and within the fc approximation.

The CBS+CV+fT+fQ equilibrium rotational constants were then corrected for the vibrational contribution. The evaluation of vibration-rotation interaction constants requires a full cubic force field (Puzzarini et al., 2010), which was computed at the ae-CCSD(T)/cc-pwCVQZ level. For the CNNCH⁺ ion, a smaller basis set, cc-pwCVTZ, was instead employed. For the (E)-NC − CH = CH − CNH⁺ species, the computational approach was slightly different because of its size. The equilibrium rotational constants were obtained from the minimization of the “CBS+CV” energy gradient, thus retaining only the first three terms on the right-hand side of Eq. 5 in conjunction with the same basis sets mentioned above. In addition, the anharmonic computation required for the vibrational corrections was carried out at the revDSD/junTZ level of theory.

For all the species considered, the VPT2 calculations performed on the computed anharmonic force field allowed for obtaining other spectroscopic parameters, such as quartic centrifugal distortion constants of the S-reduced Watson Hamiltonian (D_J, D_JK, D_K, d₁, d₂) (Papousek and Aliev, 1982) and fundamental vibrational frequencies (ν). These latter can be obtained because the computations of the cubic force constants give access, without any other additional cost, to the semi-diagonal quartic force constants. For protonated 1,2-dicyanoethene, quartic centrifugal distortion constants were computed at a different level of theory than the vibrational corrections to rotational constants, indeed considering the harmonic force field at the fc-CCSD(T)/jun-cc-pVTZ level.

On top of the best geometry, the electric dipole moment (μ) and the nuclear quadrupole-coupling constants (NQCCs, χ_ii, with i = a, b, c), which are both first-order properties, were obtained. While the electric dipole moment is required to predict the type and intensity of rotational transitions, the NQCCs shape the rotational spectrum (hyperfine structure). NQC is the interaction between the quadrupole moment of a nucleus and the electric gradient at the same nucleus. This determines a splitting of the rotational energy levels, thus causing a splitting of the rotational transitions, i.e. the so-called hypefine structure. Noted is that a nucleus with I > 1/2 has a non-zero quadrupole moment. The hyperfine structure can be useful for the assignment and analysis of rotational lines in the laboratory and in astronomical observations (Turner, 1974; Puzzarini et al., 2010; Melosso et al., 2021); thus, spectral simulations strongly benefit from the introduction of such interaction terms. Both NQCCs and electric dipole moment components were obtained at the ae-CCSD(T)/aug-cc-pwCVQZ level for protonated diatomics, while for CNNCH⁺ the aug-cc-pwCVTZ basis set was used. For protonated (E)-1,2-dicyanoethene, the fc-CCSD(T)/jun-cc-pVTZ level of theory was employed. Additionally, the vibrational corrections to the NQCCs and the electric dipole moment were also considered and obtained from the anharmonic treatment mentioned above.

The computational strategy outlined above was also employed to describe three molecules that are experimentally well-characterized in the literature. These are N₂H⁺ and HCO⁺, which offer a good reference for the new protonated diatomics, and NCCNH⁺ for the direct comparison with CNNCH⁺.

All DFT and MP2 computations were carried out using the Gaussian16 suite of programs (Frisch et al., 2016), while the CFour quantum-chemistry program (Stanton et al., 2016; Matthews et al., 2020) was used for all computations based on the CC theory. Concerning VPT2 calculations, its generalized version (GVPT2) (Schiff, 1955; Piccardo et al., 2015) was employed whenever using the Gaussian16 package.

3 Results

As starting point, the work investigated the stable form of Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and ${{\text{HSi}}_2}^{+}$. For this purpose, equilibrium structures ranging from linear to cyclic were considered and the most stable forms of each species are shown in Figure 1, where their electronic state is also reported. The same figure also provides the equilibrium structure of CNNCH⁺ and (E)-NC − CH = CH − CNH⁺, which were somewhat more straightforward to derive since bent structures would be very unstable.

FIGURE 1. Equilibrium structure, point group symmetry and electronic ground state of (1) Cl₂H⁺, (2) ${{\text{HP}}_2}^{+}$, (3) ${{\text{HSi}}_2}^{+}$, (4) CNNCH⁺, and (5) protonated (E)-1,2-dicyanoethene.

The main molecule driving protonation in the ISM being ${{\text{H}}_3}^{+}$, the first mandatory step is to investigate whether the species considered in this work can be produced by reaction 1. The first screening is based on PA. If PA of P₂, Si₂, Cl₂, CNNC, and NCCHCHCN is greater than that of H₂, the protonation occurs and the reaction proceeds at the Langevin rate (Yamamoto, 2017). The latter is typically around 10^–9 cm³s⁻¹ and is rather independent from the temperature. The data required for this comparison are reported in Table 1, which also includes the data for N₂, CO, and NCCN, whose protonated forms have already been observed in the ISM.

TABLE 1. Computed proton affinities (kJ mol⁻¹) for the molecules considered in this work. Electronic energies and ZPE are in Hartrees.

Table 1 indicates a PA value of 418 kJ⋅mol⁻¹ for H₂, in good agreement with what reported by Yamamoto (2017) and Hunter and Lias (1998), i.e. 422.61 kJ⋅mol⁻¹. A similar agreement is noted for N₂ and CO, with values of 493.5 kJ⋅mol⁻¹ and 591.1 kJ⋅mol⁻¹, respectively, to be compared with 493.8 kJ⋅mol⁻¹ and 594 kJ⋅mol⁻¹ given in Yamamoto (2017). The PA values reported in this work refer to 0 K and are more accurate than those available in the literature. The differences noted can be ascribed to the (small) thermal correction, the literature values referring to a temperature of 298.15 K.

From Table 1, it is noted that a PA value greater than that of H₂ has been obtained for Cl₂, thus suggesting the possible formation of Cl₂H⁺ in the ISM. A further support to this is provided by the fact that PA(Cl₂) is intermediate between PA(N₂) and PA(CO), with both N₂H⁺ and HCO⁺ detected in the ISM. In particular, our PA(Cl₂) value of 547.6 kJ⋅mol⁻¹ is in good agreement with that established experimentally by Cacace et al. (1998), i.e. 548.1 ± 13 kJ⋅mol⁻¹.

The PA value of P₂ is 669.9 kJ⋅mol⁻¹ and it is very close to that of NCCN, i.e., 652.7 kJ⋅mol⁻¹, thus confirming the plausible formation of the protonated form of P₂. Finally, large PA values have been obtained for CNNC, (E)-NCCHCHCN, and Si₂. For the first two species, the PA values are 717.8 kJ⋅mol⁻¹ and 735.5 kJ⋅mol⁻¹, respectively, which are comparable with those of molecules like H₂O, H₂S, and H₂CO (Yamamoto, 2017). The Si₂ species shows the largest PA, 841.0 kJ⋅mol⁻¹, this result being likely due to its triplet electronic ground state. Indeed, a similar value (869 kJ⋅mol⁻¹) is obtained for the CCS radical in the same electronic state (Yamamoto, 2017). Notably, HCCS⁺ as well as H₃O⁺ and H₂COH⁺ have already been observed in the ISM (Cabezas, C. et al., 2022; Wootten et al., 1991; Bacmann et al., 2016; Ohishi et al., 1996)

3.1 Formation routes

While the discussion above on PA values suggests that all the species under consideration could be formed in the ISM, the next step is the investigation of their formation routes. For this purpose, three test cases have been chosen: Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and NCCNH⁺. According to the literature (Klippenstein et al., 2010; Zhu et al., 2019), the reaction between any molecule and ${{\text{H}}_3}^{+}$ is barrierless and proceeds through the formation of a stable intermediate characterized by a weak interaction between the protonated species and the H₂ molecule. This also applies to the reactive PESs considered in this work, which are shown in Figure 2. While a full account on their energetics is provided in the Supporting Material (SM), in the following, only the junChS+ZPE energies will be considered for our discussion.

FIGURE 2. Reactive PES for: the Cl₂ + ${{\text{H}}_3}^{+}$ reaction (panel A), the P₂ + ${{\text{H}}_3}^{+}$ reaction (panel B), and the CNNC + ${{\text{H}}_3}^{+}$ reaction (panel C). The full list of stationary-point energies is reported in the SM.

Let us start by considering the formation of Cl₂H⁺ from the reaction between ${{\text{H}}_3}^{+}$ and Cl₂ (panel A in Figure 2). The first intermediate formed, denoted as MIN1, is located at −143.2 kJ⋅mol⁻¹ with respect to the energy of the reactants. This minimum can interconvert to its spatial equivalent structure via TS1 (−85.7 kJ⋅mol⁻¹) and form Cl₂H⁺ and H₂ without any exiting barrier. A similar reaction mechanism is also observed for P₂, panel B of Figure 2. However, two different products can be formed in this case. Starting from the reactants the open routes are: a barrierless approach to MIN1 (−237.7 kJ⋅mol⁻¹) and a barrierless approach to MIN2 (−192.4 kJ⋅mol⁻¹). These two intermediates are connected by the transition state TS2, located at −191.8 kJ⋅ mol⁻¹. In addition, MIN1 can also form its spatial equivalent structure through TS1. MIN1 leads to product Prod1 without any exiting barrier, with Prod1 being the most stable stationary point of the ${{\text{HP}}_2}^{+}$ + H₂ reactive PES. Instead, MIN2 forms Prod2, i.e. the linear isomer of ${{\text{HP}}_2}^{+}$, which might be potentially present in the ISM, even though it lies ∼50 kJ⋅mol⁻¹ higher in energy than the cyclic form.

The last panel of Figure 2 (panel C) shows the reaction mechanism for the protonation of CNNC. In this case, there are three open barrierless approaches that lead to the formation of MIN1, MIN2 and MIN3. The first two intermediates are very close in energy: MIN1 lies at −305.5 kJ⋅mol⁻¹ and MIN2 at −301.8 kJ⋅mol⁻¹. They can interconvert through TS3, which differs from MIN2 for the orientation of the H₂ moiety with respect to the linear backbone. The motion of H₂ is a large amplitude one, TS3 being indeed characterized by an imaginary frequency of only 14 cm⁻¹. Consequently, the two structures are very close in energy, i.e. they differ by 0.1 kJ⋅mol⁻¹ at both the revDSD/junTZ and junChS levels of theory augmented by corresponding ZPEs. From Figure 2, it is noted that both MIN1 and MIN2 can form–without going through any barrier–Prod1, i.e. linear CNNCH⁺ and H₂. TS2, lying at −46.4 kJ⋅mol⁻¹, represents the barrier to be overcome for distorting the linear chain and forming MIN3, which can then interconvert with its equivalent structure via TS1 or form Prod2, i.e. a bent form of CNNCH⁺ + H₂.

The three protonation reactions of Figure 2 demonstrate that Cl₂H⁺, ${{\text{HP}}_2}^{+}$ and CNNCH⁺ can be formed in the ISM. The energetic characterization there reported is expected to have an accuracy of 1–5 kJ⋅mol⁻¹ (Alessandrini et al., 2019; 2021b) and can thus be used to guide the kinetic analysis of these reactions. While there are no doubts about the isomeric form produced in the reaction between Cl₂ and ${{\text{H}}_3}^{+}$, for the other two reactions, two possible products are present. However, the stability of MIN1 and the small energy difference between MIN1 and Prod1 should favor the formation of cyclic ${{\text{HP}}_2}^{+}$ and linear CNNCH⁺. Therefore, only these forms have been considered in the following spectroscopic characterization.

3.2 Spectroscopic characterization

To enable the experimental detection of the protonated species considered in this work, a computational spectroscopic characterization involving fundamental vibrational frequencies and rotational spectroscopy parameters has been carried out. In the following, the results are presented according to the class of species they belong: (i) protonated diatomics, where the experimental reference is offered by N₂H⁺ and HCO⁺; (ii) the linear CNNCH⁺ cation compared with the NCCNH⁺ counterpart; (iii) the protonated form of (E)-1,2-dicyanoethene. These species, with the only exception of the linear CNNCH⁺ molecule, are planar asymmetric rotors. Their planarity leads to a null dipole moment component along the c axis (μ_c = 0, with ab thus being the molecular plane). For ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$, the μ_a component is equal to zero as well because of symmetry reasons. For the linear rotor CNNCH⁺, the only non-null dipole moment component is that along the molecular axis, and it is simply denoted as μ.

3.2.1 Protonated diatomics

The rotational parameters of Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and HSi₂⁺ are reported in Table 2, while those of N₂H⁺ and HCO⁺ are provided in the SM. It should be pointed out that protonated CO and N₂ are linear rotors, and thus they are characterized by a single rotational constant and one quartic centrifugal distortion term. Differently, Cl₂H⁺, ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$ are asymmetric rotors, very close to the prolate limit (i.e., Ray’s asymmetry parameter, κ, being 0.999). Therefore, these species have three rotational constants and five quartic centrifugal distortion parameters.

TABLE 2. Rotational parameters (in MHz) of Watson’s S reduced Hamiltonian (I^r representation) for Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and ${{\text{HSi}}_2}^{+}$. Dipole moment components (in debye) are also reported.

The “CBS+CV+fT+fQ” level of theory augmented by vibrational corrections at the ae-CCSD(T)/cc-pwCVQZ level of theory, provides the rotational constant for N₂H⁺ with an error as low as 0.01%, as shown in Table S1. For HCO⁺, the error is slightly larger, this being 0.07%. However, both values are within the uncertainty predicted by benchmark studies for the “CBS+CV+fT+fQ” approach (Puzzarini et al., 2008). Even though N₂H⁺ and HCO⁺ are linear molecules, the accuracy reached for such species should hold also for the rotational constants of the bent Cl₂H⁺ cation and the cyclic ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$ species, at least when the mean error of the three rotational constants is considered.

To accurately simulate rotational spectra, other parameters are required. These are the NQCCs of the N and Cl atoms together with quartic centrifugal distortion constants for all species. As mentioned in Computational details section, the latter parameters have been computed at the ae-CCSD(T)/cc-pwCVQZ level of theory, while NQCCs have been obtained, on top of the “CBS+CV+fT+fQ” geometry, by computing the electric-field gradient at the quadrupolar nuclei at the ae-CCSD(T)/aug-cc-pwCVQZ level. The NQCCs have also been corrected for the vibrational contribution evaluated at the ae-CCSD(T)/cc-pwCVQZ level.

For N₂H⁺, the experimental D_J quartic centrifugal distortion constant is reproduced with a discrepancy of about 3%, while for HCO⁺ the error is 2%. Similar uncertainties are expected for the quartic centrifugal distortion constants of the protonated species of Table 2. Moving to NQCCs, for N₂H⁺, the discrepancy with respect to experiment is on average 1% and, thus, a similar accuracy is expected for the constants of Cl₂H⁺, as found in other works in the literature (Cazzoli et al., 2006; Puzzarini et al., 2010).

Finally, to simulate the rotational spectrum, the electric dipole moment components were computed at the same level of theory as NQCCs. The obtained values indicate that ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$ have a (very) weak b-type rotational spectrum, the only non-null dipole moment component (μ_b) being equal to 0.06 D for ${{\text{HSi}}_2}^{+}$ and to 0.13 for HP₂. In the case of ${{\text{HP}}_2}^{+}$, the intensity of the rotational lines is also affected by spin statistics. Indeed, ³¹P has a nuclear spin I of 1/2 (fermion) and the total wave function has to be anti-symmetric with respect to the exchange of the two equivalent P nuclei. The electronic and vibrational ground states being total symmetric, the overall symmetry of the total wave function is determined by the rotational and nuclear terms. Three ortho (symmetric) and one para (anti-symmetric) states arise for the nuclear wave function. Instead, the symmetry of the rotational contribution is determined by the K_a + K_c sum¹ because the nuclei exchange occurs by rotation around the b axis. As a consequence, the rotational transitions associated with odd K_a + K_c values are three times more intense than those with even K_a + K_c values.

The Cl₂H⁺ species has both a-type and b-type spectra, the corresponding dipole moment components being 2.28 D and 1.57 D, respectively. Thus, a-type transitions should have an intensity comparable to that of N₂H⁺ and HCO⁺, while b-type transitions are expected to be weaker. An example is shown in Figure 3, where the blue spectrum reports a μ_a transition of Cl₂H⁺ at around 98307 MHz and a weaker μ_b transition predicted at 98352 MHz. These are compared with the J = 1 ← 0 transition of N₂H⁺, for which both the theoretical simulation (this work, in red) and the experimental data (Cazzoli et al. (2012), in black) are considered. The experimental transition occurs at 93173 MHz and the theoretical simulation lies close in frequency, indeed being centered at 93181 MHz. In relative terms, this discrepancy should also apply to the simulation of Cl₂H⁺. Conservatively, based on benchmark studies (Alessandrini et al., 2018), an uncertainty of 0.1% can be expected for each rotational constant of Cl₂H⁺. This would locate the $J_{K_a,K_c}=7_{\mathrm{1,7}}\leftarrow 6_{\mathrm{1,6}}$ transition of Figure 3 in the range 98307.6 ± 18 MHz (3σ deviation). Lastly, the two insets of Figure 3 show the hyperfine splittings based on the computed NQCCs. For N₂H⁺, the structure reproduces very well the experimental one in terms of both intensity and shift with respect to the unperturbed frequency. The same should apply to Cl₂H⁺.

FIGURE 3. Theoretical (red) and experimental (black) J =1←0 transition of N₂H⁺ compared with the simulated spectrum of Cl₂H⁺ (blue) in the nearby frequency region. For the latter, both a-type $(J_{K_a,K_c}=7_{\mathrm{1,7}}\leftarrow 6_{\mathrm{1,6}})$ and b-type $(J_{K_a,K_c}=13_{\mathrm{1,13}}\leftarrow 14_{\mathrm{0,14}})$ transitions are visible. The two insets show the hyperfine structures due to nitrogen and chlorine nuclei for N₂H⁺ and Cl₂H⁺, respectively.

For Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and ${{\text{HSi}}_2}^{+}$, vibrational frequencies have been obtained at the ae-CCSD(T)/cc-pwCVQZ level. The results are reported in Table 3, while those of N₂H⁺ and HCO⁺ are provided in the SM, where they are also compared with the experimental counterparts. To understand the accuracy of the predictions for the protonated species addressed in this work, we start by analyzing the results for N₂H⁺ and HCO⁺. For the former cation, the lowest fundamental frequency is due to the bending mode. This is predicted at 700.94 cm⁻¹ within the harmonic approximation and about 9 cm⁻¹ lower in energy at the anharmonic level. The experiment locates this band at 698.6353 (14) cm⁻¹ (Foster and McKellar, 1984). The other two bands of N₂H⁺ were observed at 3233.96085 cm⁻¹ (Nakanaga et al., 1990) and 2257.8667 (13) cm⁻¹ (Owrutsky et al., 1986). Both of them are well predicted by the VPT2 treatment, with an error of only 8 cm⁻¹. For HCO⁺, the fundamental bands are predicted (incorporating anharmonicity) to lie at ν₁ = 3090.31 cm⁻¹, ν₂ = 2187.67 cm⁻¹, and ν₃ = 830.05 cm⁻¹. These are in good agreement with the experimental data, with discrepancies of about 2–4 cm⁻¹. In fact, the experimental observed frequencies are: ν₁ = 3088.7951 (31) cm⁻¹ (Amano, 1983), ν₂ = 2183.9496 (6) cm⁻¹ (Foster and McKellar, 1984), and ν₃ = 828.2305 (9) cm⁻¹ (Kawaguchi et al., 1985).

TABLE 3. Harmonic and anharmonic vibrational frequencies (cm⁻¹) for Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and ${{\text{HSi}}_2}^{+}$. Anharmonic intensities in km mol⁻¹. All quantities are at the ae-CCSD(T)/cc-pwCVQZ level of theory.

A similar accuracy is expected for the other protonated diatomics. The lowest fundamental band is observed for the stretching involving third-row atoms (ν₂ band for ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$, ν₃ band for Cl₂H⁺). This vibrational mode lies at 518.0 cm⁻¹ for Cl₂H⁺, 728.05 cm⁻¹ for ${{\text{HP}}_2}^{+}$, and 459.89 cm⁻¹ for ${{\text{HSi}}_2}^{+}$. The ν₃ band for ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$, and the ν₂ band for Cl₂H⁺ correspond to the asymmetric stretching involving the H atom, which is predicted at 949.02 cm⁻¹, 971.54 cm⁻¹, and 858.1 cm⁻¹ for ${{\text{HSi}}_2}^{+}$, ${{\text{HP}}_2}^{+}$, and Cl₂H⁺, respectively. The fundamental frequencies for symmetric stretching (ν₁) are, in the same order as before, 1435.82 cm⁻¹, 1659.00 cm⁻¹, and 2619.23 cm⁻¹. For the Cl₂H⁺ species, the most intense band is ν₁, while for ${{\text{HP}}_2}^{+}$ and ${{\text{HSi}}_2}^{+}$ it is ν₃. In conclusion, in view of their accuracy (Puzzarini et al., 2019; Barone et al., 2021)—which is conservatively expected to be better than 10 cm⁻¹—the data obtained in this work can allow the investigation of the vibrational spectrum of these species as well as the ro-vibrational ones².

3.2.2 CNNCH⁺ and NCCNH⁺

The rotational parameters of CNNCH⁺ are compared with those of NCCNH⁺ in Table 4. The error associated with our computational strategy can be derived by considering the experimental values of NCCNH⁺ (Gottlieb et al., 2000). For B₀ of NCCNH⁺, the discrepancy between experiment and theory is about 4 MHz, which means 0.09% in relative terms; thus, the value of the rotational constant of CNNCH⁺ lies within the interval 5225.59 MHz ± 4.7 MHz. For D_J of NCCNH⁺, the deviation of the computed value (0.480 × 10⁻³ MHz) from experiment is 9.3% and the same confidence range is expected to apply to CNNCH⁺. Moving to NQCCs, in both molecules, we have two nitrogen atoms which are denoted as the outer N and the inner N; while such denotation is rather straightforward for NCCNH⁺, in the case of CNNCH⁺, the outer N is the furthest from H, i.e., N2 in Figure 1. In NCCNH⁺, the NQCC of the outer N is well reproduced, the error being 0.2%, while a discrepancy of about 12% is noted for the inner NQCC, the experimental value being 0.250 MHz and the computed one 0.221 MHz. However, it has to be noted that the absolute discrepancy is nearly the same: ∼50 kHz for χ(outer) and ∼30 kHz for χ(inner). Therefore, the hyperfine splittings are well reproduced as shown in the left inset in Figure 4, where the experimental and theoretical hyperfine structures of the J = 3 ← 2 transition of NCCNH⁺ are depicted. Focusing on the unperturbed transition, the experimental frequency of NCCNH⁺ is 24 MHz lower than the simulated one. This deviation can be considered as the uncertainty affecting the J = 3 ← 2 transition of CNNCH⁺, which is shown in Figure 4 together with its hyperfine structure. Concerning the latter, an analogous accuracy to that observed for NCCNH⁺ is expected. Indeed, in absolute terms, the uncertainty affecting NQCCs should be of the order of 50 kHz. Since the dipole moment of CNNCH⁺ is large (about 6 D) and only 0.4 D smaller than that of NCCNH⁺, the detection of this new protonated form is reasonable and future works will benefit from the present data.

TABLE 4. Rotational parameters and fundamental vibrational frequencies of CNNCH⁺ and NCCNH⁺.

FIGURE 4. Theoretical (red) and experimental (black) J =3←2 transition of NCCNH⁺ compared with the simulation of the same transition of CNNCH⁺ (blue). The insets show in detail the corresponding hyperfine structures.

Table 4 also reports the fundamental frequencies and intensities of the vibrational modes for both NCCNH⁺ and CNNCH⁺. In this case, no experimental reference is available in the literature, but vibrational frequencies computed at the ae-CCSD(T)/cc-pwCVTZ level are sufficiently accurate to guide the assignment of the (ro-)vibrational spectrum for both species (Puzzarini et al., 2019; Barone et al., 2021). The most intense vibrational band is ν₃ for CNNCH⁺ and ν₁ for NCCNH+.

3.2.3 NCCHCHCNH⁺

The last molecule considered in this work is the E isomer of NCCHCHCNH⁺. According to the computed electric dipole moment components, a-type transitions are the most prominent features of the rotational spectrum because the associated dipole component is as large as 10.1 D. Based on the literature on this topic (Puzzarini et al., 2008; Puzzarini, 2016), in relative terms, the uncertainty affecting the equilibrium “CBS+CV” rotational constants is expected to be 0.6%, which drastically reduces to about 0.1% when vibrational corrections are incorporated. Actually, this uncertainty can be as small as 0.05%, as demonstrated, for example, by Melosso et al. (2022); Alessandrini et al. (2021a). For an accurate simulation of the rotational spectrum of this species, the NQCCs of both nitrogen atoms have been computed, on top of the “CBS+CV” geometry, at the fc-CCSD(T)/jun-cc-pVTZ level. Quartic centrifugal distortion constants have been obtained from a harmonic force field at the same level of theory. All these data, collected in Table 5, can guide future experimental measurements as well as astronomical search of this molecule in view of its large dipole moment (in particular that associated with a-type transitions) and the small uncertainties affecting our spectroscopic parameters. Notably, the observation of this species in the ISM would indirectly confirm the feasibility of the reaction between the CN radical and vinyl cyanide, as suggested in Marchione et al. (2022).

TABLE 5. Computed spectroscopic parameters of (E)-NCCHCHCNH⁺.

In Table 5, the vibrational fundamental frequencies of (E)-NCCHCHCNH⁺ are reported. It is noted that the strongest band is the one associated with the stretching of the terminal H atom. The weakest band has an intensity of 0.6 km⋅mol⁻¹ and it is the out-of-plane bending of the CN and CNH terminal groups. According to the literature on this topic (Barone et al., 2015), double-hybrid functionals should reproduce the experimental fundamental frequencies with a mean error of about 10.4 cm⁻¹.

4 Conclusion

In the ISM conditions, protonation can occur easily thanks to the abundance of ${{\text{H}}_3}^{+}$. For this reason, protonated forms of “invisible” (to rotational spectroscopy and thus to radioastronomy) apolar species are considered tracers of the latter and can be used to derive an estimate of their abundance. In this work, we suggest the use of the protonated species ${{\text{HP}}_2}^{+}$, Cl₂H⁺, ${{\text{HSi}}_2}^{+}$, CNNCH⁺, and (E)-NCCHCHCNH⁺ as proxies for the corresponding apolar form (P₂, Cl₂, Si₂, CNNC, NCCHCHCN). The most stable structures of these species have been preliminary investigated at the revDSD/junTZ level of theory. This pointed out that Cl₂H⁺ has a bent structure, while ${{\text{HP}}_2}^{+}$ (in contrast with N₂H⁺ which is linear) and ${{\text{HSi}}_2}^{+}$ have a cyclic structure. CNNCH⁺ is linear in analogy to the NCCNH⁺ isomer, which has already been characterized and detected in the ISM. Finally, a planar structure has been found for (E)-NCCHCHCNH⁺.

State-of-the-art computational methodologies have been used to assess the PA values of Cl₂, P₂, Si₂, CNNC, and (E)-NCCHCHCN to be compared with that of H₂. Indeed, a PA greater than that of H₂ implies that protonation due to ${{\text{H}}_3}^{+}$ is favored. This is the case for all the species considered in this study, and the computed PA values have been found in good agreement with experimental and previous theoretical estimates available in the literature.

To provide more details on the protonation mechanisms, the formation routes leading to Cl₂H⁺, ${{\text{HP}}_2}^{+}$, and CNNCH⁺ have been explored. The junChS+ZPE energies, evaluated on top of revDSD/junTZ geometries, pointed out the presence of simple reaction paths starting from the corresponding apolar species and ${{\text{H}}_3}^{+}$. All these mechanisms are dominated by a barrierless step from reactants to adducts, but also from these latter to products. In all the three cases, the minimum energy path leads to the most stable protonated forms, namely bent Cl₂H⁺, cyclic ${{\text{HP}}_2}^{+}$, and linear CNNCH⁺. The first two formation routes are similar to those producing N₂H⁺ and HCO⁺; therefore, it can be safely assumed that a similar mechanism leads to ${{\text{HSi}}_2}^{+}$. We also expect that the formation of (E)-NCCHCHCNH⁺ presents similarity with that of CNNCH⁺.

Since the investigated protonated species are potentially present in the ISM, an accurate computational spectroscopic characterization has also been carried out. For protonated diatomics, the rotational constants are expected to reproduce the experimental counterparts with an error of 0.1%, or even smaller as in the case of HCO⁺ and N₂H⁺. For CNNCH⁺ and (E)-NCCHCHCNH⁺, because of the reduction in the level of theory employed, a larger uncertainty has been obtained, but still on the order of 0.1%. For all the species investigated, a full set of rotational parameters, also including quartic centrifugal distortion constants and NQCCs, has been obtained and used to estimate low-J rotational transitions. For CNNCH⁺, the J = 3 ← 2 rotational transition is predicted to be centered at around 30353.5 MHz with an uncertainty of about 24 MHz. The spectroscopic characterization of (E)-NCCHCHCNH⁺ is of particular interest because it results from protonation of the E isomer of 1,2-dicyanoethene, which in turn is the main product of the reaction between the CN radical and acrylonitrile (Marchione et al., 2022). This reaction is of great astrochemical relevance, but the only possible option for detecting (E)-NCCHCHCN is via its protonated proxy. This work reports–for the first time–the rotational constants of this species and points out that protonation leads to a large dipole moment of about 10 D, which further suggests the detectability of (E)-NCCHCHCNH⁺. Therefore, our simulations will be useful to support and complement experimental measurements as well as to guide possible radioastronomical observation.

For all the aforementioned molecules, the fundamental vibrational frequencies have accurately been obtained by exploiting a methodology that incorporates the effects of anharmonicity. Therefore, the expected accuracy of our predictions is 10 cm⁻¹ or even better, which is largely sufficient for the assignment of experimental infrared bands.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

SA worked on conceptualization of the work as well as investigation. SA, LB, and MM carried out the formal analysis and wrote a first draft of the manuscript. SA and CP reviewed the manuscript and provided the funding. Final editing was carried out by CP, LB, MM, and SA.

Funding

This work has been supported by MIUR (PRIN Grant Number 202082CE3T) and by the University of Bologna (RFO funds). The COST Action CA21101 “COSY—Confined molecular systems: from a new generation of materials to the stars” is also acknowledged. This work was partially supported by the Open Access Publishing Fund of the Scuola Normale Superiore.

Acknowledgments

The SMART@SNS Laboratory (http://smart.sns.it) is acknowledged for providing high-performance computing facilities.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2023.1128896/full#supplementary-material

Footnotes

¹The asymmetric rotor energy levels are labeled using (in addition to the J quantum number) two “pseudo” quantum numbers: K_a and K_c, which correspond to the K values of the two limiting prolate (referred to the a symmetry axis) and oblate (referred to the c symmetry axis) cases of the symmetric rotor, respectively.

²For those interested in the prediction of the rotational constants of vibrational excited states, the required ${\alpha }_r^{\gamma }$ values are available in Supplementary Table S2.

References