The Role of State-of-the-Art Quantum-Chemical Calculations in Astrochemistry: Formation Route and Spectroscopy of Ethanimine as a Paradigmatic Case.

The gas-phase formation and spectroscopic characteristics of ethanimine have been re-investigated as a paradigmatic case illustrating the accuracy of state-of-the-art quantum-chemical (QC) methodologies in the field of astrochemistry. According to our computations, the reaction between the amidogen, NH, and ethyl, C2H5, radicals is very fast, close to the gas-kinetics limit. Although the main reaction channel under conditions typical of the interstellar medium leads to methanimine and the methyl radical, the predicted amount of the two E,Z stereoisomers of ethanimine is around 10%. State-of-the-art QC and kinetic models lead to a [E-CH3CHNH]/[Z-CH3CHNH] ratio of ca. 1.4, slightly higher than the previous computations, but still far from the value determined from astronomical observations (ca. 3). An accurate computational characterization of the molecular structure, energetics, and spectroscopic properties of the E and Z isomers of ethanimine combined with millimeter-wave measurements up to 300 GHz, allows for predicting the rotational spectrum of both isomers up to 500 GHz, thus opening the way toward new astronomical observations.

1. Introduction

Until the latter part of the 20th century, it was believed that the vastness of interstellar space consisted mainly of hydrogen atoms. Apart from molecular hydrogen and a few other simple diatomic species, the harsh conditions of the interstellar medium (ISM) were thought to be incompatible with polyatomic molecules exhibiting even a small degree of complexity. This idea began to be questioned roughly fifty years ago with the discovery of the first ‘complex’ molecules (formaldehyde in 1969 [1], methanol in 1970 [2], and formic acid in 1971 [3]). Since these first hints, the pace of molecular detections in space has accelerated. The original paradigm alluded to above has now been completely erased by the discovery of more than two hundred molecular species in the ISM and circumstellar shells [4]. Molecules detected in space range from simple hydrides (such as H2, H2O, and NH3), to hydrocarbon chain species (such as the cyanopolyynes, HCN), to simple organics (such as alcohols and aldehydes), to highly reactive species such as ions and radicals (e.g., H2COH+, HC3NH+, and CH3O), to polycyclic aromatic hydrocarbon molecules (PAHs) [5]. The latter are the most abundant polyatomic species in space, although no individual PAH has yet been identified. This short and generic list makes it clear that interstellar chemistry is rich and diverse. Focusing on atomic elements, the contribution of hydrogen and helium amounts to about 98%; therefore, that of heavier elements (such as carbon, nitrogen, and oxygen) is only about 2%. Nevertheless, as mentioned above and in the following, this small fraction of heavy elements makes possible a great variety of chemical compounds.

Among the molecules identified in the ISM, in the present context, the interest is on the so-called “interstellar complex organic molecules” (iCOMs), namely C-bearing molecules with at least six atoms [6], such as glycolaldehyde [7] and acetamide [8]. In this respect, molecular spectroscopy plays a crucial role in finding out whether a particular molecule is present in an astronomical environment and, if so, in what abundance (e.g., References [9,10,11,12,13,14]). Indeed, most of the gas-phase molecules mentioned above have been discovered thanks to radioastronomy, and thus by means of their rotational signatures. However, still unassigned features appear in radio-astronomical spectra, thus pointing out that we are far from a complete census of the interstellar molecules. Even though iCOMs were first detected decades ago, the processes that lead to their production are still a matter of debate. While it is nowadays accepted that they can be synthesized either on dust-grain surfaces or by means of gas-phase chemistry [15,16,17,18,19,20], the available literature reveals that the understanding of chemical evolution is at a rather primitive stage and current chemical models still show large deficiencies. In particular, in the last two decades, we have witnessed the growth of the grain-surface chemistry and the progressive abandonment of the gas-phase reactivity. However, recent detections of iCOMs in cold regions [17,21], where the surface chemistry cannot be efficient, have pointed out that gas-phase chemistry has been overlooked. The idea at the basis of grain-surface chemistry is that radical species trapped in icy mantles move and encounter a reaction partner, thus giving rise to rich chemistry [22]. This mechanism is called the Langmuir–Hinshelwood mechanism. A second possibility is offered by gaseous atoms and molecules that directly hit atoms and molecules on dust grains to form products. This is called the Eley-Rideal mechanism, which is –however– less important in interstellar clouds [22]. It is thus evident that surface mobility is an important factor in grain-surface reactions. While at temperatures above 30 K all molecular species have a reasonable mobility, at lower temperatures the mechanisms described above cannot be efficient (in particular, at T = 10 K only hydrogen atoms have a non-negligible mobility).

On the one hand, the detection of iCOMs reveals an undisputed molecular complexity in the universe, which is—however—not entirely disclosed. On the other hand, there is still much to be understood about how iCOMs, and more generally molecules, can be formed in the harsh conditions typical of the ISM and how they can evolve toward more complex species. It is thus clear that we are facing two major challenges in astrochemistry, namely the detection of ‘new’ iCOMs and the understanding of chemical reactivity, and that they are strongly interconnected [23]. For this reason, we present an integrated experimental and theoretical strategy that address both challenges [24]: state-of-the-art computational approaches are combined with laboratory spectroscopy. Ethanimine, a prebiotic iCOM (as will be explained later), has been chosen as case study to illustrate this strategy.

Computational chemistry offers an invaluable aid in providing two main pieces of information, namely:

The investigation of reactive potential energy surfaces (PESs) from both energetic (thermochemistry) and kinetic points of view. Two possibilities can be actually envisaged: (i) starting from purposely chosen precursors the formation route of the sought product (i.e., a molecule already identified in the ISM) is derived; (ii) starting from small reactive species the possible pathways are elaborated. In both cases, the harsh conditions of the ISM (extremely cold—down to 10 K—regions, extremely low density—even 10 particles/cm) are used as constraints.

Referring to point (1.ii), the accurate prediction of the spectroscopic parameters of those products that can be of interest and for which such information is still missing is carried out.

The second piece of information (point (2)) is then used to guide and support laboratory measurements that, in the field of rotational spectroscopy, are mandatory to further proceed toward astronomical searches. Finally, astronomical data can confirm the plausibility of the developed formation pathway. For example, the identification (in astronomical surveys) of the sought product of the reaction investigated can provide an indirect confirmation. Astronomical observations can also be used to support the effectiveness of a gas-phase formation route [25]. Furthermore, as in the case of ethanimine, that is, in the presence of two different isomers, the computed branching ratios can be compared with the relative astronomical abundance.

As briefly mentioned above, ethanimine has been chosen as a prototypical molecular system for presenting and discussing our strategy. Its astrophysical relevance is related to its prebiotic potential as a possible precursor of amino acids [26], like alanine, by reaction with HCN and H2O, or with formic acid [27,28,29]. Both isomers of ethanimine have been detected in Sagittarius B2 North, SgrB2 (N), as a result of the GBT-PRIMOS project (this acronym standing for Green Bank Telescope–PRebiotic Interstellar MOlecule Survey) [29]. Furthermore, this molecule can also be of interest in the study of the nitrogen-based Titan atmosphere [30]. As is well known, Titan is considered a good approximation of the primitive Earth [31,32], thus explaining its relevance in the study of prebiotic chemistry in space.

1.1. The Thermochemistry-Kinetics-Spectroscopy Strategy

As explained above, our strategy can start from a potential iCOM or from a potential gas phase reaction between two molecules already detected in the ISM. The example addressed in this contribution belongs to the first case. In any event, the starting point is the design of a feasible and accessible reactive PES leading to the species of interest in the first case or starting from the chosen reactants in the latter option. This preliminary investigation is performed at a relatively low computational level, relying on density functional theory (DFT). Usually, different pathways can be derived and only those that are feasible in the conditions typical of the astronomical environment under consideration are further investigated at a higher level [33]. Therefore, the more promising pathways are re-investigated in order to improve the structural determination of the stationary points as well as to compute the energetics at the state of the art. To conclude, kinetic calculations are carried out, these providing the conclusive information on the feasibility of the suggested mechanisms and on the branching ratios of the reaction products [34,35]. On the spectroscopic side, the target species are better characterized in terms of structural, molecular and spectroscopic properties in order to provide reliable predictions for the laboratory spectroscopy studies [24,36]. Finally, astronomical observations can be performed and might be used to support the outcomes of the theoretical investigation (see e.g., Reference [35,37]).

In the case of ethanimine, several pieces of information and reliable data were already available prior to our investigation. This was indeed the reason why it was chosen as case study to develop and test our strategy. Ethanimine was already detected in the ISM, but the knowledge on its rotational spectrum was limited. Possible gas-phase mechanisms were also available in the literature, even if not entirely satisfactory from an astronomical point of view.

1.1.1. Formation Pathway

As far as the formation mechanism of ethanimine in the ISM is concerned, among the gas-phase reactions considered by Quan et al. [26], the NH + C2H5 reaction is certainly the best candidate, because the analogous NH + CH3 reaction was demonstrated to be the dominant formation route of methanimine in the model proposed by Suzuki et al. [38]. Recently, Balucani and coworkers [39] performed an exhaustive search of different possible reaction paths, concluding that only the following three channels are open in the ISM conditions:

While the computational level of ref. [39] might not be sufficient to obtain quantitative results, the computed energy barriers governing the other considered paths are so high that they can be safely excluded in the present study, with the only exception of channel 3:

The reaction of CH2 with CH2NH was also proposed [40]. However, the computational level employed to analyze the mechanism was quite low and only the singlet form of CH2 leads to submerged barriers, thus requiring a spin-orbit coupling mechanism from the more stable triplet form. Furthermore, the first step of this reaction is the abstraction of one hydrogen atom from CH2NH by CH2. Therefore, the effective reaction involves the CH3 and CHNH radicals:

While the CH3 radical has been detected in several regions of the ISM, only the CHNH+ cation has been identified so far, the corresponding neutral radical remaining still unobserved. As a consequence we have not investigated in detail this alternative mechanism.

The rate coefficients employed by Quan et al. [26] in their model are 8.25 Å cm/s for NH + C2H5→E−CH3CHNH + H (R2a) and 2.75 Å cm/s for NH + C2H5→ Z−CH3CHNH + H (R2b) with no temperature dependence. However, in their model, they did not take into consideration that the channel (R1) consumes a fraction of the reactants, thus reducing the total production rate of ethanimine. The lack of temperature dependence being quite surprising, the computations of ref. [39] indeed pointed out that the rate constants for channels R2a and R2b are very similar to one another and both of them nearly double when increasing the temperature from 10 K to 200 K (i.e., the interval considered in Reference [26]).

1.1.2. Spectroscopic Characterization

As clear from the preceding section and the literature cited, ethanimine exists in two isomeric forms, E and Z (see Figure 1), with the former being the most stable. As already mentioned, both isomers have been detected in the ISM [29]. However, the experimental works on ethanimine rotational spectrum that allowed such an identification were limited to low frequencies, that is, below 140 GHz [29,41,42]. Since extrapolations from low-frequency laboratory measurements usually provide inaccurate predictions for higher frequencies and, in view of the extended astronomical observatory facilities provided by the Atacama Large Millimeter/submillimeter Array (ALMA; working frequency range: 84–950 GHz), the extension of the investigation of the rotational spectra for both ethanimine isomers was warranted. To improve the spectroscopic characterization of both isomeric forms of CH3CHNH, an accurate computational study of the energetics as well as structural and spectroscopic parameters was performed [43], thereby focusing on the description of the methyl internal rotation and on the centrifugal distortion terms.

Figure 1

Structures and selected geometrical parameters of all stationary points on the ethanimine PES optimized at the B2 level. Bond lengths are in Å, whereas HNC (), HCH () or HNH () valence angles and CCNH () dihedral angles are in degrees.

In Reference [43], an accurate characterization of the vibrational (infrared, IR) spectrum was also provided. Indeed, in view of predicting the gas-phase IR spectrum, the computed fundamental frequencies of Reference [43] can be considered more accurate and reliable than those from the low-resolution gas-phase study of Reference [44] and the matrix investigation reported in Reference [45]. In fact, both studies provided either incomplete or contradictory results. Therefore, the computational work on the IR spectra of both isomers of ethanimine might be of great help in future observation of planetary atmospheres (e.g., that of Titan) relying on, for instance, the James Webb Space Telescope to be launched in 2021.

1.2. Organization of the Paper

In the following section the computational methodology at the basis of both challenges, namely the investigation of the formation pathway and the spectroscopic characterization, will be described in detail. On these grounds, the purpose of the present paper is to provide reliable and accurate rate constants for the reaction mechanisms considered by exploiting state-of-the-art quantum-chemical (QC) and kinetic computations. At the same time, we reconsider the spectroscopic parameters of ethanimine, which were previously reported in Reference [43]. In the subsequent section, the thermochemical, kinetic and spectroscopic results will be presented and discussed. Finally, concluding remarks will address the effectiveness of our strategy.

2. Computational Methodology

Our QC strategy relies on a composite scheme in order to obtain high accuracy at an affordable computational cost. This combines coupled-cluster (CC) techniques with extrapolation to the complete basis set (CBS) limit and consideration of core-valence (CV) correlation contributions. To account for vibrational effects, also including anharmonicity, we resort to second-order vibrational perturbation theory (VPT2) [46]. For the evaluation of kinetic parameters, a master equation (ME) approach based on ab initio transition state theory (AITSTME) [47] is used.

2.1. Electronic Structure Calculations

The QC calculations described in the following have been performed with Gaussian 16 [48] and CFOUR [49] quantum-chemical packages. The former code has been employed for DFT and VPT2 computations, while the latter for MP2 and CC calculations. Furthermore, the MRCC code [50] interfaced to CFOUR has been used to perform the CC calculations including quadruple excitations.

2.1.1. Reactive Potential Energy Surface

A preliminary characterization of the reactive PES has been carried out at the B3LYP-D3/aug-cc-pVTZ level [51,52,53,54] (hereafter shortly denoted as B3), with geometry optimizations combined with harmonic force-field calculations in order to characterize the stationary points. More generally, for this introductory analysis of the reactive PES, a double-zeta quality basis set can be safely employed in conjunction with the hybrid B3LYP functional. The results obtained at the B3 level were found already in agreement with those reported in literature [39].

Since the double-hybrid B2PLYP functional [55] represents a good compromise between accuracy and computational cost for both structural and spectroscopic properties [56,57,58] and is able to provide a reliable description of reactive and non-reactive PESs [34,35,56,57,59,60,61,62], our strategy relies on it in order to provide the final description of the stationary points from a structural point of view. For this purpose, B2PLYP is usually employed in conjunction with a triple-zeta basis sets incorporating diffuse functions. In the present case, the aug-cc-pVnZ (n = T,Q) basis sets have been employed. In addition, different basis sets have been tested. However, these as well as the aug-cc-VQZ set have not pointed out any substantial modification of either the PES or geometrical parameters with respect to the combination of B2PLYP and aug-cc-pVTZ. Furthermore, to account for dispersion interactions, the DFT-D3 scheme [63] has been used coupled to the Becke-Johnson (BJ) damping function [64]. Overall, in the following as well as in Figure 1, we only refer to the B2PLYP-D3(BJ)/aug-cc-pVTZ level of theory, hereafter simply denoted as B2. The interested reader can find a detailed structural analysis in the Supplementary Information (SI) material. In passing, it should be noted that the NH species has been considered only in the most stable triplet state () [65].

The subsequent step of our strategy requires the accurate estimates of reaction energies and barriers. For this purpose, we rely on composite schemes based on coupled-cluster (CC) theory and exploited on top of B2 optimized geometries. In particular, two different types of approach have been employed: (i) schemes entirely based on CC calculations (denoted as “CBS+CV”, “CBS+CV+fT+pQ”, and “HEAT-like”), and (ii) the so-called “cheap” model, which involves an important reduction of the computational cost, while keeping high accuracy. All these approaches are detailed in Section 2.1.3.

Finally, the inclusion of zero-point vibrational energies (ZPVEs) is required. These can be computed both at harmonic and anharmonic levels. In this respect, as briefly mentioned above, the nature of the stationary points has been checked by the evaluation of analytical Hessians (harmonic force fields), which –in turn– provide harmonic ZPVEs. Anharmonic force fields have also been calculated at the B2 level, thus allowing the evaluation of anharmonic ZPVEs by means of the generalized vibrational-perturbation theory (GVPT2). With respect to standard VPT2, GVPT2 excludes resonant terms from the perturbative treatment and reintroduces them through a reduced-dimensionality variational calculation [36,46,66,67,68].

2.1.2. Spectroscopic Parameters

As mentioned in the Introduction, most of the molecular species detected in the ISM have been identified thanks to their rotational fingerprints. Indeed, the rotational spectrum of a given molecule is so characteristic that the observation of its rotational transitions provides an unequivocal proof of the presence of that species in the environment under consideration. For this reason, rotational spectroscopy is the spectroscopic technique taken into consideration in our strategy.

The leading terms of rotational spectroscopy are the rotational constants, the rotational transition frequencies mainly depending on their values [69]. In general, to provide accurate predictions of rotational transitions, it is mandatory to go beyond the rigid-rotor approximation, thus accounting for centrifugal-distortion and vibrational effects. To incorporate the latter in the prediction of rotational constants, vibrational corrections () are added to the equilibrium rotational constants (). Within VPT2, the vibrational ground-state rotational constants () are given by the expression [70]: with i denoting the inertial axis (i = a, b, c; that is, for instance, = ). The summation runs over all vibrational normal modes and the ’s are the vibration-rotation interaction constants.

It is important to note that the vibrational corrections provide a small contribution to ’s, indeed accounting only for ∼1–3% of the total value. However, as pointed out in refs. [71,72], the inclusion of vibrational corrections is fundamental for quantitative predictions. From a computational point of view, the equilibrium rotational constants are straightforwardly derived from the equilibrium structure, while vibrational corrections require anharmonic force-field calculations [73]. In view of the fact that provides the largest contribution to , the major computational effort is put on geometry optimizations. To fulfil the accuracy requirements, composite schemes are thus employed. In particular, the CBS+CV approach mentioned above (detailed in Section 2.1.3) has been chosen. Vibrational corrections have been determined at the B2 level (however employing a smaller basis set, namely a modified version of maug-cc-pVTZ [74,75]).

To account for the effect of centrifugal distortion on rotational transitions, the computation of the corresponding parameters is required. The quartic centrifugal-distortion constants are obtained as a by product of harmonic force-field computations [73], while the sextic ones are derived from the knowledge of the cubic force field. Therefore, the latter terms have been evaluated at the B2 level. In Reference [43], the harmonic force field has been computed at a higher level of theory, thereby employing the CCSD(T) method in conjunction with the cc-pCVQZ basis set, with all electrons being correlated (all). The acronym stands for single and double excitations and a perturbative treatment of triple excitations [76].

In the specific case of ethanimine, to complete the characterization of the rotational spectrum, nitrogen quadrupole coupling constants and the spectroscopic parameters related to the methyl internal rotation are needed. To evaluate the former constants, the electric field gradient (EFG) tensor is required and it has been evaluated, for both isomers, at the all-CCSD(T)/cc-pCVQZ level using the CBS+CV reference geometry. These have been further augmented by vibrational corrections at the B2 level. For all details related to the conversion of the EFG elements to the nuclear quadrupole coupling constants, interested readers are –for example– referred to Reference [73]. To model the methyl internal rotation, the barrier for the hindered rotation has been derived as energy difference between the E-/Z-CH3CHNH minima and the corresponding transition state (see Figure 2 and Reference [43] for details).

Figure 2

Methyl internal rotation in E-ethanimine

2.1.3. Composite Schemes

The composite schemes described in the following rely on or start from the CCSD(T) method. The coupled cluster diagnostic [77] has been found to be smaller than 0.02 for closed-shell molecules and smaller than 0.04 for open-shell species (see Table S3 in the Supplementary Materials), thus confirming that the non-dynamical correlation is negligible for all the systems investigated.

The CC-Based Approaches

The composite scheme denoted as CBS+CV accounts for the extrapolation to the CBS limit and CV corrections at the CCSD(T) level. Within this approach, electronic energies are obtained as follows:

The terms on the right-hand side correspond, respectively, to: (1) the Hartree-Fock (HF-SCF) electronic energy extrapolated to the CBS limit using the exponential form proposed by Feller [78] combined with the cc-pVnZ basis sets [54] (n = T,Q,5); (2) the frozen-core (fc) CCSD(T) correlation energy extrapolated to the CBS limit employing the two-point formula proposed by Helgaker et al. [79] in conjunction with the cc-pVnZ sets, with either the n = T,Q or n = Q,5 combination of bases; (3) the CV correlation correction is calculated as energy difference between all and fc CCSD(T), both with the same basis set. The cc-pCVTZ and cc-pCVQZ [80] bases have been employed depending on whether the extrapolation to the CBS limit has been done using the n = T,Q or n = Q,5 combinations of basis sets, respectively.

To further improve the energetics, a HEAT-like approach can be (and indeed has been) exploited, the reference of this scheme being the HEAT protocol [81]. The latter is known to provide sub-kJ mol accuracy. The HEAT-like model is obtained by incorporating additional terms into Equation (2): where is the CBS+CV energy (Equation (2)). Similarly to the CV contribution, corrections due to the full treatment of triple excitations, , and the perturbative description of quadruples, , are computed as energy differences between CCSDT [82,83,84] and CCSD(T) and between CCSDT(Q) [85,86,87] and CCSDT (all of them within the fc approximation) employing the cc-pVTZ and cc-pVDZ basis sets, respectively. Inclusion of fT and pQ contributions leads to the definition of the CBS+CV+fT+pQ scheme. The diagonal Born-Oppenheimer correction [88,89,90,91], , and the scalar relativistic contribution to the energy, , [92,93] are also included in the HEAT-like approach. The former corrections are usually computed at the HF-SCF/aug-cc-pVnZ level [53], with n = D,T, and the relativistic corrections are obtained at the all-CCSD(T)/aug-cc-pCVnZ level (n = D,T) and include the (one-electron) Darwin and mass-velocity terms. In the present study, the last two corrections have been obtained using double-zeta basis sets, even if their convergence with respect to contributions calculated with triple-zeta basis sets has been checked for a few stationary points.

To exploit the high accuracy of the CBS+CV energy approach (Equation (2)) for molecular structure determinations, the energy-gradient composite scheme introduced in refs. [94,95] can be employed. The resulting energy gradient to be employed in geometry optimizations is given by the following expression: where the terms on the right-hand side are defined as in Equation (2). Analogously to Equation (3), the energy gradient of Equation (4) can be further improved by including the contributions due to full treatment of triples and quadruple excitations.

The ChS Approach

In the present work, electronic energies have also been computed by resorting to a computationally less demanding strategy, the so-called “cheap” scheme (hereafter denoted as ChS). Our strategy relies on this model when larger systems are involved. In the present context, it has been employed to further demonstrate its high accuracy and reliability, which was already pointed out in References [24,33,61].

This approach has been originally developed for the determination of accurate structural parameters [96,97] and then extended to energy evaluations [98]. In this scheme, the fc-CCSD(T)/cc-pVTZ energy is taken as the starting point and augmented by contributions accounting for the extrapolation to the CBS limit and CV correlation. To keep the computational cost limited, contrary to the CBS+CV approach, these terms are computed using second-order Møller-Plesset perturbation theory (MP2) [99]:

The second term on the right-hand side is the contribution due to the extrapolation to the CBS limit, which is performed in two steps (i.e., by separating the HF-SCF and correlation contributions) as in the CBS+CV scheme, and employing the same extrapolative formulas. MP2 calculations are carried out in conjunction with the cc-pVTZ and cc-pVQZ basis sets. The third term is evaluated, as in Equation (2), as ‘all-fc’ energy difference at the MP2/cc-pCVTZ level.

2.2. Kinetic Models

Phenomenological rate coefficients have been obtained within the AITSTME approach using the Rice-Ramsperger-Kassel-Marcus (RRKM) 1D master equation system solver (MESS) code [100]. For channels possessing a distinct saddle point, rate coefficients have been determined by conventional transition state theory (TST) within the rigid-rotor harmonic-oscillator (RRHO) approximation. Instead, rate constants for barrierless elementary reactions have been computed employing phase space theory (PST) [101,102], again within the RRHO assumption, and including tunneling as well as non classical reflection effects by using the Eckart model [103].

For what concerns PST, it provides a useful, and easily implemented, reference theory for barrierless reactions, and it is largely used in computational kinetics applied to astrochemistry (see, for example, References [35,37,104,105]). Its assumption is that the interaction between the two reacting fragments is isotropic and does not affect the internal fragment motions. This approximation is especially valid if the dynamical bottleneck lies at large separations where the interacting fragments have free rotations and unperturbed vibrations. This is generally true for low temperature phenomena, as is the case for the ISM. The isotropic attractive potential is assumed to be described by the functional form , where the coefficient C is obtained by fitting the electronic energies obtained at various long-range distances of fragments. To be more precise, we performed a relaxed scan of the HN−CH2CH3 distance and then fitted the corresponding minimum energy path to a function, thus obtaining a C value of Eh.

3. Results and Discussion

The thermochemical aspects of the gas-phase reaction between the ethyl and imidogen radicals are first discussed, and then the attention is focused on the results of RRKM kinetic simulations. Subsequently, the spectroscopic characterization of ethanimine is addressed.

3.1. The NH + C2H5 Reaction: Thermochemistry

As described in the Introduction, we have characterized the most favorable reaction channels following the addition of NH to C2H5, namely (R1), (R2a), (R2b) and (R3). The resulting pathways are presented in Figure 3.

Figure 3

Reaction mechanisms for the NH + C2H5 reaction: complete basis set (CBS)+core-valence (CV) energies augmented by B2 anharmonic zero-point vibrational energies (ZPVEs).

The initial association of reactants leading to the M1 intermediate is highly exothermic, with M1 easily interconverting to M2 through a small energy barrier (TS1). Since the channel leading to CH3 + CH2NH requires the smallest number of steps and the lowest barrier (TS2; breaking of the C−C bond), thermodynamic considerations should favor this process. The evolution of M1 to E−CH2CHNH + H and of M2 to Z−CH2CHNH + H requires overcoming the transition states TS3 and TS4, respectively, which imply the breaking of one C−H bond and are less than 10 mol−1 higher in energy than TS2. This means that multiple competitive pathways can come into play. Since the E ↔ Z isomerization involves a high energy transition state, that is, TS-isom, once ethanmine is formed, the direct interconversion of the two isomers in the gas phase is very unlikely. As expected, the predicted stability of E−CH3CHNH by mol−1 with respect to the Z isomer is in agreement with the value of mol−1 by Melli et al. [43], the level of theory being nearly the same. Analogously, the E ⇄ Z isomerization barrier of mol−1 (calculated with respect to E−CH2CHNH) matches that of ref. [43], while the CCSD(T)/aug-cc-pVTZ value reported in ref. [39] results to be overestimated. Besides leading to methanimine or ethanimine, M1 could potentially undergo multiple hydrogen transfers, according to the path reported in green in Figure 3, to end up with the release of NH2 and C2H4. However, the multiple rearrangements and the relatively high barrier associated to the M3 - M4 conversion make this process the less feasible in the harsh conditions of the ISM.

The electronic energies of all the species involved in the paths shown in Figure 3, referred to the isolated reactants, are collected in Table 1, where harmonic and anharmonic ZPVE corrections are also given. In addition to CBS+CV, CBS+CV+fT+pQ, HEAT-like and ChS energies, Table 1 also reports B2 results and, for the sake of comparison, the results from Reference [39].

Table 1

Relative electronic energies together with ZPVE corrections. Values in mol−1.

	B2	CBS+CV a	CBS+CV+fT+pQ b	HEAT-Like b	ChS	Harm-ZPE c	Anharm-ZPE c	ΔH00 d	ΔH00 e
NH + C2H5	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
M1	−345.16	−352.10	−351.92	−352.09	−352.82	28.40	28.13	−323.97	−311
		(−352.01)	(−351.83)	(−352.00)					(-323)
M2	−342.24	−349.43	−349.27	−349.44	−350.28	28.94	28.59	−320.84	−308
		(−349.33)	(−349.18)	(−349.35)					(−320)
M3	−382.90	−388.91	-	-	−390.12	31.07	30.08	−358.83	−340
		(−388.89)							(−357)
M4	−336.75	−345.98	-	-	−347.39	28.74	28.30	−317.68	−301
		(−345.71)							-
TS1	−338.10	−345.27	−345.11	−345.28	−346.14	28.51	28.17	−317.10	−304
		(−345.20)	(−345.04)	(−345.21)					(−316)
TS2	−223.85	−226.32	−229.11	−229.22	−227.56	18.47	18.43	−207.89	−197
		(−226.58)	(−229.37)	(−229.48)					(−211)
TS3	−195.51	−202.17	−203.93	−204.09	−202.88	9.31	8.85	−193.32	−180
		(−202.41)	(−204.18)	(−204.34)					(−197)
TS4	−192.92	−199.72	−201.47	−201.62	−200.42	9.48	9.07	−190.65	−177
		(−199.95)	(−201.70)	(−201.85)					(−194)
TS-isom	−94.55	−98.47	−98.17	−98.40	−100.55	−4.13	−4.00	−102.47	−84
		(−98.81)	(−98.51)	(−98.74)					(−101)
TS5	−190.43	−197.33	-	-	−199.01	18.79	18.19	−179.14	−165
		(−197.27)							-
TS6	−176.42	−182.49	-	-	−184.40	17.70	16.87	−165.62	−146
		(−182.51)							-
TS7	−236.01	−237.83	-	-	−238.77	20.67	20.63	−217.20	−207
		(−238.27)							-
CH2NH + CH3	−245.85	−251.18	−251.36	−251.34	−253.69	7.24	7.47	-243.71	-230
		(−251.38)	(−251.57)	(−251.55)					(−243)
E−CH3CHNH + H	−213.82	−222.30	−221.99	−221.92	−224.10	4.00	4.14	−218.16	−202
		(−222.65)	(−222.34)	(−222.27)					(−217)
Z−CH3CHNH + H	−210.94	−219.52	−219.22	−219.13	−221.31	4.06	4.20	−215.32	−199
		(−219.85)	(−219.55)	(−219.46)					(−214)
C2H4 + NH2	−242.67	−248.70	-	-	−250.94	8.26	8.64	-240.06	−228
		(−248.78)							(−241)

CBS+CV scheme: the n = T,Q set of bases for the extrapolation to the CBS limit and cc-pCVTZ for the the CV contribution. Within parentheses, the results for the n = Q,5 set (CBS) and cc-pCVQZ (CV). CBS+CV: n = T,Q set (CBS) and cc-pCVTZ (CV). Within parentheses, n = Q,5 (CBS) and cc-pCVQZ (CV). Relative ZPVE corrections at the B2 level. CBS+CV electronic energies augmented by anharmonic ZPVE corrections. Data from ref. [39]: values at the CCSD(T)/aug-cc-pVTZ level, within parentheses W1 values. ZPVE corrected values.

From the inspection of Table 1, it is evident that the CBS+CV energies obtained by extrapolating the correlation energy using either cc-pVTZ/cc-pVQZ or cc-pVQZ/cc-pV5Z are in close agreement, the largest difference being around mol−1 for TS7. Incorporation of the fT and pQ contributions provides negligible effects, the only exceptions being TS2, TS3 and TS4, for which deviations of about 2– 3 mol−1 are noted. As expected from our previous experience [24,33,61], the results obtained by using the ChS model are accurate. In fact, they are within mol−1 from the corresponding CBS+CV values, and within 2 mol−1 with respect to CBS+CV+fT+pQ results. Finally, by incorporating the DBOC and relativistic corrections, the HEAT-like values have been obtained. These negligibly differ from the corresponding CBS+CV+fT+pQ results, the sum of the two contributions being on the order of 0.1– mol−1.

In Table 1, the values (second last column) correspond to the CBS+CV energies augmented by anharmonic B2 ZPVEs. It is noted that they are in agreement with those calculated in Reference [39] at the W1 level (last column), even if the latter are available for a reduced number of species. On the other hand, the CCSD(T)/aug-cc-pVTZ level (last column), employed in Reference [39] for all stationary points, shows deviations as large as 19 mol−1 with respect to our CBS+CV results, with an average overestimation of about mol−1. It is particularly interesting that the B2 model chemistry performs better than CCSD(T)/aug-cc-pVTZ computations, indeed the largest deviation of the former from CBS+CV data is 9 mol−1, systematically overestimating the latter results.

3.2. The NH + C2H5 Reaction: Rate Coefficients

As described in Section 2.2, for the reaction paths shown in Figure 3, channel specific rate constants have been computed using AITSTME. For these calculations, the CBS+CV energies, corrected by anharmonic ZPVEs, have employed for all intermediate and transition states of the PES. The rate coefficients have been evaluated in the 10–300 K temperature range and at pressure of 1 × 10−12 atm, the results being collected in Table 2. The corresponding temperature dependence plots are shown in Figure 4.

Table 2

Product-formation rate constants (in cm molecule s) at 1 × 10−12 atm as a function of the temperature.

T (K)	CH3 + CH2NH	E−CH3CHNH	Z−CH3CHNH	C2H4 + NH2
10	1.13 × 10−10	8.87 × 10−12	6.29 × 10−12	7.19 × 10−13
20	1.28 × 10−10	8.45 × 10−12	6.35 × 10−12	5.98 × 10−13
30	1.38 × 10−10	8.64 × 10−12	6.61 × 10−12	5.90 × 10−13
40	1.45 × 10−10	8.91 × 10−12	6.85 × 10−12	6.00 × 10−13
50	1.50 × 10−10	9.18 × 10−12	7.08 × 10−12	6.14 × 10−13
60	1.55 × 10−10	9.41 × 10−12	7.27 × 10−12	6.27 × 10−13
70	1.59 × 10−10	9.67 × 10−12	7.47 × 10−12	6.45 × 10−13
80	1.63 × 10−10	9.90 × 10−12	7.64 × 10−12	6.61 × 10−13
90	1.66 × 10−10	1.01 × 10−11	7.80 × 10−12	6.76 × 10−13
100	1.69 × 10−10	1.03 × 10−11	7.95 × 10−12	6.90 × 10−13
110	1.71 × 10−10	1.05 × 10−11	8.08 × 10−12	7.03 × 10−13
120	1.74 × 10−10	1.07 × 10−11	8.21 × 10−12	7.16 × 10−13
130	1.76 × 10−10	1.08 × 10−11	8.33 × 10−12	7.28 × 10−13
140	1.78 × 10−10	1.10 × 10−11	8.44 × 10−12	7.39 × 10−13
150	1.80 × 10−10	1.11 × 10−11	8.55 × 10−12	7.50 × 10−13
160	1.82 × 10−10	1.12 × 10−11	8.65 × 10−12	7.60 × 10−13
170	1.84 × 10−10	1.14 × 10−11	8.75 × 10−12	7.70 × 10−13
180	1.86 × 10−10	1.15 × 10−11	8.84 × 10−12	7.80 × 10−13
190	1.87 × 10−10	1.16 × 10−11	8.93 × 10−12	7.90 × 10−13
200	1.89 × 10−10	1.18 × 10−11	9.02 × 10−12	7.99 × 10−13
210	1.90 × 10−10	1.19 × 10−11	9.11 × 10−12	8.09 × 10−13
220	1.92 × 10−10	1.20 × 10−11	9.19 × 10−12	8.18 × 10−13
230	1.93 × 10−10	1.21 × 10−11	9.27 × 10−12	8.26 × 10−13
240	1.94 × 10−10	1.22 × 10−11	9.35 × 10−12	8.35 × 10−13
250	1.96 × 10−10	1.23 × 10−11	9.42 × 10−12	8.44 × 10−13
260	1.97 × 10−10	1.24 × 10−11	9.49 × 10−12	8.52 × 10−13
270	1.98 × 10−10	1.25 × 10−11	9.57 × 10−12	8.61 × 10−13
280	1.99 × 10−10	1.26 × 10−11	9.64 × 10−12	8.69 × 10−13
290	2.00 × 10−10	1.27 × 10−11	9.71 × 10−12	8.78 × 10−13
300	2.01 × 10−10	1.28 × 10−11	9.77 × 10−12	8.86 × 10−13

Figure 4

Rate coefficients as a function of temperature for the four products.

The kinetic analysis shows that, in agreement with previous computations, only the channels leading to methanimine (R1), E-ethanimine (R2a), and Z-ethanimine (R2b) are significant. The contribution of the channel leading to C2H4 + NH2 is indeed negligible because of the limitations explained above. Focusing on the temperature dependence of the rate constant, an interesting feature can be noted at very low temperature (i.e., 10 K). In fact, at such low temperatures, the trend of k for the channel leading to C2H4 + NH2 shows a larger variation with temperature than the other paths. This can be explained as a quantum tunneling effect because the TSs of this path are higher in energy and the rearrangements involve H transfer.

As mentioned in Section 2.1.1, the two reactants, CH3CH2 and NH, are in a spin-doublet and a spin-triplet state, respectively. As a consequence, this results in a total of six spin-microstates, which can be classified in reactive (overall: spin-doublet state) and non-reactive (overall: spin-quartet state). Since the number of the reactive spin-microstates is 2 and that of the non-reactive ones is 4, a statistical factor of 1/3 must be applied to the total reaction rate constant in order to account for this ratio. This has been indeed taken into account in our calculations.

The branching ratios at selected temperatures are shown in Table 3. Is clear that, at all the temperatures considered, the channel leading to CH3 + CH2NH is by far the dominant one, indeed accounting for more than 87% of the total yield. The ratio between the E- and Z- isomers of ethanimine spans from 1.3 to 1.4 in the entire analyzed range. Even if this value is greater than that obtained by Balucani et al. [39], it is still underestimated with respect to the factor of ca. 3 derived from the observation of Loomis et al. [29]. While our theoretical estimate is difficult to improve (because already at the state of the art), this discrepancy deserves to be addressed in some details. First of all, the astronomical value is affected by a large uncertainty mostly due to the fact that, for the analysis of E-ethanimine features, the authors had to fix its excitation temperature to that derived for Z-CH3CHNH, this resulting in an approximate value of the column density of the E species (see Reference [29] for a thorough account). Second, a non negligible contribution of dust-grain chemistry to the production of ethanimine cannot be excluded.

Table 3

Product branching ratios at various temperatures.

Branching Ratios	CH3 + CH2NH	E−CH3CHNH + H	Z−CH3CHNH + H	C2H4 + NH2
10 K	87.7%	6.9%	4.9%	0.6%
100 K	89.9%	5.5%	4.2%	0.4%
300 K	89.6%	5.7%	4.3%	0.4%

Finally, kinetic calculations have also been performed using the CBS+CV+fT+pQ energies, corrected by anharmonic B2 ZPVEs. In these calculations, the reaction channel R3 has been excluded, since it has been demonstrated that it gives non relevant contribution to the total reactive flux. The results are almost identical to those obtained exploiting CBS+CV energies. For the sake of completeness, they are only reported in the Supplementary Materials.

3.3. E-/Z-CH3CHNH: Spectroscopic Characterization

In Reference [43], an integrated experimental-theoretical strategy for the spectroscopic characterization required for guiding astronomical searches of molecular signatures in space has been presented together with its application to ethanimine. To exploit the interplay of experiment and theory, we have relied on the VMS-ROT software [106].

Despite the fact that state-of-the-art QC computations are able to provide predictions for rotational transitions with an accuracy better than 0.1% (see, e.g., Reference [71] and Table 4), this is not generally sufficient for guiding astronomical searches and/or assignments, thus requiring experimental determinations of the corresponding spectroscopic parameters. In fact, when dealing with broadband unbiased astronomical surveys, rotational transition frequencies need to be known with an accuracy preferably better than 100 kHz, an accuracy that can be easily obtained from experimental studies.

Table 4

Computed and experimental rotational parameters (values in MHz) of ethanimine.

	E-CH3CHNH			Z-CH3CHNH
	Theory		Experiment ^a	Theory		Experiment ^a
	Best Estimates ^b	B2 ^c		Best Estimates ^b	B2 ^c
A_0	53,178.26	53,398.97	53,120.561(30)	50,002.63	50,288.21	49,964.87(93)
B_0	9780.14	9744.02	9782.7720(47)	9831.98	9781.22	9832.4823(96)
C_0	8702.82	8679.22	8697.0263(46)	8652.81	8621.21	8646.0305(94)
Δ_J	6.48×10^−3	6.39×10^−3	6.4641(49)×10^−3	6.99×10^−3	6.96×10^−3	6.938(13)×10^−3
Δ_K	0.568	0.575	0.5763(34)	0.468	0.480	0.468 ^d
Δ_JK	−0.0165	−0.0165	−0.01403(21)	−0.0163	−0.0182	−0.01219(23)
δ_J	1.09×10^−3	1.06×10^−3	1.1033(19)×10^−3	1.25×10^−3	1.24×10^−3	1.2657(65)×10^−3
δ_K	−0.0535	−0.0518	−0.06709(59)	−0.0522	−0.0464	−0.0642(19)
V_3^e	563.1	488.8	566.37(20)	523.3	446.4	517.41(33)

Watson A-reduction. Values in parenthesis denote one standard deviation and apply to the last digits of the constants. Equilibrium CBS+CV rotational constants augmented by B2 vibrational corrections. Quartic centrifugal distortion constants at the CCSD(T)/cc-pCVQZ level. For details, see [43]. This work. Fixed at the best-estimated theoretical value. values in cm−1.

The spectroscopic parameters of the two isomers of ethanimine have been computed as explained in Section 2.1.2 and, in more details, in Reference [43]. Table 4 collects a selection of computed spectroscopic parameters, which are compared with their experimental counterparts (see ref. [43] for an extended version of this table). In order to improve the prediction of the rotational spectrum, the computed rotational constants have been replaced by those available from Reference [42]. The resulting simulated spectra, shown in Figure 5, have then used to guide rotational spectrum measurements, thus leading to a straightforward assignment of the rotational transitions up to 300 GHz [43].

Figure 5

Simulation of the rotational spectra of E- (top panel) and Z-CH3CHNH (bottom panel) at T = 100 K obtained with the VMS-ROT software. The absolute values of the dipole moment components are also reported. The red circles highlight the frequency ranges for which measurements were available prior to the investigation of Reference [43].

Back to Table 4, there is excellent agreement between the computed and experimental rotational constants, with the maximum discrepancy (i.e., 0.1%) noted for the constant. For quartic centrifugal distortion constants, a quantitative agreement between theory and experiment is apparent as well. In Table 4, the determination of the term is also reported. As briefly mentioned in Section 2.1.2, this is due to the internal rotation of the −CH3 group (see Figure 2 and Figure 5), which leads to a periodic potential energy along the torsional angle, with being the energy barrier (the reader is referred to Reference [43] for a detailed account). This parameter is important in the prediction and analysis of the rotational spectrum because, the barrier being finite, a tunneling effect occurs. This leads to a splitting of the threefold degeneracy into two levels, a nondegenerate A level and a doubly degenerate E level. As a consequence, rotational transitions in the A and E torsional sublevels are observed, and their separation depends on the value of .

In Table 4, for comparison purposes, the results at the B2 level are also reported. Despite the reduced computational cost, the B2 parameters are accurate: the discrepancies for the rotational constants range, in relative terms, from 0.2% to 0.5%. The worsening in the prediction of the centrifugal distortion constants (with respect to the all-CCSD(T)/cc-pCVQZ level) is nearly negligible and, for the barrier, the disagreement is well within 1 mol−1. These results are particularly significant for the study of larger systems, which are currently not amenable to the most accurate computations reported in the present study.

4. Conclusions

Astrochemistry is an interdisciplinary field involving chemistry, physics, and astronomy, and is strongly interconnected to astrobiology. The main aim of astrochemistry is to understand the chemical evolution of the universe, thus investigating the formation and destruction of molecules in space, but also their interaction with radiation and their feedback on physics of the environments. In the present contribution, the paradigmatic case of ethanimine has been used to drive the readers in a fascinating journey among the different computational tools presently available to investigate two major challenges in astrochemistry—the discovery of new molecules in space and the understanding of how they can be formed. Ethanimine has been chosen because of its prebiotic potential, thus being a first step toward the demonstration of the presence of biological building-block molecules in the ISM.

Because of difficulties in experimentally mimicking the extreme conditions that characterize the interstellar medium, accurate state-of-the-art computational approaches play a fundamental role. We have demonstrated how state-of-the-art QC methodologies are able to accurately describe reactive PESs, to be complemented by kinetic calculations, and predict the spectroscopic parameters of the reaction products. Even if our state-of-the-art investigation on the gas-phase formation of ethanimine is not able to perfectly reproduce the [E−CH3CHNH]/[Z−CH3CHNH] ratio of ca. 3 issuing from astronomical observations, the result obtained (1.4) is of the correct order of magnitude and it can be considered as a satisfactory value in view of the complexity of the chemistry in space, involving several reactions at the same time both in the gas phase and on dust-grains. From a spectroscopic point of view, it has been shown how state-of-the-art QC methodologies can guide, support and complement rotational spectroscopy studies, the latter providing the rest frequencies required by astronomical searches.