Solving Schrodinger: WKB, instantons, and the double well

We numerically test the predictions of the ‘instanton’ path integral predictions for tunneling in a symmetric double well. A symmetric double well is typically described by a quartic polynomial, V0 (y/Q0-1)^2 (y/Q0+1)^2, with minima at +-(Subscript[Q, 0]).

V[y_] := ...;

We set the potential strength to make the oscillation frequency in the right well equal to ω. We use the rms harmonic oscillator (Subscript[a, 0]) as usual to set the scale.

Clear[V0]
Solve[(D[V[y],{y,2}] /. y->Q0) ==..., V0]

V0 = ... ;
a0= ...;

We can calculate the WKB factor for this well, if (Subscript[Q, 0])= n (Subscript[a, 0]).

S0 = Integrate[..., Assumptions->{m>0,ω>0, Q0>0}]

S0/ .../. Q0->...

(Subscript[Q, 0])=5 (Subscript[a, 0]) works pretty well numerically...

Q0 = ...;

Specify the constants.

ℏ = ...
ω = ...
m = ...

Real-space lattice of Np points x, from -L/2 to L/2.

...
Plot[V[x],{x,-L/2,L/2}]
x = Table[...];

Time steps as a fraction of the period. Unitary spatial evolution UpotDtOver2.

...
UpotDtOver2 = ...;

Kinetic energy unitary operator UkinTildeDt:

...
UkinTildeDt = ...;

Initial wavefunction ψ[0], Gaussian centered at -Q0.

Evolve including the double-well potential part of the evolution: Upot(dt/2) applied before and after the kinetic term Ukin(dt).

ψ[0] = ...;
ψMax = Max[Abs[ψ[0]]];
T = ...;
Nt = ...;
times = ...;
For[...,
ψ[n] = ...;

Plotψ2[n_] := ListPlot[Transpose[{x, Abs[ψ[n]]^2}], Joined->True, PlotRange->{0,ψMax^2}, PlotLabel->"ω t = "+ω times[[n+1]]]

How long does it take for the initial Gaussian on the left to oscillate to the right and return? We include ω t in the label (title) for the plot, for your convenience.

ListAnimate[Table[Plotψ2[n],{n,0,Nt-1,100}]]

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

RHS=Table[Sum[... dx, {m,1,Np/2}],{n,0,Nt-1}];
rhsPlot = ListPlot[Transpose[{times,RHS}], Joined->True]

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 (Subscript[H, TLS])={{0,-Δ},{-Δ,0}}, where {1,0} is the state localized in the left well and {0,1} is the state on the right. We solve for the time evolution, and then evaluate the probability that a state that starts in the left well is in the left well after time t:

Htls = ...;
Eigenvalues[...]
Utls = MatrixExp[-I ...]
Abs[{1,0}. ... .{1,0}]...

We now fit the data to this TLS model. First estimate Δ roughly; it’s less than a percent of ℏ ω. You’ll need to start Δ within perhaps 30% of the correct value to avoid getting stuck in local minima.

model = ...;
ΔFit = Δ /. FindFit[Transpose[{times,RHS}], model,{{Δ,...ℏ ω/ ...}}, t]; 
fitPlot = Plot[model /. Δ-> ΔFit, {t,0,times[[-1]]}];
Show[fitPlot, rhsPlot]
EigensplittingRatio = ...

Both WKB and the instanton method tell us that Δ = ℏ(Subscript[ω, 0]) (e^(-(Subscript[S, 0])/ℏ)), with (Subscript[S, 0]) = ∫ Sqrt[2mV(y)] dy, integrated between the bottom of the two wells. The prefactor (Subscript[ℏω, 0]) is complicated to calculate, but it should be of order ℏω. We can use our numerical method to estimate (Subscript[ω, 0]) / ω:

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 was (Sqrt[((6(Subscript[S, 0]))/π)]) ℏω. How close is your answer to theirs? Does it seem likely that the difference is due to the inaccuracy of the WKB/instanton approximation?

GPPrefactorRatio = ...