Quantum Monte Carlo tunneling from quantum chemistry to quantum annealing

TL;DR

This paper demonstrates that quantum Monte Carlo simulations can efficiently recover tunneling rates in high-dimensional models, with implications for quantum chemistry and quantum annealing.

Contribution

It extends previous results to continuous variable models, showing QMC's effectiveness in capturing tunneling phenomena in complex systems.

Findings

QMC simulations recover ground state tunneling rates via instanton paths.

QMC results have implications for quantum chemical reactions.

Quantum tunneling can be a resource in quantum annealing.

Abstract

Quantum Tunneling is ubiquitous across different fields, from quantum chemical reactions, and magnetic materials to quantum simulators and quantum computers. While simulating the real-time quantum dynamics of tunneling is infeasible for high-dimensional systems, quantum tunneling also shows up in quantum Monte Carlo (QMC) simulations that scale polynomially with system size. Here we extend a recent results obtained for quantum spin models {[{Phys. Rev. Lett.} {\bf 117}, 180402 (2016)]}, and study high-dimensional continuos variable models for proton transfer reactions. We demonstrate that QMC simulations efficiently recover ground state tunneling rates due to the existence of an instanton path, which always connects the reactant state with the product. We discuss the implications of our results in the context of quantum chemical reactions and quantum annealing, where quantum tunneling is expected to be a valuable resource for solving combinatorial optimization problems.

Equations

$$Z=\int dx⟨x|e^{-\beta H}|x⟩$$

$$Z\propto \int d{x}_{1}d{x}_2\cdots d{x}_P\exp \sum _{i=1}^P{S}_i\text{leavevmode},$$

$$\exp (-S[{\mathbf{x}}^{\ast }(\tau )])\propto \Delta \text{(}instanton),$$

$$\exp (-S[{\mathbf{x}}^{\ast \ast }(\tau )])\propto {\Delta }^2\text{(}double\text{leavevmode}instanton),$$

$$V(x)=\{\begin{array}{ll}\lambda (x-x_0{)}^4-(x-x_0{)}^2,& x\ge x_0\\ 0,& -x_0\le x\le x_0\\ \lambda (x+x_0{)}^4-(x+x_0{)}^2,& x\le -x_0\end{array}$$

$$H=\Theta +V_A\text{leavevmode},$$

$$\Theta =-\frac{g^2}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)\text{leavevmode},$$

$$V_A(x,y)=\{\begin{array}{ll}\frac{1}{2}(x+1{)}^2+\frac{1}{2}{\omega }_y^2(y+y_0{)}^2,& x<0\\ \frac{1}{2}(x-1{)}^2+\frac{1}{2}{\omega }_y^2(y-y_0{)}^2,& x\ge 0\text{leavevmode},\end{array}$$

$$V_S(x,y)=\frac{1}{8}(x-1{)}^2(x+1{)}^2+\frac{{\omega }_y^2}{2}{\left(y+\frac{\gamma }{{\omega }_y^2}(x^2+1)\right)}^2\text{leavevmode}.$$

$$S$$

$$\simeq S[{\mathbf{x}}^{\ast }]+\frac{1}{2}{\int }_0^{\beta }d\tau \text{leavevmode}\text{leavevmode}\mathbf{y}(\tau )\hat{G}[{\mathbf{x}}^{\ast }]\mathbf{y}(\tau ),$$

$${\hat{G}}^{-1}[{\mathbf{x}}^{\ast }]({\tau }_1,{\tau }_2)=⟨\mathbf{y}({\tau }_1)\mathbf{y}({\tau }_2){⟩}_{{\mathbf{x}}^{\ast }},$$

Full text

I Introduction

Quantum mechanical tunneling (QMT) plays a fundamental role in a broad range of disciplines, from chemistry and physics to quantum computing. QMT can be observed in chemical reactionsNagel and Klinman (2006); Bell (2013); Zuev et al. (2003); Nakamura and Mil’nikov (2013) and affects the description of water and related aqueous system at room temperature.Ceriotti et al. (2016); Richardson et al. (2016) It is essential for understanding – even at the qualitative level – the phase diagrams of correlated materials, such as dense hydrogen, which is the simplest condensed matter system.Bonev et al. (2004); Pickard and Needs (2007); Dalladay-Simpson et al. (2016); Mazzola and Sorella (2017)

QMT can also been engineered in quantum annealers,Johnson et al. (2011); Bunyk et al. (2014) to solve optimization problems using quantum effects.Farhi et al. (2000); Kadowaki and Nishimori (1998); Farhi et al. (2001); Das and Chakrabarti (2008); Crosson and Harrow (2016) Here, quantum tunneling could provide a large advantage,Boixo et al. (2016) in particular when the energy landscapes display tall but thin barriers, which are easier to tunnel through quantum-mechanically, rather than to climb over by means of thermally activated rare events, whose frequency is exponentially suppressed as the height of the barrier increases.

In general, simulating real-time quantum dynamics requires the direct integration of the time dependent Schrödinger equation. This is a formidable task as the Hilbert space of the systems grows exponentially with the number of constituents, which makes the unitary evolution of a quantum system only possible for fairly small problem sizes, on the order of $40$ to 50 spins. The characterization of quantum dynamics simplifies when it is dominated by tunneling events. In this case, the useful quantities we want to predict is the transition rate between initial and final state (e.g. reactants and product in chemical reactions), and the pathway of the transition.

For simplicity let us first consider tunneling in a deep double-well system, well described by the lowest two eigenstates of the unperturbed tunneling system.$|\psi_{0}\rangle$, $|\psi_{1}\rangle$ which can be expressed as linear combinations of the degenerate states $|\psi_{L}\rangle$ and $|\psi_{R}\rangle$, localized respectively in the left and right well (see Fig. 1). The isolated system exhibits characteristic oscillatory behavior between the unperturbed states, $|\psi_{L}\rangle$ and $|\psi_{R}\rangle$, under the action of the Hamiltonian $H$, with frequency proportional to the tunneling matrix element $\langle\psi_{L}|H|\psi_{R}\rangle=\Delta/2$.

Coherence is easily destroyed by the presence of external noise, as is the case in the proton transfer reactions and in QA. Coupling to an environment can then stop the oscillatory and the transition rate is given by the incoherent tunneling rate, proportional to $\Delta^{2}$.Weiss et al. (1987) This is also the relevant tunneling rate in the adiabatic evolution of quantum annealing (QA), where the annealing time must scale as $\Delta^{-2}$, in order to avoid Landau-Zener diabatic transitions from the ground state to the first excited state.Farhi et al. (2000); Kadowaki and Nishimori (1998)

QMT also appears in QMC simulations, which can be efficient and scale polynomially with system size for quantum many-body problems without a sign problem ( i.e. that the system should obey bosonic statistics or distinguishable particles). Path integral Monte Carlo (PIMC) has been successfully applied to a broad range of continuum and lattice models. In particular PIMC simulationsBerne and Thirumalai (1986); Ceperley (1995a) have addressed problems in which QMT is important, such as proton delocalization in waterCeriotti et al. (2013); Morrone and Car (2008), hydrogenChen et al. (2013) and QA.Santoro et al. (2002); Heim et al. (2015)

PIMC is based on the path integral formalism of quantum mechanics and samples the density matrix corresponding to the quantum Hamiltonian $H$ by means of a classical Hamiltonian $H_{cl}$ on an extended system having an additional dimension, the imaginary time direction. The original quantum system is thus mapped into a classical one, which can be simulated by standard Monte Carlo sampling.

Although QMC techniques are rigorously derived to describe equilibrium properties, we here show that equilibrium PIMC simulations also provide important dynamical quantities, and in particular the the quantum tunneling rate. In Ref. Isakov et al., 2016 we have studied tunneling events in a ferromagnetic Ising model. The Ising ferromagnet can be described by an effective double well model, with the total scalar magnetization as reaction coordinate. We have numerically demonstrated that PIMC tunneling events occur with a rate $k$ which scales, to leading exponential order, as $\Delta^{2}$ – identical to the physical dynamics. We have also seen that with open boundary conditions (OBC) in imaginary time, the tunneling rate can becomes $\Delta$, thus providing a quadratic speed-up.

In this paper we investigate the scaling relation between the PIMC tunneling rate and $\Delta$ for a broader class of problems, of paradigmatic importance in quantum chemistry. We explore models where the effective one-dimensional picture of tunneling should break down.Takada and Nakamura (1994) Our results for continous variables extend the ones for the Ising modelIsakov et al. (2016) and we find that the QMC tunneling rate always follow the $\Delta^{2}$ scaling (or better with OBC). We argue that this is a manifestation of a general phenomenon, that QMC can efficiently simulate the tunneling splitting of the ground state energy levels in multidimensional systems.

II Instantons and QMC

II.1 Path Integral Monte Carlo

PIMC and path integral molecular dynamic (PIMD) techniques directly arise from the Feynman path integral formulation of quantum mechanics and are used to simulate thermodynamic equilibrium. To briefly introduce this approach for continuous space we start from the expression for the partition function $Z$:

[TABLE]

where $x$ is the particle position (the generalization to arbitrary dimensions is straightforward), $\beta=1/k_{B}T$ is the inverse temperature and $H$ is the Hamiltonian of the system. Typical real space Hamiltonians are sums of two non-commuting operators $H=\Theta+V$, where $\Theta=1/2m\leavevmode\nobreak\ \partial^{2}/\partial x^{2}$ is the kinetic operator ($m$ being the particle mass), and $V(x)$ is the potential energy. We first notice that the operator $e^{-\beta H}$ corresponds to an evolution in imaginary time $\beta$. We use the Trotter-Suzuki approximation $e^{-\delta_{\tau}(\Theta+V)}\approx e^{-\delta_{\tau}\Theta}e^{-\delta_{\tau}V}$ for small $\delta_{\tau}$.Ceperley (1995a)

Splitting the imaginary time evolution into $P$ small time steps of length $\delta_{\tau}=\beta/P$, the path integral expression for Eq. (1) then becomes

[TABLE]

where $S_{i}=K_{i}+U_{i}$ is the action of each step. $K_{i}=(x_{i-1}-x_{i})^{2}/(2\delta_{\tau}/m)$ is the kinetic part and $U_{i}=\delta_{\tau}/2(V(x_{i-1}-x_{i})$, in the so-called primitive approximation. Notice that $x_{1}=x_{P}$ (closed boundary conditions in imaginary time), for evaluating the trace of the density operator.

This provides an analogy between a quantum system and a classical system with an additional dimension: Eq. (2) is a classical configurational integral and the multidimensional object $(x_{1},\cdots,x_{P-1})\equiv\mathbf{x}(\tau)$ can be viewed as a ring-polymer, whose elements are connected by springs. Each element is labeled by its position along the imaginary time axis, with $0\leq\tau<\beta$. We refer to the Ref. Ceperley, 1995b for a detailed review of path-integrals. An essential feature of Eq. (2) is that the integrand is positive, and hence the distribution $\exp{\sum_{i=1}^{P}S_{i}}$ can be sampled by means of Metropolis Monte Carlo methods or Molecular Dynamics (MD) simulations. The main difference between a pure Monte Carlo vs a MD approach is that the latter samples from the canonical distribution by evolving an appropriate equation of motion, whereas the former uses stochastic Monte Carlo dynamics.

II.2 Instantons in PIMC

Connections between exact quantum dynamics and PIMD approaches, such as Centroid Molecular DynamicsCao and Voth (1993) and Ring Polymer Molecular DynamicsCraig and Manolopoulos (2004) have been discussed Jang et al. (2014); Braams and Manolopoulos (2006); Hele et al. (2015) in the context of real space simulations. Here we follow an alternative approach and summarize the picture of Refs. Isakov et al., 2016; Jiang et al., 2017 based on the instanton theory of tunneling through energy barriers.

In a PIMC or PIMD simulation one samples paths $\mathbf{x}(\tau,t)$ at each update along the simulation time axis $t$, and these paths are distributed according the functional $S(\mathbf{x}(\tau))$ as in Eq. (2). We can define an underlying pseudodynamics used to sample the paths to be given by a first order Langevin dynamics, $\partial\mathbf{x}(\tau,t)/\partial t=-\delta S/\delta\mathbf{x}(\tau,t)+\eta(\tau,t)$. In this case the analogy between quantum statistic and classical statistical mechanics has already been worked out in the stochastic quantization approach in the context of quantum field theory.Parisi et al. (1981) Here, the velocity of the (deformations of) path $\partial\mathbf{x}(\tau,t)/\partial t$ is linked to the generalized force $\delta S/\delta\mathbf{x}(\tau,t)$ and a Gaussian white noise $\eta(\tau,t)$ satisfying the obvious fluctuation-dissipation relation. We can numerically integrate the discretized version of the equation of motion (with time-step $\delta_{t}$), $\mathbf{x}(\tau,t+\delta_{t})=\mathbf{x}(\tau,t)-\delta_{t}\leavevmode\nobreak\ \delta S/\delta\mathbf{x}(\tau,t)+\sqrt{2\delta_{t}}{\bf z}(\tau,t)$, where ${\bf z}(\tau,t)$ is a deformation path, which, after a Trotter discretization is a vector of uniformly random distributed number in the range $[-1,1]$. This defines a Markov chain whose fixed point is the desired distribution, in the $\delta_{t}\rightarrow 0$ limit.

If the system displays two degenerate minima, then the transition state of the pseudodynamics is given by the point $\mathbf{x}_{TS}(\tau)$ satisfying $\delta S(\mathbf{x}_{TS}(\tau))/\delta\mathbf{x}(\tau)=0$ with the condition that $\mathbf{x}_{TS}(\tau)$ is not entirely contained in one of the attraction basins corresponding to the two minimaParisi et al. (1981); Sega et al. (2007); Autieri et al. (2009); Mazzola et al. (2011).

Finding this transition state is generally very complicated, but in the case of a double well potential $V(x)$ it can be done analytically. Here the dominant contribution to the integral comes from the stationary action path $\mathbf{x}^{**}(\tau)$ (determined exactly by the condition $\delta S(\mathbf{x}(\tau))/\delta\mathbf{x}(\tau)=0$) which is called instantonColeman (1977); Forkel (2000); Chudnovsky and Tejada (1998). This trajectory in imaginary time corresponds to a particle moving in the inverted potential $-V(x)$ (see Fig. 2). Following Ref. Isakov et al., 2016 it is possible to evaluate the action $S$ at this point and the amplitude is given by

[TABLE]

where $\mathbf{x}^{*}(\tau)$ is the open trajectory which connects the two classical turning points under the barrier, near the minima. Notice that, when computing the (diagonal) density matrix $\rho(x)$ periodic boundary conditions (PBC) in imaginary time are required. Now the integral over the closed paths it is dominated by the imaginary time trajectory $\mathbf{x}^{**}(\tau)$ that moves under the barrier starting, reaches the turning point, and returns. Therefore, the saddle point estimation of the integral gives a squared tunneling amplitude

[TABLE]

due to the cost of creating an instanton and an anti-instanton (see Fig. 2). Returning to the PIMD pseudodynamics, according to Kramers theoryHänggi et al. (1990), the escape rate is $k\propto e^{-S(\mathbf{x}_{TS})}$, and therefore $k\propto\Delta^{2}$ if standard closed path integrals are used, whereas $k\propto\Delta$ if the paths are opened. In Sect. III we extend the study of Ref. Isakov et al., 2016 and demonstrate that the quadratic speedup in the tunneling rate in the case of open boundary path integrals holds also in multi-dimensional continuous space problems.

III One-dimensional double well potential

Let us consider the following one-dimensional double well potential,

[TABLE]

with $\lambda,x_{0}>0$. We can separately tune the width and the height of the barrier, varying $\lambda$ and $x_{0}$. The height of energy barrier is $\Delta V=1/4\lambda$, and the distance between the two minima is $d=2(x_{0}+\sqrt{1/2\lambda})$ (see inset of Fig. 3) Decreasing $\lambda$ reduces the energy splitting $\Delta$, as the two wells become deeper and more separated. The parameter $x_{0}$ only increases the well separation but doesn’t change the potential energy barrier height. Moreover, a variation of $x_{0}$ leaves the characteristic frequency of the potential wells unchanged, i.e. the kinetic energy associated to the localized states $|\psi_{L}\rangle$ and $|\psi_{R}\rangle$.

Following Ref. Isakov et al., 2016 we measure the mean first tunneling time (MFTT), defined as the number of updates required to find the system in the right well, if the particle has been localized in the left one at the beginning of the simulation. From Fig. 3 we see that the MFTT scales as $1/\Delta^{2}$ when PBC are used, whereas it scales as $1/\Delta$ for OBC, as the parameters $x_{0}$ and $\lambda$ change. The gap $\Delta$ is obtained using a discrete variable representation (DVR) techniqueColbert and Miller (1992). This scaling relation holds for PIMC with local Metropolis updates and PIMD (using both first and second order Langevin thermostats), at large $\beta$, and in the limit of small time steps $\delta_{\tau},\delta_{t}\textrm{(for PIMD)}\rightarrow 0$ limit. This means that the scaling of tunneling rate in a double well model $k\propto\Delta^{2}$ is correctly reproduced.Weiss et al. (1987).

Why do the open paths tunnel faster from the point of view of PIMC pseudodynamics? To answer this question we first observe that, for sufficiently low temperatures, the center of the open-path $\mathbf{x}^{*}(\tau\approx\beta/2)$ sample from the ground state distribution $|\psi_{0}(x)|^{2}$, whereas the tails, $\mathbf{x}^{*}(\tau\approx 0)$ and $\mathbf{x}^{*}(\tau\approx\beta)$, sample from the ground state distribution $\psi_{0}(x)$. Therefore the tails spend more time inside the barrier (see Fig. 4) compared to the center, which follows instead the more localized $\psi_{0}^{2}$ distribution. Once that one of the two tails crosses the barrier, then the rest of the open-polymer may easily follow, so that the whole polymer ”tunnel” faster compare to its PBC counterpart. This also means that, with OBC, it is possible to sample from the equilibrium distribution $\rho(x)\approx|\psi_{0}|^{2}$, using the center of the path, while having a considerable speed-up in the sampling. We notice that this feature is not surprising as the OBC technique is closely related to the so called path integral ground stateSarsa et al. (2000) (PIGS) technique. Indeed, in the PIGSSarsa et al. (2000) approach, sampling from the tails gives the mixed distribution $\psi_{0}(x)\psi_{T}(x)$, but in our case the trial wavefunction is $\psi_{T}(x)=1$. Therefore, we propose that OBC should be used not only in the context of quantum annealing but much more broadly also in material simulations, as far as low-temperature conditions are investigated.

IV Multidimensional tunneling

The double well model provided in Sect. III is a prototypical example of one-dimensional tunneling. One could argue that, despite having many-spins degrees of freedom, the spin models investigated in Ref. Isakov et al., 2016 are also effectively one-dimensional models, as the relevant reaction coordinate is the total scalar magnetization $M$. Indeed the instantonic nature of the transition state can be seen if we plot $M(\tau)$ as a function of the imaginary time parameter $\tau$.

It is much more straightforward to devise models that require multidimensional tunneling in continuos space, rather than spin models.Stella et al. (2006a); Inack and Pilati (2015) To this end we borrow insights from quantum chemistry, where simplified model for characterizing proton tunneling have been devised.Auerbach and Kivelson (1985); Makri and Miller (1987); Takada and Nakamura (1994); Nakamura and Mil’nikov (2013) In particular in Ref. Takada and Nakamura, 1994 a semiclassical theory of multidimensional tunneling is formulated, unraveling its qualitative differences compared to one-dimensional tunneling. It was found that in multidimensional tunneling two regimes can be identified: the pure tunneling case, which is effectively one-dimensional, where the tunneling path can be defined uniquely, and the mixed tunneling regime when tunneling occurs very broadly, i.e. where a set of dominant semiclassical paths $\{{\bf x}_{TS}\}$ is not defined. In the first case, the action which defines the semiclassical wavefunction is purely imaginary, whereas in the latter the action is complex. We refer the interested reader to Ref. Takada and Nakamura, 1994 for the analytical details.

Investigating QMC simulations for such mixed tunneling models, where the QMC scaling relation with the exact QMT rate might be expected to break down, we instead find that the quantum tunneling rate given by QMC scales as the adiabatic quantum evolution also in this case.

IV.1 QMC tunneling rate scaling

We first consider the simple shifted parabola bidimensional model of Ref. Takada and Nakamura, 1994, which is a minimal model for the antisymmetric mode coupling mechanism for proton tunneling in malonaldehyde, a well-studied molecular test case. The Hamiltonian reads

[TABLE]

with

[TABLE]

where $g>0$ is a dimensionless parameter which sets the strength of the quantum fluctuations. The potential is

[TABLE]

where $y_{0}>0$ and $\omega_{y}^{2}>0$ are dimensionless harmonic potential parameters. This potential represents two parabolas, located respectively in the half-planes $x<0$ and $x>0$, with centers shifted along the $y$-axis by an amount $2y_{0}$. In the case of malonaldehyde molecule, the coordinate $x$ represents the motion of transferring the proton, while $y$ represents the C–O stretching mode. The sign of $a^{2}=y_{0}^{2}-g/\omega_{y}$, distinguishes between the two ground state QMT cases: pure tunneling for $a^{2}>0$ and mixed tunneling for $a^{2}<0$.

In Fig. 5 we present results of PIMC simulations with local updates, using PBC, at large $\beta$ (very low temperature), and in the converged time step $\delta_{\tau}\rightarrow 0$ limit, to describe faithfully ground state tunneling. The path deformations are obtained by displacing each bead at a time by an amout $(dx,dy)$. The displacements are Gaussian distributed with zero mean and the variance is tuned in order to obtain a Metropolis acceptance probability of $\approx 40\%$.

Again we study the MFTT obtained with PIMC simulations as a function of the parameter $\omega_{y}$, in the range $[0.05,2]$ and for three different choices of $g=0.2,0.3,0.4$, and for two shifting values $y_{0}=0$ and $1$. Following Ref. Isakov et al., 2016 we define the MFTT as the number of PIMC updates required to observe an instantonic state. In turn, we algorithmically define an instanton path, as spending approximately the same fraction of imaginary time in either well.

With these parameters ranges111The actual value would be $g\approx 0.1$, here we artificially enhance $\hbar$ in order to increase the observed tunneling rate. we can roughly mimic proton transfer reactions in malonaldehyde.Makri and Miller (1987) For this molecule, it is found that, if the tunneling is described only by a one-dimensional process, the tunneling rate is reduced by two orders of magnitude compared to experimental and recent theoretical valuesBosch et al. (1990); Cvitaš and Althorpe (2016). Furthermore, it was argued in Ref. Takada and Nakamura, 1994 that deviations from the one dimensional picture leads to a mixed tunneling regime where no well defined tunneling path exists. Therefore it could be possible that QMC underestimates the exact tunneling rate. In the context of quantum annealing problem, this might mean that the performances of QA and its simulated version through QMC could be very different, under these “mixed tunneling” conditions.

We first perform tests for $y_{0}=0$, where we are always in the regime $a^{2}<0$ , i.e. the mixed tunneling regime. The gap $\Delta$ is constant as a function of $\omega_{y}$, and we observe the same in QMC, where the MFTT remains constant. Its precise value depends on the parameter $g$. We use these data to fix the proportionality constant $\alpha(g)$, which we use later to compare the MFTT to the value $\alpha(g)\Delta^{-2}$. Notice that in the $\omega_{y}\rightarrow 0$ limit the two wells become parallel and indefinitely extended along the $y$ direction. In this limit, we observe an infinite number of tunneling paths that connects reactants on left well with the products on the right.

Next, we set $y_{0}=1$ and repeat the simulations. This time the scaling of $\Delta^{-2}$ as a function of $\omega_{y}$ is non-trivial. Nevertheless it approaches the previous value – for any given $g$ – in the $\omega_{y}\rightarrow 0$ limit, as the two wells are infinitely long and the shift given by $y_{0}$ becomes irrelevant. In this case, we cross the transition point $a^{2}=0$ between the the two regimes of tunneling, when $\omega_{y}=g$. We still observe a satisfactorily agreement between the QMC MFTT data series and the $\alpha(g)\Delta^{-2}$ functions. While a residual difference between the QMC data and the expected $\Delta^{-2}$ behaviour can now be appreciated, this difference is small over the broad $\omega_{y}$-range investigated, i.e. nowhere near the 2 order of magnitude worst case scenario, reported in Ref. Bosch et al., 1990. We note that in Ref. Bosch et al., 1990 the full multidimensional problem is reduced to a one-dimensional calculation. This is the origin of the observed large deviation from the theoretical value. Notice also that, once we fix the constant $\alpha(g)$ the QMC tunneling time is always smaller than $\alpha(g)\Delta^{-2}$, so QMC tunnels slightly more efficiently than QA, even with PBC.

IV.2 QMC reaction pathways and fluctuations around the instanton solution

In this Section we explicitly track the QMC pseudo-dynamics transition states and compare to the instantonic trajectory computed by minimizing the action $S$. Let us consider the symmetric mode coupling potential,

[TABLE]

This potential energy surface is continuos and has been widely used as a model for proton tunneling. In the typical example of malonaldehyde, the coordinate $x$ represents the motion of the proton transferring between the Oxigen atoms, while $y$ gives the scissors-like motion of the O-C-C-C-O frameTakada and Nakamura (1994). We use the dimensionless potential parameters $(\omega_{y},\gamma,g)=(0.48,0.39,0.10)$ to fit the model to the ab-initio potential energy surfaceBosch et al. (1990). In this way we can directly compare the transition paths given by the PIMC simulation with other techniques, such as the Ring Polymer Instanton (RPI) method,Richardson and Althorpe (2011); Richardson et al. (2011); Cvitaš and Althorpe (2016); Richardson et al. (2016) recently introduced to calculate energy splitting. In the RPI framework, one first needs to locate the saddle point of the action $S$ (the instanton), and then evaluates the splitting energy by computing the functional integral up to second order in the fluctuations around the dominant contribution. This approach employs neither PIMC sampling nor PIMD, as the instanton path is obtained via action’s minimization and the initial guess is an OBC path that already connects the reactant to product state.

In Fig. 6 we plot a sample of transition paths produced by the PIMC pseudodynamics and recognize their instantonic character. We compare some OBC transition paths sampled with our PIMC simulation, against the RPI solution recently published in Ref. Cvitaš and Althorpe, 2016. We see that these instanton paths form a bundle around the RPI saddle point solution, and are qualitatively distant from the minimum energy path (MEP), which would be typical of a classical thermally activated processCvitaš and Althorpe (2016). It is remarkable that a simple PIMC simulation obtains the instanton path, which is otherwise computed only by a complex minimization procedure as in the RPI scheme.

Another advantage of PIMC is that we can directly sample the statistical fluctuations around the dominant solution $\mathbf{x}^{*}(\tau)$. To second order the action can be expanded as Parisi et al. (1981)

[TABLE]

where $\mathbf{y}(\tau)$ is a fluctuation path, satisfying $\mathbf{y}(0)=\mathbf{y}(\beta)=\mathbf{0}$, and $\hat{G}=-\frac{d^{2}}{d\tau^{2}}+V^{\prime\prime}[\mathbf{x}^{*}]$ is the fluctuation operator, or Hessian, adopting the notation of Ref. Richardson and Althorpe, 2011. Here $V^{\prime\prime}[\mathbf{x}^{*}]$ is the second derivative of the potential computed along the instanton path $\mathbf{x}^{*}(\tau)$.

In practice one always deals with discretized trajectories in imaginary time. Therefore also the operator $\hat{G}$ is discretized using finite differences and then diagonalized to obtain the normal modes and frequency of the fluctuations. The resulting product of Gaussian integrals allows one to evaluate Eq. IV.2.

However, it could be cumbersome to evaluate $G$ for ab-initio potentials (as they require evaluation of second derivaties of the potential), or in the case of rugged energy landscapes, where local curvature at the saddle point $V^{\prime\prime}[\mathbf{x}^{*}]$ does not correctly represent the actual amplitude of the quantum fluctuationsFaccioli ; Mazzola et al. (2011). On the other hand, the inverse operator $G^{-1}$ can be computed stochastically with PIMC, using the relation

[TABLE]

where the right hand side denotes the statistical average of the fluctuations, around a given path $\mathbf{x}^{*}$. This approach gives a more effective and fast estimate of the curvature of the potential surface in the above cases.

We note that PIMC sampling techniques have already been used to compute tunneling splittings in molecular and condensed matter systemsKuki and Wolynes (1987); Ceperley and Jacucci (1987); Alexandrou and Negele (1988); Mátyus et al. (2016). Here we propose an alternative and simple way to compute ratios of quantum mechanical rate constants. Suppose that the potential displays several minima, i.e. that we have one reactant state $R$ and two, or more, possible product states $P_{1},P_{2}$. By computing the ratio of the average PIMC tunneling times, with OBC, required to perform the transitions $R\rightarrow P_{1}$ and $R\rightarrow P_{2}$ respectively, we can estimate the ratio of the tunneling splittings $\Delta_{R,P_{2}}/\Delta_{R,P_{1}}$ corresponding to the two quantum mechanical transitions, provided that enough statistics of instanton events can be gathered in a reasonable amount of time. This approach is predictive, as the instantons are generated by the PIMC pseudodynamics, without any a-priori knowledge of the final product state.

This problem is closely connected with quantum annealing, where, starting from some high energy “reactant” states $R$’s, one would like to reach the “product” state $P_{i}$ having the lowest possible potential energy, after a sequence of tunneling eventsKnysh (2016). The relative probability of finding this state, compared to other metastable ones, is well described by PIMC simulations.

V Conclusions

We have studied the tunneling of path integral based equilibrium simulations in continuos space models, generalizing a previous studyIsakov et al. (2016) on ferromagnetic spin models. We demonstrates that the PIMC tunneling rate scales as a $\Delta^{2}$ if periodic boundary conditions in imaginary time are used, while it scales as $\Delta$ with open boundary conditions. These scaling relations seem to be a general property of path integral methods, as long as reasonable semi-local updates are employed during the Markov chain pseudodynamics (see Sect. II.2). In this case, in double well potentials, it is possible to directly identify the transition state of the path integral pseudodynamics –a purely classical process– and therefore compute its classical reaction rate using Kramers theory.

This transition state is the instanton path, and we remark here that this trajectory is sampled by the PIMC pseudodynamics using local updates, i.e. we don’t need to engineer such kind of global update moves as in Refs. Stella et al., 2006b; Alexandrou and Negele, 1988. Indeed the latter approach would invalidate the premises and discussions presented in this paper and artificially increase the reaction rate observed with PIMC. On the other hand, building in instantonic updates in the Metropolis procedure requires a full knowledge of the system, i.e. knowing in advance the transition states. Having this knowledge one would solve beforehand the quantum annealing problem, for example, without even running any PIMC simulation.

The quadratic speed-up in tunneling efficiency is a robust feature of OBC simulations for tunneling through individual barriers. In the context of simulations, therefore, we propose that open path integral simulations should be used instead of PBC and will accelerate the sampling, whenever ground state properties are desired.

We also turned our attention to simplified models for proton transfer, where multidimensional tunneling is deemed to be important, and a semiclassical description of tunneling as an effective one-dimensional process has been seen to fail. Nevertheless, the scaling relation of the PIMD transition rate, compared to the exact incoherent tunneling rate $\Delta^{2}$ holds also in this case.

The above finding is very interesting because often in a multidimensional potential the smooth tunneling path (instanton) connecting the minima of the potential does not exists due to the effects of a so-called dynamical tunneling Huang et al. (1990); Takada and Nakamura (1994); Zamastil (2005); Razavy (2013). The lack of an instanton path was also studied, e.g., in the case of the two-dimensional shifted parabola model Ref. Takada and Nakamura, 1994 considered in our study. On the other hand, as explained above in Sec.  II.2 based on the Kramers theory arguments, the existence of the instanton path is a key requirement that leads to the identical QMT and PIMC scaling laws.

To explain this conundrum we observe that the difficulty with instanton description in the case of a tunneling in a multidimensional potential usually occurs when one needs to match the solutions given by Wentzel–Kramers–Brillouin (WKB) theory in classically allowed and forbidden regions at the boundary formed by caustics. Caustics result in the complex (oscillatory) behavior of the wavefunction under the barrier Huang et al. (1990); Takada and Nakamura (1994); Zamastil (2005). This oscillatory behavior results in a phase problem in QMC.

We argue that at zero temperature this situation does not occur generically, because a classically allowed region in configuration space “collapses” into the point corresponding to the minimum of the potential. As usually, to study tunneling one should consider the wavefunction under the barrier that nearly coincides with the ground state wavefunction near one of the minima of $V$, exponentially decaying away from it. The mechanical action $S[{\bf x}(\tau)]=\int_{\infty}^{\tau}(m\dot{\bf x}^{2}(\tau_{1})+V({\bf x}(\tau_{1})))d\tau_{1}$ for the wavefunction is associated with the unstable Lagrangian manifold Guckenheimer and Holmes (2013) formed by real-valued trajectories $({\bf x}(\tau),{\bf p}(\tau))$ in the phase space moving in the imaginary time $\tau$ in the inverted potential $-V$ (above ${\bf p}(\tau)$ is a system momentum). The trajectories emanate at time $t=-\infty$ from the corresponding maximum of the $-V$. In general, projections of the Lagrangian manifold onto the coordinate space ${\bf x}$ can have caustics, cusps (and more complex singularities in the dimensions higher then twoArnold (1989)) in classically forbidden region. These singularities lead to multi-valuedness of the action surface and some of its branches become complex. However the minimum action surface is real- and single-valued. It possesses lines where the surface gradient is discontinuous (see Fig. 7) They correspond to the so-called switching lines in configuration space Dykman et al. (1994). Points at different sides of the switching are reached by a topologically different imaginary-time paths as shown in Fig. 7(b). Therefore any point ${\bf x}$ can be reached by the most-probable path that provides the minimum of the action and never crosses a switching line. An instanton is a particular member of the minimum-action family of paths that connects the two maxima of the potential $-V$. It corresponds to the heteroclinic orbit (${\bf x}^{*}(\tau)$ ,${\bf p}^{*}(\tau)$) contained in the unstable Lagrangian manifold shared by the two maxima Ref. Sharpee et al., 2002. This explains why ground state tunneling splitting for the particle in a multidimensional potential is always described by the instanton with a real-valued action and therefore can be simulated efficiently by PIMC.

This confirms that PIMC simulation of QMT in the ground state can be done without any loss of efficiency compared to what real system would do. The fact that PIMC simulations have the same scaling with the problem size as physical quantum annealing was recently experimentally confirmed, again on a spin system on chimera graphDenchev et al. (2016). In this context, it is unlikely that QA to find a ground state of optimization problem can achieve an exponential speedup over classical computation, only by using QMT as a computational resource.

We remark here that these conclusions hold only when so-called stoquastic Hamiltonians are used, i.e. Hamiltonians which allow PIMC simulations. This is the case of the Hamiltonians used in this paper. In most reaction simulations protons are assumed to be distinguishable particles, and, in QA the standard transverse field Ising Hamiltonian is also stoquastic.Johnson et al. (2011); Bunyk et al. (2014) This provides additional evidences that QA machines should implement non-stoquastic Hamiltonians which display sign-problem,Seki and Nishimori (2012); Mazzola and Troyer (2017); Hormozi et al. (2016) to avoid efficient simulations by QMC methods.

Acknowledgements. GM acknowledges useful discussions with P. Faccioli, S. Sorella and G.E. Santoro. MT acknowledges hospitality of the Aspen Center for Physics, supported by NSF grant PHY-1066293. Simulations were performed on resources provided by the Swiss National Supercomputing Centre (CSCS). This work was supported by the European Research Council through ERC Advanced Grant SIMCOFE by the Swiss National Science Foundation through NCCR QSIT and NCCR Marvel. This paper is based upon work supported in part by ODNI, IARPA via MIT Lincoln Laboratory Air Force Contract No. FA8721-05-C-0002. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purpose not-withstanding any copyright annotation thereon.