# Tunneling problems by quantum Monte Carlo

## Abstract

We develop a new numerical scheme which allows precise solution of coherent tunneling problems, i.e., problems with exponentially small transition amplitudes between quasidegenerate states. We explain how this method works for the single-particle (tunneling in the double-well potential) and many-body systems (e.g., vacuum-to-vacuum transitions), and gives directly the instanton shape and tunneling amplitude. Most importantly, transition amplitudes may be calculated to arbitrary accuracy (being limited solely by statistical errors) no matter how small are their absolute values.

Tunneling phenomena are among the most intriguing consequences of quantum theory. They are of fundamental importance both for the high-energy and condensed matter physics, and the list of systems which behavior is governed by tunneling transitions ranges from quantum chromodynamics  to Josephson junctions and defects in crystals (see, e.g., Ref. ).

Precise analytic treatment of tunneling in complex systems is very hard, if not impossible. The most crucial simplification is in reducing the original problem to the semiclassical study of the effective action for some collective variable $`𝐑`$, with the assumption that all the other degrees of freedom adjust adiabatically to the motion of $`𝐑`$ . In certain cases, one may also include dissipative effects due to “slow” modes other than the selected collective variable which do not follow the dynamics of $`𝐑`$ adiabatically . Typically, the parameters of the effective action can not be found analytically (although one may relate some of them to the linear response coefficients) and have to be deduced from experiments. As far as we are aware, at present there are no tools to address the tunneling problem numerically. “Exact” diagonalization works only for relatively small systems, and, even in small systems, its accuracy is not sufficient to resolve very small energy splittings $`\mathrm{\Delta }E`$, say, when $`\mathrm{\Delta }E/E10^{10}`$, due to round-off errors (unless specific no-round-off arithmetics is used).

In this letter we develop a quantum Monte Carlo (MC) approach which allows precise calculations of tunneling amplitudes (and instanton shapes) no matter how small are their absolute values. The MC scheme contains no systematic errors, and its accuracy is limited only by statistical noise. The key point of our approach is in simulating imaginary-time dependence of transition amplitudes $`A_{ji}(\tau )`$ between selected reference states $`|\eta _1`$ and $`|\eta _2`$ (see below). In doing so we have to solve the problem of collecting reliable statistics in a case when $`A_{ji}(\tau )`$ varies, say, over hundreds of orders of magnitude (!) between different points in time. First, we describe in detail how to evaluate tunneling splitting and instanton shape in the double-well potential. We proceed then to the many-body problem of vacuum-to-vacuum transitions by considering the case of 1D quantum antiferromagnet with exchange anysotropy. Finally, we discuss the generality of the method suggested and demonstrate our numeric results for tunneling in the double-well potential (with comparison to the exact-diagonalization data where possible).

Consider the standard problem of particle motion in external potential:

$$H=m\dot{x}^2/2+U(x),$$
(1)   
which has two minima at points $`x=\eta _1`$ and $`x=\eta _2`$, and large tunneling action $`S=_{\eta _1}^{\eta _2}p𝑑x=_{\eta _1}^{\eta _2}[2mU(x)]^{1/2}𝑑x1`$. These minima are supposed to be near-degenerate, i.e., the lowest eigenstates of the Hamiltonian (1), $`H\mathrm{\Psi }_\alpha =E_\alpha \mathrm{\Psi }_\alpha `$, form a doublet with

$$E_2E_1e^S\omega _i\omega _i,$$
(2)   
where $`\omega _i`$ are the classical vibration frequencies in the potential minima $`\omega _i=[U^{\prime \prime }(\eta _i)m]^{1/2}`$ (we set $`\mathrm{}=1`$). The lattice analog of the Hamiltonian (1) reads

$$H=t\underset{<ll^{}>}{}d_l^{}d_l^{}+\underset{l}{}n_lU_l,n_l=d_l^{}d_l,$$
(3)   
where $`d_l^{}`$ creates a particle on the site number $`l`$, and the first sum is over nearest-neighbor sites.

The transition amplitude from the state $`|\eta _1`$ to the state $`|\eta _2`$, where $`|\eta _{1,2}=\delta (x\eta _{1,2})`$, in imaginary time $`\tau `$ is given by

$$A_{ji}(\tau )=\eta _2|e^{H\tau }|\eta _1\underset{\alpha }{}\alpha |\eta _1\eta _2|\alpha e^{E_\alpha \tau }.$$
(4)   
We now make use of the inequality (2) to define the asymptotic regime $`E_2E_1\tau ^1\omega _i`$, see, e.g., Refs. ,

$$A_{ji}(\tau )e^{\overline{E}\tau }\underset{\alpha =1,2}{}\alpha |\eta _1\eta _2|\alpha [1(E_\alpha \overline{E})\tau ].$$
(5)   
where $`\overline{E}=(E_2+E_1)/2`$.

It is convenient to split the double-well potential in two terms $`U(x)=U^{(1)}(x)+U^{(2)}(x)`$, where $`U^{(1,2)}(x)`$ is identical to $`U(x)`$ to the left/right of the barrier maximum point and remains constant afterwards. Introducing system ground states in each minimum as $`H^{(i)}\mathrm{\Psi }_G^{(i)}=E_G^{(i)}\mathrm{\Psi }_G^{(i)}`$, where $`H^{(i)}=m\dot{x}^2/2+U^{(i)}(x)`$, we may rewrite

$`\mathrm{\Psi }_1`$ $`=`$ $`u\mathrm{\Psi }_G^{(1)}+v\mathrm{\Psi }_G^{(2)},`$ (6)   
$`\mathrm{\Psi }_2`$ $`=`$ $`v\mathrm{\Psi }_G^{(1)}u\mathrm{\Psi }_G^{(2)},`$ (7)   
where $`(u^2,v^2)=1/2\pm \xi /2E`$, $`E^2=\mathrm{\Delta }^2+\xi ^2`$, with obvious identification of the energy splitting $`2E=E_2E_1`$ and bias energy $`\xi =(E_G^{(1)}E_G^{(2)})/2`$. Here $`\mathrm{\Delta }`$ is the tunneling amplitude, which defines energy splitting $`E_2E_1=2\mathrm{\Delta }`$ in the degenerate case. Substituting Eq. (7) into Eq. (5) we finally obtain

$`A_{ii}(\tau )`$ $``$ $`e^{\overline{E}\tau }Z_i^2,`$ (8)   
$`A_{ji}(\tau )`$ $``$ $`\tau e^{\overline{E}\tau }Z_iZ_j\mathrm{\Delta }.`$ (9)   
Here $`Z_i=\mathrm{\Psi }_G^{(i)}|\eta _i`$ projects the reference states $`|\eta _i`$ on the corresponding ground states in each minimum. All corrections to Eq. (9), e.g., the neglect of the overlap integrals $`\mathrm{\Psi }_G^{(1)}|\eta _2`$ and $`\mathrm{\Psi }_G^{(2)}|\eta _1`$, are small in parameter $`e^S`$. Note also that Eq. (9) does not depend on the bias $`\xi \omega _i`$, and is insensitive to the behavior of $`\mathrm{\Psi }_G^{(i)}`$ in the deep underbarrier region. This (rather standard) consideration relates transition amplitudes to the tunneling amplitude through the asymptotic analysis of $`A_{ij}(\tau )`$ in time.

MC simulation of the transition amplitude is almost identical to the standard simulation of quantum statistics \[the partition function at a given temperature $`T=1/\beta `$ may be written as $`𝒵(\tau =\beta )=\mathrm{Tr}_\eta A(\eta ,\eta ,\tau )`$\]. Fixed boundary conditions, as opposed to the trace over closed trajectories (configurations), are trivial to deal within any scheme. The real difference is in sampling different time-scales - thermodynamic calculations are typically done with $`\beta =const`$. Now we are forced to consider trajectories with different values of $`\tau `$ and to treat imaginary time as “dynamic” (in MC sense) variable. The idea of utilizing time dependencies of trajectories in MC simulations was extensively discussed in connection with the Worm algorithm  and polaron Green function .

We now turn to the problem of normalization. This problem might seem intractable in view of close analogy between $`𝒵(\tau )`$ and $`A(\tau )`$. Formally, in the limit $`\tau 0`$, the amplitude is trivial to find in most cases, e.g., for the Hamiltonians (1) and (3) it is given by the free particle propagation

$$A_{ji}(\tau 0)\{\begin{array}{cc}e^{m(\eta _i\eta _j)^2/2\tau }\sqrt{\frac{m}{2\pi \tau }}\hfill & \mathrm{continuous}\hfill \\ (t\tau )^{|\eta _i\eta _j|}/|\eta _i\eta _j|!\hfill & \mathrm{discrete},\hfill \end{array}$$
(10)   
and this knowledge may be used to normalize MC statistics for $`A_{ji}(\tau )`$ (in the discrete case $`|\eta _i\eta _j|=`$integer). However the absolute values of $`A_{ji}(\tau )`$ at short and long times will typically differ by orders and orders of magnitude and simply none MC statistics will be available at short times.

The solution to the puzzle lays in the possibility to use an arbitrary fictitious potential $`A_{\text{fic}}(\tau )`$ in Metropolis-type updates  in time-domain

$$\frac{P_{acc}(𝒜)}{P_{acc}(𝒜)}=e^{S_𝒜+S_{}}\frac{A_{\text{fic}}(\tau ^{})}{A_{\text{fic}}(\tau )}\frac{W_𝒜(\tau )}{W_{}(\tau ^{})},$$
(11)   
where $`P_{acc}(𝒜)/P_{acc}(𝒜)`$ is the acceptance ratio for the update transforming initial trajectory $``$ (having duration $`\tau `$) to the trajectory $`𝒜`$ (having duration $`\tau ^{}`$), and $`e^S_{}`$ is the statistical weight of the trajectory $``$. The normalized distribution functions $`W_{}(\tau ^{})`$, according to which a new value of $`\tau `$ is seeded, are also arbitrary; the best choice of $`W`$’s follows from the conditions of (i) optimal acceptance ratio (as close to unity as possible), and (ii) simple analytic form allowing trivial solution of the equation $`_{\tau _a}^\tau ^{}W(\tau ^{})𝑑\tau ^{}=r_{\tau _a}^{\tau _b}W(\tau ^{})𝑑\tau ^{}`$ on the time interval $`(\tau _a,\tau _b)`$ , where $`0<r<1`$ is the random number . Each trajectory adds a contribution $`=1/A_{\text{fic}}(\tau )`$ to the time-histogram of $`A_{ij}(\tau )`$.

One may use fictitious potential to enhance statistics of trajectories with certain values of $`\tau `$ “by hand”, e.g., by making $`A_{\text{fic}}`$ zero outside some time-window. To get a reliable and properly weighted statistics both at short and long times we need $`A_{\text{fic}}(\tau )1/A_{ij}(\tau )`$ to compensate completely for the severe variation of the transition amplitude between different time-scales. This goal is achieved as follows. The initial stage of the calculation, called thermolization, prepares the fictitious potential using recursive self-adjusting scheme - starting from $`A_{\text{fic}}(\tau )=1`$ in a given time-window $`(\tau _{min},\tau _{max})`$ and zero otherwise, we collect statistics for $`A_{ij}(\tau )`$ to the temporary time-histogram and after $`M>10^6÷10^7`$ updates we renew the fictitious potential as

$$A_{\text{fic}}(\tau )=\{\begin{array}{cc}A_{ij}(\tau _0)/A_{ij}(\tau ),\hfill & \tau _1<\tau <\tau _2\hfill \\ A_{ij}(\tau _0)/A_{ij}(\tau _1),\hfill & \tau _{min}<\tau <\tau _1\hfill \\ A_{ij}(\tau _0)/A_{ij}(\tau _2),\hfill & \tau _2<\tau <\tau _{max}\hfill \end{array}$$
(12)   
where $`\tau _1`$ and $`\tau _2`$ are the points (to the left and to the right of some reference point $`\tau _0`$) where temporary statistics becomes unreliable and has large fluctuations. It makes sense to select $`\tau _0`$ close to the maximum of $`A_{ij}(\tau )`$ (this point may be tuned a posteriory), and points $`\tau _1`$ and $`\tau _2`$ are formally defined as the first points in the histogram where smooth variation of $`A_{ij}(\tau )`$ ends: $`A_{ij}(\tau _{1,2}+\mathrm{\Delta }\tau )/A_{ij}(\tau _{1,2})<\delta `$ (here $`\mathrm{\Delta }\tau `$ is the difference between the nearest points in time-histogram, and $`\delta `$ is a small number, say 0.01). The thermolization stage continues until $`\tau _1=\tau _{min}`$ and $`\tau _2=\tau _{max}`$, and fictitious potential stops changing (withing a factor of two). After that the actual calculation starts with a fixed $`A_{\text{fic}}(\tau )`$, and a new histogram for $`A_{ij}(\tau )`$ is collected.

The idea of using fictitious potential proportional to the inverse of the transition amplitude is clear - it allows to collect reliable statistics on different time-scales with comparable relative accuracy. With this tool at hand one can easily normalize $`A_{ij}(\tau )`$ using known analytic results for short times, and deduce tunneling amplitude and Z-factors from the analysis of the long-time asymptotic, Eq. (9). We can not but note that fictitious potential in time-domain very much resembles the so-called “guiding wavefunction” in the Green-function MC methods , with essential difference that here it is used to reach exponentially rare configurations.

Instanton shape calculation is a much easier task since it can be done by considering trajectories for the transition amplitude $`A_{ij}(\tau )`$ with fixed but sufficiently long $`\tau 1/\omega _i`$; again, to ensure that $`A_{ij}`$ is dominated by just one instanton trajectory we need $`\tau 1/\mathrm{\Delta }`$. For any given MC trajectory $`x(\tau )`$ one has first to define the instanton center position in time, and to recount all times from this center-point $`\tau _c`$ \[instanton center statistics is almost uniform in $`(0,\tau )`$ (except near the ends of the time interval) due to the generic “zero-mode” present in instanton solutions \]. This can be done by looking at the average time

$$\tau _c=\frac{_B\tau 𝑑\tau }{_B𝑑\tau },$$
(13)   
where the integral is taken over the barrier region between the wells $`U(x(\tau ))>E_G`$. The instanton shape is obtained then by collecting statistics of $`x(\tau \tau _c)`$ to the time histogram. In this simple example when the notion of collective coordinate is not necessary (or formally, $`𝐑=x`$), we do not need to define separately the estimator for $`𝐑`$ (see the opposite example below).

The case of vacuum-to-vacuum transition in the many-body problem is formulated in precisely the same manner, and Eqs. (4-9), (11-12) holds true once the identification of the states $`|\eta _{1,2}`$ is done and the short-time limit, Eq. (10), is calculated. Consider as a typical example a 1D spin chain with $`2L`$ sites and antiferromagnetic (AF) couplings between the nearest-neighbor spins

$$H=\underset{<ij>}{}\left[J\stackrel{}{S}_i\stackrel{}{S}_j+J^{}S_i^zS_j^z\right],$$
(14)   
which has near degenerate ground state with the lowest doublet separated from the rest of the system spectrum by finite gap for $`J^{}>0`$. The natural choice of $`|\eta _{1,2}`$ is then an ordered AF state with $`S_i^z|\eta _{1,2}=(\pm 1)^i|\eta _{1,2}`$. Note, that for a large system Z-factors $`\mathrm{\Psi }_G^{(i)}|\eta _i`$ will be also exponentially small, and even diagonal amplitudes $`A_{ii}`$ have to be calculated with the use of $`A_{\text{fic}}`$ \[in the single-particle case one may obtain Z-factors by ignoring $`A_{\text{fic}}`$-trick\]. The short-time behavior is given by

$`A_{ii}(\tau )`$ $``$ $`1`$ (15)   
$`A_{ij}(\tau )`$ $``$ $`\left(J\tau /2\right)^L\{\begin{array}{cc}2\hfill & (\mathrm{ring})\hfill \\ 1\hfill & (\mathrm{open}\mathrm{chain}),\hfill \end{array}`$ (18)   
To decipher the instanton we need now some a priori knowledge about the relevant collective variable (if such knowledge is not available one has to study different possibilities). For example, if tunneling proceeds via two domain walls well-separated from each other (thin-wall approximation ), then the collective variable $`𝐑`$ is the distance between the walls, and the underbarrier region in Eq. (13) is related to the existence of two separated walls. These definitions are not too specific and work only approximately; this is however the generic difficulty of dealing with collective variables which are meaningful only in the macroscopic limit. Obviously, the knowledge of the instanton shape does not allow precise evaluation of $`\mathrm{\Delta }`$, and gives only rough estimation of $`\mathrm{ln}\mathrm{\Delta }`$.

To test the proposed scheme and to compare results to the exact diagonalization (ED) data we have applied our algorithm to the lattice model (3) with $`U(x)=U_0[(x/\eta )^21]^2`$. In what follows we measure all energies in units of the hopping amplitude $`t`$ and count them from the potential minimum. We set $`U_0=1`$ and consider two interwell separations: $`\eta =10`$ and $`\eta =40`$. For $`\eta =10`$ the ED data for the ground state energy, Z-factor, and tunneling amplitude are: $`E_G=(E_1+E_2)/2=0.1923`$, $`Z^2=0.2465`$, $`\mathrm{\Delta }=(E_2E_1)/2=3.6078\times 10^6`$. Our MC data give $`E_G=0.192(2)`$, $`Z^2=0.246(2)`$, and $`\mathrm{\Delta }=3.61(3)\times 10^6`$. The case of large $`\eta `$ is more subtle since only $`E_G=0.0495`$ and $`Z^2=0.1255`$ may be tested against ED - one may use textbook semiclassical analysis of the corresponding continuous model to see that $`\mathrm{\Delta }(\eta =40)10^{24}÷10^{23}`$, that is far beyond the standard computer round-off errors.

In Fig. (1) and Fig. (2) we present our MC data for the diagonal and off-diagonal amplitudes for $`\eta =40`$, and fits to the expected long-time, Eq. (9), and short-time behavior, Eq. (10) \[we also include the lowest-order correction for the potential energy at $`\tau 0`$ which tells that $`A_{ij}e^{\overline{U}\tau }`$ where $`\overline{U}=_\eta ^\eta U(x)𝑑x`$\]. The variation of $`A_{ij}`$ in Fig. 2 is about two-hundred orders! From Fig. (1) we deduced $`E_G=0.0494(2)`$ and $`Z^2=0.126(1)`$ in agreement with ED. From Fig. (2) we then obtain $`\mathrm{\Delta }=1.7(4)\times 10^{23}`$. MC simulation of the above figures took about 5 days each on PII-266 processor.

In Fig. (3) we present our results for the instanton trajectory $`x(\tau )`$ (at $`\eta =40`$) along with the known semiclassical results for the continuous model (1) . The accuracy of the data is self-evident, although we argue that MC data rather represent the trajectory with finite energy $`E=E_G`$ while analytic results correspond to $`E=0`$.

We note that the present technique is very hard to implement in the discrete-time schemes with finite Trotter parameter $`\mathrm{\Delta }\tau `$. On one hand, long-time asymptotic regime (see Fig. 2) requires to consider $`\tau `$ as long as $`640`$. On another hand at short times, the requirement of smooth variation of the amplitude $`\mathrm{\Delta }\tau d\mathrm{ln}[A_{ij}(\tau )]/d\tau 1`$, after substituting Eq. (10), means $`\mathrm{\Delta }\tau \tau /(2\eta )10^3`$ for $`\tau =0.1`$. To avoid large systematic errors due to time-discretization at short times one has to use $`\mathrm{\Delta }\tau =10^4`$(!!). Apart from enormous memory usage (there will be about $`10^7`$ time slices) that small Trotter parameter severely slows down the efficiency of the code, in fact, we are not aware of any MC simulation with $`\mathrm{\Delta }\tau 10^4`$.

It is worth mentioning that similar technique makes it possible to study directly $`𝒵(\tau )`$ over different time-scales - normalization of the partition function does not matter since none physical quantity depend on it. Thus one can obtain temperature dependences of the free energy, entropy etc. in a single MC run.

This work was supported by the RFBR Grants 98-02-16262 and 97-02-16548 (Russia), IR-97-2124 (European Community).