Tagging effects on the mid-infrared spectrum of microsolvated protonated methane

Abstract

Although bare protonated methane is by now essentially understood at the level of intramolecular large-amplitude motion, scrambling dynamics and broadband vibrational spectra, the microsolvated species still offer plenty of challenges. One aspect is the effect of the attached solvent molecules on the infrared absorption spectra of microsolvated CH₅⁺ complexes compared to the bare parent molecule. In this study we analyze, based on ab initio molecular dynamics simulations, protonated methane molecules that have been microsolvated with up to three hydrogen molecules, i.e. CH₅⁺·(H₂)_n. In particular, upon introducing a novel multi-channel maximum entropy methodology described herein, we are able to decompose the infrared spectra of these weakly-bound complexes in the frequency window from 1000 to 4500 cm⁻¹ into additive single mode contributions. Detailed comparisons to the bare CH₅⁺ parent reveal that these perturbed modes encode distinct features that depend on the exact microsolvation pattern. Beyond the specific case, such understanding is relevant to assess tagging artifacts in vibrational spectra of parent molecules based on messenger predissociation action spectroscopy.

---

I. Introduction

Since its discovery, published in 1952 by Tal'roze and Lyubimova,¹ protonated methane (CH₅⁺) has challenged the traditional understanding of what a molecular structure is. Consequently, this molecule has been intensely studied in order to understand its intriguing behavior and, over time, became the prototype of ”fluxional” (or ”floppy”) molecules.^2–12 It is the smallest possible carbocation and forms a three-center two-electron bond^13–15 where its constituents continuously exchange their sites, resulting in all five hydrogen atoms being dynamically equivalent even in the quantum ground state.^14,15 This peculiar phenomenon is dubbed ”hydrogen scrambling” (or just ”scrambling”) and is caused by intramolecular large-amplitude motion on an extremely shallow potential energy surface (PES); see ref. 11 for a comprehensive recent review.

Bare protonated methane has been thoroughly studied by theoretical means which goes back to the early days of computational quantum chemistry. This includes pioneering semiempirical work followed by accurate studies of the structures of the ground state and various isomers, their energy differences and the barriers involved in scrambling as well as harmonic vibrational frequencies.^13,16–23

The converged static picture²³ is that the global minimum of the PES, the e – C_s symmetric structure, is separated from two energetically higher-lying s – C_s and C_2v transition states by only about 0.1 and 1 kcal mol⁻¹, respectively. This scenario led early on to the conjecture that the molecule should indeed easily undergo scrambling at room temperature,¹³ which has been explicitly confirmed two decades later not only at room temperature,²⁴ but even at low temperatures close to the nuclear quantum ground state¹⁴ based on ab initio path integral simulations in conjunction with topological bonding analyses;¹⁵ these key findings have been confirmed by several other groups using different methods as comprehensively reviewed.¹¹

In addition to a wealth of static quantum chemistry studies, ab initio molecular dynamics (AIMD) simulations²⁵ were used to compute, disentangle and assign the fully anharmonic infrared (IR) spectra of CH₅⁺ and all its H/D isotopologues;^9,11,26–28 the overriding impact of isotopic substitution on the properties of the CH₅⁺ isotopologues has been predicted early.²⁹ The decomposition and assignment of these broadband IR spectra in terms of modes was obtained from time-autocorrelation and cross-correlation function approaches based on the underlying ab initio trajectories;^27,30,31 see ref. 11 for a comprehensive review of these (ab initio) molecular dynamics based approaches to theoretical spectroscopy. It has been shown that AIMD simulations at a temperature of 300 K, thus relying on classical Newtonian dynamics of the nuclei, emulate the quantum effects of the nuclei quite well and also reproduce the experimental spectrum of CH₅⁺, while at lower temperatures of around 100 K the peaks in the IR spectrum become much sharper as a result of frozen scrambling dynamics.^9,26

In addition to these early ab initio path integral and AIMD simulations,^{5,14,15,24,26,29} diffusion Monte Carlo and vibrational configuration interaction calculations based on parameterized PESs^32–39 not only confirmed the previous findings^{5,14,15,24,26,29} that have been obtained from “on the fly” quantum simulations, but they also contributed significantly to our current knowledge of bare protonated methane. Most notable in this respect are the pioneering IR spectra that have been reported for several H/D isotopologues of protonated methane.^35,36 All of the available information taken together leads to the conclusion that bare protonated methane including its isotopologues is now essentially understood in terms of its peculiar intramolecular large-amplitude motion, scrambling dynamics, and broadband IR spectra.¹¹ Yet, largely unexplained remain the available high-resolution vibrational spectra^{5,12,37,40,41} of bare CH₅⁺, which continue to challenge state-of-the-art theoretical rovibrational spectroscopy.^42–46

This favorable situation for bare CH₅⁺ is in stark contrast to our current understanding of microsolvated CH₅⁺ species, which still offer many puzzles concerning the interplay of scrambling dynamics and microsolvation effects. Tagging of protonated methane using hydrogen molecules is possible due to the rather weak but still significant interactions between these so-called tags and the CH₅⁺ core. This approach has been pioneered experimentally by the group of Y. T. Lee^47–50 upon introducing an ingenious technique that is now known as messenger predissociation action spectroscopy; it is noted in passing that the accompanying simulations in ref. 48 must be considered with care due to technical shortcomings. Using this flavor of action spectroscopy, Lee and coworkers were able to measure vibrational action spectra of CH₅⁺·(H₂)_n complexes in the range of n = 1 up to 6, thus excluding the limiting bare case n = 0. The experiment accessed two rather narrow frequency windows from about 2700 to 3200 cm⁻¹ and 4060–4140 cm⁻¹, thus observing the C–H bands of microsolvated CH₅⁺ complexes and the H₂ stretching modes, respectively, in the mid-IR regime. These action spectra^47–50 remain so far the only existing experimental IR spectra of microsolvated CH₅⁺. We note in passing that the influence of tags on quasi-rigid molecules is generally small,^51,52 whereas studying the potentially significant changes of the IR spectra of highly fluxional molecules such as protonated methane is the challenge.

At variance with bare CH₅⁺, theoretical investigations of CH₅⁺·(H₂)_n species are rather scarce.^48,53–58 Moreover, structure optimizations of CH₅⁺·(H₂)_n lead to a wealth of minimum energy structures of differently tagged species. From the set of optimized structures⁵⁶ it is clearly seen that the tagging molecules prefer specific sites depending on n, see Fig. 1 where the global minima for n = 1, 2 and 3H₂ tags are displayed in panels b, e and h, respectively, compared to the bare CH₅⁺ reference molecule in panel a. Our atom/site labeling scheme (defining the permutations, see below) is the same that has been introduced for bare protonated methane,²⁷ where site 1 is the “outer moiety atom”, site 2 the “inner moiety atom”, site 3 denotes the “in-plane atom” while the sites 4 and 5 are located ”out-of-plane”. Thus, atoms 1 and 2 form the so-called H₂ “moiety” while atoms 3, 4 and 5 constitute the CH₃ “tripod”, the carbon nucleus of which forms a three-center two-electron bond^13–15 together with the two protons of the moiety. The singly tagged CH₅⁺·H₂ species has been found to have three close-lying low-energy minimum structures, while the doubly tagged case has three such structures and the n = 3 case features only a single low-energy minimum structure.⁵⁶

Fig. 1 Optimized low-lying minimum energy structures of H₂-microsolvated protonated methane based on ref. 56. The structures are sorted according to increasing energy.⁵⁶ As a reference, the bare protonated methane is shown in its e – C_s ground-state structure in panel a where the site labeling scheme (defining our standard ordering²⁷ for permutations) is also introduced. The singly tagged species has three distinct minimum structures that differ in the location of the tagging molecules as indicated: species b is the energetically lowest-lying minimum, whereas species c is ΔE = 0.35 kcal mol⁻¹ and species d is ΔE = 1.77 kcal mol⁻¹ above the global minimum energy. The doubly tagged species also has three distinct minimum structures where the two tag molecules are located at the positions 1 and 2 (species e being the global minimum), at 2 and 3 (species f with ΔE = 1.19 kcal mol⁻¹), and at 1 and 3 (species g with ΔE = 1.67 kcal mol⁻¹). Finally, the protons 3, 2 and 1 are the preferred sites when three tagging H₂ molecules are attached to CH₅⁺ as depicted in panel h.

Analyzing the optimized structures of tagged CH₅⁺ species, such as those depicted in Fig. 1, exclusively provides a static perspective in terms of equilibrium structures. But what about the impact of tagging on large-amplitude motion and ultimately on the scrambling dynamics of the CH₅⁺ core? This question is highly relevant to messenger spectroscopy, where ideally the tags should not (significantly) perturb the microsolvated parent species. In stark contrast, in microsolvating protonated methane, even the weak interactions of the tags with the bare CH₅⁺ species apparently do significantly impact on the vibrational spectra.^47–49 This finding, in turn, has been exploited to investigate experimentally the extend to which the vivid scrambling dynamics of CH₅⁺ itself might be affected by perturbations due to the tagging molecules.

Starting with the case of bare CH₅⁺, it is known that so-called “full scrambling” requires two processes as sketched in Fig. 2 (obviously, the following description rests on a quasi-classical picture that involves distinguishable protons whereas quantum effects are discussed in ref. 58). The first one shown in panel a consists of breaking and making of the H₂ moiety as a result of (relative) motion of the inner moiety proton toward the in-plane proton via the C_2v transition state and reformation of an accordingly displaced moiety, i.e. the two out-of-plane protons are only spectators. The resulting pseudorotational motion involves only three out of five protons, exchanges only one of the protons involved in the e – C_s equilibrium structure, and thus leads to so-called “partial scrambling”.⁵⁸ The second large-amplitude motion, see Fig. 2b is an internal rotation, where the two moiety atoms are rotating, via the s – C_s transition state, w.r.t. the pseudo C₃-axis of the tripod. This is an easy process in view of its very small activation energy, on the order of 0.1 kcal mol⁻¹ (cf.ref. 23) on the PES, that leads to an exchange of the inner and outer moiety protons within the e – C_s structure. When combined with partial scrambling, internal rotation also involves the out-of-plane protons and eventually leads to full scrambling of all five protons. The low-resolution IR spectra resulting from this complex scrambling dynamics have been computed for CH₅⁺ and all its isotopologues and used to assign the experimental broadband spectra.^26,27

Fig. 2 Schematic illustration of the pseudorotations and internal rotations in panels (a and b), respectively, that enable hydrogen scrambling in bare CH₅⁺. Note that the atom numbering is merely used here to visualize the scrambling which results in permutations of the atoms, whereas all such permutations are properly taken into account when performing the spectral analysis (see text).

Having understood the IR spectrum of bare CH₅⁺ and its relation to hydrogen scrambling, the question that remains to be answered is about the tagging effect on IR spectra when using messenger vibrational spectroscopy techniques. At the level of only structure, it has been recently demonstrated by low-temperature path integral simulations⁵⁸ that microsolvation by H₂ influences scrambling in distinct ways depending on the number of tagging molecules. The effects of one and three tagging hydrogen molecules were shown to be very similar, while two tagging H₂ lead to a distinctly different behavior. Internal rotation (cf.Fig. 2b) is found to be frozen in the case of one and three tagging molecules, while in the n = 2 case internal rotation is still active but partial scrambling is reduced drastically. In contrast, partial scrambling via pseudorotation (cf.Fig. 2a) is still possible for the n = 1 and 3 species of CH₅⁺·(H₂)_n. Such subtle yet significant perturbation of the scrambling dynamics of the CH₅⁺ core due to tagging raises the question as to the impact of the H₂ tags on the IR spectrum of protonated methane depending on their number n.

In this study, we computed the IR spectra of microsolvated protonated methane with up to n = 3 tagging hydrogen molecules in the mid-IR frequency window from 1000 to 4500 cm⁻¹ based on our AIMD simulation approach to theoretical spectroscopy.¹¹ Using classical nuclear dynamics, the scrambling motion is frozen at the chosen simulation temperature, which allows us to precisely address tagging effects on the bending and stretching bands of protonated methane at mid-IR frequencies, whereas the prominent signals due to hydrogen scrambling dynamics are known to be located in the far-IR regime.^9,26 At variance with the usual approach relying on Fourier transforms of dipole autocorrelation functions we introduce here a multi-channel signal processing approach that belongs to the family of maximum entropy (ME) methods (sometimes also dubbed MaxEnt). The resulting spectra are smooth and can subsequently be decomposed into so-called mode-specific IR spectra after introducing generalized normal modes. Having these modes at hand, it is possible to dissect the mid-IR spectrum of H₂-microsolvated CH₅⁺ and thus to reveal the effect of tagging when using n = 1, 2, or 3 tags.

II. Computational approach

A. Ab initio molecular dynamics simulations

All calculations have been performed using Kohn–Sham density functional theory in the generalized gradient approximation (GGA) within the CPMD program package⁵⁹ as described in ref. 57 (see also ref. 11 and 58); these available AIMD trajectories underly the present analyses. In order to represent faithfully the relatively weak interactions between the CH₅⁺ core and the tagging H₂ molecules as well as their mutual interactions, a dispersion correction⁶⁰ has been added to the particular GGA functional (LDA + B) used. This efficient “LDA + B + D” electronic structure treatment of CH₅⁺·(H₂)_n has been carefully validated by benchmarking⁵⁷ it with respect to coupled cluster and many-body perturbation theory data as published in the comprehensive ESI to ref. 58. The box length was set to 25 a.u. with a plane wave cutoff of 35 Ry and using a 0.34 fs time step for AIMD. In order to prevent the hydrogen molecules from drifting too far away from the CH₅⁺ core, leading eventually to physically meaningful but undesired dissociation events, a soft restraining potential has been applied during the AIMD simulations akin to ref. 58. After sufficient equilibration in the NVT ensemble at 110 K, a total of 10 NVE trajectories per system (with a length of 10 ps each) has been obtained starting from initial conditions that have been sampled from the respective NVT run; note that full scrambling is frozen at 110 K when using classical nuclei as done in these AIMD simulations. All reported data were obtained from averaging the respective properties over all trajectories for each number n = 0, …, 3 of tagging H₂ molecules.

B. Total IR spectra from single-channel maximum entropy method

Assuming classical nuclear dynamics, the total IR absorption cross section α(ω) is obtained from the well-established method of calculating the Fourier transform (FT) of the autocorrelation function of the total dipole moment of the sample (see e.g.ref. 11 and references given therein). Using the corresponding dipole velocities, α(ω) is given bywithwhere Ṁ(t) is the total dipole moment velocity of the sample in a simulation box of volume V at a temperature T and n(ω) is the refractive index (which is unity in vacuum); we note in passing that using this charge-current formulation is mathematically equivalent to using dipoles (after considering the resulting 1/ω² prefactor due to the Fourier transform of the time derivatives) and moreover required for charged (sub-)systems. Furthermore, the key computational quantity, I(ω), is directly proportional to the experimental observable α(ω) when n(ω) = 1. Here, it is assumed that the medium is isotropic, hence the orientational averaging over all directions with the factor of 1/3 and the scalar product of the dipole moment velocities in I(ω). It is mentioned in passing that the specific prefactor corresponds to what is often called the “harmonic quantum correction factor”, which follows from the quantum/classical correspondence in statistical mechanics and has been demonstrated to perform well in molecular IR spectroscopy;⁶¹ recall that the ω² factor resulting from the harmonic quantum correction factor is exactly canceled by the 1/ω² term stemming from the use of dipole derivatives.

Instead of directly evaluating eqn (2) with the time discretized autocorrelation function, a fundamentally different approach is followed by the maximum entropy (ME) method. It is based on the concept of computing the IR spectrum by using a set of model parameters that are determined from the available data (in our case these data are the x, y, z components of the total dipole velocity vector Ṁ as obtained from our trajectories). In the FT approach, the problem arises of how the edges of the available data are treated. The ME method deals with this problem by stating that the best assumption to the data outside the available range is to treat it as unbiased as possible, leading to the equivalent statement of maximizing the information theoretic entropy, being a measure of ignorance in the sense of Shannon and Jaynes, under the constraint that the autocorrelation function agrees with the existing data; we refer to standard literature^62–64 for derivation and discussion.

The so-called single-channel ME approach directly yields I(ω) and thus the absorption cross section as

Here, σ_(p)² denotes the residual variance, p is the order to which the model parameters are estimated, γ^(p)_k are the model coefficients for the chosen order and Δt is the time interval between sampling points along the trajectory. It is noted in passing that “single-channel” in the context of ME approaches refers to the analysis of scalar data, e.g. the x, y, z components of the dipole velocity vector in the present case, being treated as separate statistically independent scalar time series with the same spectrum due to the assumption that no preferred Cartesian direction exists in view of the isotropy of the medium.

The model coefficients and residual variance can be determined by using the Burg algorithm,^65–67 which takes advantage of a recursive scheme for efficient computation. Using this algorithm, the residual variance is obtained from

σ_(p)² = σ_(p−1)²(1 − (γ^(p)_p)²) (4)

andwhere the sum over i denotes the summation over the different directions (i.e. x, y, z) and trajectories, while N_i is the corresponding number of available frames. In the simplest case of just a single dipole velocity trajectory, i just indexes the three cardinal directions and N_i is the length of this trajectory. Note that we assume <Ṁ(t)> = 0. All model parameters can be determined with the help of a subset of model parameters (the special case of k = p), the so-called reflection coefficients,whereHere a^(p)_i(t) is the forward prediction error and b^(p)_i(t) the backward prediction error defined byandbut determined in the algorithm recursively via

a^(p)_i(t) = a^(p−1)_i(t + 1) − γ^(p)_pb^(p−1)_i(t), (11)

b^(p)_i(t) = b^(p−1)_i(t) − γ^(p)_pa^(p−1)_i(t + 1). (12)

Note that a⁽⁰⁾_i(t) = b⁽⁰⁾_i (t) = Ṁ_i(t) for 1 ≤ t ≤ N_i, and γ^(p)₀ = −1 for all p. The model coefficients of order p, i.e. γ^(p)_k, are determined by the coefficients of order p − 1 and the reflection coefficient for order p,

γ^(p)_k = γ^(p−1)_k − γ^(p)_pγ^(p−1)_p−k, (13)

with 0 < k < p. With this information, the spectrum can be readily evaluated using the equations above.

To summarize this concise description, we initialize the algorithm with σ₍₀₎² as the variance of the data series and a⁽⁰⁾_i(t) = b⁽⁰⁾_i(t) with the data set, i.e. dipole velocity trajectories {Ṁ_i(t)}, where the x, y, z components of the vectors Ṁ(t) are treated isotropically as separate scalar trajectories. To obtain the new sets of model parameters, we use the recursive scheme by using eqn (11)–(13), where γ^(p)_p is determined through eqn (6) with the help of eqn (7) and (8). Each set of model parameters is determined for a specific order p, where an upper bound to the order must be determined. We used the combined information criterion⁶⁸ to determine the optimal order p in the single-channel ME approach. At the end the model parameters γ^(p)_k, the reflection coefficient γ^(p)_p and the residual variance σ_(p)², obtained viaeqn (4), are used in eqn (3) to compute the spectrum.

In stark contrast to using the Fourier approach in conjunction with any smoothing, the maximum entropy spectra are computed here by minimizing the assumptions on the available data series, thus implementing the ‘maximum parsimony’ approach or ‘Occam's razor’ principle. Furthermore, based on the algorithmic selection as implemented by us, the number of parameters is autonomously determined by the algorithm itself and thus without any user interference. The resulting spectrum looks similar to a smoothed Fourier spectrum but does not require any specifications and thus (arbitrary) choices of window functions and window lengths. The reported spectra are thus in a sense ‘optimally smoothed’ spectra that are independent from the specific researcher.

C. Multi-channel maximum entropy method

In order to obtain a consistent positive semi-definite cross-correlation matrix spectrum as required for the subsequent decomposition of the total spectrum into mode-specific IR spectra, α_k(ω), according to Sec. II D, a multi-channel ME method is needed where we analyse vector data in contrast to scalar data in the single-channel method. For every center j and every component ζ = x, y, z of center j we assign a separate channel, thus resulting in 3n_centers channels and a 3n_centers by 3n_centers matrix of (cross-)correlation spectra, where n_centers is the number of centers. In each channel we have the product v⁽ⁱ⁾_jζ·q_j, where v⁽ⁱ⁾_jζ is the velocity of center j and component ζ in the i-th trajectory and q_j is the charge of center j. These channels are then simultaneously fitted to obtain model coefficient matrices Γ^(p)_k, analogously to the γ^(p)_k in the single-channel method.

We based our particular implementation on the Vieira-Morf algorithm⁶⁹ as described in ref. 70, where we took advantage of the time reversibility of the molecular dynamics propagation. This leads to the following set of equations

Σ^(p) = Σ^(p−1) − Γ^(p)_pΣ^(p−1)Γ^(p)†_p, (14)

Γ^(p)_p = (Σ^(p−1))^1/2(Θ_aa^(p))^−1/2Θ^(p)_ab·(Θ^(p)_aa)^−1/2(Σ^(p−1))^−1/2, (15)

Γ^(p)_k = Γ^(p−1)_k − Γ^(p)_pΓ^(p−1)_p−k, (16)

a^(p)_i(t) = a^(p−1)_i(t + 1) − Γ^(p)_pb^(p−1)_i(t), (17)

b^(p)_i(t) = b^(p−1)_i(t) − Γ^(p)_pa^(p−1)_i(t + 1). (18)

These equations are analogous to the single-channel method, where the Γ^(p)_k are the model coefficient matrices, Σ^(p) the covariance matrix of the residual errors, and a^(p)_i(t) and b^(p)_i(t) the forward and backward prediction errors of the different trajectories indexed by i. Note that similarly to the single-channel case, is the negative unit matrix for all p and the data sets from the molecular dynamics, i.e. the dipole velocity trajectories {Ṁ_i(t)} enter via

(a⁽⁰⁾_i(t))_jζ = (b⁽⁰⁾_i(t))_jζ = v⁽ⁱ⁾_jζ(t)·q_j, (19)

where N_i is again the length of the i-th trajectory.

Finally,

allows one to compute the dipole velocity cross-correlation matrix from eqn (25), needed in the next section to generate the generalized normal modes and thus the mode-specific spectra IR spectra α_k(ω), via

I(ω) = A_p(ω)Σ^(p)A^†_p(ω), (23)

whereand thus the IR absorption cross-correlation matrix in Sec. II D.

D. Mode-specific IR spectra and generalized normal modes from multi-channel maximum entropy cross-correlation analysis

In our previous work, mass-weighted velocities of the nuclei have been used to construct the so-called generalized normal coordinates,^30,31 building upon a scheme that uses internal coordinates. This approach, yielding only the modes which are encoded in the vibrational density of states, has been extended subsequently⁷¹ to yield mode-specific IR spectra, instead of mode-specific power spectra as before, by using Cartesian coordinates in conjunction with weighting the particle velocities by the effective charge q_i. In contrast to our previous methods,^30,31 the velocities of both nuclei and orbitals are considered here,⁷¹ where the latter are described in terms of the time derivatives of the moving positions of the centers of the corresponding maximally localized Wannier functions (see e.g. Sec. 7.2 in ref. 25) that have been determined along the AIMD trajectory from the underlying electronic structure. Including the electronic degrees of freedom in this way⁷¹ allows us to assign IR intensities to the individual mode-specific absorptions (instead of just the intensities according to the vibrational density of states when using only the nuclear degrees of freedom). Our previous approach⁷¹ has been further improved in the current work (i) by using the ME methodology of a time-reversible algorithm to enforce the proper symmetry of the correlation functions in conjunction with (ii) an automatic order-selection algorithm that ensures the optimal determination of the number of required parameters and thus ‘optimal smoothing’ independently from any user choice.

Within our latest decomposition approach,⁷¹ the dipole velocity cross-correlation matrix is given by

where v_iζ(t) is the ζ-th component of the velocity of the i-th center (being either a nucleus or a Wannier function center) at time t rotated into the molecular reference frame (which is given by a suitable reference structure as described in Sec. II E herein). The off-diagonal elements of this cross-correlation matrix, I_iζ,jξ(ω), are minimized by approximate (Jakobi) diagonalization,

The thereby obtained diagonal elements represent the modes while the remaining off-diagonal elements are a measure of the residual correlations.

Finally, the mode-specific IR absorption cross sections are given by the functions α_k,

α_k(ω) = Λ_k,k(ω)(d_k·d_k), (27)

which together with the total off-diagonal cross term,add up to the usual total IR spectrum,here, d_k are effectively the dipole direction vectors of the k^th generalized normal mode as determined from the unitary transformation matrix D.

As a result of this systematic spectral decomposition, the total IR absorption cross section can be “re-synthesized” to a good approximation as a sum of mode-specific contributions α_k(ω),

that carry proper mode-resolved IR intensities. This allows one to analyze the contributions of specific modes to the total IR spectrum by quantifying their relative importance in determining the overall IR lineshape.

E. Reference structures

In order to compute generalized normal modes as described in the previous section, it is essential to construct reference structures with respect to which the particle displacements and thus these modes can be defined. The methodology described in ref. 30 can only be applied to cases where a single reference structure is sufficient, which is relevant only to quasi-rigid molecules with a well-defined unique equilibrium structure (such as CH₄ for instance). In the case of CH₅⁺·(H₂)_n, trying to construct a single reference structure is obviously not a useful concept in view of the existence of several energetically close-lying isomers as depicted in Fig. 1. It is therefore expected that several such structures, possibly even all of them, are visited during the dynamics of these weakly-bound complexes. For instance, there are three such dynamically accessible structures for both n = 1 and n = 2 according to Fig. 1, being distinct in real space, whereas only a single energetically low-lying structure is relevant for the triply tagged species. Thus, each molecular dynamics trajectory must be split into well-defined segments which are associated to a specific reference structure (constructed as explained below).

However, each permutation of the atom labeling due to a successful hydrogen scrambling event within a distinct real space structure formally creates another reference structure when computing and diagonalizing the cross-correlation matrix defined viaeqn (25); cf. the molecular symmetry group concept⁴³ where Longuet-Higgins’ permutation-inversion groups are a special case. This permutation sampling would chop up each trajectory segment that is associated to a specific real space structure into pieces that are associated to a distinct atom labeling pattern, thus creating additional reference structures for mode analysis when using eqn (25). This trivial effect of only changing the permutation state of reference structures can be eliminated by transforming each distinct atom labeling pattern back to our standard ordering²⁷ as defined in Fig. 1a.

In the present case, the starting point of the computation of the reference structures were the respective minimum energy structures of the systems of interest, all of which are shown in Fig. 1. The trajectories were analyzed and split into segments that are close in structure to one of these energetic minimum structures, thus considering not only the global but also all local minima. To this end, all frames were transformed into standard order (which is defined according to the labeling in Fig. 1a and thus includes permutational ordering) and rotated using root-mean-square distance criteria to match as closely as possible any of the available energetic minimum structures. The center of mass motion has been taken out by transformation into the barycenter coordinate system, i.e. the square-root mass-weighted center (all this has been carried out in Cartesian coordinates).

This procedure generates trajectory segments that are associated to the available minimum energy structures. Subsequently, the configurations of these trajectory segments were averaged to create what we call the reference structure of the respective trajectory segment. They thus define structural motifs also in the case of fluxional molecules, where different minima are regularly populated due to large-amplitude motion. This procedure led in all cases to very symmetric reference structures. Moreover, these reference structures turn out to be, apart from slight changes in bond lengths, identical to the corresponding energetic minimum structures and are therefore not shown here. These reference structures and the associated trajectory segments are then used in order to compute the generalized normal modes that describe best the dynamics in these trajectory segments and thus the IR spectra corresponding to the respective structural motif.

For bare protonated methane at 110 K, a single reference structure turned out to be sufficient (after transforming the instantaneous permutation pattern to standard ordering in each trajectory frame) which was essentially the global minimum energy structure, see Fig. 1a. For the singly tagged system, i.e. CH₅⁺·(H₂)_n with n = 1, two reference structures were required. Analysis of the trajectories showed that roughly 70% of the generated configurations could be assigned to the minimum energy structure shown in Fig. 1b and ≈30% to the one depicted in panel c, whereas the structure in panel d was not observed.

In the case of the doubly tagged n = 2 system, about 97% of the trajectory frames could be assigned to the minimum energy structure shown in panel e, which corresponds to the global minimum. The remaining 3% of the frames were assigned to the structures shown in panel f and panel g. Given its overwhelming population, only the global minimum structure e has been considered in the mode decomposition of the IR spectrum. Finally, the triply tagged system has only one reference structure as displayed in Fig. 1h.

III. Results and discussion

In the following, the total IR spectra α(ω) of the investigated microsolvated systems, i.e. CH₅⁺·(H₂)_n for n = 1 up to 3, are analyzed and assigned in terms of the additive mode-specific absorption cross sections α_k(ω) in the frequency range of 1000–4500 cm⁻¹. As already alluded to in the Introduction, the present analysis focuses specifically on tagging effects on the stretching and bending modes of protonated methane which are located in the mid-IR frequency regime. In addition, the limiting n = 0 case corresponding to the bare CH₅⁺ is analyzed as well using the identical methodology in order to provide a consistent reference to discuss changes of the generalized modes and thus of the total lineshape due to microsolvation as a function of n. Our frequency window includes the bending and stretching modes of protonated methane and its H₂-tagged species, which allows us to discuss how H₂-tagging affects the bending and stretching modes of bare CH₅⁺. In the subsequent presentation, we will address bare protonated methane first followed by discussing the tagged species. For the latter, we analyze the n = 3 case first since that turned out to be most easily understood, followed by n = 2 and finally by n = 1.

A. Bare protonated methane

We first of all compare the mid-IR spectrum of bare CH₅⁺ as obtained using the single-channel maximum entropy method to the one obtained from the Fourier transform at two temperatures, 110 and 300 K, see Fig. 3. It is clearly seen that the spectra computed using the ME and FT methods agree well with each other. In addition to reproducing all features of the FT spectra (without introducing additional modulations) the ME method is seen to reduce the noise level drastically using the same statistics, thus resulting in smoother curves. Although similar results could be obtained by applying an adequate windowing function when computing the FT spectrum, we want to stress that the presented ME spectra are not based on such arbitrarily chosen smoothing procedures, but rather follow from an algorithmic procedure that selects the parameters based on an automatic minimization procedure that results in an ‘optimally smoothed’ spectrum. Furthermore, although the mid-IR spectrum at 110 K shows more structure as expected, it is qualitatively comparable to the one at 300 K.

Fig. 3 Mid-IR spectra of bare protonated methane at 300 and 110 K in the top and bottom panel, respectively, computed from AIMD simulations with classical nuclei using the single-channel ME method (blue and circles), see Sec. II B, and the usual FT approach (black). The orders for the ME spectra computations were 591 for 300 K and 700 for 110 K. Note that the FT spectra are off-set for the sake of presentation.

In addition to computing the total IR spectra, we also computed the individual modes using our new multi-channel ME approach. The obtained generalized normal modes are shown in detail in the ESI† and agree well with the previously extracted modes.²⁷

In the bare case, five stretching modes above 2000 cm⁻¹ are extracted with our novel procedure, see Fig. 3b and Table 1 for frequencies, in accordance with earlier results²⁷ that were, however, simply based on power spectra. The extracted modes are, in order of decreasing frequency, the asymmetric as well as symmetric stretching modes of the tripod hydrogen atoms (atom 4 and atom 5) abbreviated as t^a and t^s, the stretching mode of the hydrogen atom at the in-plane position (atom 3) called tⁱ, and the symmetric as well as asymmetric stretching modes of the moiety hydrogen atoms (atom 1 and atom 2) called m^s and m^a, see Fig. 1a for the atom labels. The tripod stretching modes are of similar intensity and shape while the m^s mode is significantly wider and of lower intensity compared to the m^a mode. These modes are visualized in the ESI† by means of dynamical animations.

Table 1 Generalized normal mode frequencies in cm⁻¹ for bare CH₅⁺ and the three microsolvated species, see Fig. 6 and 7 for the mode labelings. The n = 1 case can only be described using two reference structures, see text and Fig. 1b and c, which are labeled here as 1–I and 1–II, respectively

Mode	CH5^+	CH5^+(H2)1–I	CH5^+(H2)1–II	CH5^+(H2)2	CH5^+(H2)3
H^i2	—	—	—	—	4374
H^om2	—	—	4336	4343	4343
H^im2	—	4327	—	4334	4340
t^a	3309	3291	3294	3294	3304
t^s	3218	3198	3203	3200	3203
t^i	3080	3081	3099	3098	3058
m^s	2706	2710	2718	2673	2660
m^a	2429	2288	2293	2242	2286
b^w	1500	1501	1509	1513	1498
b^t	1477	1466	1473	1472	1473
b^s	1430	1442	1444	1446	1439
b^r	1276	1284	1309	1308	1324
b^b	1176	1207	1161	1176	1171

---

All bending modes of bare CH₅⁺ are located between 1000 and 1800 cm⁻¹. We extracted five bending modes in this frequency range, see Fig. 3b and Table 1 for frequencies, which also agree with the generalized normal modes that have been obtained earlier²⁷ but using standard FT instead of multi-channel ME techniques. Thus, the following modes were extracted and will be investigated as to their response to microsolvation of CH₅⁺ with n = 1 to 3H₂ molecules: wagging bend (b^w), twisting bend (b^t), scissoring bend (b^s), rocking bend (b^r), and butterfly bend (b^b). These modes are also animated in the ESI† for the sake of direct comparison to the corresponding modes in case of the n > 0 species.

At variance with the previous approach,²⁷ however, the current mode spectra are not just projected vibrational density of states, but carry instead the proper IR intensity of each distinct mode, i.e. they are the mode-specific IR absorption cross sections α_k(ω) according to eqn (27). Since the remaining cross-correlations according to eqn (28) are small, these mode-specific spectra approximately add up to the corresponding total IR spectrum, α(ω). It is demonstrated in Fig. 4 that, indeed, this approach allows one to re-synthesize the total IR spectrum “mode by mode”, which is shown from bottom to top starting with the highest-frequency mode and progressing to the lowest one according to Table 1.

Fig. 4 Mode-by-mode re-construction of the total IR absorption cross section α(ω) of bare CH₅⁺ at 110 K based on successively adding up all mode-specific IR absorption cross sections α_k(ω) according to eqn (27) from bottom to top starting with the highest-frequency modes, see text. The modes were obtained by diagonalization of the cross-correlation matrix (cf.eqn (23)) computed by the multi-channel ME method. The topmost spectrum is the total IR spectrum obtained from the single-channel ME method (where total rotations and all index permutations have been taken out whereas they are included in Fig. 3), the second spectrum from top (thick line) is the re-synthesized total IR spectrum according to eqn (30), whereas the third curve from top (dotted line) quantifies the residual cross-correlations according to eqn (28). Note that there are minor deviations between the single-channel and the total spectrum that result from the sum of all modes obtained by the multi-channel approach. Deviations between the two methods are expected when using finite time-series since the multi-channel method has access to more data in the form of the different channels and thus can better resolve fine structure compared to the single-channel approach (a close-up of the bending region is shown in the ESI†). Moreover, the more detailed multi-channel analysis also leads to a different number of parameters as a result of the ME methodology as determined autonomously by the algorithm given the larger database.

Thus, the decomposition technique introduced here allows one to work out – at a fully quantitative level – which modes contribute to what extent to peaks, shoulders or wings of the total lineshape, which appears to be an interesting theoretical spectroscopic analysis tool beyond the specific case discussed herein.

B. Tagged protonated methane: CH₅⁺(H₂)_n

The mid-IR spectra of the three microsolvated systems compared to the bare case are shown in Fig. 5. Adding hydrogen tags to bare CH₅⁺ (cf. bottom-most red curve with circles) clearly affects the mid-IR spectrum in systematic ways. Obviously, each added hydrogen tag leads to an additional signal at about 4500 cm⁻¹ that is caused by the respective hydrogen stretching mode. This mode is IR forbidden for an isolated H₂ molecule, but it is seen here due to the interaction of the tag with the CH₅⁺ core, thus breaking the symmetry of the attached H₂via electronic polarization effects. Besides this striking but trivial change compared to bare protonated methane, already upon adding the first hydrogen tag the asymmetric moiety stretch signal m^a, which is at about 2400 cm⁻¹ in the bare case, shows a considerable red-shift by more than 100 cm⁻¹ together with a significant gain in intensity compared to all other bands. Both, frequency shift and intensity increase persist upon tagging CH₅⁺ with more than a single H₂ molecule and are most pronounced for the CH₅⁺·(H₂)₂ species. The other stretching mode of the H₂ moiety in CH₅⁺, being the symmetric stretch m^s at around 2700 cm⁻¹ in the n = 0 case, is also found to gain systematically more intensity, which is most significant for the CH₅⁺·(H₂)₂ species, and red-shifts up to 46 cm⁻¹ in the case of n = 3. This is easily understood at the qualitative level since the hydrogen molecules stay attached at the moiety (similar to the structures depicted in Fig. 1) during the dynamics at 110 K and thus perturb mostly the moiety dynamics; hence a pronounced change in the bands caused by the moiety itself is expected.

Fig. 5 IR spectra of microsolvated protonated methane (with three, two and one H₂ molecules attached from top to bottom) compared to bare protonated methane all at 110 K and computed with the multi-channel ME method. See Sec. II C.

The tripod stretches show continuous changes when more tags are added. In the bare system, the three tripod stretches are distinct from each other and this distinction is conserved upon tagging. Only subtle changes within the tripod stretching modes are observed that mainly involve intensities and minor frequency shifts. Since the first tagging hydrogen molecule can be situated at two different positions we present two spectra in the case of using one tag. They show subtle differences in the tripod stretching signals, in particular the intensity of the t^a mode shows significantly different intensities depending on the distinct tagging site. Similarly, the m^s mode has a higher intensity in the situation where the tag resides at site 2. Further details will be discussed in later sections when focusing on the specific mode absorptions. Adding a second messenger hydrogen molecule leads to a better distinction of the in-plane tⁱ and asymmetric tripod t^a stretching modes while the symmetric tripod stretching mode t^s is almost unchanged. Finally, adding the third tagging hydrogen molecule results in a very much broadened signal of the in-plane tripod stretching mode tⁱ while the other two resonances, t^a and t^s, remain as sharp as they were before.

How can these observations be qualitatively understood? The first hydrogen molecule is expected to be attached at the inner moiety hydrogen atom (being atom 2 in Fig. 1a) according to the energetic global minimum structure of CH₅⁺·H₂ shown in panel b. Due to the intermolecular interaction a change of the symmetric and asymmetric moiety stretches is expected. Furthermore, a change of the in-plane tripod stretching mode is expected since this atom is now better localized in space, which is reflected in the sharper peak. Since pseudorotation is still possible in this situation, the tⁱ and the m^a modes are broader compared to the case with two hydrogen tags. The second hydrogen tag would then be located at the remaining moiety hydrogen atom 1 according to the global minimum structure of the CH₅⁺·(H₂)₂ species, see panel e in Fig. 1. Thus, internal rotation as sketched in Fig. 2b will be easily possible in view of the symmetric tagging scenario of the moiety. At the same time, the in-plane hydrogen atom 3 will be better localized because atom displacements as illustrated in Fig. 2a, and thus pseudorotation leading to partial scrambling, will be hindered resulting in a sharp resonance of tⁱ. Finally, the third H₂ molecule would be added to the in-plane hydrogen atom 3, see panel h in Fig. 1. Thus, once three tags are added, all three in-plane atoms 1, 2, and 3 in CH₅⁺·(H₂)₃ are interacting equally with H₂ molecules, which is a scenario that is equivalent to that in the bare CH₅⁺ limit where all three in-plane atoms are untagged, compare panel h to panel a in Fig. 1. This approximate equivalence in the presence of microsolvation, which is distinctly different from all other tagging patterns that are collected in Fig. 1, allows for easy pseudorotational motion according to the process sketched in Fig. 2a and eventually enables “partial scrambling” in the n = 3 case. This topological argument suggests that partial scrambling should be more easily possible for the n = 3 case compared to both n = 1 and 2, the onset of which in terms of small-amplitude motion leads to a broader signal of the tⁱ and the m^a modes.

Addressing finally the perturbations of the bending region, it appears at first sight as if these resonances feature only some mild changes due to tagging. However, it is difficult to discern at the resolution level of Fig. 5 the much more intricate changes of the bending band of CH₅⁺ due to its interaction with the tags. Instead of discussing this here, we refer to Fig. 7 and to the corresponding detailed analysis of the mode-specific IR spectra to be presented together with the mode assignments.

At this stage we conclude, based on visual inspection of the computed total IR spectra, that tagging of protonated methane by a few H₂ molecules has a significant effect on all parts of its mid-IR spectrum, from the characteristic triple-peaked tripod stretching region at highest frequencies 3000–3500 cm⁻¹ to the two moiety stretches around 2500 cm⁻¹ down to the bending regime at 1100–1600 cm⁻¹. However, the assessment of spectral changes due to tagging is still qualitative as it is simply based on directly comparing the peaks of the microsolvated cases to those of the bare reference, CH₅⁺. Thus, it remains to be seen if our spectral decomposition and mode assignment techniques are also able to disentangle the IR spectrum of protonated methane once it is microsolvated.

C. Mode assignment of the microsolvated cases

The analysis was carried out as described in the methods section leading to the modes shown in Fig. 6 and 7 for the stretching and bending bands, respectively. The corresponding frequencies are compiled in Table 1 whereas the atomic displacements that generate these localized modes are visualized in the ESI† in terms of dynamical animations.

Fig. 6 Mode decomposition and assignment of microsolvated protonated methane with labeled stretching resonances as follows. Hⁱ₂: H₂ stretch at in-plane site, H^im₂: H₂ stretch at inner moiety site, H^om₂: H₂ stretch at outer moiety site, t^a: antisymmetric out-of-plane stretch in tripod, t^s: symmetric out-of-plane stretch in tripod, tⁱ: in-plane stretch in tripod, m^s: symmetric stretch in moiety, m^a: antisymmetric stretch in moiety; the vertical dotted lines mark the corresponding resonance frequencies of bare CH₅⁺ shown in panel e. The panels a–d show the microsolvated species from top with 3 tags to bottom with 1 tag. The nuclear displacements that underlie these generalized normal modes are visualized in the ESI† by dynamical animations.

Fig. 7 Mode decomposition and assignment of microsolvated protonated methane with labeled bending resonances as follows. b^w: wagging bend b^t: twisting bend b^s: scissoring bend b^r: rocking bend b^b: butterfly bend; the vertical dotted lines mark the corresponding resonance frequencies of bare CH₅⁺ shown in panel e. The panels a–d show the microsolvated species from top with 3 tags to bottom with 1 tag. Note that some mode-specific spectra have been scaled to allow for a better visualization using scaling factors as follows. b^w: 10, b^t: 100, b^s: 10, b^r: 1, and b^b: 1. The nuclear displacements that underly these generalized normal modes are visualized in the ESI† by dynamical animations.

In the following, the stretching and bending modes will be discussed separately starting with the former ones. Within each mode class we proceed from the simplest case, which is the one with three messenger species attached being closest to bare CH₅⁺, to the most difficult one. The latter is the situation when a single H₂ molecule interacts with bare CH₅⁺ since dynamical site-exchange events become significant. Bare protonated methane, which serves as our reference, has already been reviewed in terms of its mode assignment in Sec. III A to which we now refer.

1. Stretching modes.
Three tags. After adding three messenger H₂ species to bare CH₅⁺, we obtain three modes in the region of 4500 cm⁻¹ caused by the stretching motion of the three hydrogen molecules. These species, located at the atoms 1, 2 and 3, are clearly distinguished from each other after having disentangled the total spectrum, namely the stretching mode due to the H₂ molecules attached to the in-plane position (Hⁱ₂) and to the inner (H^im₂) and outer (H^om₂) moiety sites of CH₅⁺, being atoms 1, 2, and 3, respectively, according to the scheme defined in Fig. 1a. The H₂ stretching motion from the tag bound to site 2, i.e. H^im₂, has the lowest quasi-intramolecular frequency, indicating that the interaction between the inner moiety atom and the tag is the strongest.

Next lower in frequency are the tripod stretching modes t^a (3304 cm⁻¹), t^s (3203 cm⁻¹) and tⁱ (3058 cm⁻¹), which are red-shifted w.r.t. bare CH₅⁺. The intensity of the tⁱ mode is higher w.r.t. the bare reference system and the other two tripod stretches, which reflects the additional interaction of the tripod atoms with the hydrogen tagging molecules. The increased intensity of the tⁱ mode indicates a better defined stretching mode compared to the bare case.

The symmetric and asymmetric stretching modes of the moiety, m^s (2660 cm⁻¹) and m^a (2286 cm⁻¹), respectively, show a red-shift w.r.t. the bare case (2706 and 2429 cm⁻¹, respectively), which is more pronounced for the asymmetric case. The intensities of both modes match better, compared to the previous cases, due to a significant reduction in the peak height of the m^a mode. This is in contrast to the bare and the other microsolvated cases where the m^s mode always has a significantly lower intensity. This shows that these modes are now more pronounced indicating that the structure must be better defined overall. We note that for the sake of simplicity and given the mostly rather sharp resonances, we use the peak height of the mode-specific absorption cross section as a proxy for intensity, although the integral of the full peak should be used in principle.

Two tags. The set of stretching modes obtained for CH₅⁺·(H₂)₂ is shown in Fig. 6b. We obtained the expected two hydrogen stretching modes from the two H₂ tags, labeled H^om₂ and H^im₂, and they also show the same trend as before, i.e. the H^om₂ mode (at 4343 cm⁻¹) is located at a higher frequency than the H^im₂ mode (4334 cm⁻¹).

When referring to what we have seen for the CH₅⁺·(H₂)₃ complex, the removal of the tag located at the in-plane position, i.e. atom 3, results in distinct changes in the C–H stretching modes of the CH₅⁺ core for the n = 2 species. We observe that now the t^s and t^a modes are similar in intensity to each other, while the t^s mode remains similar to the triply tagged case. It is furthermore seen that the tⁱ mode (3098 cm⁻¹) is blue-shifted in contrast to the t^a (3294 cm⁻¹) and t^s (3200 cm⁻¹) modes that are, as in the n = 3 case, red-shifted w.r.t. the bare reference molecule. This indicates that protonated methane is differently affected upon microsolvating it when using two or three solvent species.

Moreover, also the lineshape of the moiety stretching modes is very different from what is observed in the triply tagged case, whereas it is similar to the situation found for the bare system, i.e. the m^s mode is of significantly lower intensity than the m^a mode. Although different in their lineshape, the modes of the n = 2 species show the same frequency shifting pattern as those of the triply tagged case, both m^s (2673 cm⁻¹) and m^a (2242 cm⁻¹) modes being red-shifted.

Since sites 1 and 2 are tagged in the n = 2 case whereas site 3 is not, partial scrambling will be less likely to happen since this process would result in an energetically unfavorable situation given that the resulting structure would be 1.19 kcal mol⁻¹ higher in electronic energy (as provided in the caption of Fig. 1). Assuming partial scrambling to be inactive, we would expect from this observation the tⁱ mode to feature a very narrow and more harmonic type of C–H stretching vibration. Internal rotation, on the other hand, will still be allowed since the structure resulting after this process will be identical to the one before this rotation after a slight readjustment of the two weakly attached H₂ molecules. This situation might explain the different behavior of the tⁱ mode shifting as well as the lineshape changes of the moiety stretching modes and the t^a mode.

One tag. In the most intricate case of the CH₅⁺·(H₂) complex, the single tagging molecule can be attached to sites 1, 2, and 3 as depicted in Fig. 1 in panels b, c, and d, respectively. It turned out that essentially only the two moiety sites 1 and 2 were found to be populated during the AIMD simulations (for about 70 and 30% of the total trajectory length, respectively), whereas the H₂ molecule was not observed to be attached to site 3, belonging to the tripod. This finding is consistent with the energetic ordering of the optimized structures,⁵⁶cf. caption to Fig. 1, where isomer c involving site 1 is only 0.35 kcal mol⁻¹ above the global minimum, whereas attachment of H₂ at the tripod site 3 is as much as 1.77 kcal mol⁻¹ higher in energy. Thus, only the two reference structures corresponding to panels b and c, respectively, reflecting directly the two lowest-energy static structures, must be used in the following dynamical vibrational analysis.

By using two reference structures, two sets of modes are obtained as shown in Fig. 6d for reference structure I, where the tag is located at site 2, and in Fig. 6c for reference structure II involving the H₂ attachment to site 1. In both cases, a single stretching mode corresponding to the hydrogen molecule is extracted. The H^im₂ mode in Fig. 6d is located at slightly lower frequency (4327 cm⁻¹) than the H^om₂ mode (4336 cm⁻¹) in Fig. 6c. This is consistent with the weaker interaction of the H₂ tag when attached to site 1 compared to site 2, which is also indicated by a shorter intermolecular bond length when involving site 2 (1.87 Å) versus site 1 (1.89 Å).

In both reference situations, the modes look very similar at first glance, but subtle differences are still discernable upon more careful inspection. The t^a stretching mode observed in Fig. 6c due to reference structure II shows a long tail on its red-wing that overlaps with the blue wing of the t^s mode. In contrast, in reference structure I shown in Fig. 6d all three tripod stretching modes are well separated. The red-shift of the t^a (3291 and 3294 cm⁻¹) and t^s (3198 and 3203 cm⁻¹) modes is comparable to each other. Comparing the tⁱ mode in both sets reveals distinct differences regarding shifting and lineshape. In reference structure I, the tⁱ mode is with 3081 cm⁻¹ essentially unshifted whereas in the case of reference structure II this mode is blue-shifted (3099 cm⁻¹), which is similar to what has been observed in case of the doubly tagged case (3098 cm⁻¹).

Also the moiety stretching modes again look quite similar to each other as they both show a similar ratio as the bare system, where the m^s mode has a considerably lower intensity than the m^a mode. Note, however, that the m^a mode is subject to a significantly increased linewidth and that the m^s mode shows only a minor intensity increase compared to the bare system in panel e of Fig. 6. Both moiety stretches are red-shifted, where the m^a mode shows a more pronounced red-shift than the m^s mode. Additionally, the m^a mode shows the smallest red-shift in reference structure II (2293 cm⁻¹), whereas this mode is stronger red-shifted in the doubly tagged case (2242 cm⁻¹).

Qualitative trends. To sum up the discussion of the stretching modes which determine the overall shape of the high frequency part of the IR spectrum of CH₅⁺·(H₂)_n, we conclude that tagging affects the associated peaks in distinct ways as follows. First of all, the t^s mode is relatively uninfluenced by adding the first and second tagging molecule to the bare species, while the t^a mode changes readily upon microsolvation. Upon adding the third tagging H₂ molecule, both the t^a and t^s mode show the same behavior of reduced intensity and red-shifting w.r.t. bare CH₅⁺. The tⁱ mode, on the other hand, displays a distinctly different behavior as a function of the number of tagging molecules: one finds essentially no shift once the bare system is tagged at its so-called inner moiety site or a red-shift if it is triply tagged, whereas a blue-shift is observed if it is microsolvated by a single hydrogen molecule which is located at the outer moiety site, or alternatively in the case of adding two tags. Focusing on the moiety stretching modes, it becomes clear that the m^a mode shows the most significant change upon microsolvation in form of a strong red-shift. This can be understood when considering the intricate three-center electronic bonding situation in (bare) protonated methane^13–15 as summarized in the Introduction. The tripod stretching modes involve exclusively the usual two-center two-electron bonding for each of the three C–H covalent bonds involved. These bonds are rather strong¹⁵ and, therefore, are not that much affected as a result of non-covalent interactions with the tagging hydrogen molecule (at the in-plane location). In stark contrast, when considering the moiety stretching modes, it is a three-center two-electron bond that binds the two hydrogens involved to the carbon site. This bonding type is clearly weaker than the two-center one¹⁵ and, therefore, these softer C–H bonds are more strongly influenced when forming non-covalent bonds with the tag, which explains a more pronounced frequency red shift when they are tagged. The impact on the m^s mode is different for the singly tagged species as a blue shift is observed, which is explained by the asymmetric nature of the tagging in this limiting case, which hinders the full displacement of the tagged hydrogen moiety atom compared to the bare case. The moiety atoms are again symmetrically tagged when adding more than one hydrogen molecule, hence the expected red-shift of the modes is observed as in the m^a case.

These features herald at the level of mid-IR spectra the purely structural changes upon microsolvation that have been obtained recently from time-independent ab initio path integral simulations at low temperature,⁵⁸ namely that the quantum localization of CH₅⁺ due to microsolvation is very similar in the case of adding either one and three hydrogen molecules, while the case of two H₂ is distinctly different. This suggests, already at this stage, that high-resolution messenger vibrational spectroscopy should be a suitable tool^47–50 to uncover even subtle modifications of hydrogen scrambling due to microsolvation of protonated methane.

2. Bending modes. The microsolvation-induced changes of the bending dynamics of protonated methane, which appeared rather minor at the level of the full spectrum, cf.Fig. 5, turn out to be particularly rich when analyzed on the proper frequency and intensity scales as done in Fig. 7.
Three tags. In the case of attaching n = 3H₂ tags to protonated methane, the b^w and b^t modes show a strongly decreased intensity together with a red-shift (to 1498 and 1473 cm⁻¹, respectively) w.r.t. the parent system, CH₅⁺. Furthermore, these resonances are changed from relatively sharp to broad signals, which clearly reflects the significant impact of the tagging molecules on the intramolecular bending dynamics of the microsolvated CH₅⁺ core molecule.

The b^s mode is blue-shifted (to 1439 cm⁻¹) with a similarly low intensity as in the bare system. Comparison with the other microsolvated cases studied shows that both in the bare and the n = 3 case the b^s mode is of lowest intensity while in the CH₅⁺·(H₂)₂ species it gains intensity and turns out to be strongest in the n = 1 limit. The b^r mode, on the other hand, is blue-shifted (by 48 cm⁻¹) but comparable in lineshape to the bare system. The b^b mode is red-shifted (to 1171 cm⁻¹) and mainly conserved in shape with an increase in intensity.

Two tags. Upon removing the tagging molecule located at atom 3, the b^t mode of the doubly tagged species has a similar lineshape as the one in the singly tagged species in reference structure II including a red-shift (1472 cm⁻¹). The lineshape changed from that of a rather broad feature in the triply tagged case to a narrow signal for n = 2. No change in the lineshape is observed for the b^w mode, which shifts to the blue (1513 cm⁻¹) compared to bare system. Also the b^r mode has a similar structure but shows a blue-shift (1308 cm⁻¹). While the b^s mode is of low intensity in the bare and triply tagged systems, it is found to gain intensity in the n = 2 case.
One tag. As in the case of the stretching band, two sets of bending modes are obtained depending on where the single H₂ molecules is attached at the moiety, i.e. at site 2 (yielding reference structure I) or at site 1 (reference structure II).

The b^w mode shown in Fig. 7c for reference structure II is broad and of low intensity, which is strikingly different from the one observed for reference structure I in Fig. 7d, where it is found to become an extremely strong signal that is even more intense than in bare protonated methane. The b^t mode of reference structure II is similar to the one obtained from the doubly tagged case with a similar shift (1472 cm⁻¹) and lineshape, which underlines the similarities between those two tagging cases akin to what has been found when analyzing the stretching bands. On the other hand, the b^t mode obtained from reference structure I is very low in intensity. For both tagging sites, the b^s mode has gained intensity w.r.t. the bare case. While the b^r mode is essentially similar in all tagging cases in terms of intensity, it varies only in the magnitude of the blue-shift compared to the corresponding CH₅⁺ resonance. The shift is least pronounced in the case of reference structure I.

IV. Conclusions and outlook

The effect of microsolvation by one, two and three H₂ molecules on the bending and stretching bands in the mid-infrared absorption spectrum of tagged protonated methane, CH₅⁺·(H₂)_n, has been investigated systematically in comparison to bare CH₅⁺ serving as our internal reference. The spectra in the frequency window from 1000 to 4500 cm⁻¹ have been computed based on ab initio molecular dynamics trajectories obtained at 110 K with classical nuclei using an advanced maximum entropy methodology. The maximum entropy spectra are computed here by strictly minimizing the assumptions on the available data series, thus implementing the ‘maximum parsimony’ approach or ‘Occam's razor’ principle at the level of the algorithm. The reported spectra are thus in a sense ‘optimally smoothed’ spectra that are independent from the user, in contrast to standard Fourier smoothing procedures that rely on specifications and thus (arbitrary) choices of window functions and window lengths. Moreover, given the same trajectory data, this method is shown to provide smoother spectra than those obtained via Fourier transforming the dipole autocorrelation function, while preserving the correct intensities. This technique is able to process multi-channel data which enables us to extract generalized normal coordinates as well as the associated mode-specific IR absorption coefficients from molecular dynamics simulations. The full spectra can be “re-synthesized” to a good approximation simply as a sum of these mode-specific partial spectra, which allows one to quantify their contributions to the total IR lineshape in a mode-by-mode fashion.

The computed total mid-IR spectra of CH₅⁺·(H₂)_n feature subtle but significant changes w.r.t. bare CH₅⁺, which depend crucially on the number of tagging molecules, n. The case with three tags turns out to be most easily understood since only a single tagging motif, where all three H₂ molecules remain attached to specific proton sites, is found to contribute to the dynamics. In stark contrast, there exist three microsolvation motifs in the n = 2 and n = 1 cases as defined by optimized minimum energy structures. For CH₅⁺·(H₂)₂, however, only the lowest-energy structure is found to contribute to the dynamics to an overwhelming extent. Most difficult in terms of both spectral decomposition and assignment is the n = 1 species, where the two lowest-energy structures are significantly populated. Moreover, they are found to dynamically interchange along the trajectories and produce distinct spectral fingerprints.

The impact of microsolvation on the mid-IR spectrum of bare CH₅⁺ is complex and affects the peak positions and intensities of the well separated stretching and bending regions in characteristic ways. The stretching frequency regime is shaped by three and two modes at higher and lower frequency that can be associated to stretching dynamics involving mainly the tripod and moiety atoms, respectively, of the microsolvated CH₅⁺ core molecule. Also the bending band consists of five modes with vastly different intensities in this case. Quite interestingly, the two microsolvation motifs that are dynamically visited in the n = 1 case lead to pronounced differences in the bending modes in terms of their intensities and frequency shifts.

An important qualitative finding based on our systematic decomposition of the bending and stretching spectra of microsolvated protonated methane is that the pronounced frequency shift of the asymmetric moiety stretching resonance serves as a unique mid-IR marker to distinguish the tagged species from bare protonated methane. The strong frequency response of this moiety stretch due to even weak non-covalent intermolecular interactions with tagging species has been traced back to the three-center two-electron bonding nature which binds two hydrogens to the carbon site, in stark contrast to the three standard two-center two-electron C–H bonds that form the tripod. This observation is particularly useful since the two moiety stretching signals are, even at higher temperatures, well separated features in the mid-IR spectrum of protonated methane.

We close by stressing that multi-channel maximum entropy analysis of (ab initio) molecular dynamics trajectories is a powerful tool that can be applied broadly to microsolvated molecules in order to understand their IR spectra. In particular, the approach allows one to quantify and assign the distinct influence of the microsolvating species on the vibrational dynamics of the parent molecule.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We gratefully thank Alex Witt and Sergei Ivanov for supplying the ab initio trajectories of the microsolvated species that have been used for the present analyses. This work was supported by the Cluster of Excellence RESOLV (EXC 1069) funded by Deutsche Forschungsgemeinschaft and computational resources were provided by HPC@ZEMOS, HPC-RESOLV, RV-NRW, and BOVILAB@RUB.

References

1. V. L. Tal'roze and A. K. Lyubimova, Dokl. Akad. Nauk SSSR, 1952, 86, 909.
2. R. J. Saykally, Science, 1988, 239, 157.
3. T. Oka, Philos. Trans. R. Soc. London, Ser. A, 1988, 324, 81.
4. G. E. Scuseria, Nature, 1993, 366, 512.
5. D. Marx and M. Parrinello, Science, 1999, 284, 59.
6. P. R. Schreiner, Angew. Chem., Int. Ed., 2000, 39, 3239.
7. D. Gerlich, Phys. Chem. Chem. Phys., 2005, 7, 1583.
8. S. Borman, Chem. Eng. News, 2005, 83, 45.
9. P. Kumar and D. Marx, Phys. Chem. Chem. Phys., 2006, 8, 573.
10. A. B. McCoy, Int. Rev. Phys. Chem., 2006, 25, 77.
11. S. D. Ivanov, A. Witt and D. Marx, Phys. Chem. Chem. Phys., 2013, 15, 10270.
12. T. Oka, Science, 2015, 347, 1313.
13. V. Dyczmons and W. Kutzelnigg, Theor. Chim. Acta, 1974, 33, 239.
14. D. Marx and M. Parrinello, Nature, 1995, 375, 216.
15. D. Marx and A. Savin, Angew. Chem., Int. Ed. Engl., 1997, 36, 2077.
16. A. Gamba, G. Morosi and M. Simonetta, Chem. Phys. Lett., 1969, 3, 20.
17. G. A. Olah, G. Klopman and R. H. Schlosberg, J. Am. Chem. Soc., 1969, 91, 3261.
18. S. Ehrenson, Chem. Phys. Lett., 1969, 3, 585.
19. V. Dyczmons, V. Staemmler and W. Kutzelnigg, Chem. Phys. Lett., 1970, 5, 361.
20. P. C. Hariharan, W. A. Lathan and J. A. Pople, Chem. Phys. Lett., 1972, 14, 385.
21. P. R. Schreiner, S.-J. Kim, H. F. Schaefer III and P. v. R. Schleyer, J. Chem. Phys., 1993, 99, 3716.
22. M. Kolbuszewski and P. R. Bunker, J. Chem. Phys., 1996, 105, 3649.
23. H. Müller, J. Noga, W. Klopper and W. Kutzelnigg, J. Chem. Phys., 1997, 106, 1863.
24. D. Marx and M. Parrinello, Z. Phys. D, 1997, 41, 253.
25. D. Marx and J. Hutter, Ab Initio Molecular Dynamics, Cambridge University Press, Cambridge, 2009.
26. O. Asvany, P. Kumar, P. B. Redlich, I. Hegemann, S. Schlemmer and D. Marx, Science, 2005, 309, 1219.
27. S. D. Ivanov, O. Asvany, A. Witt, E. Hugo, G. Mathias, B. Redlich, D. Marx and S. Schlemmer, Nat. Chem., 2010, 2, 298.
28. A. Witt, S. D. Ivanov, G. Mathias and D. Marx, J. Phys. Chem. Lett., 2011, 2, 1377.
29. D. Marx and M. Parrinello, Science, 1999, 286, 1051.
30. G. Mathias and M. D. Baer, J. Chem. Theory Comput., 2011, 7, 2028.
31. G. Mathias, S. D. Ivanov, A. Witt, M. D. Baer and D. Marx, J. Chem. Theory Comput., 2012, 8, 224.
32. A. Brown, B. J. Braams, K. Christoffel, Z. Jin and J. M. Bowman, J. Chem. Phys., 2003, 119, 8790.
33. A. Brown, A. B. McCoy, B. J. Braams, Z. Jin and J. M. Bowman, J. Chem. Phys., 2004, 121, 4105.
34. A. B. McCoy, B. J. Braams, A. Brown, X. Huang, Z. Jin and J. M. Bowman, J. Phys. Chem. A, 2004, 108, 4991.
35. X. Huang, L. M. Johnson, J. M. Bowman and A. B. McCoy, J. Am. Chem. Soc., 2006, 128, 3478.
36. L. M. Johnson and A. B. McCoy, J. Phys. Chem. A, 2006, 110, 8213.
37. X. Huang, A. B. McCoy, J. M. Bowman, L. M. Johnson, C. Savage, F. Dong and D. J. Nesbitt, Science, 2006, 311, 60.
38. X.-G. Wang and T. Carrington Jr, J. Chem. Phys., 2008, 129, 234102.
39. C. E. Hinkle, A. S. Petit and A. B. McCoy, J. Mol. Spectrosc., 2011, 268, 189.
40. E. T. White, J. Tang and T. Oka, Science, 1999, 284, 135.
41. O. Asvany, K. M. T. Yamada, S. Brünken, A. Potapov and S. Schlemmer, Science, 2015, 347, 1346.
42. R. Wodraszka and U. Manthe, J. Phys. Chem. Lett., 2015, 6, 4229.
43. H. Schmiedt, S. Schlemmer and P. Jensen, J. Chem. Phys., 2015, 143, 154302.
44. X.-G. Wang and T. Carrington Jr, J. Chem. Phys., 2016, 144, 204304.
45. H. Schmiedt, P. Jensen and S. Schlemmer, Phys. Rev. Lett., 2016, 117, 223002.
46. H. Schmiedt, P. Jensen and S. Schlemmer, Chem. Phys. Lett., 2017, 672, 34.
47. D. W. Boo and Y. T. Lee, Chem. Phys. Lett., 1993, 211, 358.
48. D. W. Boo, Z. F. Liu, A. G. Suits, J. S. Tse and Y. T. Lee, Science, 1995, 269, 57.
49. D. W. Boo and Y. T. Lee, J. Chem. Phys., 1995, 103, 520.
50. D. W. Boo and Y. T. Lee, Int. J. Mass Spectrom., 1996, 159, 209.
51. P. J. Kelleher, C. J. Johnson, J. A. Fournier, M. A. Johnson and A. B. McCoy, J. Phys. Chem. A, 2015, 119, 4170.
52. R. Navrtil, J. Jak and J. Roithov, J. Mol. Spectrosc., 2017, 332, 52–58.
53. E. Fois, A. Gamba and M. Simonetta, Can. J. Chem., 1985, 63, 1468.
54. S. J. Kim, P. R. Schreiner, P. v. R. Schleyer and H. F. Schaefer III, J. Phys. Chem., 1993, 97, 12232.
55. S. Roszak and J. Leszczynski, Chem. Phys. Lett., 2000, 323, 278.
56. A. Witt, S. D. Ivanov, H. Forbert and D. Marx, J. Phys. Chem. A, 2008, 112, 12510.
57. A. Witt, PhD thesis, Ruhr-Universität Bochum, 2010.
58. A. Witt, S. D. Ivanov and D. Marx, Phys. Rev. Lett., 2013, 110, 083003.
59. J. Hutter, et al., CPMD Program Package, http://www.cpmd.org.
60. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
61. R. Ramirez, T. Lopez-Ciudad, P. Kumar and D. Marx, J. Chem. Phys., 2004, 121, 3973.
62. T. J. Ulrych and T. N. Bishop, Rev. Geophys. Space Phys., 1975, 13, 183.
63. R. W. Herring, CRC Report No 1330, Communications Research Centre, Ottawa, Canada, 1980.
64. N. Wu, The Maximum Entropy Method, Springer-Verlag, Berlin Heidelberg, 1997.
65. J. P. Burg, Proc. 37th Meet. Soc. Exploration Geophysicists, 1967.
66. J. P. Burg, Paper presented at Advanced Study Institute on Signal Processing, NATO, Enschede, Netherlands, 1968.
67. J. P. Burg, PhD thesis, Stanford University, 1975.
68. P. Broersen, IFAC Proc., 1997, 30, 245.
69. M. Morf, A. Vieira, D. T. L. Lee and T. Kailath, IEEE Trans. Geosci. Electron., 1978, GE-16, 85.
70. P. J. Brockwell, R. Dahlhaus and A. A. Trindade, Statistica Sinica, 2005, 15, 197.
71. J. Sun, G. Niehues, H. Forbert, D. Decka, G. Schwaab, D. Marx and M. Havenith, J. Am. Chem. Soc., 2014, 136, 5031.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sc04040g

---