Instanton theory of tunnelling in molecules with asymmetric isotopic substitutions

We consider quantum tunnelling in asymmetric double-well systems for which the local minima in the two wells have the same energy, but the frequencies differ slightly. We derive a generalization of instanton theory for these asymmetric systems, leading to a semiclassical expression for the tunnelling matrix element and hence the energy level splitting. We benchmark the method using a set of one- and two-dimensional models, for which the results compare favourably with numerically exact quantum calculations. Using the ring-polymer instanton approach, we apply the method to compute the level splittings in various isotopic-substituted versions of malonaldehyde in full dimensionality and analyse the relative contributions from the zero-point energy difference and tunnelling effects.


I. INTRODUCTION
The phenomenon of quantum tunnelling in a symmetric double well is a well-known concept. 1,2 The earliest molecular example studied was the splitting of the ground-state energy levels of the ammonia molecule due to inversion tunnelling. 3 However, the prototypical system of intramolecular hydrogen-atom tunnelling is malonaldehyde, 4,5 whose ground-state tunnelling splitting has been measured to high accuracy 6 and studied with a wide variety of theoretical methods. [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] Instanton theory, [25][26][27][28] which is based on a semiclassical approximation to the path-integral formulation of quantum mechanics, 29,30 provides a particularly powerful approach for computing tunnelling splittings, especially when reformulated into the ring-polymer instanton method 20,31,32 or related approaches. [22][23][24]33 These methods have been applied to calculate the ground-state tunnelling splitting in numerous molecular systems in full dimensionality. [19][20][21][22][23][24][34][35][36][37][38][39][40][41][42] In contrast, there has been much less theoretical work on tunnelling between asymmetric wells. 43 Some exceptions include studies of tunnelling between similar (but not equivalent) minima in the water hexamer cage. 44,45 However, in this paper, we shall study a subtly different problem, one in which both well minima have exactly the same potential energy, but different frequencies. This can occur (within the Born-Oppenheimer approximation) in molecules or clusters with certain asymmetric isotopic substitutions. Previous studies of tunnelling switching in meta-D-phenol 46 and the energy-level structure in the dimer HF -DF 47 provide examples which exhibit this behaviour. In both these cases, however, the effect of tunnelling between the ground vibrational states of each well was predicted to be much smaller than the zero-point energy (ZPE) difference between the two wells.
In this work, we study cases where the asymmetry is weak enough such that tunnelling may still have a significant effect on the ground-state energy levels of the a) These authors contributed equally b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch molecule. We develop a theory which can be applied using ring-polymer instanton methodology 20,31,32 to compute the effect of tunnelling on the level splitting in such molecular systems. This theory is based on minimumaction tunnelling pathways, called instantons, which connect the bottom of the wells. This is in contrast to alternative approaches which go by various names such as WKB theory, 1,48-50 semiclassical path integral, 51 Bohr-Sommerfield quantization 52,53 or periodic-orbit theory, 54 for which the tunnelling pathways travel between two turning points defined by the energy of a localized vibrational state. These approaches may give good results for model systems, but unlike instanton theory, they cannot be easily extended to multidimensional problems (without a priori choices of path such as that of Ref. 55).
We will apply our new instanton approach to benchmark one-and two-dimensional model systems in various parameter regimes and thereby test its accuracy. Then, we apply it in full dimensionality to compute the level splitting in a set of isotopomers of malonaldehyde, including both symmetric and asymmetric substitutions and compare with experimental results where available. 4,5

II. THEORY
We study tunnelling in a system defined by the Hamil-tonianĤ =p 2 /2m + V (x). This expression is valid for both one-dimensional and multidimensional systems, as in the case of the latter we treat position, x, and momentum, p, as vectors which have been mass-weighted such that all degrees of freedom have the same effective mass, m. We will consider potentials, V (x), which exhibit a special type of asymmetric double-well structure, as illustrated in Fig. 1. In particular, there are two minima, at x and x r , which have exactly the same potential energy and which we choose to set to zero, i.e. V (x ) = V (x r ) = 0. However, the frequencies, ω /r = ∇ 2 V (x /r )/m, in the two wells may differ leading to different ZPE values, E 1 2 ω and E r 1 2 ω r (within the harmonic approximation). 56 For such a system, we can define an effective two-level Hamiltonian using a basis of states localized in the left arXiv:2007.07540v1 [physics.chem-ph] 15 Jul 2020 FIG. 1: Schematic representation of the potential, V (x), for an asymmetric system (solid curve) with non-tunnelling zero-point energy levels shown with horizontal dashed lines. For comparison the potential of a symmetric system is shown with the grey dotted curve with its non-tunnelling zero-point energy levels E 0 (horizontal grey dotted lines). The solid horizontal lines indicate the lowest two levels of the asymmetric system including tunnelling. and right wells, where d = 1 2 (E r − E ) and E 0 = 1 2 (E + E r ). Here, − Ω is the coupling due to tunnelling between the left and right wells, which will be evaluated by the instanton theory developed in this paper. By introducing this effective Hamiltonian, we have implicitly assumed that the asymmetry is weak enough such that the ground vibrational state of the left well couples to the ground vibrational state of right well and we have neglected couplings to other levels. 57 This ensures that the tunnelling will remain coherent, which is the regime of interest as an instanton theory for incoherent tunnelling (i.e. rate constants) is of course already well known. 25,26,28,31,32,[58][59][60][61] Simple quantum-mechanical principles based on the effective Hamiltonian of Eq. (1) allow for a reduced description of the dynamics and spectrum of the full system. The solutions are well known as this is an equivalent problem for all two-level quantum systems, 62 e.g. molecular orbitals defined in terms of a linear combination of atomic orbitals in a heteronuclear diatomic or a spin-1 2 in a magnetic field. The effective Hamiltonian has eigenvalues E ± = E 0 ± d 2 + ( Ω) 2 , and thus for an energyresolved spectrum the observable of interest is the level splitting, ∆ = 2 d 2 + ( Ω) 2 , which has contributions both from the tunnelling effect as well as the ZPE difference. The ground-and excited-state eigenvectors of the Hamiltonian matrix can be written as ψ − = (cos φ 2 , sin φ 2 ) and ψ + = (− sin φ 2 , cos φ 2 ), where φ = tan −1 ( Ω/d) is the mixing angle (0 ≤ φ ≤ 180 • ). This angle is a measure of localization of the stationary states and is a function only of the ratio d/ Ω. That is, for mixing angle φ = 90 • , the eigenstates are maximally delocalized (with a population of 1 2 in each well), but if the angle approaches 0 or 180 • , the ground-and excited-state wavefunctions are each localized to one well. One can also obtain the time-dependence of the populations from the effective Hamiltonian. For example if the system starts in the left well, the time-dependent population in the right well is sin 2 φ sin 2 (∆t/2 ) and thus the mixing angle, φ, also determines the amplitude of population transfer. This analysis shows that one requires only the values of d and Ω to make predictions for all these observables of interest. It is easy to obtain a good estimate of d via the ZPEs in each well evaluated using the harmonic approximation as described above. However, the tunnelling contribution, Ω, cannot be obtained from a harmonic approximation 63 and requires a more involved treatment, to which we dedicate the rest of this section.
As has been done for the symmetric case, 20,28,32 we can start our derivation by considering the partition function of the effective Hamiltonian [Eq. (1)], which is Z = 2 e −βE0 cosh(β d 2 + ( Ω) 2 ), where β is the reciprocal temperature. However, we must generalize the standard approach in order to have an appropriate reference partition function for the non-tunnelling case. Here, we define Z g = 2 √ Z Z r in terms of the geometric average of Z = e −β(E0−d) and Z r = e −β(E0+d) , the partition functions for the left/right well respectively (i.e. in the absence of tunnelling, where Ω = 0), such that Z g = 2 e −βE0 . The ratio of the total partition function to this reference is which for later convenience can be equivalently written as a Taylor series in Ω: The effective Hamiltonian only defines the two lowest states of the system. Thus, writing the partition function in terms of this two-level Hamiltonian employs an approximation which is only valid if the temperature is chosen such that the thermal energy (β −1 ) is much smaller than the spacing between the vibrational states of a single well ( ω /r ). In Sec. II A, we will therefore take the low-temperature, β → ∞, limit of the full partition function in order that it corresponds with the expressions based on H eff .

A. Instanton Theory
To derive an method for calculating Ω and hence ∆ within instanton theory, we require expressions for Z and Z g based on the full Hamiltonian,Ĥ, which we can compare with the previous analysis based on H eff . The partition function Z = Tr[e −βĤ ] can be defined within the path-integral formalism as a sum over all closed paths in the space of x. 29 Some of these paths are located completely in one of the wells, whereas others pass back and forth between the wells n times. Each pass is called a "kink" and in order for the path to be closed in this simple double-well case, n must be an even number. We define Z n as the contribution to the partition function from paths with n kinks such that 20,28,32 We will evaluate semiclassical approximations to Z n for each value of n, which are rigorously defined by taking the leading term of the asymptotic expansion 64 in and for which we use the notation A B to mean A/B → 1 in the limit → 0 and, where appropriate, also low temperature or long imaginary time.
Consider first the Z 0 term, which includes contributions from all non-tunnelling paths completely located either in the left or right wells. We define K 0 (x a , x b , τ ) as the zero-kink contribution to the semiclassical approximation of x a |e −τĤ/ |x b . 32 Then we define Z by taking the integral only over the domain of the left well, Γ , 65 and equivalently for Z r , where K (τ ) = K 0 (x , x , τ ) and K r (τ ) = K 0 (x r , x r , τ ) are the kernels for a path starting and ending at the bottom of one of the wells without crossing the barrier, whereas Ξ and Ξ r are terms resulting from the steepest-descent integration over the end points. Explicit formulas are given in Appendix A which can be used to show that according to the semiclassical approximation, these partition functions are as could have been anticipated from a simple harmonic analysis in the low-temperature limit. Using Z 0 = Z +Z r and Z g = 2 √ Z Z r , the first term in Eq. (4) gives Comparing this with the first term in Eq. (3), we can identify d = 1 4 (ω r − ω ) as expected. In this way our path-integral analysis of the nontunnelling pathways has presented us with a formula for calculating d. There are of course simpler methods which give the same result, but the advantage of using the pathintegral formulation will come when considering the remaining (n > 0) terms in the series.
Generalizing the approach of Ref. 32, we can split the paths which contribute to Z n into n kinks. To define the contribution from a single kink, we first consider the semiclassical limit of the propagator x |e −τ kĤ / |x r . The paths which dominate this matrix element are single kinks of length τ k which start in the left well and end in the right. We will focus on a particular dominant path which is centred such that it passes through a point of maximum potential at imaginary time τ k /2. We call this dominant path an "instanton"; it is a minimum-action pathway and is thus equivalent to a classical trajectory in imaginary time. 26,66 Including also the non-zero fluctuations around this path leads to an instanton contribution which we call K 1 (x , x r , τ k ), where the prime indicates the restriction that the path be centred as described above. This quantity is defined explicitly in Eq. (90) of Ref. 32, and can be evaluated within the ring-polymer instanton formalism using Eq. (15).
We are now in a position to evaluate contributions to the two-kink term. Here, it is again necessary to extend the standard theory so as to keep track of the amount of imaginary time spent by the path trapped in each of the wells, which are no longer equivalent. The path is broken up into four pieces with two kinks of length τ k centred around imaginary times τ 1 and τ 2 as depicted in Fig. 2. Therefore, where we have defined K 1 (τ k ) = K 1 (x , x r , τ k ) = K 1 (x r , x , τ k ), which are equal due to time-reversal symmetry. Finally, using the formulas in Appendix A, we obtain Schematic of a 2-kink path starting and finishing in the left well in total imaginary time β . The kinks are located within a relatively short width τ k and are centred at imaginary times τ 1 and τ 2 .
where we have defined The formulas from Appendix A are only valid in the longtime limit and thus we must choose τ k large enough that this expression converges. However, the value of β is assumed to be much larger again such that the kinks are widely spaced and interactions between the kinks can be neglected.
In order to obtain the full expression for Z 2 , one needs to integrate over the imaginary times at which the kinks occur, which is simple because θ(τ k ) does not depend on these variables: 67 Therefore the required ratio is with τ k → ∞ implied. In order that this matches the second term of Eq. (3) we must again identify d = 1 4 (ω r − ω ) as before. However, we also gain an extra piece of information which implies that Ω lim and thereby defines the instanton approximation for the tunnelling matrix element. With these definitions, the agreement between the two series continues if one computes Z 4 /Z g etc., as of course it should due to the fact that the semiclassical approximation recovers the asymptotic ( → 0) limit. This can also be understood by considering the perturbation-theory expansion for Z = Tr[e −βH eff ], which generates Eq. (3) directly and has terms of a similar structure to those in Eq. (9) and its n ≥ 2 generalizations. 68 Alternatively, one can see that the correct expansion must be recovered because of the well-known isomorphism between a two-level system and an infinite periodic one-dimensional Ising model 69 as well as the less-well-known isomorphism between the latter and the non-interacting kinks. 70 Finally, we can check that our assumption of noninteracting kinks remains valid for the total partition function. In the standard instanton theory for symmetric systems, 26 one can predict that the average number of kinks which will appear is β Ω, which is typically low enough such that the intervals between them (Ω −1 ) is much larger than the width of the kinks (≈ 2/ω 0 , where ω 0 is the frequency of the wells). 30 It can be shown by similar arguments that the average number of kinks in an asymmetric system is less than or equal to that in the symmetric case such that the proof holds also for our case. 71 In the case of a symmetric system, , and thus one recovers the standard instanton result for the tunnelling splitting. 20,26,28,32 It therefore turns out that, although it was necessary to derive the more general asymmetric theory in a subtly different way, the only significant change to the final expression for the tunnelling matrix element is simply that the denominator should be replaced by the geometric mean of the two non-tunnelling kernels.

B. Ring-Polymer Instanton Theory
The instanton tunnelling pathways are uniquely defined as minimum-action paths connecting the two wells 26 and are equivalent to Newtonian trajectories travelling on the upside-down potential-energy surface. 66 However, in general multidimensional systems, these pathways are not known analytically. We thus employ the ring-polymer instanton approach 20,31,32 in order to obtain a practical method which can be used to locate the pathways and evaluate the level splitting in fulldimensional molecular systems.
In order to compute θ(τ k ), we only need to consider the contribution from a single kink or paths with no kinks and thus do not require a cyclic ring polymer but only a segment of it of length τ k . The path-integral representation of the propagator is discretized into N imaginarytime intervals using N − 1 beads, x i , which each correspond to a replica of the f -dimensional system. The imaginary-time action (also known as the Euclidean action) associated with this path is then 29 where τ N = τ k /N , x = {x 1 , · · · , x N −1 }, and x 0 and x N are the fixed end points, equal to x or x r as appropriate.
Performing the path integral within the steepestdescent approximation gives the single-kink contribution: 20,32 where J = τ N m ∇ 2 S(x), ∇ is a derivative with respect to x, and the primed determinant signifies that the zerofrequency mode associated with the permutational invariance of the instanton is not included in the product of eigenvalues. 20,32 The instanton path,x, is optimized so as to minimize the action, S(x). Thus, the implementation of this method is similar to that of ring-polymer instanton rate theory, 25,31,32,[58][59][60]72 except for the main difference that the tunnelling pathway is not a saddle point but the minimum of the action.
Similarly, the contribution from the non-tunnelling paths starting and ending in the left well is where the fluctuation matrix is , and x is a collapsed path which has all its beads located at the bottom of the left well. The expressions for K r (τ k ) and J r are equivalently defined in terms of a path collapsed in the right well. 73 Taking the ratio of kernels as defined in Eq. (10) thus gives and we have indicated explicitly that these expressions must be converged with respect to the number of beads. For a calculation converged with respect to τ k and N , the action will be independent of the location of the centre of the kink. Therefore, as for a symmetric system, J still has a zero mode, which is integrated out in the standard way. 20,26,28,32 However, unlike for a symmetric system, Φ will depend weakly on the location of the kink as the fluctuations in the two wells are different. To be consistent with the derivation presented above, we will thus require that the kink is centred such that it reaches its maximum potential energy at imaginary time τ k /2 (at least to the nearest bead). In our experience, this is not difficult to achieve. For instance, after an optimization, one can easily correct a path by removing beads from one end and adding them to the other. 74 In order to apply this formulation to molecular systems in full-dimensionality, translational and rotational degrees of freedom are treated in the usual way 32 by removing their corresponding modes from the product of the eigenvalues when computing the determinants. Finally, one does not typically know the relative orientation of x and x r . 20 In this case, the whole path including its end points, {x 0 , . . . , x N }, must be optimized, for instance using the l-BFGS algorithm 75 or recently developed specialized methods based on extensions of the the nudgedelastic band approach. [22][23][24]39,40,76 One could also avoid the calculation of eigenvalues of J by solving a related asymmetric Riccati equation. 24

III. RESULTS
We will apply the new theory to one-and twodimensional double-well models to benchmark the approach against exact quantum-mechanical results. Finally, we apply the method to a set of isotopically substituted malonaldehyde molecules in full dimensionality and compare with experimental measurements.

A. 1D Asymmetric Double Well
In order to benchmark our new theory in the simplest asymmetric case, we modify the well-known quartic double-well potential, 27,28,70 which has previously been used to test the standard ring-polymer instanton theory. 20 To create an asymmetric double-well potential, the symmetric version can be multiplied by a polynomial P (x) = ax 2 /x 2 0 + bx/x 0 + c. The overall potential will remain a double well if the the parameters are chosen such that P (x) has only imaginary roots. We also choose c = 1 such that it does not significantly change the barrier height. We therefore define the one-dimensional potential as This model describes a system with two minima of equal potential but with one well sightly broader than the other 77 and is depicted in Fig. 1. In particular, the minima are located at x = ±x 0 with V (±x 0 ) = 0 and the harmonic frequencies are ω = ω 0 √ 1 + a − b and Hence, E 0 and d are easily found.
Note that in the case that both a 1 and b 1, then d ≈ 1 4 b ω 0 , such that we can identify b as a parameter which controls the asymmetry. In order to ensure that the model has only two minima, we additionally require that a > b 2 /4, although the parameter a does not otherwise have a significant qualitative effect on the shape of the potential. In our tests, we therefore use a = b, although this choice has little bearing on our conclusions, and vary the value of b to achieve various degrees of asymmetry. The remaining parameters are chosen to match  √ V ‡ such that ω 0 is fixed and hence E 0 ≈ 0.283 in each case. Table I presents the numerical results of the new asymmetric instanton theory and compares them with benchmark quantum-mechanical results obtained by a numerical solution of the Schrödinger equation using the discrete variable representation (DVR). 78,79 The ring-polymer instanton calculations were found to be converged in each case with τ k = 30 and N = 128. We compare three regimes, one where b is chosen to describe a weak asymmetry (d Ω), one with medium asymmetry (d ∼ Ω) and one with high asymmetry (d Ω). Ring-polymer instanton theory for the fully symmetric case (with a = b = 0 and hence d = 0) has already been tested in Ref. 20 and can even be evaluated analytically. 27 It was shown that the instanton approximation is more accurate for higher barriers, where Ω is small, and less accurate for lower barriers, where Ω approaches the size of ω 0 . This is expected of a semiclassical method which gives only the leading exponential behaviour 30 and is due to the accuracy of the steepest-descent approximation, which improves as the ratio between and fluctuations of the action decreases.
In the case of very weak asymmetry, the predictions for the level splitting are almost identical to those of the symmetric case, 20 a pattern which is reproduced in the quantum results. This is because ∆ is dominated by the tunnelling contribution Ω, which is barely affected by the slight asymmetry added to the system. On the other hand, for high asymmetry, the level splitting is dominated by d and again gives good agreement with the quantum mechanics. 80 However, for a medium asymmetry such that d and Ω are of comparable size, the tunnelling matrix element, Ω, decreases with increasing barrier height as before, but now there is also a significant contribution from d to ∆. Importantly, in each case our asymmetric instanton calculations still give accurate predictions for the level splitting with errors which are similar to those found in the symmetric system. Like them, the error decreases with increasing barrier height and due to the properties of the semiclassical approximation would become exact in the → 0 limit.

B. 2D Asymmetric Mode-Coupling Potential
A commonly used simple model of two-dimensional tunnelling is provided by the symmetric mode-coupling potential. 28 In this work, we generate an asymmetric mode-coupling potential by following the approach used in Sec. III A, 81 where ω y is the frequency in the y direction and α is a measure of coupling strength. As before we choose c = 1 and a = b but vary the value of b to generate systems in different regimes. The remaining parameters were chosen following Ref. 82 as m = α = 1, ω y = 0.8 and = 0.04. The AMC potential and its symmetrized version are shown in Fig. 3. It should be noted that in order to better illustrate the asymmetry of the potential-energy surface to the reader, we have chosen such a large value for a = b that the level splitting would be completely dominated by the zero-point energy difference. Numerical calculations were however performed with smaller values of b to demonstrate the most interesting regimes.
We present the level splittings calculated for varying degrees of asymmetry in Table II. The instanton results were converged using N = 1024 and τ k = 100. The trends for the 2D case mirror those of the 1D case. In particular, the prediction for the case of weak asymmetry is almost identical to the fully symmetric case. For medium asymmetry, both d and Ω make significant contributions to ∆. Finally for the most asymmetric case where d Ω, the level splitting is completely dominated by the zero-point energy and ∆ ≈ 2d as expected. It can be seen that the agreement between the instanton approach and the quantum-mechanical benchmark is excellent, with errors under 10% in all cases.

C. Isotopic Substitutions in Malonaldehyde
The proton transfer in malonaldehyde has been used in many studies as a benchmark system for FIG. 3: Contour plots of the (a) symmetric and (b) exaggerated asymmetric mode-coupling potentials. In each case the ring-polymer instanton path is shown (by the blue curve connected by the ring-polymer beads represented by the blue circles). Also indicated is the norm of the quantum flux or current density (shown by yellow shading). 83 The instanton path follows the regions of high flux density and thus clearly gives a good description of the dominant tunnelling process in both the symmetric and asymmetric case.  The majority of these studies, however, do not examine the effects of isotopic substitution, or at best, only deuterate one hydrogen atom (called D 6 in the notation employed below). Experimental studies by Baughcum et al. 4,5 examine isotopic substitution in more detail, and in this work we shall compare our results to their measurements. We thus evaluate the level splitting for a number of isotopically symmetric and asymmetric species of malonaldehyde in full dimensionality. Following Ref. 4, we shall label the atoms as indicated in Fig. 4 and denote the various iso-topomers by stating only the substituted atoms labelled with a subscript.

Instanton calculations
The molecular system is described by the ground-state Born-Oppenheimer potential-energy surface (PES). We have employed two PESs for these calculations; one constructed by Mizukami et al. 13 (referred to as PES1) and one by Wang et al. 14 (referred to as PES2). Both are fitted to high-level electronic-structure calculations based on coupled-cluster theory. The more recent PES1 also employs explicit correlation using F12* theory 87 and reports a lower fitting error and is thus expected to be the more accurate of the two, although the differences in this case are not particularly dramatic.
Instanton optimizations were carried out by minimizing the action in Cartesian coordinates as outlined in Refs. 20 and 32. In addition, we must analyse the minimum-energy configurations on both sides of the barrier to obtain their harmonic frequencies. With this information we can compute the necessary quantities required to evaluate θ(τ k ) and hence predict Ω. Calculations were performed with increasing values of N and τ k until θ(τ k ) was converged. This was achieved at N = 1024 and an effective "temperature" of /(k B τ k ) = 30 K (i.e. τ k ≈ 10 000 a.u.). Finally, d and Ω were combined to predict ∆.
The symmetric parent molecule has been used widely as a benchmark to test the accuracy of quantum dynamics methods and potential-energy surfaces and thus results obtained with various methods are available. Two different MCTDH calculations 9,10 using PES2 obtain values for ∆ of 23.4 and 23.8 cm −1 for the parent configuration. In relatively good agreement with this, a calculation based on combining fixed-node DMC and MULTIMODE results 14 on PES2 obtained a tunnelling splitting in the range 21-22 cm −1 for the parent molecule and 2-3 cm −1 for the symmetric (D 6 ) isotopomer (both with statistical uncertainties of between 2 and 3 cm −1 ). Using the same PES, our ring-polymer instanton results give values of 24.9 and 3.29 cm −1 for these cases. Note that, a similar calculation based on a different implementation of instanton theory 22 reports tunnelling splittings within 0.1 cm −1 of those calculated in this work. Results generated using the newer PES1 will not necessarily be the same due to slight differences in the level of electronic-structure employed as well as the fitting method when constructing the PESs. Using PES1, the DMC results of Mizukami et al. 13  This figure also defines our labelling scheme, which follows that used by Baughcum et al. 4 dimensional models. Our predictions using the two different PESs are noticeably different, with one tending to overpredict and the other underpredict the experimental splittings. The fact that the results of different PESs vary on the same order of magnitude as the discrepancy with experiment suggests that the inherent error of the instanton approximation is not any larger than the errors introduced by the potential-energy surfaces. It will thus be impossible to quantitatively calculate the tunnelling splittings to high accuracy with these state-of-theart PESs, no matter what quantum-dynamics method is used. However, with instanton theory, one may still be able to explain and predict the trends obtained upon isotopic substitutions, as we will now show. The calculated level splittings for various isotopic substitutions of malonaldehyde are presented in Table III along with experimental measurements where available. 5 The instanton pathway shows a similar behaviour in each case and is not significantly affected by whether the isotopomer is symmetric or asymmetric. A representative instanton pathway is presented in Fig. 4 which provides an intuitive picture of the tunnelling mechanism in malonaldehyde. For instance, it can be observed that the H 6 atom is the only one directly involved in the proton transfer process and makes the largest contribution to the tunnelling path. We can quantify this in the same manner as we have done in our previous work 88 by evaluating the squared mass-weighted path length, which is in this case proportional to the action, S. 89 From this, we find that atom 6 contributes about 70% of the total action, whether it is deuterated or not, although the overall path length (and hence the action) does of course increase significantly when a D 6 substitution is introduced. Tun-nelling is thus strongly dependent on the isotope of this particular atom, as our results shall demonstrate.
This atom continues to dominate the action for the instantons associated with every other isotopic substitution under study. While this might seem to imply that a reduced-dimensionality approach based on the dynamics of this single atom would yield decent results, previous work 34,90 has shown that in general it is necessary to include all degrees of freedom (even those not involved in the tunnelling pathway) to obtain an accurate value for Ω and thus the level splitting ∆.

Discussion of the Results
We first discuss the results for the symmetric isotopomers, which are presented in the top half of Table III. In these cases d = 0 and the level splitting, ∆ = 2 Ω, can be called a tunnelling splitting. From the results, one can note that for as long as H 6 is not substituted with deuterium, ∆ remains similar to its value in the parent molecule. However, upon making the isotopic substitution D 6 , ∆ decreases dramatically, as expected. The experimental results also reveal that introducing 13 C or 18 O isotopes causes a reduction in the tunnelling splitting. This trend is correctly captured by the instanton method due to the fact that the action increases slightly with the substitution of heavier isotopes into the molecular skeleton, which thus results in a slight decrease in the predicted tunnelling splitting. The presence of 18 O atoms appears to have a greater effect than 13 C, mostly because they are more strongly coupled to the transferred H 6 /D 6 atom and contribute about 20% of the path action, whereas the C atoms contribute about 10%. 91 However, isotopic substitution of the other three H atoms (i.e. D 7 , D 8 and D 9 ) appears to be more complex because they are only very weakly coupled to the tunnelling pathway and contribute less than 1% to the path action. Substituting atoms for heavier isotopes is guaranteed to increase the action. However, although they cause a tiny increase in the action, there is also a larger effect in the contribution to Φ from the fluctuations around the tunnelling path. Making the substitution D 8 appears to consistently decrease the splitting slightly. On the other hand, according to some of the calculations, adding D 7 D 9 may increase the splitting. This behaviour is also observed in the reported experimental splittings for isotopomers containing D 6 D 7 D 9 . Nonetheless, this sort of competing effect is hard to predict reliably with any quantum-dynamics approach, and we see that the two different PESs give opposite trends for D 7 D 9 compared with the parent molecule. Next, we discuss the level splittings for various asymmetric isotopic substitutions of malonaldehyde, which are presented in the bottom half of Table III. These calculations were carried out using the generalized version of ring-polymer instanton theory developed in this paper. Here, three different regimes are demonstrated: |d| Ω for isotopomers ( 13 C 2/4 ), (D 6 D 8 13 C 2/4 ) and (D 6 D 8 13 C 2/4 13 C 3 ); |d| ∼ Ω for isotopomers (D 7/9 ), (D 7/9 D 8 ) and (D 6 D 8 18 O 1/5 ); and |d| Ω for isotopomers (D 6 D 7/9 ) and (D 6 D 7/9 D 8 ).
For the weakly asymmetric regime (|d| Ω), it can be observed that the level splitting is not significantly affected by the introduction of asymmetry, as was expected. This results in a level splitting for ( 13 C 2/4 ) which is only slightly smaller than in the parent molecule. This behaviour is confirmed by the available experimental observations and can be attributed to the small skeletal rearrangement of the carbon atoms, which increase the action slightly on isotopic substitution, just as in the symmetric case. For the case of (D 6 D 8 13 C 2/4 ) and (D 6 D 8 13 C 2/4 13 C 3 ), the situation is not quite so simple as there are competing effects from an increased |d| but a decreased Ω, making it difficult to reliably predict the size of the level splittings relative to (D 6 D 8 ). It is however clear that (D 6 D 8 13 C 2/4 ) and (D 6 D 8 13 C 2/4 13 C 3 ) have slightly larger splittings than (D 6 D 8 13 C 2 13 C 3 13 C 4 ) and (D 6 D 8 13 C 2 13 C 3 13 C 4 ) respectively due to their lower mass.
In the opposite extreme, |d| Ω, the level splitting is dominated by the asymmetry. The experiments 4 found no evidence of mixing between the two wells for (D 6 D 7/9 ) or (D 6 D 7/9 D 8 ). Our calculations confirm this scenario as the predicted mixing angles are about 5 • (PES1) or 7 • (PES2) in each case, which corresponds to a population ratio of more than 99:1.
For the intermediate cases (|d| ∼ Ω), both asymmetry and tunnelling play a significant role in determining the level splitting. Just like with the symmetric isotopomers, here it can be noted that when atom 6 is deuterated (i.e. isotopomers including D 6 ), a significant decrease in tunnelling (quantified by Ω) can be observed. Unfortunately no direct experimental values for the level splitting are reported for this set of isotopomers, although these are probably the most interesting cases. However, Baughcum et al. 4 do measure dipole moments and thereby roughly estimate the amplitude of the groundstate wavefunction in the two wells for (D 7/9 ). Their estimate is equivalent to a mixing angle of φ ≈ 41 • . Our instanton calculations predict a value of 31 • on PES1 or 44 • on PES2. This result is obviously quite sensitive to the PES, and the experimental estimate itself is only a rough approximation. However, it is nonetheless in agreement with the interpretation that the state is partially mixed, and demonstrates that our method is valid for this most difficult intermediate regime.
Bosch et al. 84,85 attempted to explain the experimental results using quantum-mechanical calculations based on a two-dimensional model system. Although the agreement with experiment is by no means quantitative [e.g. the predicted tunnelling splittings are 8.2 cm −1 for the parent molecule and 0.3 cm −1 for (D 6 )], 85 they were able to roughly predict some of the trends, most notably the strong decrease in tunnelling splitting upon substituting with the D 6 atom. Recent two-dimensional calculations 86 using an improved surface based on PES2 obtain more accurate results, but are still based on the reduced-dimensionality approximation, which is known not to be reliable in general. 34 Bosch et al. 85 also study the asymmetric (D 6 D 7/9 D 8 ) isotopomer and find that the stationary states remain practically unmixed (φ ≈ 1 • ) due to strong asymmetry, which is in agreement with both our work and the experiment. However, for the (D 7/9 D 8 ) isotopomer, they predict a mixing angle of only 16 • , whereas our calculation gives a value at least twice as large (30 • for PES1 and 43 • for PES2). This is quite a significant difference as it suggests that the ground-state population of the higher-energy well could be up to 6 times larger according to our theory. The main cause of the discrepancy is probably due to the fact that our PESs are based on much more accurate electronic structure than was available to Bosch et al. However, another important factor is the error introduced by their reduced-dimensionality approximation, which we are able to avoid by using a full-dimensional approach. Consequently, we conclude that the stationary states are much more mixed than was previously expected and that, even in these asymmetric systems, tunnelling still has a significant effect on the level splitting (such that it is about 5-9 cm −1 larger than in the (D 6 D 7/9 D 8 ) isotopomer even though the zero-point energy differences are almost the same). In a similar way, the (D 6 D 8 18 O 1/5 ) isotopomer will be partially mixed due to its particularly large d value caused by the isotopic substitution of the strongly coupled oxygen atom. Here the calculated mixing angle is 29 • on PES1 and 34 • on PES2 and it is therefore expected to demonstrate a particularly large level splitting relative to symmetric isotopomers with a D 6 substitution. These predictions could potentially be confirmed by future experiments.
Finally, we discuss the physical interpretation of our findings. As is clear from the theory presented in Sec. II, the level splitting is determined by two quantities d and Ω, which quantify the effects of asymmetry and of tunnelling respectively. Our theoretical approach is able to calculate these two quantities independently of each other, which allows us to identify the relative magnitude of the two components. In contrast, experimental measurements of spectroscopic transitions can only provide direct information on the value of ∆ and not the individual contributions. Based on their observations, Baughcum et al. 4 argue that tunnelling is "quenched" in asymmetrically substituted isotopomers. Whether or not this statement is consistent with our findings is perhaps somewhat a matter of semantics. In particular, if one quantifies tunnelling using the magnitude of Ω, then our calculations demonstrate that Ω is approximately unaffected by asymmetric isotopic substitutions and has a significant dependence only on the D 6 substitution. We do however find that the mixing angle may vary strongly on the type of asymmetry in such a way that one can say that it is not tunnelling but mixing which is quenched. This subtly different interpretation remains consistent with all the experimental observations. We find that the asymmetric isotopomers of malonaldehyde exist in three regimes which exhibit different quantum-dynamical behaviour. To some extent, one can design the isotopic substitutions to demonstrate a "sweet spot" for which |d| ∼ Ω, and it would be interesting to test our predictions for the level splittings in this regime with new experiments. We hope that our interpretation of the dynamics will be useful in explaining experimental observations and propose that the new method we have presented is used to help analyse and predict the level splittings in new molecules.

IV. CONCLUSIONS
In this work we have generalized instanton theory to describe tunnelling in a certain type of asymmetric system. Although the final result for the level splitting is expressed by a trivial formula defined in terms of the zero-point energy difference, d, and the tunnelling matrix element, Ω, an instanton approach to compute the latter did not previously exist. Note that even if the system is almost symmetric such that Ω dominates the splitting, one could not have simply used the standard formula from Ref. 20 as this requires perfect symmetry in order to converge.
Nonetheless, the final result shows that the new approach can be implemented in a very similar way to the standard ring-polymer instanton method and can thus be applied to complex molecular systems in full dimensionality in an efficient manner. The standard approach is currently available in i-PI 92 and MOLPRO 93 and can thus be used to study tunnelling in asymmetric systems too with very minor modifications to the code. One could also employ machine-learning techniques as has been done with other instanton approaches 88,94,95 to allow efficient application with ab initio electronicstructure calculations.
We have applied the new methodology to benchmark one-and two-dimensional systems and shown that the method can predict level splittings in asymmetric systems just as accurately as standard instanton theory applied to symmetric systems. Finally, we have calculated the level splittings in various asymmetric isotopomers of malonaldehyde. For this molecule, we find that the value of Ω is only weakly dependent on the isotopes of all atoms except the transferred proton (H 6 /D 6 ) and to a first approximation can be approximated by the value of the most similar symmetric isotopomer. The zero-point energy difference, d, is the only other factor which determines the level splitting and is controlled by the presence of isotopic substitutions in asymmetric positions. According to our definitions, the tunnelling effect and mixing angle can be studied independently. We find that although tunnelling (defined by Ω) is not strongly affected by asymmetric substitutions, the mixing angle can be. The most interesting regime (0 φ 90 • ) occurs when |d| ∼ Ω, which is exhibited by asymmetric isotopomers with a transferring H 6 atom and a D 7/9 substitution or by those with a D 6 atom and a 18 O 1/5 substitution. These simple rules would help assign and interpret any new experiments performed on similar molecules, and new calculations could easily be performed for other systems of interest.
One could envisage further extensions of this methodology, for instance to describe multi-well problems such as those occurring in water clusters. 24,[35][36][37][38][39][40] To go beyond the semiclassical approximation, it may be possible to apply path-integral sampling approaches 96,97 to these asymmetric tunnelling problems as well as anharmonic calculations of the zero-point energy.
Finally, we note that alongside its early development for the study of tunnelling in molecular systems, 98,99 instanton theory was more widely employed in quantum field theory. 26,27,70 There may thus also be uses in statistical and particle physics for this asymmetric instanton theory.