Evidence for tunnelling motion in the rotational spectrum of the aminomethyl radical CH₂NH₂

Abstract

The aminomethyl radical CH₂NH₂ has been proposed to be one of the key radical intermediates in the astrophysical and atmospheric chemistry of methylamine. Its spectroscopy, however, has never been characterized by rotationally resolved experiments, impeding its potential detection in the interstellar medium by radio astronomy. We report the first rotational spectrum of CH₂NH₂ recorded and modelled between 300 and 750 GHz. The broadband acquisition of the pure rotational transitions reveals an unambiguous inversion tunnelling motion of the –CH₂ and –NH₂ groups in the radical that splits the ground vibrational state of CH₂NH₂ into two substates coupled by a Coriolis-type interaction. The careful analysis of the spectrum opens the possibility of searching for CH₂NH₂ in the interstellar medium.

---

1 Introduction

Amines are important N-bearing molecules in the Earth atmosphere, the interstellar medium (ISM), and possibly the atmospheres of exoplanets. Methylamine (CH₃NH₂), the simplest primary amine, has been detected in several high-mass star-forming regions^1–6 in the ISM, as well as in Titan tholins⁷ and in comets.⁸ Methylamine has a unique importance in astrochemistry as it is proposed to be the precursor of interstellar glycine, the simplest amino acid.⁹

The chemistry of methylamine formation and degradation in terrestrial atmosphere and in the ISM has been widely examined.^10,11 Evidences have been brought to indicate the existence of two key intermediates, the aminomethyl (CH₂NH₂) and methyl amino (CH₃NH) radicals. In addition to the crucial role of both radicals in the formation and destruction processes of methylamine, several astrochemical models have investigated their role in the formation of glycine.^12–15 For example, quantum chemistry calculations proposed that at a temperature below 55 K, a favoured gas phase formation pathway of glycine involves a barrier-less radical coupling between CH₂NH₂ and HCO.¹⁶ In interstellar environments, in turn, several formation pathways of CH₂NH₂ have been proposed, including successive hydrogenation of HCN on ice grains,^12,15,17 and the H-abstraction reaction from methylamine initiated by OH, CN or NH₂ radicals both on ice grains and in gas-phase.^14,18 On Earth, the atmospheric oxidation of methylamine undergoes similar H-abstraction pathways leading to the formation of aminomethyl radicals: CH₃NH₂ + OH → CH₂NH₂/CH₃NH + H₂O.^19,20

Despite its important role in methylamine chemistry, the spectroscopic properties of CH₂NH₂ is poorly characterized in the literature. The UV photoelectron spectrum of CH₂NH₂ was reported,^21,22 and its reaction kinetics with O₂ was studied using photoionization mass spectrometry.²³ Its ionization energy was calculated,²⁴ and its 0 K heat of formation was derived from dissociative photoionization of ethylenediamine.²⁵ The most comprehensive theoretical investigation of CH₂NH₂ so-far is provided by Puzzarini et al.,¹⁸ where the authors performed high-level quantum chemical calculations on the structure and spectroscopic parameters of both CH₂NH₂ and CH₃NH. Recently, Joshi and Lee observed the infrared spectrum of CH₂NH₂ in a solid para-H₂ matrix at 3.2 K during the reaction between methylamine and H atoms, providing the first direct spectroscopic evidence of CH₂NH₂ formation in solid phase at the low temperature of the ISM.²⁶ To our best knowledge, no further experimental information on CH₂NH₂ is available.

On the aspect of spectroscopy, CH₂NH₂ is also an interesting molecule whose –CH₂ and –NH₂ groups are subject to large amplitude motions (LAMs), which refer to the internal motions in the molecule that deviates significantly from the harmonic oscillator approximation.²⁷ Its internal dynamics resembles the closed-shell molecule hydrazine (N₂H₄). In hydrazine, two inversion motions and one torsional motion between the two –NH₂ groups around the central N–N bond create eight non-superimposable equivalent equilibrium structures and a three-dimensional potential energy surface (PES).^28–30 In CH₂NH₂, the equilibrium structure has C_s symmetry, and the –CH₂ and –NH₂ groups are no longer equivalent in comparison with hydrazine. Therefore, the torsional double-well potential is lifted, and the number of non-superimposable equivalent equilibrium structures in CH₂NH₂ is reduced to four. The inversion of the –CH₂ and –NH₂ groups remains plausible. So far, no theoretical investigations has been done concerning the barrier and PES of this inversion motion in CH₂NH₂. In addition to the LAM itself, the spin angular momentum from the unpaired electron couples with the molecule's rotational motion and the LAM, which makes the spectroscopy even more complex.^31–34 Nevertheless, no rotational spectral information is available for CH₂NH₂ up to date, preventing a further understanding of its internal dynamics.

In this work, we report the first rotational spectrum of the CH₂NH₂ radical, in the submillimetre-wave range between 300 and 750 GHz. The broadband spectral recordings made possible by the Faraday rotation method not only provide the critical rotational constants of the radical, allowing its search in the interstellar medium, but also clearly reveal the tunnelling caused by the inversion motion of the –CH₂ and –NH₂ groups in the radical.

2 Experimental and theoretical methods

2.1 Experimental methods

In this work, CH₂NH₂ was synthesized in an H-abstraction reaction cell using F-atoms.^31,35 A mixture of fluorine heavily diluted in helium (5% concentration) is subject to a microwave discharge (2450 MHz, 50 W), upstream the reaction cell, to produce F atoms. The resulting mixture is introduced into the reaction cell in which gaseous methylamine is flowing. Methylamine (Sigma-Aldrich in Germany, 99%) was obtained as a pure sample (liquefied compressed gas) in 2020. The overall flow is ensured by a roots blower, with pressure measured above that pump. The optimum partial pressures to form CH₂NH₂ under these conditions are approximately 10 µbar of methylamine and 30 µbar of the F₂:He mixture, for a total pressure of 40 µbar. The H-abstraction reaction between methylamine and F-atoms is expected to produce both CH₂NH₂ and CH₃NH radicals. CH₂NH₂ is probably produced at larger quantity because it is more energetically stable than CH₃NH.¹⁸ However, we are unable to prove or deny the presence of CH₃NH in our spectrum because only target searches around spectral predictions for CH₂NH₂ were performed in this present work. Prediction of CH₃NH transition frequencies is unreliable due to the difficulty in modelling the coupling between the internal rotation of the CH₃ group and the electronic spin of the radical. Therefore, CH₃NH is outside the scope of the current study and has not been searched.

To record the pure rotation spectrum of CH₂NH₂, we followed a similar methodology that allowed some of us to study the dehydrogenated radicals of several large complex organic molecules.^32,33 This method associates a relatively specific synthesis method of the radical of interest with a sensitive and selective spectral acquisition scheme. This ensemble allows to record large portion of the spectrum, in which the numerous and intense transitions of diamagnetic species (precursor and other stable molecules produced by the reaction scheme) are filtered out, thus providing an immediate clear identification of transitions belonging to radical species.

The acquisition procedure rely on the application of the Faraday rotation effect to study the (sub-)millimetre-wave spectra of radical species.³² Briefly, an external alternating magnetic field—generated by an alternative current circulating in a coil surrounding the reaction cell—induces a small rotation of the polarization plane of the linearly vertically polarized input (sub-)millimetre-wave radiation at the resonant frequencies of paramagnetic species. This small rotation of the polarization is probed using a polarization grid, with its wires turned at a 45° angle and located between the output of the absorption cell and the detector. When using a lock-in detection whose reference signal is set to the frequency of this alternating field, one can selectively record the spectrum of open-shell species, mostly filtered out from diamagnetic transitions (for precursors with high magnetic susceptibility, some weak signals may be observed). Another advantage of the Faraday rotation acquisition method is the discrimination of the two-spin rotation components of the same rotational level, appearing in the spectrum with opposite phases. After lock-in demodulation, this translates into the appearance of spin-rotation doublets with opposite phases: one component with a positive signal and the other component showing a negative signal. The exact pointing direction of each spin component is subject to the polariser orientation, and effort was made to determine it a priori. In the present experiment, we found that the magnetic field modulation frequency optimizing the radical signal was equal to 9.9 kHz, leading to the generation of a magnetic field with an amplitude of about 11 Gauss.

For this spectrometer, the radiation source is a frequency multiplication chain (Virginia Diodes, Inc.) and the detector a close-circuit, cryogenic-free InSb hot electron bolometer (QMC Instruments Ltd). Spectra were recorded using 200 kHz frequency steps, and an acquisition speed of 200 ms to 400 ms per frequency point was sufficient to obtain clear spectra for strong, isolated transitions, whereas up to 10 s was adopted for a few selected weaker lines. Typical signal-to-noise ratio (SNR) is around 20 for strong lines.

2.2 Theoretical methods

For the theoretical investigations, we started from the optimized geometry of CH₂NH₂ reported by Puzzarini et al.¹⁸ We performed density functional theory (DFT) calculations with the Gaussian16 software,³⁶ and coupled cluster calculations with the Molpro2025 software.³⁷ PES scans were performed to investigate the LAMs in CH₂NH₂, harmonic frequency calculations were performed to find out the normal modes associated with the LAMs, and intrinsic reaction coordinate (IRC) scans were performed to estimate the tunnelling barrier and energy splitting.

The molecular axis system and the LAM coordinates, including atom numbers and symbols, are shown in Fig. 1. In our axis system, we fix the C–N bond on the x axis for simplicity, and we place the symmetry plane on the xz plane. The dipole moment components are aligned with the x and z axes. We introduced two dummy atoms, X₍₃₎ and X₍₄₎, to assist the proper construction of the Z-matrix. These two dummy atoms are the central point of the two H atoms in the –CH₂ and –NH₂ groups, respectively, and therefore they are also on the xz plane. In this system, we note γ₁ and γ₂ to be the inversion angles of –CH₂ and –NH₂, respectively. γ₁ is the supplementary angle of the dihedral angle of the –CH₂ plane and the xy plane, equal to 180° − ∠312, and γ₂ is the supplementary angle of the dihedral angle of the –NH₂ plane and the xy plane, equal to 180° − ∠421. Both γ₁ and γ₂ can vary between −180° and 180°, but more realistically between −90° and 90°. γ₁ = γ₂ = 0° means both –CH₂ and –NH₂ are in the xy plane. Additionally, we note τ to be the torsional angle between the –CH₂ plane and the –NH₂ plane, which corresponds to the dihedral angle ∠3-1-2-4. τ can also vary between 0° and 360°. τ = 180° is the equilibrium anti-configuration, and τ = 0° indicates the syn-configuration which is of higher energy. The values of the structure parameters are listed in Table 1.

Fig. 1 The molecular axis system of CH₂NH₂ and its inversion coordinates γ₁, γ₂ (top panel) and torsional coordinate τ (bottom panel). In the top panel, atoms (C₍₁₎, N₍₂₎, X₍₃₎, X₍₄₎, H₍₆₎, H₍₈₎) and bonds in solid lines are in front of or on the xz plane, i.e., the y component ≥0. Atoms (H₍₅₎, H₍₇₎) and bonds in dotted lines are behind the xz plane, i.e., the y component <0. The dipole moment components of the equilibrium structure along the a and c inertia axes of the molecules are also shown in the schema. They also lie in the xz plane.

Table 1 Equilibrium structure parameters of CH₂NH₂ in our molecular axis system shown in Fig. 1 based on the optimized geometry of CH₂NH₂.¹⁸γ₁ and γ₂ are at their equilibrium values. Unit of length in Å and unit of angle in °

Bond length		Bond angle
C–N	1.39168120	∠X(3) CN	152.35013902
C–H(5), C–H(6)	1.07730759	∠X(4) NC	140.37867543
N–H(7), N–H(8)	1.00700339	∠H(5) CH(6)	118.90541376
C–X(3)	0.54754093	∠H(7) NH(8)	111.58138854
N–X(4)	0.56615513	γ 1 (eq.)	27.64986098
		γ 2 (eq.)	39.62132457

---

Under the above-mentioned coordinate system, a 2-D relaxed PES scan concerning the inversion angles γ₁ and γ₂ were performed. We optimized first the structure of CH₂NH₂ at each scanned point by fixing γ₁ and γ₂ to their expected values, and let other structure parameters free to vary. The optimization was carried out under the B2BLYP-D3PJ/aug-cc-pVTZ level of theory. Following coupled cluster calculations were performed on each optimized geometry under the RCCSD(T)-F12/aug-cc-pVTZ-F12 level of theory. An additional 1-D rigid PES scan, as it is computationally less expensive, was performed under the B2BLYP-D3PJ/aug-cc-pVTZ level of theory on the torsional angle τ, in order to verify the conformation barrier between the equilibrium anti- geometry and the more energetic syn- geometry. Additionally, harmonic frequency calculations under the same DFT level of theory were performed to find the vibrational normal modes and their frequencies associated with the inversion of –CH₂ and –NH₂.

We can estimate the tunnelling splitting from the theoretical potential using the Wentzel–Kramers–Brillouin (WKB) approximation:†³⁸

In this equation, ω is the small amplitude harmonic oscillation angular frequency at the bottom of the well, a′ is the tunnelling point where the zero-point-energy (ZPE) crosses the potential, i.e., , μ is the reduced mass, and E₀ is the barrier height. The challenge here is to transfer the multi-atomic motion into this 1D potential for integration, and to find the correct x and reduced mass corresponding to this reduction. An appropriate way is to use the IRC q, because in this case, the reduced mass is implicitly included in the coordinate, as , and therefore . We first benchmarked this approximation with the classical umbrella motion of ammonia. For ammonia, the experimental barrier height E₀ is 2023 cm⁻¹ and ω = 950 cm⁻¹.²⁷ We performed the IRC scan the same way as what we did for CH₂NH₂, using the same level of theory, and we obtained theoretical E₀ = 1679 cm⁻¹ and ω = 1037 cm⁻¹. They slightly deviate from the experimental values, but agree well with other theoretical investigations.³⁹ If we adopt the theoretical E₀ and ω, and numerically integrate on the IRC potential curve to get S₀, we obtain ΔE = 3.9 cm⁻¹. If we rescale the IRC potential to match the experimental potential barrier and adopt the experimental ω, we obtain ΔE = 1.2 cm⁻¹, which is in accordance with the experimental value of 0.8 cm⁻¹.⁴⁰ This test proves that our approach to use the WKB approximation is valid.

3 Results

3.1 Theoretical investigation of the inversion tunnelling

The theoretical results, presented in Fig. 2, provide rich information on the internal dynamics of CH₂NH₂. Fig. 2(a) shows the relaxed 2-D PES with respect to γ₁ and γ₂, where the two equivalent equilibrium structures are located diagonally at the bottom left and top right corners. At the centre of the PES, where (γ₁,γ₂) = (0°,0°), a local maximum of 620 cm⁻¹ is found. The saddle points, where the transition states of the inversion are located, are however off the centre at (γ₁,γ₂) ≈ (−19.5°,5°) and (γ₁,γ₂) ≈ (19.5°,−5°). The dashed curves in Fig. 2(a) mark the tunnelling path of the inversion motion on the PES. The tunnelling effect leads to a splitting of the ground vibrational state into two substates: 0⁺ and 0⁻. The inversion is directed along the c principal axis and preserves the ac-plane of symmetry. CH₂NH₂ has dipole moment components of μ_a = 0.9 D and μ_c = 0.5 D.¹⁸a-Type transitions occur within each tunnelling substate (0⁺ ↔ 0⁺ or 0⁻ ↔ 0⁻), whereas c-type inversion-rotation-tunnelling transitions occur between two tunnelling substates (0⁺ ↔ 0⁻).

Fig. 2 PES of CH₂NH₂ with respect to (a) the two inversion angles γ₁ and γ₂; (b) the intrinsic reaction coordinate; (c) the torsional angle τ. The two white dashed lines in (a) connecting the two minima represent the tunnelling path in (b). The intervals between the solid contours in (a) are 50 cm⁻¹, and those between the dashed contours in (a) are 500 cm⁻¹ for clarity. The blue part of the IRC potential in (b) is the part used in the integration to obtain S₀ in q. (1). (a) Is performed by relaxed PES scan which first optimized the structure on each grid point at the B2BLYP-D3BJ/aug-cc-pVTZ level of theory, and then performed single point energy calculations at the RCCSD(T)-F12/aug-cc-pVTZ-F12 level of theory. (b) Is performed using the IRC scan in Gaussian16 at the B2BLYP-D3BJ/aug-cc-pVTZ level of theory. (c) Is performed as a rigid PES scan using the B2BLYP-D3BJ/aug-cc-pVTZ level of theory.

Fig. 2(b) demonstrates the reduced 1D inversion potential under the IRC system. The potential is asymmetric, in agreement of the tunnelling path marked in Fig. 2(a). We note that the IRC potential cannot completely show the “double well” shape. This is because the IRC calculation has to start at the transition state and will stop at the minimum. The potential outside this reaction coordinate range cannot be explored. Despite this limit, we can find the potential barrier to be 518.3 cm⁻¹ (B2BLYP-D3PJ value), much lower than that in hydrazine (2072 cm⁻¹)²⁹ or in ammonia (2023 cm⁻¹).⁴¹ Two harmonic vibrational modes associated with the inversion motion are found at 604.4 cm⁻¹ and 669.1 cm⁻¹, respectively. The ZPE, 318.4 cm⁻¹, is thus taken as the half of the average value of the two harmonic frequencies. The ground vibrational state of CH₂NH₂ is therefore near the top of the barrier, marked by the dashed horizontal line in Fig. 2(b). A notable tunnelling splitting is expected. Using the benchmarked WKB method, despite of the asymmetry of the IRC potential, we integrated over the part of the potential higher than the ZPE, which is the blue part in Fig. 2(b). We obtained S₀/ħ = 1.487, and thus ΔE = 49.2 cm⁻¹.

On the other hand, the rigid 1-D PES scan with respect to τ, plotted in Fig. 2(c), shows the torsional barrier is over 4300 cm⁻¹. From Fig. 2(a), it can also be seen that the syn- geometry is not a stable conformer, but rather an unstable geometry on the tunnelling path of the inversion motion. As a result, we conclude that the torsion will not lead to observable tunnelling splitting in the submillimetre-wave spectrum.

3.2 Experimental spectra and analysis

To start the search for spectral lines, our preferred targets were a-type transitions which fall into the submillimetre-wave range. The weaker c-type transitions occur between substates and are much more challenging to search for because their frequencies depend on the exact energy difference between the substates, and because the R-branch of these transitions are above 1.5 THz, far beyond the frequency range of our spectrometer. Because of the large frequency span between the intense a-type transitions of CH₂NH₂, we chose to perform targeted search in restricted spectral ranges, contrary to the case of our previous studies concerning CH₃COCH₂³² and CH₃OCH₂³³ radicals when large portions of the (sub-)millimetre-wave range were recorded to search for reproducible spectral patterns. We started the spectral search around 360 GHz based on the predicted frequency from the theoretical spectroscopic parameters,¹⁸ using an effective Watson-type A-reduced Hamiltonian and the SPFIT/SPCAT programs in the CALPGM suite⁴² to predict transition frequencies and then fit all measured lines. As mentioned before, thanks to the Faraday rotation acquisition technique, the up-down spectral line pair from the spin-rotation components offers us a unique signature to identify these components from the same rotational transition. We first identified the N = 7 ← 6 and N = 8 ← 7, K_a = 0 and K_a = 1 lines of both 0⁺ and 0⁻ substates. The spectrum of N = 8_0,8 ← 7_0,7 is shown in Fig. 3(a) as an example. This observation quickly confirmed the tunnelling motion suggested by our theoretical investigation, which re-assures us that we were detecting the sought radical species. Using these lines frequencies, a preliminary fit with all spectroscopic parameters fixed from the ab initio calculation,¹⁸ except for the B and C constants, were performed, leading to a converging fit and a refined prediction for the searches of more lines of CH₂NH₂. At this stage, we treated the 0⁺ and 0⁻ substates separately as if they were two different molecules, in order to reduce the complexity of the prediction-search iteration for more spectral lines.

Fig. 3 Sample spectra of CH₂NH₂ in the submillimetre-wave region. Experimental spectra are shown in black curves. Predictions from model are shown in red and blue vertical sticks. Blue sticks represent the 0⁺ state, and red sticks represent the 0⁻ state. Panel (a) shows the N_Ka,Kc = 8_0,8 ← 7_0,7 transition. The two lines pointing up are from the two J = N + 1/2 spin components, and the two lines pointing down are from the two J = N − 1/2 spin components. Panel (b) shows the N = 13 ← 12, K_a = 5 transitions. The four lines pointing up are from the J = N + 1/2 spin components, and the four lines pointing down are from the J = N − 1/2 spin components. Solid sticks represent the K_c = 9 ← 8 branch, and dashed sticks represent the K_c = 8 ← 7 branch. A calculated spectrum using the RAS model prediction and a Gaussian profile of full width half maximum of 1.8 MHz is shown in green below the experimental spectrum.

After successfully finding all lines below K_a = 3, we encountered difficulty due to line congestion and inaccurate predictions for higher K_a values. Especially for the upper 0⁻ state, we observed an increasing trend of perturbation from the model prediction for K_a = 3, which leads to difficulty in continuing the search for lines of higher K_a values. We analysed the reduced energy levels (E_r = E − 0.5(B + C)N″(N″ + 1)) with respect to N″ (N of the lower rotational state) of the lines that start to significantly deviate from the semi-rigid rotor model. Demonstrated in Fig. 4, the K_a = 5 levels of the 0⁺ state become close to the K_a = 2, 3 levels of the 0⁻ state, assuming a ΔE ≈ 50 cm⁻¹. This behaviour agrees with our observation of the increasing deviation from semi-rigid rotor model for these K_a levels, and it indicates the existence of a Coriolis-type perturbation between the two substates.

Fig. 4 Reduced energy (E_r = E − 0.5(B + C)N″(N″ + 1)) diagram of the two tunnelling substates of CH₂NH₂. The energies are from the J = N + 1/2 branch. Blue squares represent the 0⁺ state and red dots represent the 0⁻ state.

A combined fit of all so-far measured lines from the two substates were successful after adding the Coriolis coupling term F_ac and the energy difference ΔE between the two substates. We note that this fit is highly non-linear and its convergence requires the initial guess of the E value to be not too far away from the optimized value. Guided by this model, all lines up to K_a = 7, and a few strong K_a = 8 lines lying in our frequency range have been successfully measured and then added to the fit.

The final fit was performed using two different effective models. The first one, as described above, is to treat the two tunnelling substates as two separate vibrational states coupled by a Coriolis interaction. We refer to this model as the “separate model”. The second model is using the reduced axis system (RAS) formalism, first proposed by Pickett,⁴³ which divides the Hamiltonian into three parts, given in eqn (2)

Ĥ_RAS = Ĥ_A + Ĥ_Δ + Ĥ_C (2)

where Ĥ_A is the standard Watson-type A-reduced Hamiltonian (in I^r representation), Ĥ_C is the Coriolis coupling Hamiltonian, and Ĥ_Δ defined in eqn (3)

ΔE is the energy difference between the two tunnelling substates, and N, N_z, and N_± = N_x ± iN_y are the rotational angular momentum operators.^44,45 A unique set of rotational, centrifugal distortion, and spin-rotation coupling and its centrifugal distortion corrections (Δ^S(1)_N, Δ^S(1)_K, and Δ^S(1)_NK) are used in the RAS model, which can be viewed as the average values of the corresponding terms of the two substates from the separate model. They can then be directly compared with the theoretical values.

333 unique lines are assigned to 383 transitions of the two tunnelling substates, covering the lower state N quantum number from 6 to 14, and the lower state K_a quantum number from 0 to 8. All possible a-type transitions within this quantum number range have been observed. We kept, however, only lines sufficiently resolved in the final fit, and excluded heavily congested ones, such as several lines of K_a = 4, 5, and 6. We attributed a frequency uncertainty of 100 kHz to strong, isolated spectral lines (SNR ≥ 10), and 200 kHz (3 ≤ SNR < 10) or 400 kHz (SNR < 3) to spectral lines of weaker intensity or blended lines. The full table of the measured line list, including the transition quantum numbers, line frequencies, measurement uncertainties, and their frequency differences from the best-fit models, is included in digital form (SPFIT.lin and .fit files) in the SI and on the online data archive (see Data availability). We used a total of 27 parameters in the separate model, and 26 parameters in the RAS model, to reach reduced root-mean-square (rms) residuals of 1.01 for both models. The fit result is shown in Table 2 and is compared with the theoretical values by Puzzarini et al.¹⁸

Table 2 Fit result of CH₂NH₂ using the separate model and the RAS model^a

Parameter	Unit	Separate model			RAS model	Theoretical value^18
		v = 0^+	v = 0^−	Average^b
a Numbers in parentheses are 1σ standard errors in units of the last digit. b The averaged value of the v = 0^+ and v = 0^− substates. c These parameters are shared between the two tunnelling substates in the separate model. d Number of unique line frequencies. The sum of nl of the v = 0^+ and v = 0^− substates is larger than the nl in the RAS model by 1 because of a blended line containing both the transitions from the two substates. e Number of assigned transitions. f Number of parameters. g Values from fc-CCSD/cc-pVDZ calculation instead of from ref. 18 C. Puzzarini, private communication.
A	MHz	149 137.20(67)	147 450.08(90)	148 293.6(11)	148 294.13(52)	146 501.69
B	MHz	27 488.7033(76)	27 413.9557(85)	27 451.330(11)	27 451.3370(56)	27 393.55
C	MHz	23 601.5190(71)	23 635.5126(82)	23 618.516(11)	23 618.5115(51)	23 642.74
Δ N	kHz	49.020(10)	48.845(12)	48.933(15)	48.9341(77)	46.4208^g
Δ NK	MHz	0.32495(36)	0.31322(67)	0.31909(76)	0.31913(16)	0.265194^g
Δ K	MHz	3.68(27)	2.02(29)	2.85(39)	3.06(20)	2.13999^g
δ N	kHz	7.4766(76)	7.1867(94)	7.332(12)	7.3351(62)	6.22786^g
δ K	MHz	0.3713(31)	0.2884(26)	0.3313(40)	0.3329(20)	0.256068^g
Φ KN	kHz	—	0.0214(29)	0.0107(29)	0.0157(17)	0.00156
ε aa	MHz	−194.79(24)	−188.72(31)	−191.75(39)	−191.84(20)	−199.82
ε bb	MHz	−62.464(50)	−60.833(82)	−61.6485(96)	−61.796(43)	−58.87
ε cc	MHz	—	0.561(83)	—	0.124(33)	6.409
Δ ^S N	kHz	—	−0.86(16)	—	—	—
Δ ^S(1) N	kHz	—	—	—	−0.713(41)	—
Δ ^S(1) K	MHz	—	—	—	−0.0616(38)	—
Δ ^S(1) NK	kHz	—	—	—	−0.386(50)	—
F ac	MHz	287.64(70)		—	286.81(64)
F acN	kHz	5.29(24)		—	4.68(19)
F acK	MHz	−0.149(22)		—	−0.106(18)
ΔE^c	GHz	1541.58(52)		—	1543.52(17)
	cm^−1	51.421(17)		—	51.4862(57)
E N	MHz	—	—	—	10.1883(12)	—
E K	MHz	—	—	—	830.93(31)	—
E NN	kHz	—	—	—	−0.0907(74)	—
E NK	kHz	—	—	—	−6.84(44)	—
E 2	MHz	—	—	—	13.5851(20)	—
E 2N	kHz	—	—	—	−0.1528(55)	—
E 2K	MHz	—	—	—	−0.03344(88)	—
n l		178	156	—	333	—
n t		207	176	—	383	—
n p		27		—	26
		6–14, 0–8		—	6–14, 0–8
rms	kHz	166		—	169
Weighted rms		1.01		—	1.01

---

4 Discussion

The fit results of both models have reached satisfactory level to the experimental accuracy. The average values of the rotational constants A, B, and C, and the quartic centrifugal distortion constants of the separate model are in excellent agreement with the values from the RAS model. They are both in good agreement with the theoretical values from Puzzarini et al. as well. Nevertheless, we note that the signs of the quartic centrifugal distortion constants reported by Puzzarini et al. could be unreliable because of the choice of the coordinate system in the calculation. More reliable values from a calculation at the fc-CCSD/cc-pVDZ level‡ are therefore listed in Table 2 in place of the values from Puzzarini et al. A sextic centrifugal distortion constant Φ_KN is needed for the 0⁻ substate in the separate model to improve the fit. Other centrifugal distortion constants of sextic or higher order are found statistically insignificant. It is the same case in the RAS model where Φ_KN is needed. The diagonal spin-rotation coupling terms ε_aa and ε_bb are also in good agreement with the theoretical values, but ε_cc is one order of magnitude smaller than the theoretical value. In addition to the diagonal terms, the off-diagonal term (ε_ab + ε_ba)/2 is fixed to 0 instead of 18.73 MHz as suggested by the theoretical calculation. This off-diagonal term should be null for symmetry reason, and it is possible that the theoretical value is also not reliable for the same reason as the centrifugal distortion constants. The Coriolis coupling term and the separation energy ΔE between the two substates are also in good agreement between the separate model and the RAS model, indicating that the correlation of the these parameters with other ones are not significant.

The inversion-rotation-tunnelling spectrum would appear around 1.5 THz, which is beyond the frequency range of our spectrometer. Therefore, the energy difference ΔE between the two substates cannot be measured directly. Instead, it is inferred indirectly from the Coriolis coupling effect. The experimental fits from both models found ΔE around 51 cm⁻¹, in excellent agreement with the WKB estimation of 49.2 cm⁻¹. This agreement confirms that our indirect approach of finding the ΔE is valid and reliable. It would be valuable for further experiments to measure the c-type inversion-rotation-tunnelling transitions between the two tunnelling substates in the far-infrared region, using adapted light sources, such as the synchrotron radiation. A direct measurement of the tunnelling splitting energy would be appreciable to provide more accurate benchmark for the theory in treating the internal dynamics in this radical.

The CH₂ and NH₂ groups both have two identical H nuclei (H5, H6 for CH₂ and H7, H8 for NH₂). Four spin states, (ortho-, ortho-), (ortho-, para-), (para-, ortho-), and (para-, para-), with respective weights of 9, 3, 3, and 1, The total wavefunction Ψ_total is symmetric with respect to the feasible permutation P = (56)(78), equivalent to a π-rotation along the a axis. Ψ_total is the product of the electronic, vibrational, rotational and spin wavefunctions:

Ψ_total = Ψ_elec × Ψ_vib × Ψ_rot × Ψ_spin (4)

Ψ _elec is X̃²A″, anti-symmetric with respect to the π-rotation along the a axis. It can be seen from the shape of the highest occupied molecular orbital (HOMO) of CH₂NH₂ shown in Fig. 5. Ψ_vib is symmetric for the v = 0⁺ state and anti-symmetric for the v = 0⁻ state. As a result, for v = 0⁺, odd rotational states should be associated with symmetric spin states (ortho-, ortho-) and (para-, para-), whereas even rotational states (even K_a) should be associated with the anti-symmetric spin states (ortho-, para-) and (para-, ortho-). The weight ratio between the even and odd rotational states (odd K_a) is therefore (3 + 3) : (9 + 1) = 3 : 5. This ratio is reversed for v = 0⁻. The spin weights of even and odd rotational states were then taken into account when generating the intensities of individual transitions and calculating the partition function. An excellent agreement between the calculated and measured line intensities shown in Fig. 3 supports the above analysis of spin statistics. The spin statistics also explain the selection rule of the transitions. Since only the levels with the same spin state symmetry can be connected, a-type transitions connect levels within a substate where the parity of K_a and Ψ_vib do not vary, and c-type transitions connect levels between the two substates where the parity of K_a and Ψ_vib both flip.

Fig. 5 The highest occupied molecular orbital (HOMO) of CH₂NH₂. The red is positive part and blue is negative part.

With the fitted parameters, we were able to predict the frequencies of all a-type transitions of CH₂NH₂ in the submillimetre-wave region. The prediction is predictive within our measured frequency range, 300–750 GHz. Fig. 3(b) shows an example of the congested lines where N = 13 ← 12 and K_a = 5. In such situation of line confusion, we only kept in the fit two sufficiently isolated lines: the ones of the lowest and highest frequency. Nevertheless, by comparing the experimental spectrum with a simulated spectrum using the RAS model prediction and a Gaussian profile of full width half maximum of 1.8 MHz, we can conclude that a satisfactory match between the model and experiment is reached. The full predictions using both models in the SPCAT format can be found in the SI. The effective partition functions Q_eff, obtained by direct summation of energy levels up to N = 99, are listed in Table 3. They are slightly different between the two models due to the use of different formalism. The listed partition functions have taken into account the spin statistics of H and are (5 + 3)/2 = 4 times of the true rotational partition function Q_rot. This factor of 4 is necessary to allow SPCAT to generate correct line intensities. From a physics point of view, the nuclear spin degeneracy factor is 3 for N and 2⁴ for the four hydrogens. To generate correct line intensities in SPCAT for hyperfine predictions while maintaining the 5 : 3 and 3 : 5 statistical weights, the listed partition functions must be scaled by the total degeneracy of the spins, 3 for N-hyperfine only and 48 for full hyperfine. Note that this factor of 48 accounts for the fact that SPCAT sums over all 48 nuclear spin states while applying the weighted intensity factors. Vibrational partition function Q_vib is also listed in Table 3. Q_vib is obtained from the harmonic vibrational frequencies calculated under the B2PLYPD3/dAug-cc-pVTZ level of theory.

Table 3 Rotational (Q_eff) and vibrational (Q_vib) partition functions of CH₂NH₂. Q_eff was used by SPCAT for linelist prediction taking into account the spin statistics. Q_eff is 4 times of the true rotational partition function Q_rot

Temperature (K)	Q eff		Q vib
	Separate model	RAS model
300.000	40 375.3593	40 370.8108	1.264802
225.000	25 311.9993	25 307.9276	1.099366
150.000	12 905.3074	12 902.2905	1.017907
75.000	3894.4876	3893.1459	1.000177
37.500	1146.8040	1146.4459	1.000000
18.750	366.5367	365.9456	1.000000
9.375	132.2735	127.1029	1.000000

---

The hyperfine components from the ¹⁴N and H atoms could not be resolved due to the Doppler limited linewidth of our experimental spectra. Instead, these components added an extra broadening to the measured lines. It is, however, worthwhile to evaluate their contribution to the line frequencies. If we adopt all the theoretical hyperfine coupling constants from Puzzarini et al. into our best-fit RAS model, observable frequency shits may occur to the predicted lines. The strong hyperfine components can span larger than the measured linewidth, along with weaker components spanning even further away, especially for lines of lower N levels. The inclusion of only the ¹⁴N quadrupole coupling constants and the inclusion of all hyperfine coupling constants also lead to different results, as we observed the frequencies of the strongest components from the all hyperfine model deviate from the weighted average frequency of the three strongest ¹⁴N quadrupole hyperfine components. Here we demonstrate the worst case in Fig. 6(a), the 4_0,4 ← 3_0,3 lines, lying around 203 GHz. These lines have also been searched for and recorded by an enormous accumulation time due to their low intensities at room-temperature. Their frequencies are also provided in the SI and in the online data repository. They were not included in our spectral fit, because the predicted line frequencies from our hyperfine-free model deviate from the experimental line centre up to 1 MHz, clearly visible in Fig. 6(a). Since resolving individual hyperfine components is not possible, we use a calculated line profile summing up all predicted hyperfine components (green curve in Fig. 6) to evaluate its effect to the frequency shift of lines. In Fig. 6(a), we show that the inclusion of the ¹⁴N hyperfine splitting can result in a visible frequency shift from the hyperfine-free line position. When the hyperfine splitting of all atoms is included, the frequency centre of the summed line profile slightly shifts again from the ¹⁴N-only case. In addition, hundreds of weaker components appear and spread even outside the frequency range of the plotted spectrum. This complicated effect may explain the origin of the large “obs.-calc.” frequency discrepancy that discouraged us from including the 4_0,4 ← 3_0,3 transitions into our hyperfine-free models. This case indicates that it is not always “safe” to assume that the three hyperfine components from ¹⁴N are sufficient to reproduce Doppler-limited line profile at low J values. Hyperfine effects of all atoms must be taken into account to reproduce a line profile that matches best to the experimental spectrum, although none of the individual component can be resolved. Nevertheless, in Fig. 6(b), we demonstrate a much more coherent case at higher J values: the 8_0,8 ← 7_0,7 lines (also presented in Fig. 3). In this case, the three models predict almost identical line centres, which assures that our hyperfine-free models are nevertheless reliable within the observed frequency and J, K_a range.

Fig. 6 The effect of hyperfine splitting (hfs) on the line positions. Panel (a) represents the case of the 4_0,4 ← 3_0,3 transition, which was not included in the fit. Panel (b) represents the case of the 8_0,8 ← 7_0,7 transition, identical to Fig. 3(a). Green curves are calculated line profiles constructed from a summation of Gaussian profiles of full width half maximum of 1.8 MHz for each individual line. Pairs of vertical dashed lines of 2 MHz spacing around each experimental line centre are drawn to demonstrate the level of frequency shifts. For the 4_0,4 ← 3_0,3 transition, the line position is sensitive to the hfs. Inclusion of hfs of all atoms is necessary to match the calculated line profile to experimental line positions. For the 8_0,8 ← 7_0,7 transition, the hfs does not cause remarkable shift to the line positions.

In addition to the hyperfine effects, we suppose that another subtle coupling effect between the tunnelling motion and the electronic spin can contribute to the frequency shift, such as a similar effect recently demonstrated between the internal rotation motion and electronic spin.³⁴

The lack of experimental data below 300 GHz, especially high resolution data that can decompose the hyperfine components, prevents us from the quantitative evaluation of the above mentioned effects. All their contributions remain hypothetical. It means that the generation of an accurate frequency prediction for astronomical search of CH₂NH₂ in the microwave and in the 3 mm band will not be perfect using the hyperfine-free models that we have provided. It will be invaluable to fill this data gap, by future experiments between 50 and 300 GHz, preferably in a jet-cooled environment and with techniques that can reach spectral resolution at the kHz level.

Conclusions

In conclusion, we have measured and modelled the a-type pure rotational transitions of the aminomethyl radical CH₂NH₂ for the first time between 300 and 750 GHz, from N″ = 6 to N″ = 14, and K_a = 0 to K_a = 8. The spectrum reveals the existence of a collective inversion tunnelling motion of the –CH₂ and –NH₂ groups in the radical. The measured transitions were fitted with two effective models, a “separate model” that treats the two tunnelling substates separately, and a “RAS model” that treats the two substates using a single set of parameters under a reduced axis system. Both models adopt a Coriolis-type interaction between the two tunnelling substates, and derive the tunnelling splitting energy ΔE from this interaction. Both models reproduce the line frequencies of all the measured transitions, but extrapolation at lower frequencies for low-K_a lines have proven challenging and should be treated with caution. The rotational and centrifugal distortion constants from both models are in good agreement with the theoretical values. The ΔE is derived to be 51.421(17) cm⁻¹ and 51.4862(57) cm⁻¹, respectively, from the separate model and the RAS model. These values agree well with the theoretical estimation by the WKB method, which is 49.2 cm⁻¹. Our results provide the first analysis and frequency prediction of CH₂NH₂ in the submillimetre-wave region. Further spectroscopic measurements in the far-infrared region with brighter light sources are expected to improve the model and provide a more accurate measurement of the tunnelling splitting energy. We have also demonstrated that the hyperfine splitting effects could lead to observable frequency shifts from the our hyperfine-free models at millimetre-wave and lower frequencies, urging further spectroscopic measurements below 300 GHz at lower temperatures and higher spectral resolution. With the new model³⁴ that correctly treats the coupling between torsion and spin rotation, the study of CH₃NH may now also be possible.

Author contributions

L. Z.: conceptualization, investigation, formal analysis, visualization, writing – original draft, writing – review & editing, funding acquisition. R. C.: investigation, writing – review & editing. L. M.: investigation, project administration, writing – review & editing, funding acquisition. M.-A. M.-D.: investigation, analysis, writing – review & editing. O. P.: investigation, writing – review & editing, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data reported in this article are available at “Rotational spectrum of the aminomethyl radical CH₂NH₂” at https://doi.org/10.5281/zenodo.17765006. Specifically, the PES data to generate Fig. 2, the input and output files to reproduce the fit results in Table 2, and frequency predictions in the SPCAT format using the separate model and RAS model, are included. These files are also available in the supplementary information (SI) and in the online data repository. The predicted spectra are also available in the SPCAT format from the Lille Spectroscopic Database⁴⁶ (https://lsd.univ-lille.fr), where customized predictions can be generated using options such as intensity units, temperature, and frequency range. It provides additional flexibility in the data access. The SPFIT/SPCAT package⁴² to fit the lines and generate the frequency predictions can be found at http://info.ifpan.edu.pl/~kisiel/asym/asym.htm#pickett, or from the molecular spectroscopy site of NASA's Jet Propulsion Laboratory https://spec.jpl.nasa.gov/. The supplementary information includes the PES data to generate Fig. 2, the input and output files to reproduce the fit results in Table 2, and frequency predictions in the SPCAT format using the separate model and RAS model. See DOI: https://doi.org/10.1039/d5cp02849c.

Acknowledgements

The authors would like to thank Prof. C. Puzzarini for her helpful discussion on the signs of the quartic centrifugal distortion constants, and for providing us with a set of more reliable theoretical values. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 894508. L. Z. acknowledges also the start-up fund from the Agence Nationale de la Recherche (ANR-22-CPJ2-0030-01). The experimental work was performed thanks to financial support from the LabEx PALM (ANR-10-LABX-0039-PALM), from the Région Ile-de-France through DIM-ACAV+, and from the Agence Nationale de la Recherche (ANR-19-CE30-0017-01). The calculation work was carried out using the CALCULCO computing platform, supported by DSI/ULCO (Direction des Systèmes d'Information de l'Université du Littoral Côte d'Opale). We acknowledge support from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES, and the financial support from GDR EMIE 3533.

References

1. N. Kaifu, M. Morimoto, K. Nagane, K. Akabane, T. Iguchi and K. Takagi, Astrophys. J., 1974, 191, L135–L137.
2. N. Fourikis, K. Takagi and M. Morimoto, Astrophys. J., 1974, 191, L139.
3. S. Zeng, I. Jiménez-Serra, V. M. Rivilla, S. Martín, J. Martín-Pintado, M. A. Requena-Torres, J. Armijos-Abendaño, D. Riquelme and R. Aladro, Mon. Not. R. Astron. Soc., 2018, 478, 2962–2975.
4. E. G. Bøgelund, B. A. McGuire, M. R. Hogerheijde, E. F. van Dishoeck and N. F. W. Ligterink, Astron. Astrophys., 2019, 624, A82.
5. M. Ohishi, T. Suzuki, T. Hirota, M. Saito and N. Kaifu, Publ. Astron. Soc. Jpn., 2019, 71, 86.
6. T. Suzuki, L. Majumdar, P. F. Goldsmith, K. Tokuda, H. Minamoto, M. Ohishi, M. Saito, T. Hirota, H. Nomura and Y. Oya, Astrophys. J., 2023, 954, 189.
7. M. Cable, S. Hörst, C. He, A. Stockton, M. Mora, M. Tolbert, M. Smith and P. Willis, Earth Planet. Sci. Lett., 2014, 403, 99–107.
8. D. P. Glavin, J. P. Dworkin and S. A. Sandford, Meteorit. Planet. Sci., 2008, 43, 399–413.
9. P. D. Holtom, C. J. Bennett, Y. Osamura, N. J. Mason and R. I. Kaiser, Astrophys. J., 2005, 626, 940–952.
10. L. Marrodán, T. Pérez and M. U. Alzueta, Combust. Flame, 2024, 259, 113130.
11. M. Förstel, A. Bergantini, P. Maksyutenko, S. Góbi and R. I. Kaiser, Astrophys. J., 2017, 845, 83.
12. D. E. Woon, Astrophys. J., 2002, 571, L177.
13. S. Pilling, L. Baptista, H. M. Boechat-Roberty and D. P. P. Andrade, Astrobiology, 2011, 11, 883–893.
14. R. T. Garrod, Astrophys. J., 2013, 765, 60.
15. T. Suzuki, L. Majumdar, M. Ohishi, M. Saito, T. Hirota and V. Wakelam, Astrophys. J., 2018, 863, 51.
16. A. Sato, Y. Kitazawa, T. Ochi, M. Shoji, Y. Komatsu, M. Kayanuma, Y. Aikawa, M. Umemura and Y. Shigeta, Mol. Astrophys., 2018, 10, 11–19.
17. T. Suzuki, M. Ohishi, T. Hirota, M. Saito, L. Majumdar and V. Wakelam, Astrophys. J., 2016, 825, 79.
18. C. Puzzarini, Z. Salta, N. Tasinato, J. Lupi, C. Cavallotti and V. Barone, Mon. Not. R. Astron. Soc., 2020, 496, 4298–4310.
19. M. V. Kantak, K. S. De Manrique, R. H. Aglave and R. P. Hesketh, Combust. Flame, 1997, 108, 235–265.
20. P. Glarborg, C. S. Andreasen, H. Hashemi, R. Qian and P. Marshall, Int. J. Chem. Kinet., 2020, 52, 893–906.
21. J. Dyke, E. Lee and M. Z. Niavaran, Int. J. Mass Spectrom. Ion Processes, 1989, 94, 221–235.
22. T. Jansen, I. Trabjerg, S. Rettrup, P. Pagsberg and A. Sillesen, Acta Chem. Scand., 1999, 53, 1054–1058.
23. M. P. Rissanen, A. J. Eskola, T. L. Nguyen, J. R. Barker, J. Liu, J. Liu, E. Halme and R. S. Timonen, J. Phys. Chem. A, 2014, 118, 2176–2186.
24. T. G. Wright and T. A. Miller, J. Phys. Chem., 1996, 100, 4408–4412.
25. G. Muller, K. Voronova, B. Sztáray and G. Meloni, J. Phys. Chem. A, 2016, 120, 3906–3916.
26. P. R. Joshi and Y.-P. Lee, Commun. Chem., 2022, 5, 1–7.
27. H. V. L. Nguyen, I. Gulaczyk, M. Kreglewski and I. Kleiner, Coord. Chem. Rev., 2021, 436, 213797.
28. S. Tsunekawa, T. Kojima and J. Hougen, J. Mol. Spectrosc., 1982, 95, 133–152.
29. W. Łodyga, M. Krȩglewski and J. Makarewicz, J. Mol. Spectrosc., 1997, 183, 374–387.
30. M. Krȩglewski, J. Cosléou and G. Wlodarczak, J. Mol. Spectrosc., 2002, 216, 501–504.
31. O. Chitarra, M.-A. Martin-Drumel, B. Gans, J.-C. Loison, S. Spezzano, V. Lattanzi, H. S. P. Müller and O. Pirali, Astron. Astrophys., 2020, 644, A123.
32. R. Chahbazian, L. Juppet and O. Pirali, J. Phys. Chem. Lett., 2024, 15, 9803–9810.
33. R. Chahbazian, M.-A. Martin-Drumel and O. Pirali, Astrophys. J., 2025, 982, 162.
34. L. H. Coudert, O. Pirali, M.-A. Martin-Drumel, R. Chahbazian, L. Zou, R. A. Motiyenko and L. Margulès, Phys. Rev. Lett., 2025, 134, 173001.
35. O. Chitarra, O. Pirali, J.-T. Spaniol, T. S. Hearne, J.-C. Loison, J. F. Stanton and M.-A. Martin-Drumel, J. Phys. Chem. A, 2022, 126, 7502–7513.
36. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
37. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, I. Miller, F. Thomas, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
38. A. Garg, Am. J. Phys., 2000, 68, 430–437.
39. C. Léonard, S. Carter and N. C. Handy, Chem. Phys. Lett., 2003, 370, 360–365.
40. C. E. Cleeton and N. H. Williams, Phys. Rev., 1934, 45, 234–237.
41. D. Papoušek, J. Stone and V. Špirko, J. Mol. Spectrosc., 1973, 48, 17–37.
42. H. M. Pickett, J. Mol. Spectrosc., 1991, 148, 371–377.
43. H. M. Pickett, J. Chem. Phys., 1972, 56, 1715–1723.
44. D. Christen and H. S. P. Müller, Phys. Chem. Chem. Phys., 2003, 5, 3600–3605.
45. L. Margulès, B. A. McGuire, M. L. Senent, R. A. Motiyenko, A. Remijan and J. C. Guillemin, Astron. Astrophys., 2017, 601, A50.
46. R. A. Motiyenko and L. Margulès, Astron. Astrophys., 2025, 699, A348.

---

Footnotes

† We slightly adapted the labels from Garg³⁸ for clarity.

‡ C. Puzzarini, private communication.

---