Problems Chapter 11

Quantum Mechanics
                       K. Konishi, G. Paffuti

Problem

Write the semiclassical quantization conditions for the case with two turning points using Cauchy's theorem. Assume that p[x] has only one cut, between the two turning points, and no essential singularities in the complex plane.

Solution

Definitions

Let (Subscript[e, 1])and (Subscript[e, 2]) the two classical turning points, where E = V[x]. In general

p[x] = (Sqrt[2m(E-V[x])])

and the function p[x] has a cut between the zeros of E-V[x]. If this is the only cut

p[x] = R[x](Sqrt[(x-(Subscript[e, 1]))((Subscript[e, 2])-x)])

If there are not essential singularities in R, even at infinity, R[x] is a rational function. This will be assumed from now on.

To extend () to the complex plane we have to define the branch of the square root. We have two possibilities: or the positive determination corresponds to upper edge of the cut or to the lower edge. In the two cases the action integral

I =2(∫_(Subscript[e, 1])^(Subscript[e, 2]))p[x]dx = ∮ p[x]dx

can be expressed as one of the two complex integral below

We choose the first alternative, i.e. p[x] positive on the upper edge of the cut, we leave to the reader the task of verifying that the conclusions below can be obtained also with the other choice.

To understand clearly the conventions used it is convenient to introduce a complex variable z and to put

z-(Subscript[e, 1]) = (Subscript[ρ, 1])(e^(i (Subscript[θ, 1])));       z - (Subscript[e, 2]) = (Subscript[ρ, 2])(e^(i (Subscript[θ, 2])))

The choice of the cut correspond to

(Sqrt[z-(Subscript[e, 1])]) = (Sqrt[(Subscript[ρ, 1])])(e^(i ((Subscript[θ, 1])/2)));       (Sqrt[z - (Subscript[e, 2])]) = (Sqrt[(Subscript[ρ, 2])])(e^(i (((Subscript[θ, 2])/2)))).

and taking p[x] as the analytic continuation of

q[z] =-i R[z](Sqrt[(z-(Subscript[e, 1]))(z - (Subscript[e, 2]))])

coincides with p[x] on the upper edge of the cut, where (Subscript[θ, 1])= 0 and (Subscript[θ, 2])= π. On the lower edge (Subscript[θ, 2])= - π and

q[z]→ R[x] (Sqrt[(Subscript[ρ, 1])])(Sqrt[(Subscript[ρ, 2])])(e^(-iπ)) = - p[x]

It is important to keep track of phases along real axis. For x > (Subscript[e, 2]) we have (Subscript[θ, 1])= 0 and (Subscript[θ, 2])= 0 so

q[z]→ R[x] (Sqrt[(Subscript[ρ, 1])])(Sqrt[(Subscript[ρ, 2])]) (e^(-i (π/2)))= - i R[x] (Sqrt[(Subscript[ρ, 1])])(Sqrt[(Subscript[ρ, 2])])

For x < (Subscript[e, 2]), (Subscript[θ, 1])= π, (Subscript[θ, 2])= π or  (Subscript[θ, 1])= -π, (Subscript[θ, 2])=- π in both cases

q[z]→ R[x] (Sqrt[(Subscript[ρ, 1])])(Sqrt[(Subscript[ρ, 2])]) (e^(∓ i π-i (π/2)))= + i R[x] (Sqrt[(Subscript[ρ, 1])])(Sqrt[(Subscript[ρ, 2])])

The fact that the two determinations coincide means in fact that the function has no cut for x < (Subscript[e, 1]).

From a more simple point of view we are taking

(Subscript[e, 2])- z = (e^(- i π))(z-(Subscript[e, 2])).

Cauchy's theorem

Apart from the cut the integrand can have singularities (poles) due to R[z] or at infinity. We can exclude these singularities by drawing contours in the complex plane ℂ. With the convention used for the cut contour wrap clockwise the singularities and the situation is shown below

The function is analytic in the domain excluded by the contours then Cauchy theorem states

I  - 2π i (∑_a)Res[q[z]] - 2π i (Subscript[Res, ∞])[q[z]] = 0

Res stands for the residue, the minus sign is due to the clockwise direction of the integrals.

Examples

R[x] = 1

The simplest example is

p[x] = (Sqrt[(x-(Subscript[e, 1]))((Subscript[e, 2])-x)])

In this case the integral can be trivially done with the substitution

x → (((Subscript[e, 1])+(Subscript[e, 2]))/2)+ ξ (((Subscript[e, 2])-(Subscript[e, 1]))/2);    2(∫_(Subscript[e, 1])^(Subscript[e, 2]))p[x]dx = (((((Subscript[e, 2])-(Subscript[e, 1])))^2)/2)(∫_(-1)^(+1))(Sqrt[1-(ξ^2)])dξ = (π/4)((((Subscript[e, 2])-(Subscript[e, 1])))^2)

but let us check the result with the above general method.

The only singular point is at infinity. To compute the residue at infinity of q[z] let us make the change of variables z = 1/ζ, dz = - dζ/(ζ^2)

q[z] dz =((- dζ)/(ζ^2)) (-i)(Sqrt[((1/ζ)-(Subscript[e, 1]))((1/ζ)- (Subscript[e, 2]))])=  i(( dζ)/(ζ^3))(Sqrt[(1- ζ (Subscript[e, 1]))(1- ζ (Subscript[e, 2]))])

By definition the residue is the coefficient of ζ in the Taylor expansion of the function. Using

(Sqrt[(1- ζ (Subscript[e, 1]))(1- ζ (Subscript[e, 2]))]) ∼ 1 - (1/2)((Subscript[e, 1])+ (Subscript[e, 2]))ζ- (1/8)((((Subscript[e, 2])- (Subscript[e, 1])))^2) (ζ^2) + 𝒪(ζ)

we have

(Subscript[Res, ∞])q[z] = -(i/8)((((Subscript[e, 2])- (Subscript[e, 1])))^2)

and from ()

I =  ∮ p[x]dx = 2(∫_(Subscript[e, 1])^(Subscript[e, 2]))p[x]dx = (π/4)((((Subscript[e, 2])- (Subscript[e, 1])))^2)

Note for the Mathematica users.

In Mathematica 6.0 this integral can be computed with

Simplify[Integrate[2(Sqrt[(x-a)(b-x)]),{x,a,b},GenerateConditions-> False],{0<a<b}]

In Mathematica 5.2 it is sufficient to write

Integrate[2(Sqrt[(x-a)(b-x)]),{x,a,b},Assumptions-> {0<a<b}]

R[x] = 1/x

p[x] = (1/x)(Sqrt[(x-(Subscript[e, 1]))((Subscript[e, 2])-x)])

In this case we have a singularity at infinity and at x = 0. We are supposing 0 < (Subscript[e, 1])< (Subscript[e, 2]).

At infinity, working as in the previous example

q[z] dz =((- dζ)/(ζ^2)) (-i)ζ(Sqrt[((1/ζ)-(Subscript[e, 1]))((1/ζ)- (Subscript[e, 2]))])=  i(( dζ)/(ζ^2))(Sqrt[(1- ζ (Subscript[e, 1]))(1- ζ (Subscript[e, 2]))])

and using the expansion ()

(Subscript[Res, ∞])q[z] = -(i/2)((Subscript[e, 1])+ (Subscript[e, 1]))

For z→ 0 we have to remember the choice () for phases and write

q[z] → + (i/z)(Sqrt[(Subscript[e, 1])(Subscript[e, 2])])   ⇒ (Subscript[Res, 0])q[z] = i (Sqrt[(Subscript[e, 1])(Subscript[e, 2])])

Substituting in  ()

I =  ∮ p[x]dx = 2(∫_(Subscript[e, 1])^(Subscript[e, 2]))p[x]dx = π((((Sqrt[(Subscript[e, 2])])- (Sqrt[(Subscript[e, 1])])))^2)

Let us note that correctly the integral vanish for  (Subscript[e, 1])= (Subscript[e, 2]).

Note for the Mathematica users.

In Mathematica 6.0 this integral can be computed with

Integrate[(2/x)(Sqrt[(x-a)(b-x)]),{x,a,b},Assumptions-> {0<a<b}]

In Mathematica 5.2 the computation gives the wrong result  π((((Sqrt[(Subscript[e, 2])])+ (Sqrt[(Subscript[e, 1])])))^2).

Problem

Write the semiclassical quantization conditions for a harmonic oscillator. Generalize the results for a pure (x^N)oscillator, with N even.

Solution

Quantization

The Hamiltonian is given by

H[p,q] = (1/(2m))(p^2)+ (1/2)m (ω^2)(q^2) .

At fixed energy E

p[q] = ± (Sqrt[2m E - (m^2) (ω^2)(q^2)]).

We have an oscillatory motion, with two classical inversion points

(Subscript[q, 1])= -A;  (Subscript[q, 2])= A;       A = (Sqrt[((2E)/(m (ω^2)))]) .

The action integral is easily computed as

J = (1/(2π))∮p[q]dq = (1/π)(∫_(Subscript[q, 1])^(Subscript[q, 2]))(Sqrt[2m E - (m^2) (ω^2)(q^2)]) dq = ((2E)/(π ω))(∫_(-1)^1)(Sqrt[1-(x^2)])dx = (E/ω) .

The quantization conditions are

J = ℏ( n + (1/2)).

The relation between E and J can  be trivially inverted in this case and we have for the spectrum

E = ω J =  ℏ ω( n + (1/2)).

Which is the exact result.

Angle action variables

The angle variable φ is defined as the canonically conjugate variable of J. From () we have

H[J] = ω J.

The variable φ is cyclic and its equation of motion is trivially solved:

((dφ)/dt) = ((∂H)/(∂J))=ω ;   ⇒  φ[t] = ω t + δ.

From the known solution of the equation of motion for p and q:

q= A Sin[ω t + δ]; p = m ω A Cos[ω t + δ]

or by using Hamilton Jacobi equation, one has

q = (Sqrt[((2J)/(m ω))])Sin[φ];   p = (Sqrt[2 m ω J]) Cos[φ].

Which is the explicit canonical transformation between the couple (q, p) and (φ, J).

Potential (x^N)

Previous considerations can be easily extended to Hamiltonians of the form

H[p,q] = (1/(2m))(p^2)+ (1/2)g (q^N) .

With N even.

At fixed energy

p[q] = ± (Sqrt[2m E - m g (q^N)])

Turning points are

(Subscript[q, 1,2])= ∓ a;     a = (((((2E)/g)))^(1/N)).

The action integral is given by

J = (1/(2π))∮p[q]dq = (1/π)(∫_(-a)^a)(Sqrt[2m E - m g (x^N)]) dx = (2/π)(Sqrt[2m E ])(∫_0^a)(Sqrt[1- ((((x/a)))^N)]) dx =C (a/π)(Sqrt[2m E ]) .

The constant C is given by

C=2(∫_0^1)(Sqrt[1- (z^N)]) dz = (((Sqrt[π]) Gamma[1+(1/N)])/( Gamma[(3/2)+(1/N)]))

The quantization condition is always

J = ℏ( n + (1/2)).

We can invert the relation () to get

E = D  (J^((2N)/(N+2)));         1/D = ((((C/π)(Sqrt[2m])((((2/g)))^(1/N))))^((2N)/(N+2)))= ((((C/π)(Sqrt[2m])))^((2N)/(N+2)))((((2/g)))^(2/(N+2)))

This relation has some interesting implications.

The transformation from (q,p) to angle - action variables is a canonical transformation.  If we would start our quantum theory with an Hamiltonian like () and proceed naively we go rapidly into troubles. The commutation relations

[J,φ]= (ℏ/i)

fix he spectrum of J, as φ is an angular variable. As it is well known from the theory of angular momentum the eigenvalues of J are

J =ℏ n;    n = 0,±1,±2...

For a cyclic Hamiltonian the eigenvalues of J fix the spectrum of H, in the present case

E = D (((ℏ n))^((2N)/(N+2)))

This is the correct semiclassical spectrum, i.e. for large n, as it is apparent from () but it is not the true quantum spectrum. This is another evidence that canonical classical transformation are not an invariance of Quantum Theory: cartesian coordinates are "special" and von Neumann representation theorem assure the uniqueness of the Schrödinger representation up to unitary transformations in that case. Vice versa as the spectrum () is not the quantum spectrum this means that canonical transformations are not realized by a unitary transformation.

The semiclassical spectrum () predicts a non analytic behavior in g

E ∼ (g^(2/(N+2)))

This is indeed true in the quantum spectrum.

Problem

Write the semiclassical quantization conditions for a central field. In particular find the hydrogen spectrum using the Bohr Sommerfeld quantization conditions.

Solution

The Hamiltonian

Let V[r] the central potential. To correctly identify canonical momenta in radial problems it is convenient to start from the Lagrangian of the system

ℒ = (μ/2)(v^2) - V[r]

Writing the line element in radial coordinates

(ds^2) = (dr^2) + (r^2)((dθ)^2)+ (r^2)((Sin[θ])^2) ((dφ)^2)

one has directly the velocity in radial coordinates

(v^2)= ((OverDot[r])^2) + (r^2)((OverDot[θ])^2) + (r^2)((Sin[θ])^2)((OverDot[φ])^2);   ℒ = (μ/2)( ((OverDot[r])^2) + (r^2)((OverDot[θ])^2) + (r^2)((Sin[θ])^2)((OverDot[φ])^2)) - V[r]

θ and φ are polar and azimuthal angles respectively. From () one has the canonical momenta

(Subscript[p, r]) = ((∂ℒ)/(∂(OverDot[r])))= μ (OverDot[r]) ;  (Subscript[p, θ])= ((∂ℒ)/(∂(OverDot[θ]))) = μ (r^2)(OverDot[θ]) ; (Subscript[p, φ])= ((∂ℒ)/(∂(OverDot[φ]))) = μ (r^2)((Sin[θ])^2)(OverDot[φ])

Transforming variables from velocities to momenta and substituting in the general definition of Hamiltonian

H = p (OverDot[q]) - ℒ

one gets

H = (1/(2μ))((Subscript[p, r]^2) +(1/(r^2))(Subscript[p, θ]^2)  +(1/((r^2)((Sin[θ])^2)))(Subscript[p, φ]^2)) + V[r].

φ is a cyclic variable, then (Subscript[p, φ])is a constant of motion.

Hamilton Jacobi equation and J variables

The Hamilton Jacobi equation for the Hamiltonian () is

(((((∂S)/(∂r))))^2)+ (1/(r^2))(((((∂S)/(∂θ))))^2) +(1/((r^2)((Sin[θ])^2)))(((((∂S)/(∂φ))))^2)+2μ (V[r] - E) = 0.

This equation can be separated, i.e. admits a solution in the form

S = (Subscript[S, r])[r] + (Subscript[S, θ])[θ] + (Subscript[S, φ])[φ].

By substitution in () we see first of all that φ can appear only in the term ∂S/∂φ, then this term must be constant, let call it (Subscript[A, φ]). The θ dependence is then confined in

(((((Subscript[dS, θ])/(dθ))))^2) +(1/((Sin[θ])^2))((Subscript[A, φ])^2)

which must be also constant, call it (Subscript[A, θ]^2) . Finally for (Subscript[S, r])we have an ordinary differential equation. Summarizing we have:

((Subscript[dS, φ])/(dφ))= (Subscript[A, φ]);       (((((Subscript[dS, θ])/(dθ))))^2) +(((Subscript[A, φ])^2)/((Sin[θ])^2))= (Subscript[A, θ]^2) ;      (((((Subscript[dS, r])/dr)))^2)+ ((Subscript[A, θ]^2)/(r^2)) + 2μ (V[r]-E) = 0.

The problem has been separated and we have three integration constants, (Subscript[A, φ]), (Subscript[A, θ]), and E . The meaning of the first two constant is evident if we write them explicitly, they are the z - component of angular momentum and the modulus of angular momentum:

(Subscript[A, φ])= (Subscript[p, φ])= μ (r^2)((Sin[θ])^2) (OverDot[φ]) = (Subscript[L, z]) ;

(Subscript[A, θ])= (Sqrt[(Subscript[p, θ]^2)+ ((Subscript[p, φ]^2)/((Sin[θ])^2))])= μ r (Sqrt[(((r (OverDot[θ])))^2)+ (((r Sin[θ] (OverDot[φ])))^2)])= μ |r⋀ v| = L .

We use this notation from now on.

The action integrals are

(Subscript[J, φ])= (1/(2π))∮ (Subscript[p, φ])dφ = (Subscript[L, z]);
(Subscript[J, θ])= (1/(2π))∮ (Subscript[p, θ])dθ =(1/(2π))∮ ((Subscript[dS, θ])/(dθ))dθ = (1/(2π))∮d θ(Sqrt[(L^2)- ((Subscript[L, z]^2)/((Sin[θ])^2))]) = (L - |(Subscript[L, z])|);
(Subscript[J, r])= (1/(2π))∮ (Subscript[p, r])dr =(1/(2π))∮ ((Subscript[dS, r])/dr)dr= (1/(2π))∮dr(Sqrt[2μ(E-U)- ((L^2)/(r^2))])
        = (1/(2π))∮dr(Sqrt[2μ(E-U)- (((((Subscript[J, θ])+|(Subscript[J, φ])|))^2)/(r^2))])

WKB quantization and hydrogen atom

(Subscript[L, z])

The variable φ has no turning points then the quantization rule is simply

(Subscript[J, φ])=(Subscript[L, z])= m ℏ

L

The variable θ has two turning points, given by the solution of

Sin[θ] (L^2) - (Subscript[L, z]^2) = 0

This equation imply the obvious physical condition L ≥ (Subscript[L, z]).  The quantization condition for (Subscript[J, θ]) are

(Subscript[J, θ]) =L -|(Subscript[L, z])| =  ℏ (k'+ (1/2)) ⇒ L = ℏ(k'+ |m| + (1/2));   k'=0,1,...

The additional term 1/2 is peculiar of WKB quantization, it would be absent in old Bohr - Sommerfeld quantization.

(Subscript[J, r])

The quantization rule for (Subscript[J, r]) finally fix the spectrum. As there are in general two inversion points:

(Subscript[J, r]) = ℏ ((Subscript[n, r])+(1/2))

Let us consider in particular the hydrogen spectrum, with

V[r] = - (( (e^2))/r).

The radial action integral is

(Subscript[J, r]) = (1/π)(∫_(Subscript[r, 1])^(Subscript[r, 2]))(Sqrt[2μ(E+((e^2)/r))- ((L^2)/(r^2))])

There are two inversion points for E < 0, ϵ = |E|:

(1/(Subscript[r, 1]))= (((e^2)μ+(Sqrt[(e^4)(μ^2)- 2 (L^2)ϵ μ]))/(L^2));  (1/(Subscript[r, 2]))= (((e^2)μ-(Sqrt[(e^4)(μ^2)- 2 (L^2)ϵ μ]))/(L^2)).

We have

(Sqrt[2μ(E+((e^2)/r))- ((L^2)/(r^2))])= (Sqrt[2ϵ μ])(1/r)(Sqrt[(r-(Subscript[r, 1]))((Subscript[r, 2])-r)]);      (Subscript[J, r]) = ((Sqrt[2ϵ μ])/π)(∫_(Subscript[r, 1])^(Subscript[r, 2]))(Sqrt[(r-(Subscript[r, 1]))((Subscript[r, 2])-r)]) d r;

This type of integral has already been done in problem [1], we just report the result:

(Subscript[J, r])=(π/2)((Sqrt[2ϵ μ])/π) ((((Sqrt[(Subscript[r, 2])])-(Sqrt[(Subscript[r, 1])])))^2) = - L + (e^2)(Sqrt[(μ/(2ϵ))])

Then the WKB spectrum is given by

2 ϵ = - 2 E =  (((ℏ^2) (e^4)μ)/((((Subscript[n, r])+ (1/2)+ (k'+ |m|+(1/2))))^2))= (((ℏ^2) (e^4)μ)/((((Subscript[n, r])+ k'+ |m|+1)))^2))

with

n = (Subscript[n, r])+ k'+ |m| = 1,2,...

We see that exact hydrogen spectrum is reproduced.

Let us note that in old quantum theory both terms 1/2 in L and (Subscript[J, r]) where missing and, quite arbitrarily, (Subscript[n, r]) was taken as ≥ 1, instead of ≥ 0.

Problem

Solve  the harmonic oscillator problem in the lowest order WKB approximation.

Solution

Eigenvalues

The Schrödinger equation is

-((ℏ^2)/(2m))((d^2)/(dx^2))ψ +(1/2)m (ω^2) (x^2)ψ = E ψ;    V[x] = (1/2)m (ω^2) (x^2).

The classical momentum at fixed energy is given by

p[x] = (Sqrt[2m(E -V[x])]).

The classical turning points are given by the solutions of the equation E = V[x]. In this case

x = ± a = ± (Sqrt[((2 E)/(m (ω^2)))]).

In terms of a we can rewrite p[x] as

p[x] = m ω (Sqrt[(a^2)- (x^2)]).

The quantization condition states :

(1/(2π ℏ))∮p[x] 𝕕x = (1/(π ℏ))(∫_(-a)^(+a))p[x] 𝕕x = n + (1/2);    n = 0,1,...

Inserting eq.()

n + (1/2) = ((m ω)/(π ℏ))(∫_(-a)^(+a))  (Sqrt[(a^2)- (x^2)]) 𝕕x = ((m ω)/(π ℏ))(a^2)(∫_(-1)^(+1))  (Sqrt[1- (z^2)]) 𝕕z = ((m ω)/(2 ℏ))(a^2) = (E/(ℏ ω)).

The levels are given by

(Subscript[E, n]) = ℏ ω (n+ (1/2))

Wave functions

In the text has been shown that the n-th bound state wave function in the lowest order WKB approximation is

ψ[x] = (C/(2(Sqrt[|p[x]|])))Exp[-| σ[(Subscript[x, 1]),x]|] ;  x < (Subscript[x, 1]);
ψ[x] = (C/(Sqrt[|p[x]|])) Cos[ σ[(Subscript[x, 1]),x]- (π/4)];  (Subscript[x, 1])< x < (Subscript[x, 2]);
ψ[x] = (C/(2(Sqrt[|p[x]|])))Exp[-| σ[(Subscript[x, 2]),x]|] ;  (Subscript[x, 2]) < x;

Where (Subscript[x, 1])and (Subscript[x, 2])are the classical turning points and σ is the classical action in units of ℏ, i.e.

σ[(Subscript[x, 1]),x]= (1/ℏ)(∫_(Subscript[x, 1])^x)p[x] 𝕕x

The consistency of the definition rely on the quantization condition ():

σ[(Subscript[x, 1]),(Subscript[x, 2])]= (1/ℏ)(∫_(Subscript[x, 1])^(Subscript[x, 2]))p[x] 𝕕x = π (n+(1/2))

C is a normalization constant, which can be fixed by neglecting exponential terms and using a mean value for the Cos - part:

1 ∼ (C^2)(∫_(Subscript[x, 1])^(Subscript[x, 2]))((|ψ|)^2)𝕕x = (C^2)(∫_(Subscript[x, 1])^(Subscript[x, 2]))((𝕕x)/(p[x]))((Cos[σ-(π/4)])^2)∼ ((C^2)/2)(∫_(Subscript[x, 1])^(Subscript[x, 2]))((𝕕x)/(m v[x]))= ((C^2)/(2 m))(T/2).

C = 2 (Sqrt[(m/T)])

where T is the classical period of motion.

In the present case (Subscript[x, 1])= - a;  (Subscript[x, 2])= + a and T = 2π/ω, i.e.

C = (Sqrt[((2m ω)/π)]).

It is sufficient to compute the wave function for x > 0, as n = even-odd give rise to even-odd functions.

0 < x < a

We have, using () and  m (ω^2)(a^2)/2 = E :

p[x] =  m ω (Sqrt[(a^2)- (x^2)]);

σ[-a,x] = σ[-a,0]+ σ[0,x]= (1/2) σ[-a,a]+σ[0,x] =
 (π/2) (n+(1/2))+(1/ℏ)(∫_0^x)m ω (Sqrt[(a^2)- (x^2)])𝕕x = (π/2) (n+(1/2))+ ((2E)/(ℏ ω))(∫_0^(x/a)) (Sqrt[1- (z^2)])𝕕z =
 (((π n)/2)+(π/4))+((2E)/(ℏ ω))((x/(2a)) (Sqrt[1-((x^2)/(a^2))])+(1/2) ArcSin[(x/a)])

and (with E = ℏ ω (n+1/2) )

ψ[x] = (C/(Sqrt[|p[x]|]))Cos[((π n)/2)+(E/(ℏ ω))((x/a) (Sqrt[1-((x^2)/(a^2))])+ ArcSin[(x/a)])

a < x

In this region

p[x] =  m ω (Sqrt[(x^2)- (a^2)]);

σ[a,x] = ((2E)/(ℏ ω))(∫_1^(x/a)) (Sqrt[(z^2)-1])𝕕z = ((2E)/(ℏ ω))(1/2)((x/a) (Sqrt[((x^2)/(a^2))-1])- ArcCosh[(x/a)])

and

ψ[x] = (1/2)(C/(Sqrt[|p[x]|]))Exp[-(E/(ℏ ω))((x/a) (Sqrt[((x^2)/(a^2))-1])- ArcCosh[(x/a)])]

Asymptotically

σ[a,x](Underscript[→, x→∞]) (E/(ℏ ω))(((x^2)/( (a^2)))- Log[x])=(( m ω)/ℏ) ((x^2)/2)- (n+(1/2))Log[x]

and

ψ[x]∼(1/(Sqrt[x])) Exp[- ((m ω)/ℏ) ((x^2)/2)](x^(n+(1/2)))= (x^n)Exp[- ((m ω)/ℏ) ((x^2)/2)]

Which is the correct asymptotic behavior for  (Subscript[H, n])[x]Exp[- ((m ω)/ℏ) ((x^2)/2)].

Here is graph showing the approximation of the wave function for n=0 and n=10. We see the divergence at the turning points.

Problem

Compute the spectrum for a potential (x^N)in the  lowest order WKB approximation. N is even. Study the limit of large quantum numbers.

Solution

Change of variables

The Schrödinger equation is

-((ℏ^2)/(2m))((d^2)/(dx^2))ψ +(1/2)g (x^N)ψ = E ψ;    V[x] = (1/2)g (x^N).

We change variables x = λ z in such a way that the coefficients in the left hand side of the equation become equal. A change of variable do not change the spectrum.

-((ℏ^2)/(2m (λ^2)))((d^2)/(dz^2))ψ +(1/2)g (λ^N) (z^N)ψ = E ψ;

The coefficients are equal for

((ℏ^2)/(m (λ^2)))= g (λ^N)  ⇒ λ = (((((ℏ^2)/(m g))))^(1/(N+2))).

λ has the dimension of a length, as the reader can easily check.

With this choice the Schrödinger equation can be written as

-(1/2)((d^2)/(dz^2))ψ +(1/2) (z^N)ψ = E (m/(ℏ^2))(λ^2) ψ;

With

E(m/(ℏ^2))(λ^2) = (ϵ/2)

we have the universal form

-((d^2)/(dz^2))ψ + (z^N)ψ = ϵ ψ

Energies can be recovered using (). The energy scale is given by

(Subscript[ℰ, 0]) = ((ℏ^2)/(2 m (λ^2)))

Let us note that the rescaling is a canonical transformation then action remain invariant. In particular

((p[x])/ℏ)=(1/ℏ)(Sqrt[2m(E- (1/2)g (x^N))])= (1/ℏ)(Sqrt[2m(ϵ((ℏ^2)/(2m (λ^2)))- (1/2)g (λ^N) (z^N))])=
  =(1/ℏ)(Sqrt[2m])(ℏ/(λ(Sqrt[2m])))(Sqrt[ϵ- (z^N)]) = (1/λ)(Sqrt[ϵ- (z^N)])

then

∫((p[x])/ℏ)𝕕x = ∫k[z] 𝕕z

where k[z] is the wave number in eq.(),  k[z] = (Sqrt[ϵ- (z^N)]) .

WKB quantization

The turning points in () are (with N even)

z = ± a = ± (ϵ^(1/N))

The "wave number" is

k[z] = (Sqrt[ϵ - (z^N)])

and the quantization conditions

(1/(2π))∮k[z] 𝕕z = (n+(1/2))   ⇒   (∫_(-a)^(+a))k[z] 𝕕z = π (n+(1/2))

The integral appearing in () can be expressed in terms of Euler Beta function. With N even:

(∫_(-a)^(+a))k[z] 𝕕z =(∫_(-a)^(+a))(Sqrt[(a^N) - (z^N)]) 𝕕z =2 (a^(1+(N/2)))(∫_0^(+1))(Sqrt[1 - (y^N)]) 𝕕y = 2 (a^(1+(N/2)))(((Sqrt[π]) Γ[1+(1/N)])/(2 Γ[(3/2)+(1/N)])).

We have from ()

a = ((((n+(1/2)) (Sqrt[π])((Γ[(3/2)+(1/N)])/(Γ[1+(1/N)]))))^(2/(2+N)));   (Subscript[ϵ, n]) = ((((n+(1/2)) (Sqrt[π])((Γ[(3/2)+(1/N)])/(Γ[1+(1/N)]))))^((2N)/(2+N))).

In particular for N = 4

(Subscript[ϵ, n]) = (((n+(1/2)))^(4/3))(π^(2/3)) ((((( Γ[(7/4)])/(Γ[(5/4)]))))^(4/3))

In the text the same expression is written in a slightly different but equivalent  form,

(((n+(1/2)))^(4/3))((((π/(2I))))^(4/3))= (((n+(1/2)))^(4/3))((((π/2)8 ((Γ[(7/4)])/((Sqrt[π])Γ[(1/4)]))))^(4/3)) = (Subscript[ϵ, n])

The last equality follows from x Γ[x] = Γ[1+x], for x = 1/4.

Usual energy is

(Subscript[E, n])= (1/2)((ℏ^2)/m)(1/(λ^2))(Subscript[ϵ, n])=(1/2)((ℏ^2)/m)(((((m g)/(ℏ^2))))^(2/(N+2)))(Subscript[ϵ, n]) ≡ (Subscript[ℰ, 0])(Subscript[ϵ, n])

(Subscript[ℰ, 0]) is the characteristic energy of the problem.

Here we present a comparison between numerical computed  (Subscript[ϵ, n]) and WKB results for N = 4 and N = 6 (see notebook [NB-11.1])

({{, (Subscript[ϵ, n]), (Subscript[ϵ, n]^WKB), δϵ/ϵ}, {0, 1.0603620904688216`, 0.867145326484821`, 0.18221772139984083`}, {1, 3.799673029803774`, 3.7519199235504312`, 0.012567688292855251`}, {2, 7.455697938006631`, 7.413988252810775`, 0.0055943367800933294`}, {3, 11.64474551138099`, 11.611525345197105`, 0.002852803107755059`}, {4, 16.261826018869712`, 16.233614692705245`, 0.0017348190868437048`}, {5, 21.23837291823715`, 21.213653359057176`, 0.0011639102145507615`}})                      ({{, (Subscript[ϵ, n]), (Subscript[ϵ, n]^WKB), δϵ/ϵ}, {0, 1.1448024556389298`, 0.8008300276096234`, 0.30046443937555206`}, {1, 4.338598710984494`, 4.161234888139964`, 0.040880439667185416`}, {2, 9.07308456186628`, 8.953551900790758`, 0.013174423787243289`}, {3, 14.935169635771647`, 14.831579668416374`, 0.006935975277251719`}, {4, 21.714165422536293`, 21.62241074545983`, 0.0042255677476433755`}, {5, 29.29964593762869`, 29.216579946741845`, 0.0028350510126869764`}})

Classical parameters and semiclassical limit

The classical motion is fixed by the turning points. It is a periodic motion with period

T = 2 (∫_(-a)^(+a)) ((𝕕x)/(v[x]))= 2 m(∫_(-a)^(+a)) ((𝕕x)/(p[x]))= 2 m ((λ^2)/ℏ)(∫_(-a)^(+a)) ((𝕕z)/(k[z]))=
 4 ((m (λ^2))/ℏ)(a^(1-(N/2)))(∫_0^1)((𝕕y)/(Sqrt[1-(y^N)]))=  4 ((m (λ^2))/ℏ)(a^(1-(N/2))) (((Sqrt[π]) Γ[1+(1/N)])/(Γ[(1/2)+(1/N)]))

By expressing a as a function of E, with E = (Subscript[ℰ, 0]) (a^N) we have

T = 2 (ℏ/(Subscript[ℰ, 0]))(((((Subscript[E, n])/(Subscript[ℰ, 0]))))^((2-N)/(2N)))(((Sqrt[π]) Γ[1+(1/N)])/(Γ[(1/2)+(1/N)]))

At the semiclassical level the following relation must hold :

((∂(Subscript[E, n]))/(∂n))= ((2π)/T) ℏ

This is a check of () for the case N =4, g = 1.2. The continuous line is the left hand side of equation () computed with an interpolation, the points are classical expression for 2 π/T.

Problem

Consider a particle of mass m in a one dimensional potential V[x] = g δ(x-a) + g δ(x+a). Compute in the semiclassical approximation energy and width of the metastable states. Compare the results with the solution of the Schrödinger equation with the Gamow Siegert boundary conditions.

Solve the same problem for a "radial" barrier, i.e. defined in the semispace x > 0.

Solution

One dimensional problem

The potential

V[x] = g δ(x-a)+g δ(x+a)

is the limit case of the one shown in the figure below

Metastable states can occur, which decay through the δ - barriers. The Hamiltonian of the system is

H = ((p^2)/(2m)) + g δ(x-a)+g δ(x+a)

Using

E =(1/(2m)) (ℏ^2)(k^2);  β =(( m g)/(ℏ^2))

stationary states are given by the solution of the Schrödinger equation

ψ''[x] + (k^2)ψ[x] = 2β (δ(x-a)+ δ(x+a))ψ[x] = 2β (δ(x-a)ψ[a]+ δ(x+a)ψ[-a]).

ψ[x] is continuous at x= ± a while its derivative has a discontinuity fixed by equation (). Integrating in a small interval ϵ around a singularity and taking the limit ϵ→0 we have the jumps:

(Δψ)'[a] = 2β ψ[a];   (Δψ)'[-a] = 2β ψ[-a];            where  (Δψ)'[z] = (Underscript[Lim, x→(SuperPlus[z])]) ψ'[x] - (Underscript[Lim, x→(SuperMinus[z])]) ψ'[x]

WKB approximation

We have a well of width 2 a. The semiclassical quantization conditions are the same as those of a free particle in a well

(1/(2π ℏ))∮p[x]dx = (1/(π ℏ))(∫_(-a)^(+a))pdx = ((2a)/(π ℏ)) = n;  n = 1,2,...;       
 p = ℏ (π/(2a))n;      k = (p/ℏ) = ((n π)/(2a)) .

Note. In general the quantization conditions are of the form

(1/(2π ℏ))∮p[x]dx = (Subscript[n, 0])+(∑_(t.p.)) (μ/4)

The sum is on he turning points and μ = 1 for turning points with continuous derivative, while μ=2 for turning points with discontinuous derivative, as in this case. As (Subscript[n, 0])=0,1,... we have (Subscript[n, 0])+ 2 μ/4 = (Subscript[n, 0])+1 = 1,2,3....

The corresponding energies are

(Subscript[E, n])= ((ℏ^2)/(2m))(((n^2)(π^2))/(4 (a^2))).

The width of these states is given by

Γ = 2(ℏ/T)P ,

where P is the probability of crossing a δ - barrier and

T =  ((4a )/v)= ((4a m)/p) = ((8(a^2)m)/(n π ℏ)),

is the classical period of motion. P can be estimated using tunnel effect. Consider for example the right barrier. The approximate solutions for ψ are

x<a: (Subscript[ψ, L])= (e^(i k x)) + ℛ (e^(-i k x)) ;   x>a: (Subscript[ψ, R])= 𝒯 (e^(i k x)).

Continuity of ψ and the constraint () give

(Subscript[ψ, L])[a] = (Subscript[ψ, R])[a];  (Subscript[ψ, R]^′)[a]-(Subscript[ψ, L]^′)[a]= 2β (Subscript[ψ, R])[a].

Inserting () and solving for 𝒯 one easily find

𝒯 =(( i k)/(i k - β));         P = ((|𝒯|)^2)= ((k^2)/((k^2)+(β^2)))

In this approximation we are neglecting multiple reflections, i.e. we assume P ≪ 1, or β ≫ k. In this approximation

P ≃ ((k^2)/(β^2))

Inserting this value and the expression for T in ()

Γ = ((ℏ^2)/(2m))(((π^3) (n^3))/(8 (a^4)(β^2))).

This formula is no more valid for β ∼ k.

Gamow Siegert boundary conditions

In one dimension a divergent wave means

ψ ∼ Exp[i k |x|]

In fact, the phase grows for both limits x → ± ∞. This will be the form of the solution for |x| > a, as the potential is zero.

For -a < x < a there two linearly independent solutions, Cos[kx], Sin[kx], respectively for even and odd states. Then the solution will be:

Even states

x<-a: ψ= A (e^(i k|x|));   -a<x<a:   Cos[kx];    a< x: A (e^(i k|x|));

Odd states

x<-a: ψ=- A (e^(i k|x|));   -a<x<a:   Sin[kx];    a< x: A (e^(i k|x|));

We can limit ourselves to the point x=a and impose there the constraint () :

Even states

Cos[k a]= A (e^(i k a));    i k A (e^(i k a)) + k Sin[k a] = 2β Cos[ka];

Odd states

Sin[k a]= A (e^(i k a));    i k A (e^(i k a)) - k Cos[k a] = 2β Sin[ka];

Eliminating A

({{(ik-2β)Cos[ka] + k Sin[ka] = 0,      even states}, {(ik-2β)Sin[ka] - k Cos[ka] = 0,   odd states}})  .

The solutions of these equation are  complex k, which give rise to a complex energy

E = ((ℏ^2)/(2m))(k^2)= (Subscript[E, R])- i (Γ/2) .

Large β limit

To have a finite solution as β grows ka must be approximately an odd/even multiple of π/2 for even/odd states, to compensate the growth in β. Let us pose

even states: k =(π/(2a))(2s+1) + ((Subscript[c, 1])/((a^2)β)) + ((Subscript[c, 2])/((a^3)(β^2))) + ...
odd states:   k = (π/(2a))(2s) + ((Subscript[d, 1])/((a^2)β)) + ((Subscript[d, 2])/((a^3)(β^2))) + ...

Inserting n the previous equation and expanding in 1/β one find

(Subscript[c, 1])=- (( π)/4)(2s+1);  (Subscript[c, 2]) = (π/8)(2s+1) - i ((π^2)/16)(((2s+1))^2);
(Subscript[d, 1])= s (π/2)= (π/4)(2s);  (Subscript[d, 2])= s (π/4) - i (s^2)((π^2)/4)= (π/8)(2s) - i ((π^2)/16)(((2s))^2)

Note : As n starts from 1 even states are n=1,3,5 ... while for odd states n = 2,4,..

Both expression can be written as a correction of the form

k ≃ (π/(2a))n + (((-1))^n)(π/4)n (1/((a^2)β)) + ((π/8)n - i ((π^2)/16)(n^2))(1/((a^3)(β^2))) + 𝒪((β^(-3)))

For the energy one deduces :

E = ((ℏ^2)/(2m))(k^2)≃((ℏ^2)/(2m))( (((n^2)(π^2))/(4 (a^2)))+ (((-1))^n)(((n^2)(π^2))/(4(a^3)β)) + ((3/16)(((n^2)(π^2))/((a^4)(β^2))) - i (1/16)(((n^3)(π^3))/((a^4)(β^2))))) +  𝒪((β^(-3)))

At lowest order

E =(Subscript[E, R])- i (Γ/2) ≃((ℏ^2)/(2m))(((n^2)(π^2))/(4 (a^2))) - i ((ℏ^2)/(2m))(1/16)(((n^3)(π^3))/((a^4)(β^2)))

in agreement with ().

The radial problem

In this case the space is bounded to x > 0 and the effective Schrödinger equation is

ψ''[x] + (k^2)ψ[x] = 2β δ(x-a)ψ[x] = 2β δ(x-a)ψ[a].

It can represent the wave equation for the S - wave reduced function in a radial problem.

WKB

For energies the only difference with respect previous problem is that the well has width a instead of 2 a, then energies are

(Subscript[E, n])= ((ℏ^2)/(2m))(((n^2)(π^2))/( (a^2))).

For the width Γ we have now

Γ = (ℏ/T)P

as only one barrier is present. The period is given by

T =  ((2a )/v)= ((2a m)/p) = ((2(a^2)m)/(n π ℏ))  .

The computation the probability of tunneling P is the same as before:

𝒯 =(( i k)/(i k - β));         P = ((|𝒯|)^2)= ((k^2)/((k^2)+(β^2))) ∼ ((k^2)/(β^2)).

The resulting WKB approximation for Γ is

Γ = ((ℏ^2)/(2m))(((π^3) (n^3))/( (a^4)(β^2))).

Gamow Siegert boundary conditions

In this case only solutions vanishing at x = 0 are acceptable and matching conditions at x = 0 reduce to

Sin[k a]= A (e^(i k a));    i k A (e^(i k a)) - k Cos[k a] = 2β Sin[ka];

Eliminating A

(ik-2β)Sin[ka] - k Cos[ka] = 0  .

As expected this is the condition for odd states in the previous problem.

The large β limit can be studied by writing

k = (π/a)n + ((Subscript[d, 1])/((a^2)β))+ ((Subscript[d, 2])/((a^3)(β^2)))+...

Inserting in previous equation one finds

(Subscript[d, 1])= n (π/2); (Subscript[d, 2])= ((n π)/4) - i(((n^2)(π^2))/4).

Finally for k and E

k = (π/a)n + n (π/2)(1/((a^2)β)) +( ((n π)/4) - i(((n^2)(π^2))/4)) (1/((a^3)(β^2))) + 𝒪((β^(-3)));
E = (Subscript[E, R])-i (Γ/2) ≃ ((ℏ^2)/(2m))( (((n^2)(π^2))/( (a^2)))+ (((n^2)(π^2))/((a^3)β)) + ((3/4)(((n^2)(π^2))/((a^4)(β^2))) - i (1/2)(((n^3)(π^3))/((a^4)(β^2))))) +  𝒪((β^(-3)))

which at lowest order coincides with WKB result.

Problem

Compute the spectrum of an anharmonic oscillator 1/2( (Subscript[λ, 2]) (x^2)+ (Subscript[λ, 4]) (x^4)) in the lowest order WKB approximation. (Subscript[λ, 2]) is positive. Compare the results to that of  perturbation theory.

Solution

Change of variables

The Schrödinger equation is

-((ℏ^2)/(2m))((d^2)/(dx^2))ψ +(1/2)((Subscript[λ, 2])(x^2)+ (Subscript[λ, 4])(x^4))ψ = E ψ;

The eigenvalues will depend on (Subscript[λ, 2]), (Subscript[λ, 4])and m. A simple change of variables show that the effective parameter is only one and that eq.() depends actually only on this parameter. Take

x = (((|(Subscript[λ, 2])| m ))^(-1/4)) z;  μ = Sign[(Subscript[λ, 2])] = ± 1.

Substituting in () we get

(1/2)(((((|(Subscript[λ, 2])|)/m)))^(1/2))(- (ℏ^2) ((d^2)/(dz^2))+ μ (z^2) + g (z^4))ψ = E ψ;    g = ((Subscript[λ, 4])/(Sqrt[((|(Subscript[λ, 2])|)^3)m])).

Then we have

E[m,(Subscript[λ, 2]),(Subscript[λ, 4])] = (((((|(Subscript[λ, 2])|)/m)))^(1/2)) E[1,μ,g]

This means that, for a given sign of the quadratic term, the only real parameter is g. In this problem we consider (Subscript[λ, 2])> 0 then our reference equation is (we call again x the space variable)

(-(1/2)(ℏ^2) ((d^2)/(dx^2))+ (1/2) (x^2) + (1/2) g (x^4))ψ = E ψ;    H = ((p^2)/2)+ (1/2) (x^2) + (1/2) g (x^4).

Note

The new change of variables x = (ℏ^(1/2))y can reabsorbs also the factor ℏ :

(-(1/2) ((d^2)/(dy^2))+ (1/2) (y^2) + (1/2) (gℏ) (y^4))ψ = (E/ℏ) ψ

This shows that in fact a weak coupling expansion in g is equivalent to an ℏ expansion in this model, if energies are measured in units of ℏ. As we want to keep separate g and ℏ variables in this problem to show a particular limit with g fixed and ℏ n fixed, n is the quantum number, we do not use here this last coordinate transformation. This is implicitly used in most of the numerical notebooks when "natural units", m = 1, ℏ = 1, ω =1 are used.

WKB quantization

The classical momentum in equation () is

p[x] = (Sqrt[2 E - (x^2)- (gx^4)])

It is convenient to rewrite this expression by extracting the turning points.

With

(a^2)=(( -1 + (Sqrt[1+ 8 g E]))/(2g)) ;  (b^2) = (( 1 + (Sqrt[1+ 8 g E]))/(2g))

we have

p[x] = (Sqrt[g])(Sqrt[((a^2)-(x^2))((b^2)+(x^2))])

The turning points are clearly

x = ± a

The WKB quantization condition is

(n+(1/2))= (1/(2π ℏ))∮ p[x]𝕕x = (1/(π ℏ))(∫_(-a)^(+a))p[x]𝕕x ≡ J[E].

Inserting the expression of p[x] and using the symmetry x → - x, we have

J[E] = ((2 (Sqrt[g]))/(π ℏ))(∫_0^a)(Sqrt[((a^2)-(x^2))((b^2)+(x^2))]) 𝕕x

The integral can be expressed through complete elliptic integrals. Their definition is

fE[m] = (∫_0^(π/2))(Sqrt[1- m ((Sin[φ])^2)]) 𝕕φ ;    fK[m]= (∫_0^(π/2)) ((𝕕φ)/(Sqrt[1- m ((Sin[φ])^2)])) .

For m we follow the notation used in Mathematica , some authors use (m^2)where we write m.  In terms of these integrals

(∫_0^a)(Sqrt[((a^2)-(x^2))((b^2)+(x^2))]) 𝕕x = ((Sqrt[(a^2)+(b^2)])/3)((b^2)fK[((a^2)/((a^2)+(b^2)))] - ((b^2)-(a^2))fE[((a^2)/((a^2)+(b^2)))])

It is possible to give different forms at the result using identities among elliptic integrals.

For the action integral we have

J[E] = ((2 (Sqrt[g]))/(π ℏ))((Sqrt[(a^2)+(b^2)])/3)((b^2)fK[((a^2)/((a^2)+(b^2)))] - ((b^2)-(a^2))fE[((a^2)/((a^2)+(b^2)))]).

Here is a plot of ℏ J with g = 1 :

Inverting the relation () we get the spectrum.

Limits

It is interesting to investigate small and large g values. This can be done by a series expansion in g or 1/g. By using Mathematica or consulting a book on Elliptic functions.

Small g

For small g

J[E,g](Underscript[→, g → 0]) (1/ℏ) ( E - (3/4)g (E^2) + (35/16) (g^2) (E^3) -((1155 (E^4) (g^3))/128)+((45045 (E^5) (g^4))/1024)... );

At lowest order

(Subscript[E, 0]) =ϵ = ℏ J = ℏ (n + (1/2))

the correct result for an unperturbed oscillator of unit frequency.

Using a power expansion

E = (Subscript[E, 0]) + (∑_k) (g^k) (Subscript[E, k])

the series is easily inverted with the result

(Subscript[E, WKB]) = ϵ+((3 g (ϵ^2))/4)-((17 (g^2) (ϵ^3))/16)+((375 (g^3) (ϵ^4))/128)-((10689 (g^4) (ϵ^5))/1024) + 𝒪((g^5))

Consider now the first orders in usual perturbative expansion, see notebook [NB-9.1] in chapter 9.

(Subscript[E, pert]) = ℏ(n+(1/2))+ (3/8) g (1+2 n+2 (n^2)) (ℏ^2)-(1/32) (g^2) (21+59 n+51 (n^2)+34 (n^3)) (ℏ^3)+(3/128) (g^3) (111+347 n+472 (n^2)+250 (n^3)+125 (n^4)) (ℏ^4)-(((g^4) (30885+111697 n+160470 (n^2)+142610 (n^3)+53445 (n^4)+21378 (n^5)) (ℏ^5))/2048)

Let us note that in agreement with eq.() the quantity E/ℏ depends only on the product (ℏ g).

It is easily seen that the two expansions do not agree in general, remember that ϵ = ℏ (n+1/2). As an instance for the ground state

(Subscript[E, pert])= (1/2)ℏ + (3/8) g  (ℏ^2) + ...;      (Subscript[E, WKB]) = (1/2)ℏ+((3 g )/16)(ℏ^2) + ...

But if we take the limit ℏ→0 , n→∞, with n ℏ finite, i.e. we take only the leading order in ℏ n in both expansions they coincide:

(Subscript[E, WKB]) = ℏ n+((3 g (((ℏ n))^2))/4)-((17 (g^2) (((ℏ n))^3))/16)+((375 (g^3) (((ℏ n))^4))/128)-((10689 (g^4) (((ℏ n))^5))/1024)+ ... = (Subscript[E, pert])

This confirms what stressed in the text : the semiclassical quantization formula is accurate in the double limit ℏ→0 , n→∞, with  (nℏ) finite.

To avoid misunderstanding let us repeat the conclusions in another form:

In conclusion for several regimes WKB expansion is more general than perturbative expansion For equal order of expansion results in of WKB and perturbation theory coincide, even if making a WKB expansion can be much more difficult than perform a perturbative expansion.

Large g

For large g

J[E,g](Underscript[→, g → ∞]) (((2^(5/4))  (Sqrt[π]) ((((1/g)))^(1/4)))/(3 ((Γ[(3/4)])^2))) (E^(3/4))  ⇒  E ∼ (n^(4/3)) (g^(1/3))

This is important as it indicates a cut in the complex g plane. The result implies also a growth with n faster than in the harmonic oscillator case, where (Subscript[E, n]) ∼ n.

Classical period

The classical period is given by (here m = 1)

T = ∮ ((dx)/(p[x]))= (4/(Sqrt[g]))(∫_0^a)((dx)/(Sqrt[((a^2)-(x^2))((b^2)+(x^2))])).

Also this integral can be expressed through elliptic integrals with the result

T[E] = (4/(Sqrt[g]))(1/(Sqrt[(a^2)+(b^2)]))fK[((a^2)/((a^2)+(b^2)))].

The asymptotic limits are

(1/(2π))T[E](Underscript[→, g → 0]) 1 - (3/2)g E;         (1/(2π)) T[E](Underscript[→, g → ∞]) ((((((1/g)))^(1/4)) (Sqrt[π]))/((2^(3/4))  ((Γ[(3/4)])^2)))(E^(-1/4)) .

The reader can verify that these limits satisfy the general requirement ( correspondence principle )

(dn/dE) = (T/(2π))

Problem

Consider a generic even potential V[x] with two degenerate absolute minima. Show that for small tunneling amplitudes the splitting between ground state and first excited state can be obtained by a variational technique based on the two semiclassical solutions built on the two minima.

Solution

A prototype of this problem is the double well anharmonic oscillator, with a potential of the form

V[x] = - a (x^2) + b (x^4);   a, b > 0.

The following considerations are valid for  a generic potential with this kind of shape.

We will use a particular feature of variational calculations: to obtain first order correct results in energy it is sufficient to know the states only at zero order, this is the essential feature of stationary points.

We assume that neglecting tunneling effects the system has two degenerate levels, pictorially depicted in the above figure. Let us consider a (semiclassical) solution (Subscript[φ, c])[x] describing a bound state. This state by definition goes exponentially to zero as x→ ∞. If we denote by K the penetration barrier integral the value of (Subscript[φ, c]) at x=0 will be of the order of (K^(1/2)) as semiclassical wave functions decrease exponentially in the forbidden region and x=0 is at half way in the barrier.

The same construction holds for left well, by symmetry. Consider now a function (Subscript[φ, 1])which coincides with (Subscript[φ, c]) for x>0 and is prolonged with continuity in the x < 0 region, in an exponential decreasing way. This can be done in infinite ways of course. Do the same for the left well. Then we have

(Subscript[φ, 1])[x] = (Subscript[φ, c])[x];  x > 0;      (Subscript[φ, 1])[x] = 𝒪((Sqrt[K])) for x < 0;
(Subscript[φ, 2])[x] = (Subscript[φ, c])[-x];  x < 0;   (Subscript[φ, 2])[x] = 𝒪((Sqrt[K])) for x > 0;

It is simpler to assume that the two functions are prolonged in the same way, even if this is not essential. We assume for simplicity that  (Subscript[φ, c]) are normalized in a semi-space (this can always be done, at the end we will write formulas for arbitrary normalization):

(∫_0^∞)(Subscript[φ, c]^2)[x] = 1

By construction the two functions satisfy Schrödinger equation in the two respective semi-spaces:

(H - (Subscript[E, 0]))(Subscript[φ, 1])= (Piecewise[{{0, x>0}, { 𝒪((Sqrt[K])), x<0}}])     ;    (H - (Subscript[E, 0]))(Subscript[φ, 2])= (Piecewise[{{0, x<0}, { 𝒪((Sqrt[K])), x>0}}])

We can take these two functions as a basis for a variational computation. We know that the minimum corresponds to eigenstates of the system

(Subscript[H, ij])(Subscript[c, j]) = E (Subscript[N, ij])(Subscript[c, j])

Where H and N are respectively the Hamiltonian matrix and the scalar product matrix. Let us compute these matrices. By symmetry:

(Subscript[N, 11])= (∫_(-∞)^∞)(Subscript[φ, 1])(Subscript[φ, 1])dx = (Subscript[N, 22])= (∫_(-∞)^∞)(Subscript[φ, 2])(Subscript[φ, 2])dx = 1 + 𝒪(K)

The 1 comes form normalization of (Subscript[φ, c]), the order K comes from splitting the integral in

(∫_(-∞)^∞)(Subscript[φ, 1])(Subscript[φ, 1])dx = (∫_0^∞)(Subscript[φ, 1])(Subscript[φ, 1])dx  + (∫_(-∞)^0)(Subscript[φ, 1])(Subscript[φ, 1])dx  = (∫_0^∞)(Subscript[φ, c])(Subscript[φ, c])dx  + (∫_(-∞)^0)(Subscript[φ, 1])(Subscript[φ, 1])dx = 1 +  𝒪(K)

Off diagonal elements can be estimated in the same way. In a semispace only one of the two functions is depressed then

(Subscript[N, 12])= (Subscript[N, 21]) =  𝒪((Sqrt[K]))

The same trick can be used for the Hamiltonian matrix

(Subscript[H, 11]) =(∫_(-∞)^∞)(Subscript[φ, 1])H (Subscript[φ, 1])dx = (∫_0^∞)(Subscript[φ, 1])H (Subscript[φ, 1])dx + (∫_(-∞)^0)(Subscript[φ, 1])H (Subscript[φ, 1])dx = (∫_0^∞)(Subscript[φ, 1])(Subscript[E, 0])(Subscript[φ, 1])dx + (∫_(-∞)^0)(Subscript[φ, 1])H (Subscript[φ, 1])dx  =
(Subscript[E, 0])(Subscript[N, 11])+ (∫_(-∞)^0)(Subscript[φ, 1])(H - (Subscript[E, 0])) (Subscript[φ, 1])dx = E (Subscript[N, 11])+ 𝒪(K)

The off diagonal term is a bit more difficult

(Subscript[H, 21])= (∫_(-∞)^∞)(Subscript[φ, 2])H (Subscript[φ, 1])dx = (∫_0^∞)(Subscript[φ, 2])H (Subscript[φ, 1])dx + (∫_(-∞)^0)(Subscript[φ, 2])H (Subscript[φ, 1])dx

In the first term we can use Schrödinger equation, while in the second this is not possible as (Subscript[φ, 1])satisfy this equation only for x > 0. We can make an integration by part, in such a way H will act on (Subscript[φ, 2]) which satisfy the equation for x < 0. We neglect terms at infinity:

(Subscript[(Subscript[(∫_(-∞)^0)(Subscript[φ, 2])(((d^2)/(dx^2)) )(Subscript[φ, 1])dx = ((Subscript[φ, 2]) (Subscript[φ, 1]^′))|, –∞]^0)- (∫_(-∞)^0)(Subscript[φ, 2]^′)(Subscript[φ, 1]^′)dx = (Subscript[φ, 2])[0] (Subscript[φ, 1]^′)[0]- (Subscript[φ, 2]^′)(Subscript[φ, 1])|, -∞]^0)+(∫_(-∞)^0)(Subscript[φ, 2]^(′′))(Subscript[φ, 1])dx =
(Subscript[φ, 2])[0] (Subscript[φ, 1]^′)[0]-(Subscript[φ, 2]^′)[0](Subscript[φ, 1])[0] + (∫_(-∞)^0)(Subscript[φ, 2]^(′′))(Subscript[φ, 1])dx

As (Subscript[φ, 2]) = (Subscript[φ, c])[-x] we have (Subscript[φ, 1]^′)[0]=(Subscript[φ, c]^′)[0]and (Subscript[φ, 2]^′)[0]= - (Subscript[φ, c]^′)[0]. Then with the help of Schrödinger equation in x < 0

(Subscript[H, 21])= (Subscript[N, 21])(Subscript[E, 0]) - ((ℏ^2)/(2m))(2 (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]) = (Subscript[N, 21])(Subscript[E, 0]) - ((ℏ^2)/m) (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]≡ (Subscript[N, 21])(Subscript[E, 0]) + δ

Collecting our matrix elements and writing explicitly K factors ((Subscript[σ, i])are Pauli matrices)

H = (Subscript[E, 0])N  + (Subscript[a, 1])K + δ (Subscript[σ, 1]);    N = 1 + (Subscript[b, 1])(Sqrt[K])(Subscript[σ, 1])+ (Subscript[b, 2])K (Subscript[σ, 3])

We have to solve

det (H - N E) = 0 ⇒ det ((N^(-1))H -  E)= 0.

We have

(N^(-1)) = 1 - (Subscript[b, 1])(Sqrt[K])(Subscript[σ, 1]) + 𝒪(K);     (N^(-1))H = (Subscript[E, 0])+(N^(-1)) δ (Subscript[σ, 1]) =  (Subscript[E, 0]) + (Subscript[δσ, 1])(1+ 𝒪(K))

It is important that all corrections cancel exactly in the leading term, as the (Subscript[E, 0])part of H was proportional to N, while all corrections just modify δ by higher order terms. The eigenvalues of (N^(-1))H are well known:

E = (Subscript[E, 0]) ± δ = (Subscript[E, 0]) ∓ ((ℏ^2)/m) (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]

(Subscript[E, 0]) + δ is the eigenvalue of the symmetric state, as for a decreasing function  (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]> 0 this corresponds to the ground state, as it must be. The energy splitting between symmetric and antisymmetric states is

ΔE = 2((ℏ^2)/m) (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]

For not normalized functions obviously :

ΔE = 2((ℏ^2)/m) (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]/ (∫_0^∞)(Subscript[φ, c]^2)dx

The result coincides with semiclassical result, in the semiclassical limit, but is more accurate, as it do not depend on connection formulas. As far as (Subscript[φ, c])is accurate the only errors in () are due to higher powers in K.

NOTE: This derivation make particularly clear why we can calculate the splitting. In whatever scheme, perturbative, semiclassical etc., the single levels will have some unknown corrections, higher order in the coupling, higher order in ℏ etc.. As the potential is symmetric in the difference ΔE all these corrections cancel exactly, what is left is the non perturbative corrections, due to tunneling. In eq.() it is apparent that the accuracy of the calculation depends on the accuracy of the wave function in the asymptotic region, as x=0 is far form the local minimum of the potential. This is the only point in which the calculation must be accurate.

Lowest order WKB

For classically normalized solutions the exponential tail of the wave function is, in the semiclassical approximation

(Subscript[φ, c])[x] = (C/(2(Sqrt[|p[x]|])))Exp[(1/ℏ)(∫_a^x)|p[x]|dx]

where a>0 is the turning point,  p[x] the momentum, C the normalization constant. For symmetric potentials

p'[0]∝ V'[0] = 0

then

2((ℏ^2)/m) (Subscript[φ, c])[0](Subscript[φ, c]^′)[0]=2((ℏ^2)/m)((C^2)/4)(1/(p[0]))((p[0])/ℏ)Exp[(2/ℏ)(∫_a^0)|p[x]|dx] = ((ℏ (C^2))/(2m))Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx]

ΔE =  ((ℏ (C^2))/(2m)) Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx] ≡ ((ℏ (C^2))/(2m)) K

With the lowest order normalization

C = 2 (Sqrt[(m/T)])

ΔE =   ((ℏ 2)/T) K = ((ℏ ω)/π)K

where ω is the classical oscillation in the well. In this formula the only error is in the determination of C.

Problem

Compute in the lowest order WKB approximation the energy splitting between the two lowest-lying states in a double well.

Solution

Change of variables

In previous problem it has been shown that by a change of variables

x = (ℏ^(1/2))(((|(Subscript[λ, 2])| m ))^(-1/4)) z;  μ = Sign[(Subscript[λ, 2])] = ± 1.

the Schrödinger equation

-((ℏ^2)/(2m))((d^2)/(dx^2))ψ +(1/2)((Subscript[λ, 2])(x^2)+ (Subscript[λ, 4])(x^4))ψ = E ψ;

can be transformed in a standard one

(-(1/2) ((d^2)/(dx^2))+ (μ/2) (x^2) + (1/2) g (x^4))ψ = E ψ;    H = ((p^2)/2)+μ (1/2) (x^2) + (1/2) g (x^4).

where μ = -1 and

g = (((Subscript[λ, 4])ℏ)/(Sqrt[((|(Subscript[λ, 2])|)^3)m]));    E[ℏ,m,(Subscript[λ, 2]),(Subscript[λ, 4])] = ℏ(((((|(Subscript[λ, 2])|)/m)))^(1/2)) E[1,1,μ,g]

In the following we will use the standard form (). In the present case, a double well potential, μ = - 1. At variance with respect previous problem we absorb ℏ in the change of variables as we will not be interested in the large n limit (n is the quantum number). The potential is

V[x] = -(1/2) (x^2) + (1/2) g (x^4)

A plot of the potential is:

Condition for tunneling

The inversions points are given by the solution of the equation

2 E - (x^2) - g (x^4) = 0

and are  x = ± a, ± b where

(a^2)= ((1-(Sqrt[1+8 E g]))/(2g));    (b^2)= ((1+(Sqrt[1+8 E g]))/(2g))

Points  ± a  are there only for E < 0. For E > 0 there are only two turning points,  ± b.  The classical momentum can be written as

p[x] = (Sqrt[g]) (Sqrt[((x^2)-(a^2))((b^2)- (x^2))])

We are interested in the case in which two energy levels are negative, these will correspond to the spitting of two approximately degenerate eigenstates describing a particle in the left or right well respectively. The tunnel effect remove the degeneracy.

At the semiclassical level the energy of the bound states in each of the two well are given by the quantization condition (ℏ = 1)

J[E] = (1/( π))(∫_a^b)p[x] d x = n+ (1/2)

To have a tunneling between states confined in the two wells at least one level in each of the well must be present, with negative energy. The limit case is E = 0. For E = 0, we have a = 0, b = 1/g and the action integral becomes

J[0] = (1/( π))(∫_a^b)p[x] d x = (Sqrt[g])(∫_0^(1/(Sqrt[g]))) x(Sqrt[(1/g)- (x^2)])d x = (1/( (Sqrt[g])))(∫_0^1) (Sqrt[1- (z^2)])d z = (1/(3 g  π))

To have a bound state the following condition must hold

J[0]> (1/2)   ⇒    (1/(3 g  π)) > (1/2)   ⇒  g < (2/(3π)) ∼ 0.21

The tunneling

In the text it has been shown how at a semiclassical level the tunnel effect can provide an exponentially small splitting between the first two bound states of the system. An alternative proof was given in previous problem. The result is

ΔE =  ((ℏ (C^2))/(2m)) Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx]

where C is a normalization for semiclassical state, in lowest order, with T the classical period of motion:

C = 2 (Sqrt[(m/T)])

It follows

ΔE =  ((2 ℏ)/T)  Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx]≡ ((2 ℏ)/T) K .

K is the penetration barrier factor.

Evaluation of the integrals

We see that all physical interesting quantities are given by integrals of p[x] i.e. of square roots of quartic polynomials. All this integrals can be expressed trough elliptic integrals.  The relevant integrals are:

(∫_a^b)(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = (b/3)(((a^2)+(b^2))EllipticE[(((b^2)-(a^2))/(b^2))] - 2 (a^2)EllipticK[(((b^2)-(a^2))/(b^2))])
(∫_0^a)(Sqrt[( (a^2)- x^2)((b^2)-(x^2))])dx = (b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] -  ((b^2) - (a^2))EllipticK[((a^2)/(b^2))])
(∫_a^b)(1/(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]))dx = (1/b)EllipticK[(((b^2)-(a^2))/(b^2))]

And we have respectively :

J = ((Sqrt[g])/π)(∫_a^b)(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = ((Sqrt[g])/π)(b/3)(((a^2)+(b^2))EllipticE[(((b^2)-(a^2))/(b^2))] - 2 (a^2)EllipticK[(((b^2)-(a^2))/(b^2))]);

T = (2/(Sqrt[g]))(∫_a^b)(1/(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]))dx = (2/(Sqrt[g]))(1/b)EllipticK[(((b^2)-(a^2))/(b^2))];

K = Exp[-2 D];  D = (Sqrt[g])(∫_0^a)(Sqrt[( (a^2)- x^2)((b^2)-(x^2))])dx = (Sqrt[g]) (b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] -  ((b^2) - (a^2))EllipticK[((a^2)/(b^2))])

From these expressions ΔE can be computed, and compared with numerical evaluation of eigenvalues. This is done in notebook [NB-11.2]. We report here a plot: the line is the WKB result, points are numerical eigenvalues:

This subject is explored again in notebooks [NB-11.2] and in problems.

Note: the result is not accurate for g → 0, as the subleading term in the WKB phase is not under control in this approximation: a precise determination of the constant C in eq.() require a more refined analysis, this will be done for the anharmonic oscillator in problem [11] and in notebook [NB-11.2].

Problem

Consider the Schrödinger equation for a potential V[x]. Let us suppose that V[x] is analytic, then the equation can be extended in the complex domain. Set

ψ[x] = Exp[i ((σ[x])/ℏ)];

σ[x] in general will have branch points for classical turning points. Let us suppose that there are only two of them. The monodromy requirement on ψ implies the constraint

(1/ℏ)∮ σ[x] dx = 2π n .

Show that this lead to leading order to the Bohr - Sommerfeld quantization conditions, and allows to compute higher order corrections in ℏ in the WKB quantization procedure.

Solution

The position

ψ[x] = Exp[i ((σ[x])/ℏ)];

transforms the linear second order Schrödinger equation for ψ in a nonlinear first order equation for σ'.

The only unknown quantity is σ' and in the mathematical literature the substitution () is usually written in the form

ψ[x]= Exp[(∫^x)f[x]dx]

The resulting equation is known as a Riccati equation. We will follow the more conventional notation used in physics and will adopt the form ().

NOTE  In the following we will take the signs adapted to work in the classical allowed regions. Formulas in the classical forbidden region can be obtained by analytic continuation but to avoid problems with signs and for didactic reasons at the end of the problem we give also expressions adapted to computations in the forbidden regions.

As shown in the text by substitution of () in the Schrödinger equation

-((ℏ^2)/(2m))ψ''[x] + V[x] ψ[x] = E ψ[x]

one obtains

(1/(2m))((σ'[x])^2) - i (ℏ/m) σ''[x] + V[x] = E .

With

((p[x])^2)= 2m (E - V[x]),

the classical momentum, at fixed energy,

((σ'[x])^2) - i  σ''[x] = ((p[x])^2).

One recognize that at leading order in ℏ, σ[x] is the classical action

(Subscript[σ, 0])[x] = ±(∫^x)p[x]dx.

To make an ℏ expansion means to write

σ[x] =(∑_(k=0)^∞)((((ℏ/i)))^n)(Subscript[σ, k]) = (Subscript[σ, 0])+ (ℏ/i)(Subscript[σ, 1])- (ℏ^2)(Subscript[σ, 2]) + ...

Inserting this expansion in () one obtain at lowest order

(((Subscript[σ, 0])'[x])^2)= ((p[x])^2)

This equation has two roots, let us take for definiteness the positive root. Imposing that coefficients of (ℏ^n)cancel in () one has

(∑_(k=0)^n)(Subscript[σ, k]^′) (Subscript[σ, n-k]^′) + (Subscript[σ, n-1]^(′′)) = 0

This recursion relation is easily solved order by order as (Subscript[σ, n]^′) appears only in two terms in the sum, the first  and the last.

(Subscript[σ, 1]^′)[x] =- (1/(2(Subscript[σ, 0]^′)))(Subscript[σ, 0]^(′′));       (Subscript[σ, n]^′)[x] = - (1/(2(Subscript[σ, 0]^′)))((∑_(k=1)^(n-1))(Subscript[σ, k]^′) (Subscript[σ, n-k]^′)+ (Subscript[σ, n-1]^(′′)))   n≥ 2

Explicitly:

(Subscript[σ, 1]^′)[x] =- (1/(2(Subscript[σ, 0]^′)))(Subscript[σ, 0]^(′′))   ⇒  (Subscript[σ, 1])=- (1/2)Log[p[x]];
(Subscript[σ, 2]^′)[x] =- (1/(2(Subscript[σ, 0]^′)))((((Subscript[σ, 1]^′)[x])^2)+ (Subscript[σ, 1]^(′′))[x])⇒ (1/4)((p''[x])/((p[x])^2))- (3/8)(((p'[x])^2)/((p[x])^3))
            =- (( 1)/(2 ((p[x])^(1/2))))((d^2)/(dx^2))((p[x])^(-1/2));
(Subscript[σ, 3]^′)[x] =- (1/(2(Subscript[σ, 0]^′)))(2(Subscript[σ, 1]^′)[x](Subscript[σ, 2]^′)[x]+ (Subscript[σ, 2]^(′′))[x])⇒(Subscript[σ, 3]^′)[x]=- (1/2)(d/dx)(((Subscript[σ, 2]^′)[x])/(p[x])) ....

With p[x] real all these quantities are real, this means that only even orders appears as phase factors in the expansion of ().

Near  a classical inversion point (p^2)[x] = (a-x) R[x], where for simplicity we assume that R[a]≠0, i.e. a is a single zero. This imply that σ'[x], the square root of (p^2), has  a cut starting from x=a. These kind of cuts propagate in higher order approximations. For x∈ℝ ψ[x] must be a single valued function (i.e. must be a function) this means that we are approximating a single value function with a non single valued sequence of functions.

Let us suppose that (p^2)[x] has two zeros, (Subscript[x, 1]) and (Subscript[x, 2]), than a cut is present between these two values. Approaching the cut from above or from below must give the same value. The phase factor can be defined by usual analytic continuation along an arbitrary path.

The request of monodromy for ψ imply that the difference between the values on both side of the cut must be a multiple of 2π. By deforming the contour if necessary this imply

(1/ℏ)∮ dσ=(1/ℏ)∮ σ'[x] dx = 2π n .

We stress that this an exact constraint.

The location of the turning points, and the value of the integral, depend on the energy, than () is a quantization condition for energy. If the series in ℏ converges to exact solution than eigenvalues must satisfy this relation.

At first order () reduces to Bohr Sommerfeld condition

(1/ℏ)∮ (Subscript[σ, 0]^′)[x] dx = (1/ℏ)∮ (Subscript[p, 0])[x] dx  = (2/ℏ)(∫_(Subscript[x, 1])^(Subscript[x, 2]))p[x] dx = 2π n .

Next order is a bit special. For a parametrization

p[x]= (Sqrt[(x-(Subscript[x, 1]))((Subscript[x, 2])-x)R[x]])

we have from ()

(Subscript[σ, 1]^′)[z] =- (1/2)(d/dz)Log[p] =- (1/4)((1/(z-(Subscript[x, 1])))+(1/(z-(Subscript[x, 2])))+ ((R'[z])/(R[z])))

Last term is regular, then it does not contribute to the contour integral. The first two terms have poles, not branching point singularities, and their contribution, using Cauchy theorem is

∮ (Subscript[σ, 1]^′)[x] dx  =- (1/4)(2 (2π i)) = - i π

Then we have, up to the first order

(1/ℏ)∮ σ'[x] dx = (1/ℏ)∮ (Subscript[σ, 0]^′)[x] dx +(1/ℏ) (ℏ/i)∮ (Subscript[σ, 1]^′)[x] dx = (1/ℏ)∮ (Subscript[σ, 0]^′)[x] dx -  π

and the quantization condition becomes

(1/ℏ)∮ (Subscript[σ, 0]^′)[x] dx = 2(n+(1/2))π

i.e. the WKB quantization condition.

(Subscript[σ, 1])is special as its derivative has poles instead of cuts, all other terms do not have these type of singularities.

It is expected that, apart (Subscript[σ, 1]), all other odd terms in the expansion do not contribute to the phase and so do not participate to the quantization condition. This is indeed true as they are all total derivative and their integral on a closed loop gives zero. This is most easily seen by separating modulus and phase in ()

ψ[x] = A Exp[i (S/ℏ)];   A, S ∈ ℝ

Odd terms in the expansion of σ contribute to A, even terms to S. Substituting in the Schrödinger equation and separating real and imaginary parts

(1/(2m))(A[x] ((S'[x])^2) - (ℏ^2)A''[x]) + V[x] = E;  
S''[x] A[x]+ 2 S'[x]A'[x] = 0 ⇒ ((A'[x])/(A[x]))= -(1/2)(d/dx)Log[S'[x]].

A'/A is the derivative of the odd part of σ. Log[S']once expanded in power series of ℏ has a logarithmic branch cut, the one giving rise to (Subscript[σ, 1])while other terms come from the Taylor expansion of the logarithm, and are usual functions, with at most brunch cuts between (Subscript[x, 1])and (Subscript[x, 2]). Now the analytic continuation means exactly that the derivative along a closed circuit is zero, as the continuation is defined by an integral along a path, then all this terms do not contribute to the quantization condition.

All subsequent even terms produce a correction to quantization condition. Let us consider for example the first correction, due to (Subscript[σ, 2.])At second order we must impose

(1/ℏ)∮((Subscript[σ, 0]^′)[x] + ((((ℏ/i)))^2)(Subscript[σ, 2]^′)[x]) = (1/ℏ)∮(p[x] - (Subscript[σ, 2]^′)[x]) = 2π(n+(1/2))

The integral involving (Subscript[σ, 2])can in principle be computed on complex domain, but it is often convenient to reduce it to a manageable form on real numbers. In the following manipulations it is convenient to write

((p[x])^2)= Q[x]

in order to easily identify square roots. Then

p'[x] =(1/2) ((Q'[x])/((Q[x])^(1/2)));   p''[x] =(1/2) (1/((Q[x])^(1/2)))(Q''[x]- (1/2)(((Q'[x])^2)/(Q[x])))

From eq.()

(Subscript[σ, 2]^′)[x]= (1/4)((p''[x])/((p[x])^2))-(3/8)(((p'[x])^2)/((p[x])^3))= (1/8)((Q''[x])/((Q[x])^(3/2)))- (5/32)(((Q'[x])^2)/((Q[x])^(5/2)))

We can extract a total derivative using

(d/dx)((Q'[x])/((Q[x])^(3/2)))= ((Q''[x])/((Q[x])^(3/2)))- (3/2)(((Q'[x])^2)/((Q[x])^(5/2)))

and write

(Subscript[σ, 2]^′)[x]= (1/48)((Q''[x])/((Q[x])^(3/2)))+ (5/48)(d/dx)((Q'[x])/((Q[x])^(3/2)))

Last term do not contribute to the contour integral. First term can be transformed in a manageable way by noticing that

((Q[x])^(-3/2))= (((2m(E-V[x])))^(-3/2))=- (1/m)(∂/(∂E))((Q[x])^(-1/2));  Q''[x] = -2m V''[x].

∮ (Subscript[σ, 2]^′)[x]dx = (1/24)(∂/(∂E))∮ dz((V''[z])/(((2m(E-V[x])))^(1/2)))

and the quantization condition finally becomes

∮ p[z]dz - ((ℏ^2)/24)(∂/(∂E))∮ dz((V''[z])/(((2m(E-V[x])))^(1/2)))= 2π ℏ (n+(1/2))

In terms of action integral (between the two turning points):

J + δJ = n+ (1/2)

where

J = (1/π)(∫_(Subscript[x, 1])^(Subscript[x, 2]))p[x]dx ;   δJ = - ((ℏ^2)/(24π))(∂/(∂E))(∫_(Subscript[x, 1])^(Subscript[x, 2]))((V''[z])/(((2m(E-V[x])))^(1/2))) .

This formula will be widely used in the numerical notebooks.

The classical forbidden region

It is convenient to define

q[x] = (Sqrt[2m(V[x]-E)])

and write the semiclassical expansion in the form of a Riccati equation

ψ = Exp[(1/ℏ)(∫^x)f[x] d x];   f = (Subscript[f, 0])+ ℏ (Subscript[f, 1])+ ℏ (Subscript[f, 2])+ ...

By substitution in the Schrödinger equation one obtains

(f^2) + ℏ f' = (q^2)[x]

At the lowest order we have the two solutions

(Subscript[f, 0])= ± q[x]

Take for example the plus sign. The next equation is

2 (Subscript[f, 0])(Subscript[f, 1]) + (Subscript[f, 0]^′) = 0 ⇒ (Subscript[f, 1]) =- (1/2)(d/dx)Log[|q[x]|]  ⇒ Exp[ℏ (∫^x)(Subscript[f, 1])] = (1/(Sqrt[|q[x]|]))

Let us stress explicitly that this term do not depend on the sign of q[x] chosen in the solution (Subscript[f, 0]), i.e. it is equal for both solutions.

For second order

2 (Subscript[f, 0])(Subscript[f, 2]) + (Subscript[f, 1]^2) + (Subscript[f, 1]^′) = 0 ⇒ (Subscript[f, 2])= - (1/(2(Subscript[f, 0])))((Subscript[f, 1]^2) + (Subscript[f, 1]^′))=(1/4)((q''[x])/((q[x])^2))- (3/8)(((q'[x])^2)/((q[x])^3))

and so on.

Problem

Compute the energy splitting between the two lowest-lying states in a double well,   in the limit g → 0.

Solution

Statement of the problem

A potential with a double degenerate minimum, like

has a double degenerate ground state if tunneling is neglected. We know that in one dimension the ground state is unique and, for a reflection invariant potential, we expect an even ground state and an odd first excited state. Tunneling provide the mechanism to realize this picture as the first two states are the even and odd superposition of |+⟩ and |-⟩, the semiclassical states confined in the two wells. To compute the energy splitting we have to solve Schrödinger equation taking account for tunneling. This can be done in several ways, which will be reviewed below. In the text a semiclassical approach has been used, in problem [9]  has been presented a variational - like approximation. In both cases the result was

ΔE = ξ ((2 ℏ)/T)  Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx]≡ ξ((2 ℏ)/T) K =ξ((ℏ ω)/π)K

x = ± a are the classical turning points indicated in the graph, K is the penetration barrier factor, T is the classical period of motion inside a well, ω = 2π/T the corresponding frequency. We have explicitly written a prefactor ξ ∼ 1 to stress the fact that this formula is only "exponentially correct", i.e. the prefactor is ambiguous as it depends on the precise form of the WKB wave function in the classical forbidden region.

In this problem we want to compute exactly the splitting () in the particular case of an anharmonic oscillator (we use natural units)

V[x] = -(1/2)(x^2) + g (1/2)(x^4)

The potential has two degenerate minima at

x = ± v;   v = 1/(Sqrt[2g]);    V[v] = -1/8g.

We will often add the constant 1/8g to the potential and write

V[x] = (1/2)g ((((x^2)- (v^2)))^2)=  (1/2)g (((x- v))^2)(((x+v))^2).

Small g and ℏ → 0

Consider a generic potential which can be written in the form

V[x] → (1/g)V[(Sqrt[g]) x]

This is clearly the case for an anharmonic oscillator :

-(1/2)(x^2) + g (1/2)(x^4) = (1/(2g))(- g (x^2)+ (g^2)(x^4)).

Consider now the Schrödinger equation

-((ℏ^2)/(2m)) ((d^2)/(dx^2))ψ + (1/g)V[(Sqrt[g]) x] ψ = E ψ

With the change of variables, which do not affect the spectrum,

x = (ℏ^(1/2))(m^(-1/4)) z

it becomes

-(1/2) ((d^2)/(dz^2))ψ + (1/λ)V[(Sqrt[λ]) x] ψ = ((Sqrt[m])/ℏ) E ψ ≡ ϵ ψ;   with   λ = g ℏ (m^(-1/2))

This means

E[ℏ,m,g] = (ℏ/(Sqrt[m]))E[1,1,g ℏ (m^(-1/2))]

i.e. the only real expansion parameter is  ℏ g, the classical limit ℏ→0 is closely related to the perturbative limit g → 0. This is well known, for example, in the usual anharmonic oscillator, where

(Subscript[E, pert]) = ℏ(n+(1/2))+ (3/8) g (1+2 n+2 (n^2)) (ℏ^2)-(1/32) (g^2) (21+59 n+51 (n^2)+34 (n^3)) (ℏ^3)+...

which has indeed the form (). This observation will be crucial to understand how to compute the prefactor ξ in eq.() in the weak coupling limit.

Connection formulas

The use of connection formulas is probably the simplest way to obtain the result () and is the method presented in the text. Let us briefly review it to point out the critical points. One start from the asymptotic behavior at large x, x ≫ b, and apply repeatedly the usual connection formulas

Cos[|S|+(π/4)] ↔  (e^(+|S|));   Sin[|S|+(π/4)]↔ (1/2) (e^(-|S|))

to reach the middle point x = 0.

S is the classical action in units of ℏ. Using the symbol w[a,x] for this quantity with specified extrema (the notation used in the text) we have in all details (we omit the prefactor 1/(Sqrt[p]) as it is inessential in the following ):

(e^(- w[b,x]))→ 2 Sin[w[x,b]+(π/4)] = 2 Sin[w[a,b]-w[a,x]+(π/4)]= 
2(Sin[w[a,b]Cos[w[a,x]-(π/4)]- Cos[w[a,b]Sin[w[a,x]-(π/4)])=
2(Sin[w[a,b]Cos[w[a,x]-(π/4)]+ Cos[w[a,b]Cos[w[a,x]+(π/4)])→ 
Sin[w[a,b](e^(- w[x,a]))+ 2 Cos[w[a,b](e^( w[x,a]))=
 Sin[w[a,b](e^(- w[0,a]+w[0,x]))+ 2 Cos[w[a,b](e^( + w[0,a]-w[0,x]))

For even/odd states the combination of the exponential must reduce to Cosh[w[0,x]] or Sinh[w[0,x]] then we have for the ground state and the excited state respectively

((Cos[w[a,b]])/(Sin[w[a,b]))= (1/2)(e^(- 2w[0,a]))≡ (1/2)K;         ((Cos[w[a,b]])/(Sin[w[a,b]))=- (1/2)(e^(- 2w[0,a]))≡ -(1/2)K;

The splitting will amount to a small non perturbative (in ℏ) correction to quantization formulas. At leading order K = 0 and, for the ground state in the well w[a,b] = π/2. In general w[a,b] = π/2 - μ and we have, at first order in μ from ():

((Cos[w[a,b]])/(Sin[w[a,b]))≃ Sin[μ]≃ μ = (1/2)K.

The second possibility giving  μ = - K/2. The quantization condition becomes, in the first case:

w[a,b] = (1/ℏ)(∫_a^b)(Sqrt[2m(E-V)])dx = (π/2) - (1/2)K

and expanding with E = (Subscript[E, 0])+ δE ( (Subscript[E, 0])is the unperturbed ground state) the first order correction is

(1/2) (2m δE)(1/ℏ)(∫_a^b)(1/(Sqrt[2m((Subscript[E, 0])-V)]))dx = δE (T/(2ℏ)) = -(1/2)K  ⇒  (Subscript[δE, 1]) = -(ℏ/T)K = -((ℏω)/(2π))K

Similarly for the odd combination

(Subscript[δE, 2])= + ((ℏω)/(2π))K

and we have

Δ E = (Subscript[E, 2])- (Subscript[E, 1]) = ((ℏω)/π)K

which is the result () with ξ = 1.  What is missing in this derivation which produce the ξ factor?

Asymptotic formulas

The answer to all these problem is to build a uniform expansion in ℏ, valid up to the middle point x=0. In this way all subleading factors are under control and ξ can be computed. The key point from the mathematical point of view is again the classical theorem of Liouville quoted in the file [WKBresults.nb]

Let

y''[x] + ((λ^2) P[x] + Q[x,λ])y[x] = 0

a linear second order differential equation depending on some parameter λ.  If Q is uniformly bounded al λ → ∞ then the asymptotic solutions (in λ) of equation () are given by the solutions of the simpler equation

y''[x] + (λ^2) P[x] y[x] = 0

We can use this theorem in two different, even if almost equivalent, ways.

Uniform approximation

In file [WKBresults.nb] it has been briefly discussed the technique of uniform approximation for a problem with two turning points. The original  Schrödinger equation in ψ[x] and a new Schrödinger equation in u[z]

a)   ψ''[x] + (k^2)[x] ψ[x] = 0;      b)      u''[z] + (q^2)[z] u[z] = 0;

are equivalent, i.e. they give raise to the same solutions, with

ψ[x] = u[z]/(((z'[x]))^(1/2))

if z and x are related by

(((z'[x]))^2)(q^2)[z] = (k^2)[x] - (1/2){z;x}

where {z; x} denotes the Schwarzian derivative

{z;x}= ((z''')/(z')) - (3/2)(((((z'')/(z'))))^2)

As q and k contain, implicitly, a factor 1/ℏ, we see that if the Schwarzian derivative is uniformly bounded, it plays the role of Q in eq.() and can be neglected. For two points boundary problems

(q^2)[z] = t - (z^2)

and the mapping x → z is given by

z[x]- (Sqrt[t]) = ((∫_a)^x)k[x]𝕕x;             (π/2)t =  (∫_(-(Sqrt[t]))^(Sqrt[t]))q[z] 𝕕z = (∫_(Subscript[x, 1])^(Subscript[x, 2]))k[x]𝕕x

The general solutions of eq.() are parabolic cylinder functions. There is only one function which is normalizable for the whole z-range and, as it is well known from the theory of harmonic oscillator it is

u[z]=(Subscript[H, n])[z] Exp[-(z^2)/2]

In fact, t in () plays the role of 2E in the usual oscillator and it must be t = 2 (Subscript[E, HO])= 2n+1. This n is the one which appears in ().

We will use this kind of uniform approximation to derive the energy splitting below, in a couple of different ways.

Perturbative expansion

For potentials of the form 1/g V[(Sqrt[g])x] the Schrödinger equation takes the form

-(1/2) ((d^2)/(dx^2))ψ + (1/g)V[(Sqrt[g]) x] ψ =  E ψ

and we see immediately that the asymptotic regime is connected with the limit g→0, with (Sqrt[g]) x fixed. We note that in this limit the energy E plays the role of a subleading term, this, in this language, explain why the splitting do not depend on the details of E[g], at least at leading order, but only on the form of the potential.

We will use this technique, introduced for anharmonic oscillator in the works of Brezin, Parisi, Zinn-Justin etc., see ref.[2 , 3], to study in a systematic way the anharmonic oscillator.

Connection formulas and uniform expansion

From asymptotic form of the parabolic cylinder functions, it is possible to derive a general connection formula connecting the left hand side and the right and side of a well, see text and file  [WKBresults.nb]. In the notation of eq.() (here we consider left well):

((Exp[- w[x,-b]])/(Sqrt[|p[x]|]))(Underscript[←, x→-∞]) ψ = (2^(1/4))(((2(e/t)))^(t/4))(1/(Sqrt[z']))(Subscript[D, ((t-1)/2)])[- (Sqrt[2]) z] (Underscript[→, x→+∞])
 ((((2/π)))^(1/2))(((((e )/J)))^J)Γ[J+(1/2)]Cos[π J]((Exp[+ w[-a,x]])/(Sqrt[|p[x]|]))+ Sin[π J]((Exp[- w[-a,x]])/(Sqrt[|p[x]|]))=
 ((((2/π)))^(1/2))(((((e )/J)))^J)Γ[J+(1/2)]Cos[π J]((Exp[+ w[-a,0]+ w[0,x]])/(Sqrt[|p[x]|]))+ Sin[π J]((Exp[- w[-a,0]-w[0,x]])/(Sqrt[|p[x]|]))

J is the usual action integral. We stress that at the leading order it must be J = n+1/2, this kill the exponential growing therm in Cos[π J] in previous formula.

Let us now make the same considerations as in the case of usual WKB. The wave function must be even/odd for ground state/excited state so we must have Cosh[w[0,x] and Sinh[w[0,x]]. We consider ground state for definiteness, excited state differ just by a sign in the formulas. We have the constraint

((((2/π)))^(1/2))(((((e )/J)))^J)Γ[J+(1/2)]((Cos[π J])/(Sin[π J])) = Exp[-2  w[-a,0]] = K

With J = (n+1/2) - μ

((Cos[π J])/(Sin[π J]))≃ π μ;

So at first order in the prefactor we can use J=n+1/2 and we have

((((2/π)))^(1/2))(((((e )/(n+1/2))))^(n+1/2))n! π μ = K

Ground state

For n = 0, the ground state

((2(Sqrt[e]))/(Sqrt[π]))π μ = K;  μ = (Sqrt[(π/e)])(1/2)K

which  differ from () by a factor (Sqrt[π/e]), than we have:

Δ E =  (Subscript[E, 2])- (Subscript[E, 1]) =(Sqrt[(π/e)]) ((ℏω)/π)K

Excited states

For enough deep wells we can have doubling also for excited states. For large n we can use Stirling formula

n! ∼ (n^n)(e^(-n))(Sqrt[2π n])

to deduce for large n

((((2/π)))^(1/2))(((((e )/(n+1/2))))^(n+1/2))n!→((((2/π)))^(1/2))(e^n) (Sqrt[(1/n)]) (n^(-n))(n^n)(e^(-n))(Sqrt[2π n]) = 2

and we recover the usual WKB result ():  μ = K/2 and

Δ E =  (Subscript[E, 2])- (Subscript[E, 1]) = ((ℏω)/π)K

As expected usual WKB is always accurate for n large.

In the general case

Δ (Subscript[E, n]) =  (Subscript[E, n]^odd)- (Subscript[E, n]^even) =2 ((((π/2)))^(1/2))((((((n+1/2))/e)))^(n+(1/2)))(1/(n!)) ((ℏω)/π)K

The prefactor in () is always close to one, its range is from (Sqrt[π/e])∼1.07 for the ground state to 1 for n→∞.

The anharmonic oscillator

The potential is

V[x] = -(1/2)(x^2) + g (1/2)(x^4)

and, as has been shown at the beginning of the problem:

V[x] =  (1/2)g (((x- v))^2)(((x+v))^2)+(1/(8g));  v = (1/(Sqrt[2g]))

With x = v + ξ we see that a well has and effective potential, at lowest order

V[ξ] ∼ (1/2)g (ξ^2)4(v^2)+ (1/(8g)) = (1/2)2 (ξ^2) + (1/(8g)).

This correspond to harmonic oscillations with frequency

ω = (Sqrt[2]) ;   T = ((2π)/ω)= (Sqrt[2])π

For the anharmonic oscillator it has been shown in problem [9] that

J = ((Sqrt[g])/π)(∫_a^b)(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = ((Sqrt[g])/π)(b/3)(((a^2)+(b^2))EllipticE[(((b^2)-(a^2))/(b^2))] - 2 (a^2)EllipticK[(((b^2)-(a^2))/(b^2))]);

T = (2/(Sqrt[g]))(∫_a^b)(1/(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]))dx = (2/(Sqrt[g]))(1/b)EllipticK[(((b^2)-(a^2))/(b^2))];

K = Exp[-2 D];  D = (Sqrt[g])(∫_0^a)(Sqrt[( (a^2)- x^2)((b^2)-(x^2))])dx = 
(Sqrt[g]) (b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] -  ((b^2) - (a^2))EllipticK[((a^2)/(b^2))])

where the turning points a and b are:

(a^2)= ((1-(Sqrt[1+8 E g]))/(2g));    (b^2)= ((1+(Sqrt[1+8 E g]))/(2g)).

To make a Taylor expansion in g we have to remember that in these formulas the zero of V[x] has not been settled, then we have E = -1/8g + ϵ. With this substitution a straightforward expansion gives

D = (1/(3 (Sqrt[2]) g))-(ϵ/(2 (Sqrt[2])))-2 (Sqrt[2]) ϵ Log[2]+((ϵ Log[g])/(2 (Sqrt[2])))+((ϵ Log[4 (Sqrt[2]) (Sqrt[ϵ])])/(Sqrt[2]));
T = (Sqrt[2]) π+((3 g π ϵ)/(Sqrt[2]))

J=(ϵ/(Sqrt[2]))+((3 g (ϵ^2))/(4 (Sqrt[2])))+((35 (g^2) (ϵ^3))/(16 (Sqrt[2])))

The leading order in T is in agreement with (), and as such could be predicted without any computation. At leading order in WKB the quantization condition J = n+1/2 correctly gives

ϵ = (Sqrt[2]) (n+(1/2))

as appropriate for an oscillator with ω = (Sqrt[2]). These values are correct up to order g then once inserted in the expression of D we obtain D with the same precision. We note that the leading term in D do not depend on ϵ and this has already be anticipated. Inserting these values in D and computing K = Exp[-2 D] we get

K = (2^((7/4)+((7 n)/2))) (e^((1/2)-((Sqrt[2])/(3 g))+n)) (((g+2 g n))^(-(1/2)-n)) = (2^((7/4)+((7 n)/2) - n - (1/2)))(1/(g^(n+(1/2))))((((e/(n+(1/2)))))^((n+(1/2))))(e^(-((Sqrt[2])/(3 g))))

and, with ω = (Sqrt[2]) from ()

ΔE = (((2^((9/4)+((5 n)/2))) )/((Sqrt[π]) n!))(1/(g^(n+(1/2))))(e^(-((Sqrt[2])/(3 g))))

This result is valid only for very small g, in the general case the result () is more reliable.

In particular for the ground state

ΔE = (((2^(9/4)) )/((Sqrt[π]) ))(1/(Sqrt[g]))(e^(-((Sqrt[2])/(3 g))))

For example of the accuracy of the approach we report here the result of a numerical computation for g = 0.005. For this coupling exist at least 20 "doublets". The first column are numerical energies splitting,  the second contains expression () and the last the naive WKB, i.e. expression () without the prefactor

A warning on the notations

In the literature, eg. in the articles of Zinn-Justin, the Schrödinger equation for this problem is usually written in the simpler form

-(1/2) ((d^2)/(dz^2))ψ + (1/2)(z^2)(((1-(Sqrt[λ])z))^2) ψ = (Subscript[E, Z]) ψ

To see the connection with our formulas we have to shift x, with x = ξ - v, then our Schrödinger equation  becomes

Eψ=-(1/2) ((d^2)/((dξ)^2))ψ + (1/2)(ξ^2)g((((2/(Sqrt[2g]))-ξ))^2) ψ = -(1/2) ((d^2)/((dξ)^2))ψ + (1/2)(ξ^2)((((Sqrt[2])-(Sqrt[g])ξ))^2) ψ=
 -(1/2) ((d^2)/((dξ)^2))ψ + (ξ^2)(((1-(Sqrt[(g/2)])ξ))^2) ψ

With the change of variables ξ = z (2^(-1/4)) we obtain

(E/2)ψ =  -(1/2) ((d^2)/(dz^2))ψ +(1/2) (z^2)(((1- (Sqrt[(g/2)])(2^(-1/4)) z))^2)ψ

which has the form () with

(Subscript[E, z])= (E/2);  g = 2 (Sqrt[2]) λ

Inserting in () we have

(Subscript[ΔE, Z]) = ((2 )/((Sqrt[π]) ))(1/(Sqrt[λ]))(e^(-(1/(6λ))))

which is the form usually found in literature.

The perturbative asymptotic expansion

Let us consider the shifted form of the potential, we write always x instead of ξ:

V = (g/2)(x^2)(((x + 2v))^2) = (g/2)(x^2)(((x + (Sqrt[(2/g)])))^2) = (1/2)(x^2)((((Sqrt[2]) + x(Sqrt[g])))^2) .

This is exactly of the general form () and the asymptotic small g expansion can be performed by usual WKB methods. We have to remember that formally x(Sqrt[g]) is held fixed, in this way g plays the role of (ℏ^2)in

-(1/2) ((d^2)/(dx^2))ψ + (1/g)V[x(Sqrt[g])] = E ψ .

A crucial point is that E do not contain g, so in the leading order can be neglected, and the expansion is similar to a zero energy WKB expansion. It is convenient to write as in usual WKB

ψ = Exp[-(∫^x)f[x,g]].

This transformation transforms the linear second order Schrödinger equation in a nonlinear first order equation, known as Riccati equation:

f' - (f^2)+ 2 (V/g) =2 E .

At the leading order

(Subscript[f, 0]) = (1/(Sqrt[g]))(Sqrt[2V]) ,

as expected from the interpretation g ∼ (ℏ^2). Higher orders will be computed as power series in g, just like usual WKB. From now on we will use the explicit form of V and write

(Subscript[f, 0])= x ((Sqrt[2]) + x(Sqrt[g])).

Let us note that for x(Sqrt[g]) fixed (Subscript[f, 0])is of order 1/(Sqrt[g]). At the next order (Subscript[f, 1])will be of order (Sqrt[g]) and we have

(Subscript[f, 0]^′) - 2 (Subscript[f, 0])(Subscript[f, 1]) = (Sqrt[2])+ 2x (Sqrt[g])- 2 x ((Sqrt[2]) + x(Sqrt[g]))(Subscript[f, 1])= 2 E .

E will have an expansion

E = (ω/2)+ (∑_(k=1)^∞)(Subscript[E, k])(g^k);   ω = (Sqrt[2]).

We remember that in our notation ω = (Sqrt[2]). Then

(Subscript[f, 1])= ((Sqrt[g])/((Sqrt[2]) + x(Sqrt[g]))) .

It is not difficult to write a recursion relation for the (Subscript[f, k]), but for the moment we limit ourselves to lowest order. In this approximation

ψ ∼ Exp[-(∫^x)(x ((Sqrt[2]) + x(Sqrt[g]))+  ((Sqrt[g])/((Sqrt[2]) + x(Sqrt[g]))))]= (C/((Sqrt[2]) + x(Sqrt[g])))Exp[-(( (x^2))/2)(Sqrt[2]) - ((x^3)/3)(Sqrt[g])];

C is an inessential constant which can be taken 1.

Let us now remember the general result obtained in problem [8]: if (Subscript[φ, c]) is the semiclassical wave function relative to one of the well then

ΔE = 2((ℏ^2)/m) (Subscript[φ, c])[(Subscript[x, 0])](Subscript[φ, c]^′)[(Subscript[x, 0])]/ (∫_(Subscript[x, 0])^∞)(Subscript[φ, c]^2)dx ,

(Subscript[x, 0]) is the symmetric point of the potential, in our shifted coordinates (Subscript[x, 0])= - v = - 1/(Sqrt[2g]). Substituting the expression () for φ we have

ψ[-v]ψ'[-v]= (((e^(-((Sqrt[2])/(3 g)))) ((Sqrt[2])-4 g))/((Sqrt[2]) (Sqrt[g])))≃ (1/(Sqrt[g]))(e^(-((Sqrt[2])/(3 g)))).

For the normalization integral we have to remember that we are performing an expansion in g and write

(∫_((-1)/(Sqrt[2g]))^∞)((((1/((Sqrt[2]) + x(Sqrt[g])))))^2)Exp[-2(( (x^2))/2)(Sqrt[2]) -2 ((x^3)/3)(Sqrt[g])]∼ (∫_(-∞)^(+∞))(1/2)Exp[- (x^2)(Sqrt[2])]dx = ((Sqrt[π])/(2^(5/4))) .

Substituting in (), in natural units:

ΔE = ((2^(9/4))/(Sqrt[π]))(1/(Sqrt[g]))(e^(-((Sqrt[2])/(3 g)))) ,

which is the previous result ().

Higher orders

The procedure can be iterated in a straightforward way. We refer to the literature [3] for the whole perturbative series, here we quote just the first few orders:

ΔE =(1/(Sqrt[π]))(2^(9/4))((e^(-((Sqrt[2])/(3g))))/(Sqrt[g]))(1-((71 g)/(24 (Sqrt[2])))-((6299 (g^2))/2304))

Uniform wave functions

Perturbative based expansion

This method is just another way to put previous results. We give an account of this method as it is probably the simplest one and give an insight into some deeper aspects of this problem.

To get rid of complicated numerical factors which would spoil the simplicity of the method let us shift the variables with the origin at the left well:

((-1)/2)ψ'' + (x^2)(((1+(Sqrt[(g/2)])x))^2)ψ = E ψ

then make the rescaling

x = α ξ = (2^(-1/4))ξ;  ϵ = 2 (α^2)E = (Sqrt[2])E ; λ = (g/(2(Sqrt[2])))

In this notations the Schrödinger equation take the simpler form

-ψ'' + (ξ^2)(((1- (Sqrt[λ])ξ))^2)ψ = ϵ ψ

The symmetric point of the potential (the local maximum) is at (Subscript[ξ, 0]) = 1/(2(Sqrt[λ])), the inversion points for the left well are, at lowest order in λ, at (Subscript[x, 1])=± (Sqrt[ϵ]). At lowest order in λ, the value for the parameter ϵ (the double of the "energy" for this oscillator) is 2n+1, for the n-th state.

In the forbidden region the equation () will have a WKB solution of the form

ψ = ((Subscript[C, 1])/(Sqrt[k[ξ]]))Exp[-(∫_(Subscript[x, 1])^ξ)k[x]] + ((Subscript[C, 2])/(Sqrt[k[ξ]]))Exp[(∫_(Subscript[x, 1])^ξ)k[x]];     k[x]∼(Sqrt[(x^2)-ϵ])

The quasi degenerate states come in doublets of odd/even states. For even/odd states we must have respectively:

even:  ψ'[(Subscript[ξ, 0])] = 0;    odd: ψ[(Subscript[ξ, 0])]=0

As k'∝ V' which vanishes at the symmetric point, the previous conditions amount to

(Subscript[C, 1]) Exp[-(∫_(Subscript[x, 1])^(Subscript[ξ, 0]))k[x]] ∓ (Subscript[C, 2])Exp[(∫_(Subscript[x, 1])^(Subscript[ξ, 0]))k[x]]= 0  ⇒ ((Subscript[C, 2])/(Subscript[C, 1]))= Exp[-2 (∫_(Subscript[x, 1])^(Subscript[ξ, 0]))k[x]] ≡ ± K

This condition will be a quantization condition for energy and will give us also the non perturbative correction to energy levels. The only problem is to find an accurate solution of the Schrödinger equation inside the well which can be matched with the asymptotic WKB solution () in the intermediate region   (Subscript[x, 1])≪ ξ ≪ 1/(Sqrt[λ]).

At lowest order in λ

-ψ'' + (ξ^2)ψ = ϵ ψ

The solution of this equation which vanish for ξ → - ∞ (outside the potential) is, apart a multiplicative constant:

ψ = (Subscript[D, ((ϵ-1)/2)])(-(Sqrt[2])ξ)

This solution is bounded also for ξ→+∞ for ϵ = 2n+1, and this would give the usual energy levels. In general at large positive ξ

-i (2^(-(1/4)+(ε/4))) (e^(((i π ε)/2)-((ξ^2)/2))) (Sqrt[(1/ξ)]) (ξ^(ε/2))+(((2^((1/4)-(ε/4))) (e^((ξ^2)/2)) (Sqrt[π]) (Sqrt[(1/ξ)]) (ξ^(-ε/2)))/(Gamma[((1-ε)/2)]))

For large ξ the classical action is easily found to be:

(∫_(Sqrt[ϵ])^ξ)k[x]dx≃ ((ξ^2)/2)+(1/4) (-ϵ-2 ϵ Log[2]+2 ϵ Log[(1/ξ)]+ϵ Log[ϵ]);  and  (Sqrt[ k[ξ]])∼ (ξ^(1/2))

then

(Subscript[ψ, WKB])∼ (Subscript[C, 1])(1/(Sqrt[ξ]))(e^(-((ξ^2)/2)))(e^(ϵ/4))(2^(ϵ/2))(ϵ^(-(ϵ/4)))(ξ^(ϵ/2)) + (Subscript[C, 2])(1/(Sqrt[ξ]))(e^(+((ξ^2)/2)))(e^(-(ϵ/4)))(2^(-(ϵ/2)))(ϵ^(ϵ/4))(ξ^(-(ϵ/2)))

Comparing with ()

(Subscript[C, 1])= (2^(-(1/4))) (e^((i π (ε-1))/2))(2^(-(ϵ/4)))(e^(+(ϵ/4)))(ϵ^(+(ϵ/4))) ;  (Subscript[C, 2])= (2^(1/4))(e^(+(ϵ/4)))(2^(+(ϵ/4)))(ϵ^(-(ϵ/4)))((Sqrt[π])/(Γ[((1-ϵ)/2)]))

and the quantization condition becomes:

((Sqrt[π])/(Γ[((1-ϵ)/2)]))(2^((ϵ+1)/2))(e^(+(ϵ/2)))(ϵ^(-(ϵ/2)))(e^(-((i π (ε-1))/2)))=± K

This equation can be solved by iteration. At lowest order K = 0 and Γ must have a pole, i.e. 1-ϵ = - 2 n ⇒ ϵ = 2n+1, the usual quantization condition. To compute the first correction let us write ϵ = 2n+1 + δ. Using the relation

(1/(Γ[-x]))=- (( Sin[π x])/π)Γ[1+x] ⇒ (1/(Γ[-(n+(δ/2))]))∼- (( 1)/π)Sin[n π+π(δ/2)]Γ[n+1+(δ/2)]∼ (((-1))^(n+1))n! (δ/2)

we can write (the imaginary exponential gives  a factor (((-1))^n) )

- (Sqrt[π]) (2^(n+1))(e^((n+(1/2))))(((2n+1))^(-(n +(1/2)))) n!(δ/2) =± K    ⇒    δ =∓ (2/(Sqrt[π]))(2^(-(1/2)))(e^(-(n+(1/2))))(((n+(1/2)))^(n+(1/2)))(1/(n!))K

In this system with our units Δϵ = 2δ =  (Sqrt[2])ΔE then ΔE = (Sqrt[2]) δ and

ΔE = (2/(Sqrt[π]))(e^(-(n+(1/2))))(((n+(1/2)))^(n+(1/2)))(1/(n!))K

This formula coincides with (), we remember that the classical frequency in this system is ω =(Sqrt[2]).

We stress some points.

The last point is especially interesting in a basic course in Quantum Mechanics as it give an example of another way to compute eigenvalues, and taking into account  exponentially small effects, as we have seen. We give an idea of the method and leave to the interested reader the task of filling the details.

First one has to write the subleading WKB term, this is not a problem apart the fact that the initial point in the integral for action cannot be taken (Subscript[x, 1])as the WKB is not uniform at turning points and divergencies arise in the integrals. This can be avoided by using another starting point, b > (Subscript[x, 1]). A possibility is to perform a λ expansion and write the whole action as

(∫_(Subscript[x, 1])^ξ)(Subscript[k, 0])[x]dx + (Sqrt[λ]) (∫_b^ξ)(Subscript[k, 1])[x]dx + λ (∫_b^ξ)(Subscript[k, 2])[x]dx

(Subscript[k, 0]) is the zero order momentum, the other term represent higher order in the expansion of the square root in λ and higher order in WKB, (Subscript[σ, 2]) at this order. If F is a primitive of one integral, the result is F[x] - F[b], F[b] is just a redefinition of the constants C, and F[x] can be used for the matching.

To do the matching we have now to use a wave function correct to order λ. This can be done with a method borrowed from uniform approximation technique. We know that a model

u''[z]+ (t - (z^2))u[z] = 0;

is equivalent to the original model,

ψ'' +(ϵ- (ξ^2)(((1-(Sqrt[λ])ξ))^2))ψ = 0;     ψ[ξ] = u[z]/(((z'[ξ]))^(1/2))

if

(((z'[ξ]))^2)(t - (z^2)) = (ϵ- (ξ^2)(((1-(Sqrt[λ])ξ))^2)) - (1/2){z;ξ};

{z; ξ} being the Schwarz derivative

{z;ξ}= ((z''')/(z')) - (3/2)(((((z'')/(z'))))^2)

At leading order (λ = 0), clearly t = ϵ and z = ξ.  It easily found a polynomial solution at order λ:

z = ξ -  (Sqrt[λ])(((ξ^2)/3)+((2 ϵ)/3)) + λ (-((ξ^3)/18)-((19 ξ ϵ)/36));    t = ϵ+(1/2) (1+3 (ϵ^2))λ

The value of t is interesting: you can compute with usual perturbation theory ϵ up to order λ, with the result

ϵ = (2n+1)-2(1+3 n+3 (n^2)) λ

Substituting in the previous relation one find t = 2 n + 1. Vice versa we know that the parameter t must be 2n+1 to assure overall integrable parabolic cylinder functions, then this method provide a computation of the perturbative correction to energy! .

The rest of the calculation follow the same way as zero order: we have to write the asymptotic expansion in z of the solution, always a parabolic cylinder function, then express z as a function of x and keep terms up to order λ. Finally we can compute the coefficients (Subscript[C, 1,2]) with the matching conditions and insert them in the quantization condition ().

WKB expansion

The procedure described up to now is based on perturbation theory but it can be easily generalized to a full fledged WKB analysis. This will also dispense us from asymptotic matching conditions and make contact with connection formulas described above.

The uniform approximation technique for ψ[x] tell us that the two equations

u''[z] + (t-(z^2))u[z] = 0;    ψ''[x] + (k^2)[x]ψ[x] = 0;

have equivalent solutions

ψ[x] = ((u[z[x]])/(Sqrt[z'[x]]))

as far as

(((z'[x]))^2)(t-(z^2)) = (k^2)[x] - (1/2){z;x};     {z;x}≡ ((z''')/(z')) - (3/2)(((((z'')/(z'))))^2)

The Schwarz derivative is suppressed to order ℏ and can be neglected in a zeroth order WKB analysis, this is equivalent to zeroth order perturbation theory in the previous analysis.

In this approximation, see notebook [NB-11.9.nb], in the classical allowed region

(∫_(-(Sqrt[t]))^z)q[z]= (∫_(-(Sqrt[t]))^z)(Sqrt[t-(z^2)])dz= (∫_(Subscript[x, 1])^x)k[x]dx;    t = 2 J =(2/π) (∫_(Subscript[x, 1])^(Subscript[x, 2]))k[x]dx.

(Subscript[x, 1]) and (Subscript[x, 2])are the classical turning points in the well, which correspond to z = ∓(Sqrt[t]). Let us suppose that we are working in the left well of a double well symmetric potential. Let (Subscript[x, 0])the symmetric point (the local maximum for the anharmonic oscillator) and (Subscript[z, 0])the corresponding value for z. In the classically forbidden region eq.() reads

(∫_(Sqrt[t])^z)(Sqrt[(z^2)-t])dz= (∫_(Subscript[x, 2])^x)k[x]dx ≡ σ;

Using, with ζ = (z^2)/2,

(∫_(Sqrt[t])^z)(Sqrt[(z^2)-t])= (1/2) (z (Sqrt[(z^2)-t])+t Log[((Sqrt[t])/(z+(Sqrt[(z^2)-t])))])=ζ (Sqrt[1-(t/2)(1/ζ)])- (t/4)Log[((2ζ)/t)]- (t/2)Log[1+ (Sqrt[1-(t/(2 ζ))])]

By iteration it is easily found the asymptotic behavior for large σ

ζ = σ + (1/4)(t- t Log[(t/(8σ))]) + 𝒪(Log[σ]/σ)

The solution of the equation for u[z], finite at z→-∞

u[z] = (Subscript[D, ((t-1)/2)])[-(Sqrt[2])z]

This is accurate up to order ℏ up to z = (Subscript[z, 0]), the corresponding solution for ψ[x] is given by (). We can prolong the solution in a even/odd way in the region x>(Subscript[x, 0]) and this will be a solution of the Schrödinger equation provided that even/odd prolongations have respectively zero derivative at x=(Subscript[x, 0])and are continuous. z[x] is constructed by reflection in the right half space. At x=(Subscript[x, 0]) from eq.()

z'[(Subscript[x, 0])]∼ ((k[(Subscript[x, 0])])/z);      z'[(Subscript[x, 0])]z''[(Subscript[x, 0])]∼ (((((k^2)))')/(Subscript[z, 0]^4)) - ((k^2)/(Subscript[z, 0]^5)) ∼ ((V'[(Subscript[x, 0])])/(Subscript[z, 0]^4))- ((k^2)/(Subscript[z, 0]^5))

at the symmetric point V'[(Subscript[x, 0])]=0, then z''[(Subscript[x, 0])]∼ k/(z^4)which goes rapidly to zero, so the prefactor with z'term in ψ will not play any role in the following, we can directly work with u[z].

For large z we can use the asymptotic expansion for parabolic cylinder functions (), now t plays the role of ϵ and the variable is z:

-i (2^(-(1/4)+(t/4))) (e^(((i π t)/2)-((z^2)/2))) (Sqrt[(1/z)]) (z^(t/2))+(((2^((1/4)-(t/4))) (e^((z^2)/2)) (Sqrt[π]) (Sqrt[(1/z)]) (z^(-t/2)))/(Gamma[((1-t)/2)]))

To have a consistent WKB solution this expression must vanish at z=(Subscript[z, 0])for odd states or its first derivative must be zero for even states. The leading term in the derivative come from the exponential, then, as in the previous perturbative case the constraint for even/odd states is

(((2^((1/4)-(t/4))) (e^((z^2)/2)) (Sqrt[π])  (z^(-t/2)))/(Gamma[((1-t)/2)])) = ± (2^(-(1/4)+(t/4))) (e^(i π(( t-1)/2)-((z^2)/2)))  (z^(t/2));   for z = (Subscript[z, 0])

With t = 2 J and using

(1/(Γ[((1-t)/2)]))= (1/(Γ[(1/2)-J]))=(1/(Γ[1-((1/2)+J)]))= (1/π)Sin[π(J+(1/2))]Γ[J+(1/2)]= ((Cos[π J])/π)Γ[J+(1/2)];
(e^(i π(( t-1)/2)))= Cos[π (J-(1/2))] + i Sin[π (J-(1/2))]= Sin[πJ]-i Cos[π J];

we have

((Cos[π J])/π)Γ[J+(1/2)](Sqrt[π]) = ± (2^(J-(1/2))) (Sin[πJ]-iCos[π J]) (e^(- (Subscript[z, 0]^2))) (Subscript[z, 0]^(2J))

This equation has imaginary terms, this means that J, and by consequence the energy, have imaginary contributions. For J = n+1/2, the leading order, and the only result at zeroth order WKB, the imaginary term cancels, and for now we omit this term, and we look for real solutions for J:

((Cos[π J])/π)Γ[J+(1/2)](Sqrt[π])=± (2^(J-(1/2))) Sin[πJ](e^(- (Subscript[z, 0]^2))) (Subscript[z, 0]^(2J))

Substituting for (Subscript[z, 0])its asymptotic value in term of the action, eq.() ( σ is the action from (Subscript[x, 2])to (Subscript[x, 0])):

(z^2)= 2ζ= 2σ + (1/2)(2J- 2J Log[((2J)/(8σ))])= 2σ + J(1+Log[((4σ)/J)])

we obtain at the leading order

((Cos[π J])/π)Γ[J+(1/2)](Sqrt[π])=± (2^(J-(1/2))) Sin[πJ] (e^(- 2σ -J))((((J/(4σ))))^J) (((2σ))^J)=

± (2^(-(1/2))) Sin[πJ] (e^(- 2σ ))((((J/e)))^J)

i.e.

((((2/π)))^(1/2))(((((e )/J)))^J)Γ[J+(1/2)]((Cos[π J])/(Sin[π J]))= ± (e^(- 2σ )) ≡ ± K

which is exactly, as expected, the formula () obtained with connection formulas. J is a function of energy and we know from classical mechanics

(dE/dJ) = ω ⇒ J ≃ n+(1/2) + (dJ/dE)δE = n+(1/2) + ((δE)/ω)

Expanding () for small δE we recover the general formula ()

π(( δE)/ω) = ∓ ((((π/2)))^(1/2))(((((n+(1/2))/e)))^(n+(1/2)))(1/(n!))K;   ΔE = 2 δE

Connection with the path integral approach

The derivation of energy splitting through path integral and its connection with usual  Schrödinger equation is very clearly exposed in the paper of Coleman [1]. Here we show how our result can be translated in that language.

In the instanton approach the path integral is evaluated by a saddle point technique around classical "euclidean" solutions of the equation of motion. The result is expressed in the form

ΔE = 2(ω/(Sqrt[π]))C Exp[- 2 (Subscript[S, 0])]

where

ω is the parameter in V ∼ 1/2 (ω^2)(x^2) around a classical minimum

(Subscript[S, 0]) is the classical action, at zero energy, between the classical minimum and the symmetric point (Subscript[x, 0])

C is a constant which can be determined by computing fluctuations around classical solutions or by the asymptotic form of the integral

(∫_(Subscript[x, b])^(Subscript[x, 0]))((dx)/(Sqrt[2V[x]]))=- (( 1)/ω)Log[C(Sqrt[ω]) (Subscript[x, b])] + 𝒪((Subscript[x, b]))

The logarithmic singularity at (Subscript[x, b])= 0 (a classical minimum of V) is due to the form of the potential 1/2 (ω^2)(x^2). The factors ω have been introduced for convenience inside the logarithm, as due to integration outside. The dependence on ω is the same as that of an oscillator, as will be evident below.

The assumption behind the computation is a strict ℏ→0 limit. In this limit the quantized energy E, computed inside a single well, usually goes to 0 then it is reasonable to consider as leading term the action at zero energy. The first correction in ℏ will include the actual energy. We will consider the ground state of the system, and its partner in the doublet created by tunneling.

We note that this limit is equivalent to weak coupling limit in all potentials of the form V[g x], as discussed above. We assume to be in this framework, and this is also the approximation used in the paper of Coleman quoted above. In this approximation the energy of the ground state, without tunneling effects, is

E = ((ℏ ω)/2)

All we have to do to show the equivalence is to show that the first order expansion in E of our general expression ()

Δ E =  (Subscript[E, n]^odd)- (Subscript[E, n]^even) =2 ((((π/2)))^(1/2))((((1/(2e))))^(1/2)) (ω/π)K = (Sqrt[(π/e)])(ω/π) Exp[-2 S[E]];
 S[E] = (∫_(Subscript[x, 1])^(Subscript[x, 0]))(Sqrt[2V[x]-2E]) ≡ σ[(Subscript[x, 1]),(Subscript[x, 0]);E],

reproduces (). (Subscript[x, 0])is the symmetric point of the potential and (Subscript[x, 1])the classical inversion point.

The expansion in E is a bit delicate as S[E] is not analytic as E→0, due to divergence of its first derivative. Let us choose a fixed small value to break the integral, this factor will disappear at the end of the calculation. We can write

σ[(Subscript[x, 1]),(Subscript[x, 0]);E] = σ[(Subscript[x, 1]),(Subscript[x, b]);E]+ σ[(Subscript[x, b,])(Subscript[x, 0]);E]

The second integral can be easily expanded in E. At the first order and using the definition of the constant C in ()

σ[(Subscript[x, b,])(Subscript[x, 0]);E]= σ[(Subscript[x, b,])(Subscript[x, 0]);0] - E (∫_(Subscript[x, b])^(Subscript[x, 0]))((dx)/(Sqrt[2V[x]])) ≃ σ[(Subscript[x, b,])(Subscript[x, 0]);0] + (E/ω)Log[(1/C)ω (Subscript[x, b])]

For the first integral we can use the form of the potential at small x. The inversion classical point is (Subscript[x, 1])= (Sqrt[2E])/ω and

σ[(Subscript[x, 1]),(Subscript[x, b]);E]= (∫_(Subscript[x, 1])^(Subscript[x, b]))(Sqrt[(ω^2)(x^2)- 2 E])dx = ((2E)/ω)(∫_1^((Subscript[x, b])ω/(Sqrt[2E])))(Sqrt[(z^2)-1])dz = 
=((2E)/ω)(1/2) (α (Sqrt[-1+(α^2)])-Log[α+(Sqrt[-1+(α^2)])]);    α = (((Subscript[x, b])ω)/(Sqrt[2E]))

The large α expansion of the integral can be done by iteration, or by using Mathematica with the result, always at first order in E:

((α^2)/2)+(1/4) (-1-2 Log[2]+2 Log[(1/α)])

and for σ

σ[(Subscript[x, 1]),(Subscript[x, b]);E] = (ω^2)((Subscript[x, b]^2)/2)- (E/ω)Log[(Subscript[x, b])]+ ((2E)/ω)(1/4)(-1 -2 Log[2] + Log[2(E/(ω^2))])

The first term is just the zero energy action from (Subscript[x, 1])to (Subscript[x, b]). Summing with previous contribution the logarithmic divergence in (Subscript[x, b]) cancels and using the explicit value for E, ω/2, we have

σ[(Subscript[x, b,])(Subscript[x, 0]);E]= σ[(Subscript[x, b,])(Subscript[x, 0]);0]+(1/2)Log[(1/C)(ω^(1/2)) ]+ (1/4)(-1 -2 Log[2] + Log[(1/ω)]).

Then

K = Exp[-2 σ]= Exp[-2 (Subscript[S, 0])] C 2(Sqrt[e])

and

ΔE=(Sqrt[(π/e)])(ω/π) Exp[-2 S[E]] = 2C Exp[-2 (Subscript[S, 0])] (ω/(Sqrt[π]))

which is the instanton result ().

References

1) S. Coleman, The Uses of Instantons, Erice Lectures 1977, In The Whys of Subnuclear Physics (Plenum Publishing Company N.Y. 1979), and in Aspects of symmetry: Selected Erice Lectures, Cambridge University Press, 1988.

2) E. Brezin, G. Parisi, J. Zinn-Justin, Perturbation theory at large orders for a potential with degenerate minima, Phys. Rev. D 16, 408 (1977).

3) J. Zinn-Justin, Expansion around instantons in quantum mechanics, J. Math. Phys. 22, 511 (1981).

4) E.B. Bogomolny, Calculation of Instanton-anti-Instanton contributions in Quantum Mechanics, Phys. Lett. 91 B, 431 (1980).

Problem

Compute in the limit g → 0 the imaginary part Im[E] for a potential

V[x] = (1/2)(x^2)- (1/2)g (x^4) ;   g > 0.

Solution

Statement of the problem

The potential has the form

V[x] = (1/2)(x^2)- (1/2)g (x^4) ;   g > 0.

The classical inversion points a and b are given by

a= (Sqrt[((1-(Sqrt[1-8 E g]))/(2g))]);  b= (Sqrt[((1+(Sqrt[1-8 E g]))/(2g))]).

For small enough coupling the system can have metastable states, which decay by tunnel effect through the barrier. At the semiclassical level the width is given by

Γ = ξ (2/T) Exp[-(2/ℏ) S[a,b]]

where T is the classical period of motion inside the well, S is the classical action and ξ is a factor of order unit. Let us note the factor 2, which classically is due to the frequency of hits against the barriers: 2/T as in one period the particle hits two barriers. In this notebook we will compute exactly the prefactor ξ.

The width Γ is related to the imaginary part of the eigenvalues with

E = (Subscript[E, R]) - i (Γ/2);     Γ = - 2 Im[E];

and we will compute Γ in this form.

This problem has one important aspect. Starting from a positive coupled oscillator

V[x] = (1/2)(x^2)+ (1/2)λ (x^4) ;

we have seen in the study of large orders in perturbation theory, that the leading behavior of perturbation coefficients is fixed by the imaginary part of the eigenvalues of () at negative λ, what we have to compute in this problem (some details will be given at the end of the problem).

In the text it has been shown that in the semiclassical approximation a general formula for Im[E] can be written. Let us briefly report the argument to stress some particular points which will be important in the following.

The Schrödinger equation is:

-(( (ℏ^2))/(2m))ψ''[x] + V[x]ψ[x] = E ψ[x]

We assume that the potential is symmetric. Depending on the study of symmetric or antisymmetric states the equation () is equipped with the boundary conditions

even:  ψ'[0] = 0;    odd:    ψ[0] = 0.

If it is possible to impose the condition ψ∈ (𝕃^2) we know that eigenvalues are real. For large x the potential, in our case goes like - (x^4)then at large distances we have always an oscillatory regime (the region is a classical allowed region) and the solution cannot be normalized. There is the possibility of looking for metastable states, using the Gamow - Siegert approach, i.e. we look for states which behave as divergent waves as |x| grows, this automatically imply imaginary eigenvalues, as diverging waves are complex and for E real the solutions of () are real.

Multiplying  () by (SuperStar[ψ])and subtracting from it its complex conjugate one have

2i Im[E] ((|ψ|)^2) = - ((ℏ^2)/(2m))((SuperStar[ψ])ψ''- ψ (ψ^(′′* )))=- ((ℏ^2)/(2m)) (d/dx)((SuperStar[ψ]) ψ'- ψ (ψ^(′* )))

We can integrate from 0 to x and use boundary conditions (). In both cases we obtain

2 Im[E](∫_0^x) ((|ψ|)^2)dx =- (1/i) ((ℏ^2)/(2m)) ((SuperStar[ψ]) ψ'- ψ (ψ^(′* )))

With  ψ = |ψ| Exp[i θ] we have

2 Im[E](∫_0^x) ((|ψ|)^2)dx = -  ((ℏ^2)/(2m)) 2 ((dθ)/dx)

This is an important point : if the phase grows as x gets large then Im[E] is negative, otherwise it will be positive. In the Gamow Siegert approach metastable states are identified with diverging waves, then θ'> 0 and Im[E] < 0. The other possibility will be explored at the end of the problem.

In a problem like () we can choose b ≪ x and approximate ψ with WKB.

One write

E = (Subscript[E, 1])+ i (Subscript[E, 2]);  ψ = (Subscript[ψ, 1])+ i (Subscript[ψ, 2]) ;      (Subscript[ψ, 1]),(Subscript[ψ, 2])∈ ℝ

Neglecting tunneling (Subscript[E, 2])= 0, (Subscript[ψ, 2])= 0. In usual case the continuation of ψ in the forbidden region give rise to

(C/(2(Sqrt[|p|])))Exp[- σ(a,x)];      σ = |S|/ℏ,  |p[x]| =(Sqrt[2m(V[x]-E)])

The presence of an imaginary part produce an exponential growing term, which we rewrite in the form

ψ ∼ (C/(2(Sqrt[|p|])))Exp[-σ[a,b]+σ[r,b]] + i (D/(Sqrt[|p|]))Exp[- σ[r,b]]

For r > b the usual connections formula will give

ψ ∼ (1/(Sqrt[|p|]))(-(C/2)(e^(-σ[a,b]))Sin[w[b,r]-(π/4)] + i D 2 Cos[w[b,r]-(π/4)])

The crucial choice is

We look at diverging waves, then we require a solution Exp[i w[b,r] ]. This imply  4D = C Exp[-σ[a,b]] and

ψ ∼ i (C/(2(Sqrt[|p|])))Exp[-σ[a,b]]Exp[i (w[b,r] - (π/4))]

Substituting in () and using w'= p/ℏ one has

2 Im[E](∫_0^b)((|ψ|)^2)dx =- (1/i) ((ℏ^2)/(2m))(e^(-2 σ[a,b]))(C^2)2 i w'=- (ℏ/(4m))(C^2)(e^(-2 σ[a,b]))

We look at converging waves, then we require a solution Exp[-i w[b,r] ]. This imply  4D = -C Exp[-σ[a,b]] and with the same manipulations we have

2 Im[E] =- (1/i) ((ℏ^2)/(2m))(e^(-2 σ[a,b]))(C^2)2 (-i w')=+ (ℏ/(4m))(C^2)(e^(-2 σ[a,b]))

We see that the sign of Im[E] depends on boundary conditions at infinity.

The anharmonic oscillator

To fix the constant C we have to connect WKB approximation with an accurate solution for small x. This is easily done at weak coupling. For small g

a ≃ (Sqrt[2E])(1+ E g);  b = (1/(Sqrt[g]))(1- E g).

Then the second inversion point tends to infinity.

We make the matching in the intermediate region a ≪ x ≪ b. Using

2(-E + V) = g((b^2)-(x^2))((x^2)- (a^2))

and the small g limiting forms for a and b, we can write, in this region:

σ[a,x]= ((Sqrt[g])/ℏ)(∫_a^x)(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = ((Sqrt[g])/ℏ)(a^(3/2))(∫_1^(x/a))(Sqrt[((z^2)- 1)(((b^2)/(a^2))-(z^2))])dz ∼
((Sqrt[g])/ℏ)(a^(3/2)) (1/(Sqrt[2E g]))(∫_1^(x/(Sqrt[2E])))(Sqrt[((z^2)- 1)])dz = 2(E/ℏ)(1/2)(Subscript[|z(Sqrt[(z^2)-1])- Log[z + (Sqrt[(z^2)-1])]|, 1]^(x/(Sqrt[2E])))∼
(E/ℏ)(Subscript[((z^2)- (1/2) - Log[2z]), z=x/(Sqrt[2E])])= (E/ℏ)( ((x^2)/(2E))-(1/2)- Log[((2x)/(Sqrt[2E]))])

For the momentum, in the same region (a ≪ x ≪ b )

p[x] ∼ (Sqrt[g])(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]) ∼ x

and finally

ψ[x] ∼ (C/2)((((1/x)))^(1/2)) Exp[- ((x^2)/(2ℏ))+ (E/(2ℏ))](((((2x)/(Sqrt[2E]))))^(E/ℏ))

At lowest order in g, E/ℏ = n + 1/2 and we have, for the n-th state

(Subscript[ψ, n])[x] ∼ (C/(Sqrt[2]))(((2E))^(-1/4))(((((2x)/(Sqrt[2E]))))^n)Exp[- ((x^2)/(2ℏ))](e^((n/2)+(1/4)));      (Subscript[ψ, 0])= (C/(Sqrt[2]))Exp[- ((x^2)/(2ℏ))](e^(1/4)) .

For x ∼ 0 and small g the normalized wave function is given by (we have a harmonic oscillator with ω = 1)

(Subscript[ψ, n])[x] =(π^(-1/4)) (1/(Sqrt[(2^n)n!]))Exp[-(x^2)/2ℏ](Subscript[H, n])[x];  (Subscript[ψ, 0])[x] =(π^(-1/4))Exp[-(x^2)/2ℏ]

For the ground state we obtain, by matching () and ()

(C^2)= (2/(Sqrt[e π])).

We can insert this expression in (). In the left hand side the integral is relevant only in the allowed region, and it is equal to half the normalization, as it starts from 0, then (we put ℏ = 1, m = 1)

Im[E] = - (1/4)(2/(Sqrt[e π]))(e^(-2 σ[a,b]))= - (Sqrt[(π/e)]) (1/(2π))(e^(-2 σ[a,b])).

In our units ω = 1, then T = 2 π and in general previous formula can be written

Im[E] = - (Sqrt[(π/e)]) (ℏ/T)(e^(-2 σ[a,b]))

To compute σ[a, b] the simplest thing is to write its exact expression in terms of elliptic integrals (for a proof see notebook [NB-11.3] )

σ[a,b] = (∫_a^b)|p[x]|dx = (Sqrt[g])(1/3) b (((a^2)+(b^2)) EllipticE[1-((a^2)/(b^2))]-2 (a^2) EllipticK[1-((a^2)/(b^2))])

The small g expansion is (for the ground state)

σ ∼(1/(3 g))-(1/4)- Log[2]+(1/4) Log[g])

and substituting in Im[E]

Im[E] = - (Sqrt[(π/e)]) (1/(2π)) (e^(-(2/(3g))))(Sqrt[e]) 4 (1/(Sqrt[g])) = -(2/(Sqrt[π]))(1/(Sqrt[g]))(e^(-(2/(3g))))

and for Γ

Γ = (4/(Sqrt[π]))(1/(Sqrt[g]))(e^(-(2/(3g))))

For large n, using the definition of Hermite polynomials

(Subscript[H, n])[x]= (((-1))^n) (e^(x^2))((d^n)/(dx^n))(e^(-(x^2)))  ⇒ (2^n)(x^n)

and the reader can easily show, using Stirling approximation, that

(Subscript[C, n]^2) → (2/π) ;      Im[E]→ (1/(2π))Exp[-2σ] → (ℏ/T)Exp[-2σ]  ⇒   Γ → ((2 ℏ)/T)Exp[-2σ]

i.e. the usual WKB approximation.

Analytic continuation

Let us now discuss briefly the problem of analytic continuation in g of a usual harmonic oscillator:

V[x] = (1/2)(x^2)+ (1/2)g (x^4) ;

At large positive x the bounded solution has a behavior (from WKB or just by taking only leading terms in the equation)

ψ[x] ∼ Exp[- (∫^x)(Sqrt[2 g (x^4)])] ∼ Exp[-((x^3)/3)(Sqrt[2 g])]

With the analytic continuation g → Exp[i π] |g| the wave function transforms in

ψ ∼ Exp[-i((x^3)/3)(Sqrt[2 |g|])]

This is a convergent wave, and as discussed above this boundary condition gives a positive imaginary part for E, i.e. in the analytic continuation of the normalizable problem:

Im[E] = + (2/(Sqrt[π]))(1/(Sqrt[g]))(e^(-(2/(3g)))).

Problem

Compute exactly the reflection and the transmission coefficients for a potential V[x] = - F |x|.

Solution

Scattering below the barrier

Positive and negative energies correspond respectively to scattering above and below the barrier. Let us consider first the second case, E <0 , and put E = - ϵ. The Schrödinger equation is

((d^2)/(dx^2))ψ + ((2m)/(ℏ^2))( F |x| - ϵ)ψ = 0

with classical turning points

x = ± a = ± (ϵ/F)

It is convenient to introduce the parameter

λ = (((((2m F)/(ℏ^2))))^(1/3))

and rewrite eq.() in the form

(1/(λ^2))((d^2)/(dx^2))ψ - λ( - x+a)ψ = 0;   x > 0
(1/(λ^2))((d^2)/(dx^2))ψ - λ( x+a)ψ = 0;       x < 0

The wave number (or the momentum) is given by

((|k|)^2) = (((|p|)^2)/(ℏ^2)) = (λ^3)|x∓a|

In the classical allowed regions the behavior of ψ follows by semiclassical WKB approximation

ψ[x] ∼ (1/(Sqrt[p[x]]))Exp[± i (∫^x)p[x]𝕕x ]∼ (1/((|x∓a|)^(1/4)))Exp[± i (2/3)(λ^(3/2))((|x∓a|)^(3/2))]

Right movers waves correspond to a growing phase, i.e. to a positive sign in the exponent () for positive x. Let us note that if we put r = |x| then in both asymptotic regime the argument is written |r - a|.

We want to describe an incident wave scattered by the potential, this give rise to a transmitted wave and a reflected wave, then the solution we are looking for has the asymptotic behavior

A ((Exp[- i (2/3)(λ^(3/2))((|x+a|)^(3/2))])/((|x-a|)^(1/4)))+ B ((Exp[+ i (2/3)(λ^(3/2))((|x+a|)^(3/2))])/((|x-a|)^(1/4)))  as x→-∞
C ((Exp[ i (2/3)(λ^(3/2))((|x-a|)^(3/2))])/((|x-a|)^(1/4)))   as x → + ∞

Let us stress that in eq.() the incoming wave in the region x < 0 is the one with negative exponent, i.e. the one proportional to A. The phase grows form large negative values as x ≪ -a to 0 when x = - a. In the usual case the progressive wave is written Exp[i k x] but for x < 0 the phase is negative, as in the present case.

The transmission and reflection coefficients are

T = (((|C|)^2)/((|A|)^2));       R = (((|B|)^2)/((|A|)^2));

and must satisfy the unitarity constraint

T + R = 1.

In the domains x < 0 and x > 0 the solutions are simply Airy functions. Airy functions satisfy

u''[z] - z u = 0

then the general solution of () is

a)         (Subscript[b, 1]) Ai[-λ(x-a)] +  (Subscript[b, 2])Bi[-λ(x-a)]  ;     x > 0   
b)         (Subscript[c, 1]) Ai[λ(x+a)]    + (Subscript[c, 2])Bi[λ(x+a)]  ;        x < 0

Let us recall the asymptotic behavior of Airy functions, the reader can easily check it with Mathematica. With r = |x| for large r:

Ai[r] ≃ ((e^(-((2 (r^(3/2)))/3)))/(2 (Sqrt[π]) (r^(1/4))));    Ai[-r]≃ ((Sin[(π/4)+((2 (r^(3/2)))/3)])/((Sqrt[π]) (r^(1/4))));
Bi[r]≃ ((e^((2 (r^(3/2)))/3))/((Sqrt[π]) (r^(1/4))));         Bi[-r]≃ ((Cos[(π/4)+((2 (r^(3/2)))/3)])/((Sqrt[π]) (r^(1/4)))).

It is then convenient to introduce the combinations

(Subscript[E, ±])[x] = Bi[x] ±i Ai[x]

From () it follows

(Subscript[E, ±])[-r] =(1/((Sqrt[π]) (r^(1/4))))Exp[± i ((π/4)+((2 (r^(3/2)))/3))]     as r→ ∞ .

The solution satisfying boundary conditions is then of the form

(SubPlus[E])[-λ(x-a)] ;                                       x > 0;
A (SubMinus[E])[λ(x+a)] + B (SubPlus[E])[λ(x+a)];     x < 0

The notation is the same as in eq.() with C = 1, the reader has to remember the observation after eq.() to check the correspondence.

The coefficients A and B are fixed by continuity of ψ and its derivative at the origin, x= 0.

A (SubMinus[E])[λ a] + B (SubPlus[E])[λ a] = (SubPlus[E])[λ a];       A (SubMinus[E'])[λ a] + B (SubPlus[E'])[λ a] = - (SubPlus[E'])[λ a];

The solution can be written using Cramer's rule. The determinant of the coefficients is what is usually called the Wronskian for two functions

W[f,g] = f g' - f' g

and we have

det = W[(SubMinus[E]), (SubPlus[E])] = - 2 i W[Ai, Bi]

The Wronskian for the two solution of a linear differential equation is a constant, as can be easily checked, then can be easily computed or from asymptotic expansion or from Taylor expansion around origin. For example as x→ 0

Ai[x]∼ (Subscript[c, 1])- (Subscript[c, 2])x;  Bi[x]∼ (Sqrt[3])((Subscript[c, 1])+ (Subscript[c, 2])x);  (Subscript[c, 1])= (1/((3^(2/3))Γ[(2/3)]));  (Subscript[c, 2])=  (1/((3^(1/3))Γ[(1/3)]))

from which it follows W[Ai, Bi]  = 1/π and

det = - 2i (1/π)

The solutions of () are then:

A = i π (Subscript[((SubPlus[E]) (SubPlus[E'])), λa]) ;  B = - i (π/2)(Subscript[((SubMinus[E]) (SubPlus[E'])+ (SubPlus[E])(SubMinus[E'])), λa])

By using () it is easy to check unitarity (we have C = 1 in eq.() ):

((|B|)^2) - ((|A|)^2) = ((π^2)/4)(W^2) = - 1 ⇒  (1/((|A|)^2))+ (((|B|)^2)/((|A|)^2)) = 1.

From () we have (all functions are computed with argument λ a):

((|A|)^2) = (π^2)((Ai^2)+ (Bi^2))( (((Ai'))^2)+ (((Bi'))^2));    ((|B|)^2) = (π^2)(((Ai Ai'+ Bi Bi'))^2) .

Asymptotic form for λa ≫ 1

From asymptotic expansions () we can give an asymptotic estimate, for large λa, for these factors

(Ai^2)∼ (1/(Sqrt[λa]))((e^(-((4 (((λa))^(3/2)))/3)))/(4π)) ; (((Ai'))^2) ∼ (Sqrt[λa])((e^(-((4 (((λa))^(3/2)))/3)))/(4π));

(Bi^2)∼ (1/(Sqrt[λa]))((e^(+((4 (((λa))^(3/2)))/3)))/π) ;(((Bi'))^2) ∼ (Sqrt[λa])((e^(+((4 (((λa))^(3/2)))/3)))/π);

with

K = Exp[-((8 (((λa))^(3/2)))/3)]

we have

((|A|)^2) = (((((Sqrt[K])/4)+ (1/(Sqrt[K]))))^2) = ((1/K) + (1/2)+ ((K^2)/16)) (1 + 𝒪( (((λa))^(-3/2)) )

The notation means that each term is the leading term of an asymptotic expansion which starts with a (((λa))^(-3/2)) correction. Clearly only the first term is dominant for K ≪ 1.

K is the transmission factor for the barrier, indeed:

Exp[((- 2)/ℏ)(∫_(-a)^a) 𝕕x |p[x]| = Exp[- 4 (λ^(3/2))(∫_0^a) 𝕕x |x-a|] = Exp[-((8 (((λa))^(3/2)))/3)] = K

The leading order of the transmission coefficient is in agreement with usual lowest order WKB result:

T = (1/((|A|)^2))= (1/((1/K) + (1/2)+ ((K^2)/16)))≃ K

Asymptotic form for λa ≪ 1

This regime correspond to scattering for energies just below the barrier' s top. For small values of the argument we can use the expansion () to get:

((|A|)^2)∼ (π^2)(4 (Subscript[c, 1]^2))(4 (Subscript[c, 2]^2))= (4/3);  ((|B|)^2) ∼ (π^2)(((2 (Subscript[c, 1])(Subscript[c, 2])))^2) = (1/3).

In particular we see that the reflection coefficient is different from 0,  R = 1/4 .

Lowest order WKB

Let us see which result is obtained by a blind application of lowest order connection formulas of WKB, which are, we remember:

Cos[|S|+μ-(π/4)] ↔ ( (1/2)Cos[μ] (e^(-|S|))+ Sin[μ] (e^(+|S|)))

or

Cos[|S|+(π/4)] ↔  (e^(+|S|));   Sin[|S|+(π/4)]↔ (1/2) (e^(-|S|))

Starting with a progressive wave in the region x ≫ a and going through the two turning points:

Cos[S[a,x] + (π/4)] + i Sin[S[a,x] + (π/4)]→ (e^(S[x,a])) + (i/2)(e^(-S[x,a]))= 
Exp[S[-a,a]-S[-a,x]]+  (i/2)Exp[-S[-a,a] + S[-a,x]] → 
(1/(Sqrt[K])) 2 Sin[S[x,-a]+(π/4)] + (i/2)(Sqrt[K]) Cos[S[x,-a]+(π/4)]

With  ϕ = S[x,-a]+(π/4) we have the result

i (e^(-i ϕ)) ((1/(Sqrt[K]))+(( (Sqrt[K]))/4))- i (e^(i ϕ))((1/(Sqrt[K]))-(( (Sqrt[K]))/4))

Exp[i ϕ] is the progressive way so we have

((|A|)^2)= ((((1/(Sqrt[K]))-(( (Sqrt[K]))/4)))^2)= (1/K)+ (K/16) - (1/2); ((|B|)^2)= ((((1/(Sqrt[K]))+(( (Sqrt[K]))/4)))^2)= (1/K)+ (K/16) + (1/2);

This result is correct only at the leading order, as we can also check from the break of unitarity:

(1/((|A|)^2)) + (((|B|)^2)/((|A|)^2)) ≠ 1

This is a general feature of WKB at lowest order : only lowest order terms, which come from dominant asymptotic terms, turn out to be correct, subleading terms are affected by subdominant terms and are out of control in this approximation. We stress that the leading non exponential term in () is correct instead.

Scattering above the barrier

Almost nothing must be changed in the above solution, simply a → - a.  Now a = E/F. What changes is the asymptotic behaviors. For small energies λa → 0 and we have again R = 1/4, T = 3/4. For higher energies one has to use () but now for negative values of the argument, i.e. in the oscillator regime. We have to take next to leading terms in the expansion of Ai and Bi, with r = λ a

Ai ∼ -((5 Cos[(π/4)+((2 (r^(3/2)))/3)])/(48 (Sqrt[π]) (r^(7/4))))+((Sin[(π/4)+((2 (r^(3/2)))/3)])/((Sqrt[π]) (r^(1/4))));  Ai'∼ ((7 Cos[(π/4)-((2 (r^(3/2)))/3)])/(48 (Sqrt[π]) (r^(5/4))))-(((r^(1/4)) Sin[(π/4)-((2 (r^(3/2)))/3)])/(Sqrt[π]));

Bi∼ ((Cos[(π/4)+((2 (r^(3/2)))/3)])/((Sqrt[π]) (r^(1/4))))+((5 Sin[(π/4)+((2 (r^(3/2)))/3)])/(48 (Sqrt[π]) (r^(7/4)))); Bi'∼ (((r^(1/4)) Cos[(π/4)-((2 (r^(3/2)))/3)])/(Sqrt[π]))+((7 Sin[(π/4)-((2 (r^(3/2)))/3)])/(48 (Sqrt[π]) (r^(5/4)))).

These formulas allow a computation of ((|B|)^2), which vanish at leading order. For ((|A|)^2) we need the next order in the expansion. For simplicity we will write only reflection coefficient, and extract T from unitarity. The reader can check that this is actually true using for example Mathematica to compute the expansions to next order:

R ∼ (1/(16 (((λ a))^3)))=  (((F^2) (ℏ^2))/(32 m (ϵ^3)))= (((m^2) (ℏ^2))/(4  (p^3)))(F^2);      T = 1 - R

where p is the momentum of the particle.

We note that in the classical limit, ℏ → 0, R → 0, as expected. This coefficient cannot be computed by standard WKB transition formulas as there are not turning points.

Problem

Compute exactly the reflection and the transmission coefficients for a potential V[x] = - (1/2) β (x^2).

Solution

The boundary conditions

We want to describe the scattering process for a particle coming from x = - ∞ impinging on a parabolic barrier. We first consider the case of energies below the barrier's top.  With V = - 1/2 β (x^2) this means energies E ≤ 0. With E = - ϵ the  Schrödinger equation is

((d^2)/(dx^2))ψ + ((2m)/(ℏ^2))( (1/2)β (x^2) - ϵ)ψ = 0 .

To simplify some formal manipulations below we make the change of variables

x = λ z;    λ =((((( m β)/(ℏ^2))))^(-1/4))(1/(Sqrt[2]));    α = (ϵ/ℏ)(Sqrt[(m/β)]) .

The equation () takes the form

((d^2)/(dz^2))ψ + ( (1/4) (z^2) - α)ψ = 0 .

The effective wave number is

(q^2)[z] = ( (1/4) (z^2) - α)

For large distances the WKB solutions have the form

Exp[± i (∫^z)q[z] 𝕕z] ∼ Exp[ ± i ((z^2)/4)]

Positive exponents for z > 0 describe wave propagating to the right, i.e. scattered waves in our problems. For z < 0 the phase grows for negative exponent, i.e. the incident wave has the form Exp[- i (z^2)/4]. We look therefore to a solution with

ψ ∼ Exp [i ((z^2)/4)];  as z → ∞;  
ψ ∼ A  Exp [-i ((z^2)/4)]+ B Exp [i ((z^2)/4)]; as  z → - ∞;

We repeat that the incident wave is the one proportional to A in previous formula. Transmission and reflection coefficients are:

T = (1/((|A|)^2));       R = (((|B|)^2)/((|A|)^2));

and must satisfy the unitarity constraint

T + R = 1.

The solution

Solutions of equation () are a kind of parabolic cylinder function, a particular case of confluent hypergeometric functions. There are several possible definitions for the independent couple of solutions, in the following we use the notations of the book [1] and take as independent solutions

W[α,z] ;  W[α, - z]

W[α, z] has no definite parity transformation property, then, as eq.() is parity invariant, if W[α,z] is a solution, W[α,-z] is another, independent solution. The explicit form of W[α,z] is not very illuminating

W[α,z]= (2^(-3/4))((Sqrt[((Subscript[G, 1])/(Subscript[G, 3]))])H[((-3)/4),(α/2),(1/4)(z^2)]+ (Sqrt[((2(Subscript[G, 3]))/(Subscript[G, 1]))]) x H[((-1)/4),(α/2),(1/4)(z^2)]);

with

(Subscript[G, 1]) = Abs[Γ((1/4)+ i (α/2))];  (Subscript[G, 3]) = Abs[Γ((3/4)+ i (α/2))]; 
H[m,n,x] = Exp[-i x]Hypergeometric1F1[m+1-i n,2 m + 2, 2 i x];

These functions are real. Here is an example of their behavior for α = 1 and positive x.

What matters in this case is that through these basis functions it is possible to define functions which have a "simple" asymptotic behavior

E[α,z] = (1/(Sqrt[k]))W[α,z] + i (Sqrt[k]) W[α,-z];
(SuperStar[E])[α,z] = (1/(Sqrt[k]))W[α,z] - i (Sqrt[k]) W[α,-z];

with

k = (Sqrt[1 + (e^(2π α))]) - (e^(π α));     1/k =  (Sqrt[1 + (e^(2π α))]) +  (e^(π α)).

For large and positive z

E[α,z] = (Sqrt[(2/z)])Exp[i (((z^2)/4) - α Log[z] + (1/2)ϕ + (π/4))](1 + 𝒪(1/(z^2)))

with

ϕ = Arg[Γ((1/2)+ i α)]

These functions are analytic in the whole complex plane then from () we have at once the solution satisfying the right boundary condition:

ψ[z] = E[α,z]

This example is particularly instructive: the function is analytic so the solution is valid also for  z < 0. In this region its asymptotic expansion must contain both a progressive and a reflected wave, in particular a part with negative exponent. This part cannot be obtained by the analytic continuation of the asymptotic form (), as its leading term is a function of (z^2), and do not change for z → - z. The conclusion, which has a general validity, is that the asymptotic expansion of analytic continuation is not given by the analytic continuation of the asymptotic expansion. This is a particular case of the so called Stokes phenomenon.

Our problem is to obtain the asymptotic form of E[α,z] for large negative z, in order to extract A and B coefficients in ().

This is easily done passing through the basis functions W[α,z]. By inverting ()

W[α,z] = ((Sqrt[k])/2)(E[α,z] + (SuperStar[E])[α,z]);
W[α,-z] = (1/(2i (Sqrt[k])))(E[α,z] - (SuperStar[E])[α,z]);

From (), using these relations:

E[α,-z]=  (1/(Sqrt[k]))W[α,-z] + i (Sqrt[k]) W[α,z]=
(i/2)E[α,z](k - (1/k)) + (i/2)(SuperStar[E])[α,z](k+(1/k)) =- i (e^(π α))E[α,z] + i (Sqrt[1 + (e^(2π α))]) (SuperStar[E])[α,z].

From () we see that we have to look for a combination

E[α,-|z|] = A (SuperStar[E])[α,|z|] + B  E[α,|z|]

then we have immediately

A=   i(Sqrt[1 + (e^(2π α))])  ;   B = -i (e^(π α)).

The transmission and reflection coefficients follow

T = (1/((|A|)^2)) = (1/(1 + (e^(2π α))));       R = (((|B|)^2)/((|A|)^2)) =   ((e^(2π α))/(1 + (e^(2π α)))) .

The unitarity constraint is obviously satisfied.

For energies much below the threshold, i.e.  α →∞

T ∼ Exp[-2 π α];   R ∼ 1 -  Exp[-2 π α];

Phases

It can be useful to compute relative phases between incoming and scattered and reflected waves. To this purpose we have first of all to define these phases. A possibility is to write the asymptotic part of the wave function as

(1/(Sqrt[p]))(Exp[i S[a,x]] + ℛ Exp[-i S[a,x]]);    x → -∞ 
𝒯(1/(Sqrt[p])) Exp[i S[b,x]];    x → +∞

S is the classical action in units of ℏ, a and b the two turning points. This formula is normalized with unit incident flux. Clearly

T = ((|𝒯|)^2);    R = ((|ℛ|)^2)

Our solution, with unit incident flux, reads

(SuperStar[E])[α,|z|] + (B/A)  E[α,|z|];    z → -∞ 
(1/A)E[α,|z|];    z → +∞

The phases are readily identified computing classical actions. In our case the two inversion points are z= ± 2 (Sqrt[α])and

S[b,z]= (∫_(2(Sqrt[α]))^z)(Sqrt[((z^2)/4)-α])𝕕z = ((z^2)/4)-(α/2)-α Log[2]+α Log[((2 (Sqrt[α]))/z)] + 𝒪(1/z)

S[a,z]= (∫_(-2(Sqrt[α]))^(-|z|))(Sqrt[((z^2)/4)-α])𝕕z =-( ((z^2)/4)-(α/2)-α Log[2]+α Log[((2 (Sqrt[α]))/(|z|))] )+ 𝒪(1/z)= - S[b,x].

By comparison with  () we have

E[α,z]→ (1/(Sqrt[p])) Exp[i S[b,z]+ (i/2) (ϕ - α Log[(α/e)]+(π/2))]= (1/(Sqrt[p])) Exp[i S[b,z]+ (i/2) φ]
(SuperStar[E])[α,z]→ (1/(Sqrt[p])) Exp[-i S[b,z]- (i/2) (ϕ - α Log[(α/e)]+(π/2))]= (1/(Sqrt[p])) Exp[i S[a,z]- (i/2) φ]

By eliminating the extra phase φ and moving the origin for S[b,x] for reflected wave the solution () takes the form

(1/(Sqrt[p])) Exp[i S[a,z] + (B/A)(e^(i φ )) (1/(Sqrt[p])) Exp[-i S[a,z] ];    z → -∞

((e^(i φ ))/A)E[α,|z|];    z → +∞

We have then

𝒯 = (1/(|A|))Exp[i(ϕ - α Log[(α/e)]+(π/2))-i (π/2))] = (1/(|A|))Exp[-i δ];

ℛ =- ((|B|)/(|A|))Exp[i(ϕ - α Log[(α/e)]+(π/2))]= -i ((|B|)/(|A|))Exp[-i δ];

with

δ = α Log[(α/e)] + Arg[Γ((1/2)- i α)]

To this equation we can give an invariant form, as π α is just the integral of p/ℏ below the barrier, let us call σ this integral

δ = (σ/π)Log[(σ/(π e))] + Arg[Γ((1/2)- i (σ/(π )))]

This solves completely the problem of phase shifts.

Poles

Let us note that T has poles for

α = (n+(1/2))i  ⇒  E = - i ℏ (Sqrt[(β/m)])(n+(1/2))

The analytic continuation β → Exp[i π] β transforms the barrier in a well and the poles become go to the bound states of this harmonic well.

Lowest order WKB

As stated above it is impossible to derive A and B coefficients by analytic continuation but it is possible to apply the usual formula for transmission coefficient

T = Exp[((- 2)/ℏ)∫p[x]𝕕x ]

where the integral is performed between turning points. In z variables the turning points are z = ± 2 (Sqrt[α]) and

T ∼ Exp[-2 (∫_(-2(Sqrt[α]))^(2(Sqrt[α])))(Sqrt[α - ((z^2)/4)])𝕕z ] = Exp[-4α (∫_(-1)^(+1))(Sqrt[1-(y^2)])𝕕y] = Exp[-2π α]

Reproducing the previous asymptotic result.

Reflection above the barrier

The only change in formulas above is α → -α. The transmission and reflection coefficient are then given by

T = (1/((|A|)^2)) = (1/(1 + (e^(-2π α)))) ∼ 1 - (e^(-2π α));       R = (((|B|)^2)/((|A|)^2)) =   ((e^(-2π α))/(1 + (e^(-2π α)))) ∼ (e^(-2π α)).

In usual units

R ∼Exp[ -((2π|E|)/ℏ)(Sqrt[(m/β)])]

which goes rapidly to zero in the classical limit,  ℏ → 0.

References

1) M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover Publications Inc. (1965).

Problem

Show that the normalization condition C = 2(Sqrt[m/T])for a semiclassical bound state in one dimension  is valid also for low-lying quantum states if the quantization condition holds.

Solution

The following proof is due W.H. Furry [1].

Let us choose the origin of the coordinates x=0 as the minimum of the potential V[x]. Let (a,b) the allowed classical region, a, b are the classical inversion points, ad ψ the normalized semiclassical solution.

Let us consider first the region x>0 and u[E,x] a normalizable solution with an arbitrary energy E (in general this solution is not bounded as x→ - ∞)

-(( (ℏ^2))/(2m))((d^2)/(dx^2))u[x] + V[x]u[x] =  E u[x].

The solution is chosen in such a way that for E → (Subscript[E, 0]), u[E,x] → ψ[x], where (Subscript[E, 0])is the bound state energy. Writing the eq.() for two energies (Subscript[E, 1])and (Subscript[E, 2]) and subtracting them after multiplication by (Subscript[u, 2]) and (Subscript[u, 1]) (the two solutions)

-(( (ℏ^2))/(2m))((Subscript[u, 1])(Subscript[u, 2]^(′′)) - (Subscript[u, 2])(Subscript[u, 1]^(′′))) = -(( (ℏ^2))/(2m))(d/dx)((Subscript[u, 1])(Subscript[u, 2]^′) - (Subscript[u, 2])(Subscript[u, 1]^′))= ((Subscript[E, 2])-(Subscript[E, 1]))(Subscript[u, 2])(Subscript[u, 1]).

Integration from 0 to ∞ gives

(( (ℏ^2))/(2m))(Subscript[((Subscript[u, 1])(Subscript[u, 2]^′) - (Subscript[u, 2])(Subscript[u, 1]^′)), x=0]) = ((Subscript[E, 2])-(Subscript[E, 1]))(∫_0^∞)(Subscript[u, 2])(Subscript[u, 1])dx.

With (Subscript[E, 2])= (Subscript[E, 1])+δE the comparison of first order in δE gives (we write (Subscript[u, 1])=u):

(( (ℏ^2))/(2m))(Subscript[(u((∂u')/(∂E)) - ((∂u)/(∂E))u'), x=0]) = (∫_0^∞)(u^2)dx.

Analogously, with v[E, x] the solutions in (-∞, 0):

-(( (ℏ^2))/(2m))(Subscript[(v((∂v')/(∂E)) - ((∂v)/(∂E))v'), x=0]) = (∫_(-∞)^0)(v^2)dx.

If we sum these equations and take the limit E→(Subscript[E, 0])the right hand side is just the normalization of ψ, then it is equal to 1.

Let us compute the left hand side of eqs( , ).

Let us remember that x=0 is the minimum of V[x], then V'[0]=0and for the classical momentum

p[x]=(Sqrt[2m(E-V[x])])

we have p'[0]=0, then all derivatives of the momentum can be neglected in x=0. With p[0]=(Subscript[p, 0]) we have easily the semiclassical approximation for u and v:

v[0]=(C/(Sqrt[(Subscript[p, 0])]))Cos[(1/ℏ)(∫_a^0)pdx - (π/4)];  v'[0]=-(C/ℏ)(Sqrt[(Subscript[p, 0])])Sin[(1/ℏ)(∫_a^0)pdx - (π/4)];
u[0]=(C/(Sqrt[(Subscript[p, 0])]))Cos[(1/ℏ)(∫_0^b)pdx - (π/4)];  u'[0]=(C/ℏ)(Sqrt[(Subscript[p, 0])])Sin[(1/ℏ)(∫_0^b)pdx - (π/4)];

The sign in u'[0] is due to the fact that for u[x] the integral of the phase factor goes from x to b.

In taking the derivative with respect to E the reader can check the ∂C/∂E cancel in eqs( , ) and we get

(( (ℏ^2))/(2m))(Subscript[(v((∂v')/(∂E)) - ((∂v)/(∂E))v'), x=0])= -((C^2)/(2ℏ (Subscript[p, 0])))((∂(Subscript[p, 0]))/(∂E))Sin[(2/ℏ)(∫_a^0)pdx] - ((C^2)/(ℏ^2))(∫_a^0)((∂p)/(∂E))dx .
(( (ℏ^2))/(2m))(Subscript[(u((∂u')/(∂E)) - ((∂u)/(∂E))u'), x=0])= ((C^2)/(2ℏ (Subscript[p, 0])))((∂(Subscript[p, 0]))/(∂E))Sin[(2/ℏ)(∫_0^b)pdx]+((C^2)/(ℏ^2))(∫_0^b)((∂p)/(∂E))dx .

Using these equalities the sum of the left hand side of eqs( , ) in the limit E→(Subscript[E, 0])gives

2(m/(ℏ^2))=((C^2)/(2ℏ (Subscript[p, 0])))((∂(Subscript[p, 0]))/(∂E))Sin[(2/ℏ)(∫_a^b)pdx] +((C^2)/(ℏ^2))(∫_a^b)((∂p)/(∂E))dx .

The first integral in () vanish if the quantization condition holds

ϕ=(2/ℏ)(∫_a^b)pdx = 2((n+(1/2))π) ⇒ Sin[ϕ] = 0.

Using ∂p/∂E = m/p = 1/v we get

2m = (C^2)(∫_a^b)(1/v)dx = (C^2) (T/2) ⇒ C = 2 (Sqrt[(m/T)]) ,

where T is the classical period of motion.

References

1) W.H. Furry: Phys. Rev. 71, 360, (1947).

Problem

a) Study the transformation properties of

(Subscript[I, 1])=∮ (Subscript[p, i])d (Subscript[q, i])     and    (Subscript[I, 2])=(∑_i) ∮ (Subscript[p, i])d (Subscript[q, i])

under coordinate transformations and general canonical transformations.

b) Show that in a system with n degrees of freedom there exist at most n functionally independent integrals of motion.

c) Show that in the hypothesis of Liouville theorem for integrable systems there are n tangent vectors tangent to the torus and construct them. Show that the involution property of the integrals of motion implies the local existence of S such that ∇S = p.

d) Consider a n-1 parameter family of trajectories in an integrable system. Show that envelopes of the family concides with caustics.

Solution

a)

A canonical transformation (q,p) → (Q,P) can be defined by a generating function F[q,P] with

(Subscript[p, i])= ((∂F)/(∂(Subscript[q, i])));            (Subscript[Q, i])= ((∂F)/(∂(Subscript[P, i])));

see ref.[1]. Then we have

(∑_i)(Subscript[p, i]) (Subscript[dq, i]) = (∑_i) ((∂F)/(∂(Subscript[q, i])))(Subscript[dq, i]) = (∑_i) (((∂F)/(∂(Subscript[q, i])))(Subscript[dq, i]) + ((∂F)/(∂(Subscript[P, i])))(Subscript[dP, i])) - (∑_i) ((∂F)/(∂(Subscript[P, i])))(Subscript[dP, i]) ≡ dF - (∑_i) (Subscript[Q, i])(Subscript[dP, i])=
dF - d( (∑_i) (Subscript[Q, i])(Subscript[P, i]))+ (Subscript[P, i])(Subscript[dQ, i]).

As a closed integral of a total differential vanishes we have

(∑_i) ∮ (Subscript[p, i])d (Subscript[q, i])  =  (∑_i) ∮ (Subscript[P, i])d (Subscript[Q, i]) ;

i.e. (Subscript[I, 2])is invariant under general canonical transformations. The proof clearly does not work for (Subscript[I, 1])as the missing sum prevent the construction of the differential dF.

A point transformation, i.e. a coordinate's transformation,  is defined by

F[q,P] = (∑_j)(Subscript[f, j])[q](Subscript[P, j])

In this case, from ()

(Subscript[Q, i])= (Subscript[f, i])[q];             (Subscript[dQ, i])=(∑_j) ((∂(Subscript[f, i]))/(∂(Subscript[q, j])))(Subscript[dq, j]) ≡(∑_j) (Subscript[J, ij]) (Subscript[dq, j]);         (Subscript[p, i])=(∑_j) ((∂(Subscript[f, j]))/(∂(Subscript[q, i])))(Subscript[P, j]) ≡ (Subscript[((J^T)), ij])(Subscript[P, j])

J is the jacobian matrix for coordinate transformations. As for any Jacobian matrix (J^T)= (J^(-1)).  dQ and P transform respectively as a controvariant vector (a usual vector) and a covariant vector (a gradient). For the differential form (Subscript[I, 2])

(∑_i) (Subscript[p, i]) (Subscript[dq, i]) = (∑_i)(∑_(j,s))(Subscript[((J^T)), ij])(Subscript[P, j]) (Subscript[((J^(-1))), is]) (Subscript[dQ, s]) =  (∑_i)(∑_(j,s)) (Subscript[P, j])( (Subscript[J, ji]) (Subscript[J, is]^(-1))) (Subscript[dQ, s]) = (∑_j)  (Subscript[P, j])d (Subscript[Q, j]).

Then (Subscript[I, 2]) is invariant, as expected, while (Subscript[I, 1])is not invariant as the missing sum on the index i in the previous formula prevent the cancellation of Jacobian matrices.

b)

A system of k functions (Subscript[F, s])in (ℝ^N) is functionally independent if the k differentials (Subscript[dF, s])are linearly independent, i.e. if the Jacobian matrix has rank k. This definition comes form Dini theorem on implicit functions. The equation (Subscript[F, i])= 0 defines locally a surface in (ℝ^N). If the rank of the Jacobian matrix is k we can solve the system (Subscript[F, i])= 0 for k unknown quantities in terms of the other N-k unknown quantities. In our case N = 2n is the dimension of the phase space of the system and k = n for integrable systems.

Let us now show that at most n such F can exist. We can perform a canonical transoformation which brings (Subscript[F, 1])to a canonical momentum, let (Subscript[φ, 1])the conjugate angular variable. As (Subscript[F, 1])is an integral of motion (Subscript[dF, 1])/dt = 0, then by Hamilton equations H cannot depends on (Subscript[φ, 1]), it is a cyclic variable. Then our Hamiltonian has the form

H = H((Subscript[F, 1]),(Subscript[p, 2]),...(Subscript[p, n]),(Subscript[q, 2,...]),(Subscript[q, n]))

The existence of one integral has disposed of two variables. We can repeat this procedure at most n times, eliminating the n-couples (p,q).

c)

Let (Subscript[F, i]), i = 1...n the n integral of motion for a system with n degrees of freedom. The 2n coordinates of phase space will be ordered as (p,q). The symbol ∇ whitout suffixes will denote the gradient in this space. The n equations

(Subscript[F, i])(p,q) = (Subscript[I, i]);  i = 1… n;

define a 2 n - n = n dimensional manifold in phase space.

Consider the n vectors in (ℝ^(2n))

(Subscript[V, i]) = (- (Subscript[∇, q]) (Subscript[F, i]), (Subscript[∇, p])(Subscript[F, i]))

Each vector (Subscript[V, i])is tangent to the surface (Subscript[F, i])= (Subscript[I, i]) as it is orthogonal to the gradient vector ∇(Subscript[F, i])≡ ((Subscript[∇, p]) (Subscript[F, i]), (Subscript[∇, q])(Subscript[F, i])). As the quantities (Subscript[F, i])are in involution (Subscript[V, i])is tangent also to the other surfaces (Subscript[F, j])= (Subscript[I, j]), in fact  it is orthogonal to their gradient:

(Subscript[V, i]) · ∇ (Subscript[F, j]) = (- (Subscript[∇, q]) (Subscript[F, i]), (Subscript[∇, p])(Subscript[F, i]))·((Subscript[∇, p]) (Subscript[F, j]), (Subscript[∇, q])(Subscript[F, j]))= (Subscript[∇, p])(Subscript[F, i])(Subscript[∇, q])(Subscript[F, j]) - (Subscript[∇, q]) (Subscript[F, i])(Subscript[∇, p]) (Subscript[F, j]) = {(Subscript[F, i]),(Subscript[F, j])} = 0.

Then there exist a set of n regular vectors tangent to the surface (). Due to a theorem of Poincaré - Hoopf the only surface which admits n regular tangent vectors has the topology of a torus, this is the topological content of Liouville-Arnold theorem.

Now let us consider the implication of involution property on the integrability of the differential form p dq, i.e. on the existence of a function S such that ∇S = p.  As is well known a necessary property for local integrability is that locally, i.e. in sufficiently small open sets:

∮ p dq = 0

Now let us consider the surface (). The functional independence of the constraints allow to invert locally these equations and to express p coordinates in terms of q and I:

(Subscript[p, i]) = (Subscript[f, i])[q,I].

The condition for local validity of () is

((∂(Subscript[f, i]))/(∂(Subscript[q, j])))= ((∂(Subscript[f, j]))/(∂(Subscript[q, i]))).

Taking the derivative of constraints (), written as F[f(q,I),q], we have

((∂(Subscript[F, i]))/(∂(Subscript[p, α])))((∂(Subscript[f, α]))/(∂(Subscript[q, j])))+ ((∂(Subscript[F, i]))/(∂(Subscript[q, j])))= 0

This equation has the matrix form A B = - C, where

(Subscript[A, ij])= ((∂(Subscript[F, i]))/(∂(Subscript[p, j])));   (Subscript[B, ij])= ((∂(Subscript[f, i]))/(∂(Subscript[q, j])));  (Subscript[C, ij])= ((∂(Subscript[F, i]))/(∂(Subscript[q, j]))).

The condition () means B = (B^T) i.e.

(A^(-1))C = ((((A^(-1))C))^T)= (C^T) (A^(-1T)) ⇒ C (A^T) - A (C^T) = 0.

Using the explicit expression for the matrices

((∂(Subscript[F, i]))/(∂(Subscript[q, α])))((∂(Subscript[F, j]))/(∂(Subscript[p, α]))) - ((∂(Subscript[F, i]))/(∂(Subscript[p, α])))((∂(Subscript[F, j]))/(∂(Subscript[q, α]))) = 0 ⇒ {(Subscript[F, i]),(Subscript[F, j])} = 0.

The involution of the constraints is the integrability condition.

d)

Let us consider a family of trajectories for a system with n degrees of freedom:

(Subscript[x, i])= (Subscript[x, i])(t,(Subscript[c, 2]),(Subscript[c, 3]),... ,(Subscript[c, n])).

The envelope of the family is given by definition by the equation

(Subscript[J, c])= det(({{((∂(x^1))/(∂t)), ((∂(x^2))/(∂t)), ..., ((∂(x^n))/(∂t))}, {((∂(x^1))/(∂(Subscript[c, 2]))), ((∂(x^2))/(∂(Subscript[c, 2]))), ..., ((∂(x^n))/(∂(Subscript[c, 2])))}, {..., ..., ..., ...}, {((∂(x^1))/(∂(Subscript[c, n]))), ((∂(x^2))/(∂(Subscript[c, n]))), ..., ((∂(x^n))/(∂(Subscript[c, n])))}}))= 0.

Consider now an integrable system. A generic solution of equation of motion can be expressed throw angle variables as

(Subscript[φ, i]) = (Subscript[ω, i])t + (Subscript[δ, i]);   (Subscript[ω, i])= ((∂H[I])/(∂(Subscript[I, i]))).

One of the constants (Subscript[δ, i])can be reabsorbed by shifting time origin. For simplicity let us identify the other constants with our parameters (Subscript[c, i]).

The coordinates φ are global variables on the invariant torus, a caustic, corresponding to a singularity of the projection on configuration space, is given by

(Subscript[J, φ])=det(((∂ (Subscript[x, i]))/(∂ (Subscript[φ, j])))) = 0.

Within the above conventions

((∂(Subscript[x, i]))/(∂(Subscript[c, j])))= ((∂(Subscript[x, i]))/(∂(Subscript[φ, j])));  j ≥ 2;   ((∂(Subscript[x, i]))/(∂t))= (Subscript[ω, 1])((∂(Subscript[x, i]))/(∂(Subscript[φ, 1])))+ (∑_(k>1))(Subscript[ω, k])((∂(Subscript[x, i]))/(∂(Subscript[φ, k]))).

Then the rows 2 ... n of the two determinants (Subscript[J, φ])and (Subscript[J, c]) are identical. The first row of (Subscript[J, c]) is the sum of a multiple of the first row of (Subscript[J, φ]) and a linear combinations of higher rows, which do not contribute to the determinant, then

(Subscript[J, c]) = (Subscript[ω, 1]) (Subscript[J, φ]);

this gives an identification between caustics and envelopes.

References

1) H. Goldstein: Classical Mechanics, Addison-Wesley Publishing Company (1965)

Problem

Perform the ℏ → 0 limit in the Schrödinger equation with the substitution

ψ = A Exp[ (i/ℏ) S]

Solution

a)

For an Hamiltonian of the form

H = -((ℏ^2)/(2m))Δ + V[x]

a direct substitution of

ψ = A Exp[ (i/ℏ) S]

in the Schrödinger equation gives

((∂S)/(∂t)) + H[∇S,x] = 0;         ((∂(A^2))/(∂t))+∇((A^2)v) = 0;

with

v ≡ ((∂H)/(∂p));    p = ∇S .

For a general Hamiltonian substitution of () in the Schrödinger equation gives

i ℏ ((∂A)/(∂t)) - A ((∂S)/(∂t)) = (e^(-(i/ℏ)S))H[(ℏ/i)∇,x]A (e^((i/ℏ)S)).

Using commutation relation

(e^(-(i/ℏ)S))(ℏ/i)∇ (e^((i/ℏ)S)) = (ℏ/i)∇ + ∇S

and by a formal expansion in ℏ we have:

(e^(-(i/ℏ)S))H[(ℏ/i)∇,x] (e^((i/ℏ)S)) = H[(ℏ/i)∇ + ∇S,x] ≃  H[∇S,x]+(1/2)(((∂H)/(∂p))(ℏ/i)∇ + (ℏ/i)∇ ((∂H)/(∂p))) + 𝒪((ℏ^2))

Inserting this expression in () and using the definition () for v

i ℏ ((∂A)/(∂t)) - A ((∂S)/(∂t)) =  H[∇S,x] A + (1/2)((Subscript[v, s]) (ℏ/i)((∂A)/(∂(x^s)))+  (ℏ/i)(∂/(∂(x^s)))((Subscript[v, s]) A)) +  𝒪((ℏ^2))

Selecting order 0 and order 1 in ℏ we obtain again equation ().