Protonation of apolar species: From Cl2H+ to (E)-NCCHCHCNH+ through computational investigations

Radioastronomy is a powerful tool for the discovery of molecules in space but it requires molecular species to be polar. The observation of apolar species can be however enabled by protonation, which occurs from reaction with the abundant H 3 + ion whenever the proton affinity of the species under consideration is greater than that of H2. This property can be easily investigated by computational chemistry and, in this work, it has been used to asses the potential protonation of simple homo diatomics, such as Cl2, P2, and Si2, as well as apolar species containing two equivalent CN moieties, such as diisocyanogen (CNNC) and (E)-1,2-dicyanoethene. Quantum chemistry has also been exploited to investigate the mechanisms of three protonation reactions of H 3 + with Cl2, P2, and CNNC. To support laboratory measurements and astronomical observations of the resulting transient species, their rotational spectroscopic parameters were accurately computed together with fundamental vibrational frequencies. For this purpose, we have employed CCSD(T)-based computational methodologies, which provide equilibrium structures with errors smaller than 0.001 Å and 0.1° for bond distances and angles, respectively. Such an accuracy is expected to lead to rotational constants predicted, in relative terms, with uncertainties better than 0.2%. Instead, the expected accuracy on vibrational frequencies is ∼10 cm−1, thus being well suited to guide band assignments.


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 6 H 5 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 6 H 6 (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 H 3 + , which in turn is obtained from the reaction between H 2 and H 2 + , with the latter being formed via cosmic ray ionization (Yamamoto, 2017). The 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 2 . Based on the PA values reported by Yamamoto (2017), H 3 + is the primary source of N 2 H + and HCO + /HOC + , but it does not lead to the formation of O 2 H + because O 2 has a smaller PA than H 2 . 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 2 , Si 2 , and Cl 2 , 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 2 H + molecule-observed for the first time in 1974 by Turner (1974)-is considered a good tracer of N 2 , which is otherwise not observable at radio frequencies.
For this reason, the study of the simple HP 2 + , HSi 2 + , and Cl 2 H + ions can potentially provide useful data to guide laboratory and/or astronomical observations. While the P 2 and Si 2 precursors might be considered somewhat exotic, Cl 2 is a simple species stable under Earth conditions. To assess the accuracy of our work, the N 2 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,2dicyanoethene, which is the main product of the reaction between the CN radical and vinyl cyanide (CH 2 = CH − CN, Marchione et al. (2022)), also known as acrylonitrile. While the (Z) isomer of 1,2dicyanoethene, 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.
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.

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;, 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.

Proton affinity and formation routes
To understand if a protonated molecule can be formed in the ISM starting from H 3 + , its PA value has to be considered in comparison with that of the H 2 molecule. The PA at T = 0 K is computed as (McNaught and Wilkinson, 1997;Marino et al., 2000): where, in this study, both ΔH 0 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)/juncc-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 −3 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 2 H + , 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.

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 , corresponding to A 0 for γ = a, B 0 for γ = b, and C 0 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 is straightforwardly obtained from a geometry optimization and accounts for about 99% of the value of B γ 0 (Puzzarini et al., 2008;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 vibrationrotation interaction constants (α γ r ). In Eq. 4, d r denotes the degeneracy of the rth vibrational mode. The vibrational correction represents a minor contribution to B γ 0 , but its consideration is crucial to reach accurate computational results, comparable with experimental ones (Puzzarini et al., 2008;Alessandrini et al., 2018).
Since B γ e is the major contribution to B γ 0 , it should be computed at the highest possible level of theory when aiming at high accuracy. For this reason, for Cl 2 H + , HP 2 + , 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;, 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. dE ∞ (HF−SCF) dx and dΔE ∞ (CCSD(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 −3 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 aeand 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 vibrationrotation 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 1 , d 2 ) (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 2 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.

Results
As starting point, the work investigated the stable form of Cl 2 H + , HP 2 + , and 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.
The main molecule driving protonation in the ISM being H 3 + , the first mandatory step is to investigate whether the species considered in this work can be produced by reaction 1. The first

FIGURE 1
Equilibrium structure, point group symmetry and electronic ground state of (1) Cl 2 H + , (2) HP 2 + , (3) HSi 2 + , (4) CNNCH + , and (5) protonated screening is based on PA. If PA of P 2 , Si 2 , Cl 2 , CNNC, and NCCHCHCN is greater than that of H 2 , the protonation occurs and the reaction proceeds at the Langevin rate (Yamamoto, 2017). The latter is typically around 10 -9 cm 3 s −1 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 2 , CO, and NCCN, whose protonated forms have already been observed in the ISM. Table 1 indicates a PA value of 418 kJ⋅mol −1 for H 2 , in good agreement with what reported by Yamamoto (2017) and Hunter and Lias (1998), i.e. 422.61 kJ⋅mol −1 . A similar agreement is noted for N 2 and CO, with values of 493.5 kJ⋅mol −1 and 591.1 kJ⋅mol −1 , respectively, to be compared with 493.8 kJ⋅mol −1 and 594 kJ⋅mol −1 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 2 has been obtained for Cl 2 , thus suggesting the possible formation of Cl 2 H + in the ISM. A further support to this is provided by the fact that PA(Cl 2 ) is intermediate between PA(N 2 ) and PA(CO), with both N 2 H + and HCO + detected in the ISM. In particular, our PA(Cl 2 ) value of 547.6 kJ⋅mol −1 is in good agreement with that established experimentally by Cacace et al. (1998), i.e. 548.1 ± 13 kJ⋅mol −1 .
The PA value of P 2 is 669.9 kJ⋅mol −1 and it is very close to that of NCCN, i.e., 652.7 kJ⋅mol −1 , thus confirming the plausible formation of the protonated form of P 2 . Finally, large PA values have been obtained for CNNC, (E)-NCCHCHCN, and Si 2 . For the first two species, the PA values are 717.8 kJ⋅mol −1 and 735.5 kJ⋅mol −1 , respectively, which are comparable with those of molecules like H 2 O, H 2 S, and H 2 CO (Yamamoto, 2017). The Si 2 species shows the largest PA, 841.0 kJ⋅mol −1 , this result being likely due to its triplet electronic ground state. Indeed, a similar value (869 kJ⋅mol −1 ) is obtained for the CCS radical in the same electronic state (Yamamoto, 2017). Notably, HCCS + as well as H 3 O + and H 2 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)

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 2 H + , HP 2 + , and NCCNH + . According to the literature (Klippenstein et al., 2010;Zhu et al., 2019), the reaction between any molecule and 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 2 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. Let us start by considering the formation of Cl 2 H + from the reaction between H 3 + and Cl 2 (panel A in Figure 2). The first intermediate formed, denoted as MIN1, is located at −143.2 kJ⋅mol −1 with respect to the energy of the reactants. This minimum can interconvert to its spatial equivalent structure via TS1 (−85.7 kJ⋅mol −1 ) and form Cl 2 H + and H 2 without any exiting barrier. A similar reaction mechanism is also observed for P 2 , 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 −1 ) and a barrierless approach to MIN2 (−192.4 kJ⋅mol −1 ). These two intermediates are connected by the transition state TS2, located at −191.8 kJ⋅ mol −1 . 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 HP 2 + + H 2 reactive PES. Instead, MIN2 forms Prod2, i.e. the linear isomer of HP 2 + , which might be potentially present in the ISM, even though it lies ∼50 kJ⋅mol −1 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 −1 and MIN2 at −301.8 kJ⋅mol −1 . They can interconvert through TS3, which differs from MIN2 for the orientation of the H 2 moiety with respect to the linear backbone. The motion of H 2 is a large amplitude one, TS3 being indeed characterized by an imaginary frequency of only 14 cm −1 . Consequently, the two structures are very close in energy, i.e. they differ by 0.1 kJ⋅mol −1 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 2 . TS2, lying at −46.4 kJ⋅mol −1 , 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 2 .
The three protonation reactions of Figure 2 demonstrate that Cl 2 H + , 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 −1 (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 2 and 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 HP 2 + and linear CNNCH + . Therefore, only these forms have been considered in the following spectroscopic characterization.

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 2 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 HP 2 + and 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 μ.

Protonated diatomics
The rotational parameters of Cl 2 H + , HP 2 + , and HSi 2 + are reported in Table 2, while those of N 2 H + and HCO + are provided in the SM. It should be pointed out that protonated CO and N 2 are linear rotors, and thus they are characterized by a single rotational constant and one quartic centrifugal distortion term. Differently, Cl 2 H + , HP 2 + and 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.
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 2 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 2 H + and HCO + are linear molecules, the accuracy reached for such species should hold also for the rotational constants of the bent Cl 2 H + cation and the cyclic HP 2 + and 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 2 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 2 H + , the discrepancy with respect to experiment is on average 1% and, thus, a similar accuracy is expected for the constants of Cl 2 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 HP 2 + and 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 HSi 2 + and to 0.13 for HP 2 . In the case of HP 2 + , the intensity of the rotational lines is  also affected by spin statistics. Indeed, 31 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 1 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 2 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 2 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 2 H + at around 98307 MHz 1 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. and a weaker μ b transition predicted at 98352 MHz. These are compared with the J = 1 ← 0 transition of N 2 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 2 H + . Conservatively, based on benchmark studies (Alessandrini et al., 2018), an uncertainty of 0.1% can be expected for each rotational constant of Cl 2 H + . This would locate the J K a ,K c = 7 1,7 ← 6 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 2 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 2 H + . For Cl 2 H + , HP 2 + , and 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 2 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 2 H + and HCO + . For the former cation, the lowest fundamental frequency is due to the bending mode. This is predicted at 700.94 cm −1 within the harmonic approximation and about 9 cm −1 lower in energy at the anharmonic level. The experiment locates this band at 698.6353 (14) cm −1 (Foster and McKellar, 1984). The other two bands of N 2 H + were observed at 3233.96085 cm −1 (Nakanaga et al., 1990) and 2257.8667 (13) cm −1 (Owrutsky et al., 1986). Both of them are well predicted by the VPT2 treatment, with an error of only 8 cm −1 . For HCO + , the fundamental bands are predicted (incorporating anharmonicity) to lie at ν 1 = 3090.31 cm −1 , ν 2 = 2187.67 cm −1 , and ν 3 = 830.05 cm −1 . These are in good agreement with the experimental data, with discrepancies of about 2-4 cm −1 . In fact, the experimental observed frequencies are: ν 1 = 3088.7951(31) cm −1 (Amano, 1983), ν 2 = 2183.9496(6) cm −1 (Foster and McKellar, 1984), and ν 3 = 828.2305(9) cm −1 (Kawaguchi et al., 1985). A similar accuracy is expected for the other protonated diatomics. The lowest fundamental band is observed for the stretching involving third-row atoms (ν 2 band for HP 2 + and HSi 2 + , ν 3 band for Cl 2 H + ). This vibrational mode lies at 518.0 cm −1 for Cl 2 H + , 728.05 cm −1 for HP 2 + , and 459.89 cm −1 for HSi 2 + . The ν 3 band for HP 2 + and HSi 2 + , and the ν 2 band for Cl 2 H + correspond to the asymmetric stretching involving the H atom, which is predicted at 949.02 cm −1 , 971.54 cm −1 , and 858.1 cm −1 for HSi 2 + , HP 2 + , and Cl 2 H + , respectively. The fundamental frequencies for symmetric stretching (ν 1 ) are, in the same order as before, 1435.82 cm −1 , 1659.00 cm −1 , and 2619.23 cm −1 . For the Cl 2 H + species, the most intense band is ν 1 , while for HP 2 + and HSi 2 + it is ν 3 . In conclusion, in view of their accuracy Barone et al., 2021)-which is conservatively expected to be better than 10 cm −1 -the data obtained in this work can allow the investigation of the vibrational spectrum of these species as well as the ro-vibrational ones 2 .

FIGURE 3
Theoretical (red) and experimental (black) J = 1 ← 0 transition of N 2 H + compared with the simulated spectrum of Cl 2 H + (blue) in the nearby frequency region. For the latter, both a-type (J K a ,K c = 7 1,7 ← 6 1,6 ) and b-type (J K a ,K c = 13 1,13 ← 14 0,14 ) transitions are visible. The two insets show the hyperfine structures due to nitrogen and chlorine nuclei for N 2 H + and Cl 2 H + , respectively.

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 0 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 −3 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 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 a Equilibrium rotational constants from the "CBS+CV+fT+fQ" approach. Dipole moment and NQCCs computed, on top of the CBS+CV+fT+ fQ reference geometry, at the ae-CCSD(T)/aug-cc-pwCVTZ level. Vibrational corrections to the previous quantities and quartic centrifugal distortion constants at the ae-CCSD(T)/cc-pwCVTZ level of theory. Anharmonic frequencies and intensities at the ae-CCSD(T)/cc-pwCVTZ level of theory. b Experimental data taken from Gottlieb et al. (2000). c The symmetry of the vibrational modes is given in parentheses.

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. level of theory augmented by vibrational corrections at the revDSD/junTZ level. Quartic centrifugal distortion constants at the fc-CCSD(T)/junTZ level. Dipole moments (in debyes) and NQCCs computed at the same level of theory but on top of the "CBS+CV" geometry and vibrationally corrected at the revDSD/junTZ level. For further details, see text. Watson's S reduction in I r representation is used. b The symmetry of the normal modes is given in parentheses. Vibrational frequencies and intensity computed at the revDSD/junTZ level of theory using GVPT2.
sufficiently accurate to guide the assignment of the (ro-)vibrational spectrum for both species Barone et al., 2021). The most intense vibrational band is ν 3 for CNNCH + and ν 1 for NCCNH+.

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)/juncc-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).
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 −1 and it is the out-ofplane bending of the CN and CNH terminal groups. According to the literature on this topic , double-hybrid functionals should reproduce the experimental fundamental frequencies with a mean error of about 10.4 cm −1 .

Conclusion
In the ISM conditions, protonation can occur easily thanks to the abundance of 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 HP 2 + , Cl 2 H + , HSi 2 + , CNNCH + , and (E)-NCCHCHCNH + as proxies for the corresponding apolar form (P 2 , Cl 2 , Si 2 , 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 2 H + has a bent structure, while HP 2 + (in contrast with N 2 H + which is linear) and 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 2 , P 2 , Si 2 , CNNC, and (E)-NCCHCHCN to be compared with that of H 2 . Indeed, a PA greater than that of H 2 implies that protonation due to 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 2 H + , 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 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 2 H + , cyclic HP 2 + , and linear CNNCH + . The first two formation routes are similar to those producing N 2 H + and HCO + ; therefore, it can be safely assumed that a similar mechanism leads to 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 2 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 −1 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.