Wavepacket in a double well

By Simon Rochester, UC Berkeley Physics Dept.
http://budker.berkeley.edu/

We find the eigenstates and eigenenergies of a double well potential. The potential is symmetric about the origin, so the eigenstates all have either even or odd parity.  For simplicity, we consider only those states with energy less than the potential height between the wells.

We then expand an initial state (e.g. the ground state of the right-hand well) in terms of the double-well eigenstates and plot the evolution.

To use, choose values for "params" below, and then evaluate the notebook.

Define the potential

b is the width of each well
a is 1/2 of the spacing between the wells
v = 2m V/(ℏ^2), where m is the mass and V is the height of the potential in between the wells

params={b->1,a->1/20,v->300}

Draw a picture of the potential:

potential=({EdgeForm[Black],LightGray,Rectangle[{-a,0},{a,1}],Rectangle[{-a-b-1,0},{-a-b,5}],Rectangle[{a+b,0},{a+b+1,5}]}/.params);

Plot[0,Evaluate[{x,-a-b-.5,a+b+.5}/.params],PlotRange->{0,4},Prolog->potential,Axes->{True,False}]

Find the even-parity eigenstates with E < V

The parity condition determines the form of the wave function in each region.

Wave function in left well (k= (Sqrt[ϵ]), where ϵ=2 m E/(ℏ^2)):

(Subscript[ψ, 1,e])[x_]=A Sin[k(x+a+b)]

Wave function between wells (κ= (Sqrt[v-ϵ])=(Sqrt[v-(k^2)]), we assume v>ϵ):

(Subscript[ψ, 2,e])[x_]=B Cosh[κ x]

Wave function in right well:

(Subscript[ψ, 3,e])[x_]=-A Sin[k(x-a-b)]

We need to determine the eigenvalues (in terms of k) and the constants A and B.

To find the eigenvalues, we use continuity of the logarithmic derivative at x=a:

keven=((D[(Subscript[ψ, 2,e])[x],x])/((Subscript[ψ, 2,e])[x]))==((D[(Subscript[ψ, 3,e])[x],x])/((Subscript[ψ, 3,e])[x]))/.x->a

We plot the left and right sides as a function of k -- roots are given by the points of intersection:

List@@keven/.κ->(Sqrt[v-(k^2)])/.params
Plot[%,{k,0,(Sqrt[v])/.params},Exclusions->(Sin[b k]==0/.params),PlotPoints->100]

We find numerical values of the the roots by beginning a search at each zero crossing of the cotangent term:

evenroots=Table[FindRoot[keven/.κ->(Sqrt[v-(k^2)])/.params,{k,k0+(π/(2b)),k0,k0+π/b}/.params],Evaluate[{k0,0,(Sqrt[v]),(π/b)}/.params]]

We find B using continuity of the wave function at x=a:

beven=Solve[(Subscript[ψ, 3,e])[x]==(Subscript[ψ, 2,e])[x]/.x->a,B][[1]]

We find A by normalizing:

1==Integrate[((((Subscript[ψ, 1,e])[x]))^2),{x,-a-b,-a}]+Integrate[((((Subscript[ψ, 2,e])[x]))^2),{x,-a,a}]+Integrate[((((Subscript[ψ, 3,e])[x]))^2),{x,a,a+b}]/.beven
aeven=Solve[%,A][[1]]

Substituting in the expressions for A and B, we find the general form of the eigenfunctions:

(Subscript[ψ, even])[x_]=(Piecewise[{{(Subscript[ψ, 1,e])[x], -a-b≤x<-a}, {(Subscript[ψ, 2,e])[x], -a≤x≤a}, {(Subscript[ψ, 3,e])[x], a<x≤a+b}}])/.beven/.A->Abs[A]/.aeven/.κ->(Sqrt[v-(k^2)])

Each eigenvalue gives a particular eigenfunction:

evenstates=(Subscript[ψ, even])[x]/.params/.evenroots//Chop

Plot the even eigenstates:

Plot[evenstates,Evaluate[{x,-a-b-.5,a+b+.5}/.params],PlotRange->All,Prolog->potential,Axes->{True,False}]

Find the odd-parity eigenstates with E < V

Do the same for the odd-parity eigenstates.

(Subscript[ψ, 1,o])[x_]=A Sin[k(x+a+b)]

(Subscript[ψ, 2,o])[x_]=B Sinh[κ x]

(Subscript[ψ, 3,o])[x_]=A Sin[k(x-a-b)]

kodd=((D[(Subscript[ψ, 2,o])[x],x])/((Subscript[ψ, 2,o])[x]))==((D[(Subscript[ψ, 3,o])[x],x])/((Subscript[ψ, 3,o])[x]))/.x->a

{List@@kodd,List@@keven}/.κ->(Sqrt[v-(k^2)])/.params
Plot[%,{k,0,(Sqrt[v])/.params},Exclusions->Sin[k]==0,PlotPoints->100]

oddroots=Table[FindRoot[kodd/.κ->(Sqrt[v-(k^2)])/.params,{k,k0+(π/(2b)),k0,k0+π/b}/.params],Evaluate[{k0,0,(Sqrt[v]),(π/b)}/.params]]//Chop

bodd=Solve[(Subscript[ψ, 3,o])[x]==(Subscript[ψ, 2,o])[x]/.x->a,B][[1]]

1==Integrate[((((Subscript[ψ, 1,o])[x]))^2),{x,-a-b,-a}]+Integrate[((((Subscript[ψ, 2,o])[x]))^2),{x,-a,a}]+Integrate[((((Subscript[ψ, 3,o])[x]))^2),{x,a,a+b}]/.bodd
aodd=Solve[%,A][[1]]

(Subscript[ψ, odd])[x_]=(Piecewise[{{(Subscript[ψ, 1,o])[x], -a-b≤x<-a}, {(Subscript[ψ, 2,o])[x], -a≤x≤a}, {(Subscript[ψ, 3,o])[x], a<x≤a+b}}])/.bodd/.A->Abs[A]/.aodd/.κ->(Sqrt[v-(k^2)])

oddstates=(Subscript[ψ, odd])[x]/.params/.oddroots//Chop

Plot[oddstates,Evaluate[{x,-a-b-.5,a+b+.5}/.params],PlotRange->All,Prolog->potential,Axes->{True,False}]

Construct a wave packet and observe evolution

Choose this function for a Gaussian centered in the right-hand well.

gfunc=Exp[-(((((x-a-b/2)/σ)))^2)]/.params/.{σ->.2}

Choose this function for the quasi-ground-state of the right-hand well.

gfunc=(Piecewise[{{Sin[(x-a)π/b], a≤x≤b}}])/.params

Expand the wave function in the even and odd parity eigenstates found above, and plot the magnitude of the expansion coefficients:

Off[NIntegrate::"ncvb"]

evencoeff=NIntegrate[gfunc evenstates,Evaluate[{x,-a-b,a+b}/.params]]
ListLinePlot[Abs@%,PlotRange->{0,All}]

oddcoeff=NIntegrate[gfunc oddstates,Evaluate[{x,-a-b,a+b}/.params]]
ListLinePlot[Abs@%,PlotRange->{0,All}]

Form the wavepacket as a superposition of the eigenstates, including the time-dependent phase factors:

packet=Total[evencoeff evenstates Exp[-I t (k^2)/.evenroots]]+Total[oddcoeff oddstates Exp[-I t (k^2)/.oddroots]];

Plot the initial state.

Plot[((Abs[packet])^2)/.t->0,Evaluate[{x,-a-b-.5,a+b+.5}/.params],Filling->Automatic,PlotRange->{0,1.5},Prolog->potential,Axes->{True,False}]

Animate the wavepacket in time. It starts in the right-hand well, but can eventually be found in the left-hand well.

Animate[Plot[((Abs[packet])^2)/.t->time,Evaluate[{x,-a-b-.5,a+b+.5}/.params],Filling->Automatic,PlotRange->{0,1.5},Prolog->potential,Axes->{True,False}],{time,0,Infinity},AnimationRate->.05]