Instanton calculations of tunneling splittings in chiral molecules

We report the ground state tunneling splittings (ΔE±) of a number of axially chiral molecules using the ring‐polymer instanton (RPI) method (J. Chem. Phys., 2011, 134, 054109). The list includes isotopomers of hydrogen dichalcogenides H2X2 (X = O, S, Se, Te, and Po), hydrogen thioperoxide HSOH and dichlorodisulfane S2Cl2. Ab initio electronic‐structure calculations have been performed on the level of second‐order Møller–Plesset perturbation (MP2) theory either with split‐valance basis sets or augmented correlation‐consistent basis sets on H, O, S, and Cl atoms. Energy‐consistent pseudopotential and corresponding triple zeta basis sets of the Stuttgart group are used on Se, Te, and Po atoms. The results are further improved using single point calculations performed at the coupled cluster level with iterative singles and doubles and perturbative triples amplitudes. When available for comparison, our computed values of ΔE± are found to lie within the same order of magnitude as values reported in the literature, although RPI also provides predictions for H2Po2 and S2Cl2, which have not previously been directly calculated. Since RPI is a single‐shot method which does not require detailed prior knowledge of the optimal tunneling path, it offers an effective way for estimating the tunneling dynamics of more complex chiral molecules, and especially those with small tunneling splittings.


| INTRODUCTION
Quantum tunneling in degenerate rearrangements of molecular clusters has been a topic of high interest in many fields of physics, chemistry, and biology. [1][2][3] In recent years, tunneling phenomena have attracted particular attention in molecular systems where the splittings give insights into the nature of rearrangement dynamics, 4,5 even at low temperature at which rearrangements are otherwise highly suppressed. 1,6 In a symmetric double-well system, tunneling splittings are produced by the overlap of the tails of degenerate wavefunctions essentially localized in the wells and correspond to the energy difference between delocalized eigenstates that can be formed as superpositions of the localized states. Much work has been devoted to determining such splittings for a wide variety of systems comprising of small molecules and molecular clusters, in particular water clusters. 5,[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] The magnitude of the splittings in the molecular systems varies over many orders of magnitude and is sensitive to the nonequilibrium regions of the potential energy surface (PES). 3 Numerical determination of splittings thus can be used to investigate the interactions and dynamics of the systems and to interpret the experimental spectra. 5,11,26 There are a variety of theoretical methods available for computing tunneling splittings of molecules. One approach is the numerical determination of tunneling splittings by a direct solution of the Schrödinger equation for multi-dimensional or reduced-dimensional rovibrational tunneling dynamics. [27][28][29][30][31][32] Some of the notable studies include the calculation of lowest rovibrational states and tunneling splittings in the intermolecular bound states of (H 2 O) 2 , 14,27,33 calculations on the tunneling splittings of the vibrational ground state and low-lying excited states of malonaldehyde, 30 the construction of a global analytical PES for the electronic ground state of NH 3 from multireference configuration interaction and coupled cluster study 31 and calculation of the tunneling splittings of NH + 3 , OH + 3 , and so on. 32 Being computationally expensive due to the formal exponential scaling with the number of degrees of freedom, this is in full dimensionality only applicable to small gas-phase systems. Methods such as diffusion Monte Carlo (DMC) 20,34,35 and path-integral techniques [36][37][38][39][40][41][42] are useful for larger systems. Systems treated by these techniques include small water clusters to helium crystals, ferrous-ferric exchange and electron tunneling paths in proteins. 20,[34][35][36][38][39][40][41][42][43][44] However, in order to reduce statistical errors, these numerical methods probe the PES over a large number of independent geometries that are accessible to the nuclear wave function and hence become computationally expensive. Further, DMC is practically impossible to converge for systems with very small splittings, as in the water hexamer prism. 5 To break this bottleneck, approximations based on the semiempirical framework within Wentzel-Kramers-Brillouin (WKB) limit are reported in the literature, which although not as accurate as the Monte Carlo method, often give qualitative agreement with the experimentally observed splittings. 37,[45][46][47][48] In general, the WKB approach has the advantage of yielding the splittings directly from fluctuations around a single tunneling path. However, the paths must be assigned a priori. This requires one to develop an a priori model of the tunneling, which may be possible to do accurately for water dimer and certain other small systems [15][16][17]45,49 but is more difficult for larger clusters with complex dynamics. 4,5,23 Thus, each of these methods have their advantages and disadvantages.
Quack and coworkers presented quantitative calculations of mode-selective stereomutation tunneling of a wealth of small chiral molecules. 50 The list includes isotopomers of chiral hydrogen dichalcogenides of the type HXXH (X = O, S, Se, Te), hydrogen thioperoxide (HSOH), dichlorodisulfane (ClSSCl), and HClOH cations. [51][52][53][54][55][56] The stereomutation process in these compounds is dominated by the torsional motion. The magnitude of the tunneling splittings typically increases for higher torsional levels as compared to lower torsional states, i.e. excited state tunneling splittings are dominant over the ground state counterparts. 57,58 These torsional tunneling dynamics were investigated with a quasi-adiabatic channel quasi-harmonic reaction path Hamiltonian (RPH) approach, 55 61 This method is also based on a tunneling path which must be chosen a priori, typically along the torsional coordinate. In general, the torsional motions are treated anharmonically, whereas all remaining coordinates are only treated harmonically, although these are anharmonically coupled with reaction coordinates. The RPH approach spans all six vibrational degrees of freedom, where the influence of the nontorsional modes is taken into account approximately and allows to study mode-selective catalysis and inhibition of stereomutation by the various vibrational modes in these chiral molecules.
Among all these molecules, HOOH, HSSH, and HSOH are widely investigated both experimentally and theoretically, 29,51,52,55, whereas the studies on other hydrogen dichalcogenides are rather rare. 53,54 The latter is likely due to the fact that tunneling splittings are very small for the heavier hydrogen dichalcogenides and thus challenging to observe experimentally. Nevertheless, the magnitude of the tunneling splittings is of fundamental importance for the dynamics of chiral molecules and for detection attempts of parity violating interactions in chiral molecules. We refer the interested reader to Refs. 88-90 for a more detailed discussion of the interplay between tunneling and parity violation.
The instanton method is derived from a rigorous semiclassical approximation of the path-integral formalism for the quantum partition function. [91][92][93] Like the WKB method, it predicts the tunneling splittings from a single tunneling path, but unlike WKB, it does not require prior detailed knowledge of this path and instead locates the optimal imaginary-time trajectory connecting the wells. Initially, the instanton approach was introduced for rate theories, 94,95 but later also applied for predicting ground-state tunneling splittings, 95,96 where the tunneling path was identified as the classical trajectory that connects the two identical minima on the upside-down potential.
Milinikov and Nakamura 97 were the first to develop the instanton approach into a practical computational tool for calculating tunneling splittings. Their implementation was set up in internal molecular coordinates, which is favorable for reducing coupling terms in the potential energy but it may be challenging to apply the method to systems with more complicated structural arrangement. The ring-polymer instanton (RPI) method uses Cartesian coordinates instead, 93,98 which can be advantageous for the study of molecular clusters. It has been successfully applied to compute the tunneling splitting in a number of molecular systems, including ammonia, malonaldehyde, the formic acid dimer as well as the radicals hydroperoxyl and vinyl. 97,[99][100][101][102][103] In addition, the RPI has also been generalized to yield the splittings pattern in systems with multiple symmetric wells 23 such that it can be applied to study water clusters. 5,13,23,24 Similarly to the ring-polymer approach, 104 the tunneling path is discretized into beads that represent the molecular system at equally spaced imaginary-time steps and the instanton pathway is obtained numerically by minimization of the Euclidean action. The tunneling splitting is obtained in terms of a steepest-descent approximation around it. This is the origin of the semiclassical approximation inherent in the method. Further, in using Cartesian coordinates for the beads along the imaginary-time trajectory, one must also assume that tunneling splittings depend only weakly on the overall rotation of the system, 93 which commonly is expected to hold for clusters or large molecules.
Limitations of the instanton method thus include a harmonic treatment of the fluctuations around the instanton path and the neglect of rovibrational coupling terms. Anharmonicity along the instanton path is, however, accounted for. Hence, the instanton approach is expected to capture the dominant tunneling mechanism and to provide an estimate of the tunneling splittings of the right order of magnitude unless anharmonic effects perpendicular to the tunneling path are sufficiently large to produce qualitative changes in the tunneling dynamics. For the case of a low barrier height, some fluctuations can follow paths that access the top of the barrier and even become trapped in the other well, which are not well described by this approach. Therefore, the accuracy of the computed tunneling splittings increases when the symmetric wells are further separated and when the barrier between them rises. In contrast to DMC, the RPI method is thus best suited to computing small tunneling splittings.
In practice, for a given underlying PES, instanton tunneling splittings calculated for molecular systems are usually at least within a factor of 2 of exact quantum results and often better. 99 As well as the quantitative estimate, the instanton approach identifies the optimal tunneling path, which gives mechanistic insight into tunneling in complex systems. In the present paper, we have applied the RPI method for calculating the tunneling splittings of a wealth of chiral molecules including isotopomers of hydrogen dichalcogenides H 2 X 2 with X = O, S, Se, Te, Po and HSOH and S 2 Cl 2 which are then compared with the values reported in the literature where available. By this, we provide RPI reference data for tunneling splittings of axially chiral molecules.
Tunneling splittings (ΔE ± ) on H 2 Po 2 are studied herein to the best of our knowledge for the first time, and although tunneling in S 2 Cl 2 was estimated previously, 56 we provide the first direct calculation. The following section describes the mathematical formulation for ΔE ± within the RPI approach whereas the results from the latter are summarized in Section 3.

| Instanton theory
The RPI method has been previously described, 23,91-93,100 however, for uninitiated readers, here we briefly summarize the working principle and underlying theory of the method. The tunneling splitting ΔE ± depends on the energy difference between closely lying energy levels of opposite parity and is thus a temperature independent property.
Nonetheless, in order to connect to the formalism of path integrals, it is practical to introduce the partition function of the system in the low temperature limit. This limit ensures that the partition function is determined purely in terms of the lowest energy levels corresponding to the ground-state tunneling splitting. If E 0 is the vibrational zeropoint energy of an isolated well, the energy levels of the tunneling system is assumed to be symmetrically split into E 0 − ΔE ± /2 and E 0 + ΔE ± /2 (cf., Figure 1). If Z(β) and Z 0 (β) represent the partition functions of the system with and without tunneling, respectively, then the tunneling splittings ΔE ± can be defined by. 1 where β = 1/k B T is proportional to the inverse thermodynamic temperature, T, and k B is the Boltzmann constant. Thus, by evaluating the ratio of partition functions in the low-temperature limit within a semiclassical path-integral formalism, it will be possible to obtain an expression for the tunneling splitting where the partition function is represented in terms of ring polymers and the required integrals are computed using the method of steepest descent.
Application of steepest descent to the full system results in an infinite number of minima on the ring-polymer potential surface. 93 These ring polymers describe imaginary-time trajectories along which the particle tunnels back and forth between the minima of the degenerate wells. The trajectory can make n such passes, known as kinks, and each kink contributes a factor θ(β) to the partition function of the system. Including factors of 1/n! to avoid overcounting and noting that only even numbers contribute, the series can be summed to give cosh(θ(β)), from which one can identify an expression for ΔE ± in terms of a single kink: The problem is thus reduced to the determination of the contribution of an individual kink.
In a practical computation, the representation for a single isolated kink is obtained by introducing a linear polymer of N beads. The beads of the linear polymer are connected with harmonic springs to give the effective potential where β N = β/N and Q i,j = ffiffiffiffi m p j × x j represents the mass-weighted coordinates of the i th bead with m j being the mass of the corresponding degree of freedom defined by the Cartesian coordinate x j . The kink is F I G U R E 1 Potential energy curves for the full system (solid curves) of a double-well potential system with partition function Z(β) and tunneling free system (dotted curves) with partition function Z 0 (β). Schematic presentation for the tunneling splittings up to fourth vibrational state. The present article deals with ground state splittings defined as a local minimum of U N (Q) in which the beads travel from one well to the other. For practical purposes, a suitably large β is considered in order to yield adequate numerical convergence of ΔE ± and likewise a large but finite value of N is taken to ensure that β N is small.
Finally, the weight factor θ(β) for the contribution made by a single kink is obtained by where S = β N ℏU N (Q) is the classical action evaluated along the kink, and the determinant ratio is 23,91-93,100 η k represent the frequencies of the normal modes obtained from diagonalising the Hessian of the (fixed-ended) linear-polymer PES. 93 Note that, the zero frequency corresponding to the permutational mode is not included in the product. η polymer. 23,93 As stated in the previous section, the instanton includes anharmonicity only along the tunneling path, whereas the fluctuations are treated harmonically. Furthermore, the instanton approach for the tunneling splittings uses Cartesian coordinates for the beads, since this is the preferred way of applying the instanton approach to molecular clusters, although this inherently neglects the weak dependence of ΔE ± on the overall rotation of the system. 93 Thus, even after the numerical convergence is achieved, the results for ΔE ± still include the errors from the instanton approximations. Typically, the RPI method can be expected to give an estimate of ΔE ± that is within an order of magnitude or better for a given PES. 23,93,97,99,100,102,105 As per the working principle, generally a single kink minimum is located on the linear-polymer PES using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm, 106,107 starting from an initial arrangement of beads that are evenly spaced between the minima. For higher number of beads, a substantial reduction of the number of steps in the minimization was found by using as a starting configuration, the arrangement of minimal energy previously found from a lesser number of beads. Here, we use the instanton approach 23

| Computational details
The instanton optimization 108   We compute the quantity W along the mass-weighted path Q according to both the CCSD(T) and MP2 levels of theory. 116,117 where i indexes the beads used. Note that, E is the potential energy of the minimum according to the theory used to calculate V(Q).
3 | RESULTS AND DISCUSSION H 2 X 2 with X = O, S, Se, Te, Po and HSOH and S 2 Cl 2 serve as prototype molecules for the stereorotation dynamics of an axially chiral molecule (cf., Figure 2). Table 1 depicts the level of theory used, the cis-and trans-barrier heights for the interconversion between the enantiomers of these molecules. In all cases, the cis-barrier height is larger than the trans-barrier height, such that the tunneling splitting is dominated by the trans mechanism, as our calculations shall confirm.
Except for HSOH isotopomers, all other molecules reveal a C 2 symmetrical equilibrium structure. They represent two (P and M) enantiomers that cannot be superimposed through rotation. From  The values of N range from 32 to 2048 in a two-fold manner. From Table 2, it can be seen that for every value of β, the action S converges F I G U R E 2 Chiral C 2 -symmetric equilibrium structure of H 2 X 2 molecule very fast as we increase the number of beads whereas the convergence is relatively slow for ϕ. This does not cause a significant problem as only S appears in the exponent whereas ϕ appears in the prefactor and thus does not require a high accuracy. Note that, suggestions have been given for speeding-up the convergence of ϕ 118 and for the optimization itself. 100 Nonetheless, here numerical convergence for ϕ to within less than 1% is achieved with 512 beads in each case.  Table 3 gives the calculated ΔE ± values for H 2 X 2 (X = O, S) molecules.

| Convergence of the instanton calculations
It is seen that slightly different results are obtained from the various values of β, which shows that full convergence has not been reached.
We take the largest value of β as our best result and use the others as a guide to estimate the convergence error, which is in this case less than 10%. Bearing in mind that the instanton approximation itself typically introduces an error at least of this magnitude, 99 and that the approximate electronic-structure methodology will also introduce an error into the prediction, this convergence criterion is adequate.
Thus, based on the observation from Tables 2 and 3, we conclude that instanton optimization with 512 beads is sufficient to achieve convergence for the action S and the determinant ratio ϕ and hence the calculated ΔE ± values for these chiral molecules. In      Table 4), which shows that they are not making significant errors either. This is of particular importance as this system is probably the most difficult for the instanton method due to its low barrier. One would thus expect even smaller errors for the other systems, although the error of the electronic-structure method may start to dominate then.

| Tunneling splittings ΔE ± of isotopomers
Likewise, in case of HSOH and H 2 S 2 isotopomers, the calculated ΔE ± values from the RPI method match to within the same order of magnitude to the earlier reported values. 52,55,72,73,[78][79][80][81][82][83] In particular, the values from RPI and the quasi-adiabatic channel-RPH method 55,59,60 of Quack and coworkers are comparable to the same order of magnitude, although, in most cases, the values from the latter are slightly smaller. However, in case of HSOH isotopomers, RPH method estimates ΔE ± values that are roughly the same as the RPI ones. The differences of the values from RPI and RPH can be mostly attributed to the level of theory and basis sets used for the electronicstructure calculations, which produce different barrier heights. Due to the lack of experimental data on the tunneling splittings in H 2 X 2 [X = Se, Te] species, we here only consider the RPH values as reference ones. As expected, RPI and RPH again give similar results. To the best of our knowledge, no calculation of tunneling splittings ΔE ± for H 2 Po 2 isotopomers have been reported before, so that we present the RPI data as benchmark for future theoretical and experimental work.
From Table 4, it can be observed that H 2 X 2 shows smaller and smaller tunneling splittings ΔE ± in the series of heavier analogues of hydrogen peroxide, that is, for X = O via X = S and X = Se to X = Te.
The barrier height for X = O is particularly small, which clearly explains why H 2 O 2 has the largest tunneling splitting. However, the trend for the other molecules in the series is that the tunneling splittings typically decreases, despite the fact that the barrier height gets lower. This result is due to the fact that the effective barriers get progressively wider due to the lengthening of the H X bond, which thus sweeps out a longer arc during the torsional tunneling process. Plots of the potentials along the mass-weighted instanton pathway provide a measure of the width of the effective barrier and are shown in Figure 4. It turns out that in most of these cases, the increasing width of the barrier has a greater effect than the decreasing height. The two opposing effects almost cancel out, however, when going from H 2 Te 2 to H 2 Po 2 , giving H 2 Po 2 an unexpectedly large splittings for its position in the series. This observation can be attributed to the relatively low barrier height of H 2 Po 2 molecule with respect to H 2 Te 2 , whereas the barrier width only increases slightly (see Figure 4).
Unlike the DMC method, RPI is well-known for its better suitability to systems with very small splittings and in fact becomes even more accurate in this limit. Thus, the favorable comparison with F I G U R E 4 Barrier shapes as obtained from instanton calculations of the trans-pathways for the molecules in the H 2 X 2 series for N = 512 and β = 8000.0 E −1 h . The potential maximum along the path was centered at Q = 0 T A B L E 4 Trans-tunneling splittings ΔE ± (in cm −1 hc) of the isotopomers calculated on the MP2 level with 512 beads and imaginary-time-duration βℏ = 8000. In each case, the splittings decrease further after isotopic substitution of hydrogen by deuterium (D) or tritium (T), since they are associated with a higher effective reduced mass. Except for H 2 O 2 , which has a particularly low barrier, the tunneling splittings ΔE ± are found to decrease by $3 or 4 orders of magnitude upon D substitution, which are even further reduced by another $2 or 3 orders of magnitude after the T substitution (Table 4 ). Te, and Po) and HSOH molecules are utilized for estimating the same at CCSD(T) level of theory. This can be easily done by computing F I G U R E 5 Variation of bond lengths and bond angles with respect to the equilibrium structure for H 2 X 2 series along the transtrajectory for N = 512 and β = 8000.0 E − 1 h . The potential maximum along the path was centered at Q = 0 W (which essentially denotes the action S of the instanton path) as given in Equation (6). The advantage of using W is that it is simply a one-dimensional integral along the mass-weighted path length Q. 117 The computed value of W from MP2 and CCSD(T) level of theory are utilized for calculating new splittings ΔE ±,est.CCSD(T) defined by Equation (7). As already discussed in the computational section, 16 evenly-spaced beads in Q are chosen from the optimized instanton path obtained with 512 beads and β values of 8000.0 E −1 h . This is done only for the hydrogen isotopomers of the H 2 X 2 series and for HSOH (see Table 6). The difference between MP2 (Table 4) and CCSD(T) splittings remain on the level of 30% or less. Thus, none of the predictions changes by a significant amount. Similar values are also expected for heavier isotopomers.
The tunneling dynamics in S 2 Cl 2 deserves special discussion. Due to the existence of high-energy barrier for the trans-interconversion between the enantiomers and higher effective reduced mass, the tunneling time is larger, which ultimately results in an extremely low value for the tunneling splittings ΔE ± . The earlier 56 RPH based calculations were concerned with the determination of the very small hypothetical tunneling splittings in this molecule, for which the direct calculation of the extremely small torsional tunneling splittings was not possible due to numerical limitations. In view of this, an extrapolation procedure was used for determining the tunneling splittings in S 2 Cl 2 , wherein the ground state tunneling splittings of S 2 Cl 2 was estimated to be ≈10 −76 cm −1 hc. 56 RPI, however, is able to calculate the splittings directly from a single kink, which is no more difficult than for the other molecules. As already discussed, the optimized instanton pathway involves the motion of all four atoms in this molecule. We choose β values of 20,000.0 E −1 h , 30,000.0 E − 1 h , and 40,000.0 E −1 h for the MP2/6-311+G(3df) level instanton optimization, which is extended from 32 beads to 512 beads in a two-fold manner. The choice of large β values is obvious by virtue of the high-energy barrier and large effective reduced mass for the tunneling. Note that, this calculation does not necessarily require more beads, as the time-step βℏ/N can also be longer for these slow-moving heavy atoms. Table 7 summarizes the calculated ΔE ± values whereas the corresponding S and ϕ values are provided in the supplementary material. Instanton theory predicts ΔE ± values of the order of 10 −96 cm −1 hc, which are even smaller than those estimated from the extrapolation procedure. 56 This confirms the previous findings that parity violation dominates the dynamics rather than tunneling.

| CONCLUSIONS
We have calculated the tunneling splittings ΔE ± for several axially chiral molecules: H 2 X 2 with X = O, S, Se, Te, Po, and HSOH and S 2 Cl 2 along with their isotopomers employing the RPI approach. The instanton approach is a single-shot method which does not require detailed a priori knowledge of the tunneling path but instead automatically locates the optimal tunneling pathway on the ab initio potentialenergy surface. In this work we have evaluated the potential-energy surface "on the fly," although we note that recent developments may enable future calculations to be carried out at a reduced cost by employing machine-learning techniques. 120 The fluctuations perpendicular to the instanton are treated harmonically whereas anharmonicity along the tunneling pathway is accounted for. Also, RPI neglects the weak dependence of ΔE ± on the overall rotation of the system. One might expect that this degrades the accuracy of the instanton tunneling splittings to some extent. However, in practice, RPI has been seen to produce ΔE ± values which are usually within an order of magnitude or better for a given PES, 99 and there are no indications of a breakdown for the systems treated here.
The calculated trans-tunneling splittings ΔE ± are then compared to the reported counterparts from the literature. The cis-pathway is found to be less important in each case. For H 2 X 2 with X = O, S, Se, T A B L E 5 Cis-tunneling splittings ΔE ± (in cm −1 hc) of the isotopomers calculated on the MP2 level with 512 beads and imaginary-time-duration βℏ = 8000.0 E −1 h ℏ Te, and HSOH molecules and corresponding isotopomers, the instanton optimization is found to produce ΔE ± values to within the same order of the reported values. This implies that the reaction path chosen for the RPH approach was a good enough approximation to the instanton pathway, which was expected for these small systems in which the tunneling pathways are fairly simple. The only case with significantly larger deviation is the instanton optimization for S 2 Cl 2 , which gives ΔE ± values of the order of 10 −96 cm −1 hc, which is approximately 20 orders of magnitude less than the previous prediction based on an extrapolation procedure. In this case, the instanton method has an advantage over alternative methods because it can be directly applied to calculate very small splittings. The RPI results confirm the previous finding that the tunneling splittings in S 2 Cl 2 are much smaller than the parity violating energy differences, which is of particular interest. The calculated ΔE ± values for H 2 Po 2 isotopomers from this work may be used as reference values for other methods that are utilized for deriving tunneling splittings ΔE ± in chiral molecules and also even for conducting experiments in future.