Notebook

Solving Schrodinger's equation: WKB, instantons, and the double well¶

We numerically test the predictions of the 'instanton' path integral predictions for tunneling in a symmetric double well.

Start this notebook with "ipython notebook --pylab" (without the final 'inline').

A symmetric double well is typically described by a quartic polynomial, $V_0 (y/Q_0-1)^2 (y/Q_0+1)^2$ , with minima at $\pm Q_0$ . We set up constants and the unitary operators as usual.

In [ ]:

from scipy import fft, ifft
from matplotlib import animation

# Constants
...

# Derived constants for potential
a0 = ...
Q0 = ...
V0 = ...

# Real-space lattice
...

# Quartic well
V = ...
plot(x,V)

# Time evolution
...
psi = (1.+0.j) * zeros((Nt, Np))

# Unitary potential energy evolution operator
UpotDtOver2 = ...

# Fourier space lattice
...

# Unitary kinetic energy evolution operator
UkinTildeDt = ...

Define the initial wavefunction $\psi_0$ centered in the left well

In [ ]:

psi[0] = ...
psiMax = abs(max(psi[0]))

Evolve in small time steps dt, using the Baker-Campbell-Hausdorff formula:

$U(dt) = U_{pot}(dt/2) U_{kin}(dt) U_{pot}(dt/2)$

$\exp(- i (p^2/2m + V(x)) dt/\hbar) \sim \exp(- i V(x) dt /(2 \hbar)) \exp(-i (p^2/2m) dt/\hbar) \exp(- i V(x)dt /2\hbar).$ Store

$\psi_n(x) = \psi(x, n dt)$ .

In [ ]:

for n in range(1,len(times)):
    psi[n] = ...

Plot the probability density $|\psi(x)|^2$ at time $t = \mathrm{Period}/5$ , hence $n = Nt/10$ .

In [ ]:

plot(x, ...);

Animate the evolution, plotting $|\psi_n(x)|^2$ . Skip 100 frames between views.

Did you start the notebook without the "inline" option?

In [ ]:

fig = figure()
ax = axes( xlim=(-L/2,L/2), ylim = (0,psiMax**2) )
probDensityLine, = ax.plot(x, abs(psi[0])**2)

def animate(n):
    probDensityLine.set_data(x, abs(psi[n*100])**2)
    title("timeStep = " + str(100*n) + "; omega t = " + str(omega*times[100*n]))
    return probDensityLine

anim = animation.FuncAnimation(fig, animate, frames=Nt/100, interval=50)
show()

We can quantify this nicely by calculating the time-dependent probability of being in the right well.

(You can use 'slicing' to select the parts of psi**2 that are in the left hand well; A[0:50], for example, gives the array starting at zero and ending at 49.)

In [ ]:

rhs = [sum(...**2)*dx for n in range(Nt)]
plot(times, rhs)

Your plot should start at one, drop to zero, and then oscillate back.

This slow oscillation is due to the small excitation energy between the ground and first excited state of the double well. Since these two states are so low in energy compared to the others, we often approximate the double well as a two-level system (TLS). The Hamiltonian for a symmetric TLS is given by a 2x2 matrix $\left(\begin{smallmatrix}0&-\Delta\cr-\Delta&0\end{smallmatrix}\right)$ , where $\left(\begin{smallmatrix}1\cr0\end{smallmatrix}\right)$ is the state localized in the left well and $\left(\begin{smallmatrix}0\cr1\end{smallmatrix}\right)$ is the state on the right. You solve for the time evolution to get the probability rhsTheory that a state starting in the left well is in the left well after a time t.

Load the nonlinear least-squares package from scipy. Define a function that returns the residuals. Find the best fit, and plot along with the oscillation.

In [ ]:

import scipy.optimize
# For help, type scipy.optimize.leastsq?

def residual(Delta):
    return rhs - ...

[DeltaFit], ier = scipy.optimize.leastsq(residual, ... hbar*omega/ ...)
plot(times, rhs, times, ...)
EigensplittingRatio = ...
EigensplittingRatio

Both WKB and the instanton method tell us that $\Delta = \hbar \omega_0 \exp(-S_0/\hbar)$ , with $S_0 = \int \sqrt{2mV(y)} dy$ , integrated between the bottom of the two wells. The prefactor $\hbar\omega_0$ is complicated to calculate, but it should be of order $\hbar\omega$ . We can use our numerical method to estimate $\omega_0/\omega$ :

In [ ]:

S0 = ...
numericalPrefactorRatio = DeltaFit/...
numericalPrefactorRatio

The analytic form of the prefactor was calculated by Gildener and Patrascioiu [PRD 16, 423 (1977)] for our quartic potential. In our notation, their result is $\sqrt{6 S_0 / \pi} \hbar\omega$ . How close is your answer to theirs? Does it seem likely that the difference is due to the inaccuracy of the WKB/instanton approximation?

In [ ]:

GPPrefactorRatio = ...
GPPrefactorRatio