Bound States in a Pöschl-Teller Potential

Erek Bilgici

In this notebook we will solve the Schrödinger equation for the Pöschl-Teller potential in one dimension analytically for λ>1.

\(TraditionalForm\`V(
        x) = \(-\(\(ℏ\^2\ α\^2\ λ\ \((λ - 
                  1)\)\)\/\(\((2\ m)\)\ \(\(cosh\^2\)(α\ x)\)\)\)\)\)

We observe the behavior of the eigenfunctions by changing the two parameters λ , α and energy. Further we will discuss a superposition of m eigenfunctions.

As usual we set ℏ=m=1.
Most of the necessary functions are defined in the section "Definitions", so if you want to run the input cells you should first evaluate that section.

Definitions

needed packages

<< VQM`;
<< Graphics`FilledPlot`;

switch off some disturbing messages

Off[NIntegrate::"ncvb", NIntegrate::"slwcon", 
        General::"spell1", NIntegrate::"ploss"];

potential and eigenvalues

V[
          x_, λ_, α_] := -(((α^2) λ (λ - 1))/(2 Cosh[α x]^2));

Ept[
          n_, λ_, α_] := -(1/2) (((λ - 1 - n))^2) (α^2);

k as a funcions of E

k[e_] := Sqrt[2 e]

eigenfunctions

NormFacE[e_, λ_, α_] := 
        NIntegrate[
          Abs[Cosh[α x]^λ Hypergeometric2F1[
                (1/2) (λ + ((I k[e])/α)), 
                (1/2) (λ - ((I k[e])/α)), 
                (1/2), -Sinh[α x]^2]], {x, -20, 
            20}];

NormFacO[e_, λ_, α_] := 
        NIntegrate[
          Abs[Cosh[α x]^λ Sinh[α x] Hypergeometric2F1[
                (1/2) (λ + ((I k[e])/α)) + 
                  (1/2), (1/2) (λ - ((I k[e])/α)) + (1/2), (3/2), -Sinh[α x]^2]], {x, -20, 
            20}];

φe[x_, e_, λ_, α_] := 
        Cosh[α x]^λ Hypergeometric2F1[
              (1/2) (λ + ((I k[e])/α)), 
              (1/2) (λ - ((I k[e])/α)), 
              (1/2), -Sinh[α x]^2]/
            NormFacE[e, λ, α];

φo[x_, e_, λ_, α_] := 
        Cosh[α x]^λ Sinh[α x] Hypergeometric2F1[
              (1/2) (λ + ((I k[e])/α)) + 
                (1/2), (1/2) (λ - ((I k[e])/α)) + (1/2), (3/2), -Sinh[α x]^2]/
            NormFacO[e, λ, α];

φ[n_, x_, e_, λ_, α_] := 
        If[EvenQ[n], φe[x, e, λ, α], φo[x, 
            e, λ, α]];

Solution

As you can see the Pöschl Teller potential has similarity to the harmonic oscillator at low energies so we expect comparable eigenfunctions.

λ = 4.3;
α = 1;
Plot[V[x, λ, α], {x, -4, 4}];

The Schrödinger equation

-(((ℏ^2) )/(2 m))u''(x)+V(x, λ, α) u(x)==E u(x)

becomes

-u''(x)-\(\(α\^2\ λ\ \((λ - 1)\)\)\/\(\(cosh\^2\)(α\ x)\)\ \(u(x)\)\)==2 E u(x)

\(To\ solve\ this\ differential\ equation\ we\ substitude\ y = \(\(Cosh[α\ x]\^2\)\(.\)\)\)

y (1 - y)  u''+((1/2) - y) u'-u ((((λ - 1) λ)/(4 y)) + 
                (E/(2 (α^2))))==0

Splitting off a fitting power of y by setting

u = (y^(λ/2)) v(y)

we arrive at the hypergeometric differential equation:

y (1 - y) v''+((λ + (1/2)) - (λ + 1) y) v'-(1/4) v ((λ^2) + ((2 E)/(α^2)))==0

With the abbreviations

\(TraditionalForm\`a = \(1\/2\ \((λ + \(i\ \@\(2\ E\)\)\/α)\)\ \ \ \ \ and\ \ \ \ \ b = 
        1\/2\ \((λ - \(i\ \@\(2\ E\)\)\/α)\)\)\)

we get the complete solution of the Schrödinger equation:

u(x)=(y^(λ/2)) (B (Sqrt[1 - y]) (Subscript[, 2])(Subscript[F, 1])(a + (1/2),b + (1/2);(3/2);1 - y)+A (Subscript[, 2])(Subscript[F, 1])(a,b;(1/2);1 - y))

\(We\ shall\ pick\ out\ a\ fundametal\ system\ of\ two\ real\ standard\ solutions, \ u\_E\ and\ u\_O, \ even\ and\ odd\ in\ x . \ With\ B = 0, \ 
    A = 1\ we\ find\ the\ even\ standard\ solution, \)

(Subscript[u, E])(x)=(cosh^λ)(α x) (Subscript[, 2])(Subscript[F, 1])(a,b;(1/2);-(sinh^2)(α x))

and with A=0, B=i the odd standard solution,

(Subscript[u, O])(x)=(cosh^λ)(α x) sinh(α x) (Subscript[, 2])(Subscript[F, 1])(a + (1/2),b + (1/2);(3/2);-(sinh^2)(α x))

Bound states

\(Eigenvalues\ exist\ for\ negative\ energies\ so\ k \((E)\) = \(\@\(2\ E\)\) becomes\ \ imaginary\ and\ the\ abbreviations\ a\ and\ b\ become\ real . 
\ To\ find\ out\ the\ conditions\ \ of\ λ\ and\ α\ to\ get\ eigenvalues\ we\ first\ need\ the\ asymptotic\ behaviour\ of\ the\ even\ and\ odd\ solution\ at\ negative\ infinity\ large\ values\ of\ the\ \(\(argument\)\(.\)\)\)

\(TraditionalForm\`\(-\(sinh(α\ x)\)\) -> \(-\(1\/4\)\)\ e\^\(2\ α\ |x|\)\)

The two solutions yield

\(TraditionalForm\`\(u\_E\)(x) -> 
      2\^\(-λ\)\ e\^\(λ\ α\ |x|\)\ \(Γ(
          1\/2)\)\ \((\(2\^\(2\ a\)\ e\^\(\(-2\)\ a\ α\ |x|\)\ \(Γ(b - a)\)\)\/\(\(Γ(b)\)\ \(Γ(1\/2 - a)\)\) + \(2\^\(2\ b\)\ e\^\(\(-2\)\ b\ α\ |x|\)\ \(Γ(a - b)\)\)\/\(\(Γ(a)\)\ \(Γ(1\/2 - b)\)\))\)\)

and

\(TraditionalForm\`\(u\_O\)(x) -> 
      2\^\(-\((λ + 1)\)\)\ e\^\(\((λ + 1)\)\ α\ |x|\)\ \(Γ(
          3\/2)\)\ \((\(2\^\(2\ a + 1\)\ e\^\(\(-\((2\ a + 1)\)\)\ α\ |x|\)\ \(Γ(b - a)\)\)\/\(\(Γ(b + 1\/2)\)\ \(Γ(1 - a)\)\) + \(2\^\(2\ b + 1\)\ e\^\(\(-\((2\ b + 1)\)\)\ α\ |x|\)\ \(Γ(a - b)\)\)\/\(\(Γ(a + 1\/2)\)\ \(Γ(1 - b)\)\))\)\)

We simplify both solutions:

\(TraditionalForm\`\(u\_e\)(
        x) -> \(2\^\(2\ a - λ\)\ e\^\(\(-i\)\ \@\(2\ E\)\ |x|\)\ \@π\ \(Γ(b - a)\)\)\/\(\(Γ(b)\)\ \(Γ(1\/2 - a)\)\) + \(2\^\(2\ a - λ\)\ e\^\(i\ \@\(2\ E\)\ |x|\)\ \@π\ \(Γ(a - b)\)\)\/\(\(Γ(a)\)\ \(Γ(1\/2 - b)\)\)\)

and

\(TraditionalForm\`\(u\_O\)(
        x) -> \(2\^\(2\ a - \((λ + 1)\)\)\ e\^\(\(-i\)\ \@\(2\ E\)\ |x|\)\ \@π\ \(Γ(b - a)\)\)\/\(\(Γ(b + 1\/2)\)\ \(Γ(1 - a)\)\) + \(2\^\(2\ b - λ - 1\)\ e\^\(i\ \@\(2\ E\)\ |x|\)\ \@π\ \((\(Γ(a - b)\)\ 2\^\(2\ b + 1\)\ e\^\(\(-\((2\ b + 1)\)\)\ α\ |x|\))\)\)\/\(\(Γ(a + 1\/2)\)\ \(Γ(1 - b)\)\)\)

The exponents of all exponential terms get real while the first term behaves like

\(TraditionalForm\`\(\(e\^\(\(-i\)\ \@\(2\ E\)\ |x|\) == e\^\(\(-i\)\ \@\(\(-2\)\ |E|\)\ |x|\) == e\^\(\@\(2\ |E|\)\ |x|\)\)\(,\)\)\)

and the second

\(TraditionalForm\`e\^\(i\ \@\(2\ E\)\ |x|\) == e\^\(i\ \@\(\(-2\)\ |E|\)\ |x|\) == \(\(e\^\(\(-\@\(2\ |E|\)\)\ |x|\)\)\(.\)\)\)

A normalizable solution is therfore possible for |x| > 0 if, and only if, the factor of the first term vanishes. Since the Γ functions are now all taken for real arguments where poles exist at negative integers "-m"  (m=0,1,2,...) .

Plot[Gamma[m], {m, -5, 3}];

The eigenvalues follow from

\(TraditionalForm\`1\/2 - 
        b = \(\(1 - λ\)\/2 + \@\(2\ |E|\)\/\(2\ α\) = \(\(-m\)\ \ \ \  -> \ \ \ \ E\_m = \(-\(1\/2\)\)\ α\^2\ \((\(-2\)\ m + λ - 1)\)\^2\)\)\)

for even, and

\(TraditionalForm\`1 - 
        a = \(\(-\(λ\/2\)\) + \@\(2\ |E|\)\/\(2\ α\) + 
          1 = \(\(-m\)\ \ \ \  -> \ \ \ \ E\_m = \(-\(1\/2\)\)\ α\^2\ \((\(-\((2\ m + 1)\)\) + λ - 1)\)\^2\)\)\)

for odd eigenstates.
We combine the even and odd eigenstates. Since 2m are even integer numbers and 2m+1 are odd ones and modify:

\(TraditionalForm\`E\_n = \(-\(1\/2\)\)\ α\^2\ \((\(-n\) + λ - 1)\)\^2\ \ \ \ \ with\ \ \ \ \ n ≤ λ - 
          1\)

Eigenfunctions

Range

xLeftRight = 3.5;

Define plots

p1[n_, λ_, α_] := 
        Plot[Evaluate[φ[n, x, 
                Ept[n, λ, α], λ, α] + 
              Ept[n, λ, α]], {x, -xLeftRight, xLeftRight}, 
          PlotRange -> All, DisplayFunction -> Identity, 
          AxesOrigin -> {0, 0}];

p2[λ_, α_] := 
        FilledPlot[{V[x, λ, α], 
            V[0, λ, α]*1.02}, {x, -xLeftRight, 
            xLeftRight}, DisplayFunction -> Identity, 
          Fills -> GrayLevel[0.8]];

p3[n_, λ_, α_] := 
        Plot[Evaluate[φ[n, x, 
                Ept[n, λ, α], λ, α] + 
              Ept[n, λ, α]], {x, -xLeftRight, xLeftRight}, 
          AxesOrigin -> {0, Ept[n, λ, α]}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];

Illustration

plotEF[n_, λ_, α_] := 
        If[λ - 1 ≥ n, 
          Show[GraphicsArray[{Show[{p2[λ, α], 
                    p1[n, λ, α]}], 
                p3[n, λ, α]}, 
              DisplayFunction -> $DisplayFunction, 
              ImageSize -> {550, 250}, 
              PlotLabel -> 
                StringForm["n =`1`\nλ =`2`\nα =`3`\nE \
=`4`", PaddedForm[n, 2], PaddedForm[λ, {4, 2}], 
                  PaddedForm[α, {4, 2}], 
                  PaddedForm[N[Ept[n, λ, α]], {4, 1}]]]], 
          Print["There are only ", λ - 
              1, " (λ-1) eigenstates"]];

plotEF[3,6,1]

Parameter λ

Run the whole section!

\(TraditionalForm\`V(
        x) = \(-\(\(ℏ\^2\ α\^2\ λ\ \((λ - 
                  1)\)\)\/\(\((2\ m)\)\ \(\(cosh\^2\)(α\ x)\)\)\)\)\)

We discuss the behaviour of the λ parameter. As you can imagine this parameter adjusts the depth of the potential. In each plot all eigenfunctions are shown while λ gets bigger. The number of eigenstates goes like n≤ λ-1. At each integer step of λ there exists one more eigenstate. It is remarkable that if λ is an integer value even at energy (Subscript[E, λ - 1]) = 0 exists one eigenfunction. The special eigenfunctions at zero energy is drawn thicker.

Set Parameters

Clear[n, λ, α];
α = 1;

Range

xLeftRight = 3.5;
λStart = 2;
λEnd = 10;
λStep = 0.1;

Define Plot

doplotλ[λ_, α_] := 
        Show[{FilledPlot[{V[x, λ, α], 
                V[0, λ, α]*1.02}, {x, -xLeftRight, 
                xLeftRight}, DisplayFunction -> Identity, 
              AspectRatio -> 1.7, ImageSize -> {250, 400}, 
              Fills -> GrayLevel[0.8],  
              PlotRange -> {V[
                      0, λ, α]*1.03, -V[
                        0, λ, α]*0.1}], 
            Table[Plot[
                Evaluate[φ[n, x, 
                      Ept[n, λ, α], λ, α] + 
                    Ept[n, λ, α]], {x, -xLeftRight, 
                  xLeftRight}, DisplayFunction -> Identity, 
                AspectRatio -> 1.7, ImageSize -> {250, 400}, 
                PlotStyle -> 
                  If[Ept[n, λ, α] == 
                      0, {Hue[n/(λ - 1.2)], Thickness[0.02]}, 
                    Hue[n/(λ - 1.2)]]], {n, 0, 
                IntegerPart[λ - 1]}]}, 
          PlotLabel -> 
            StringForm["λ =`1`\nα =`2`\neigenstates:`3`", 
              PaddedForm[N[λ], {4, 2}], PaddedForm[α, {4, 2}], 
              PaddedForm[IntegerPart[λ], {2, 0}]], 
          DisplayFunction -> $DisplayFunction];

Animation

Do[
        doplotλ[λ, α], {λ, λStart, λEnd, λStep}];

Parameter α

Run the whole section!

\(TraditionalForm\`V(
        x) = \(-\(\(ℏ\^2\ α\^2\ λ\ \((λ - 
                  1)\)\)\/\(\((2\ m)\)\ \(\(cosh\^2\)(α\ x)\)\)\)\)\)

By modifying the parameter α you change depth and width of the potential. Numer of eigenfunctions keeps constant while the oscillating range gets smaller with larger α.

Set Parameters

Clear[n, λ, α];
λ = 5;

Range

xLeftRight = 3.5;
αStart = 0.4;
αEnd = 5;
αStep =  .1;
Off[NIntegrate::"ploss"]

Define Plot

doplotα[λ_, α_] := 
        Show[{FilledPlot[{V[x, λ, α], 
                V[0, λ, α]*1.02}, {x, -xLeftRight, 
                xLeftRight}, DisplayFunction -> Identity, 
              AspectRatio -> 1.7, ImageSize -> {250, 400}, 
              Fills -> GrayLevel[0.8], 
              PlotRange -> {V[
                      0, λ, α]*1.03, -V[
                        0, λ, α]*0.1}], 
            Table[Plot[
                Evaluate[φ[n, x, 
                      Ept[n, λ, α], λ, α] + 
                    Ept[n, λ, α]], {x, -xLeftRight, 
                  xLeftRight}, DisplayFunction -> Identity, 
                AspectRatio -> 1.7, 
                PlotLabel -> 
                  StringForm["α = `1`\nλ =`2`\n\
eigenstates:`3`", PaddedForm[N[α], {4, 2}], 
                    PaddedForm[N[λ], {4, 2}], 
                    PaddedForm[α, {4, 2}], 
                    PaddedForm[IntegerPart[λ], {2, 0}]], 
                ImageSize -> {250, 400}, 
                PlotStyle -> Hue[n/(λ - 1.2)]], {n, 0, 
                IntegerPart[λ - 1]}]}, 
          PlotLabel -> 
            StringForm["α =`1`\nλ =`2`\neigenstates:`3`", 
              PaddedForm[N[α], {4, 2}], PaddedForm[λ, {4, 2}], 
              PaddedForm[IntegerPart[λ], {2, 0}]], 
          DisplayFunction -> $DisplayFunction];

Animation

Do[
        doplotα[λ, α], {α, αStart, αEnd, αStep}];

Energy

Eigenfunctions have to be square integrable. This is fulfilled if the functions quickly go to zero at |x| > 0. The asymptotic behaviour occurs exactly at eigenvalues (Subscript[E, n]). If the energy differs all eigenfunctions diverge. The following animation will show this behaviour. If you change any input parameters like α, λ or n you also should change the plotrange "xLeftRight". In small steps the energy will be reduced form just a little above the n'th eignevalue to a small value below. If the energy has reached exactly an eigenvalue the function will be drawn thicker and red colored.

Set Parameters

Clear[α, λ, n];
n = 6;
λ = 8.1112;
α = 2;

Range

xLeftRight = 6;
EStart = 
        Ept[n, λ, α] +  .001;
EEnd = 
        Ept[n, λ, α] -  .001;
EStep = - .00005;

Define Plot

doplotE[e_] := 
        Plot[Evaluate[φ[n, x, 
              e, λ, α]], {x, -xLeftRight, xLeftRight}, 
          PlotRange -> All, 
          PlotLabel -> 
            StringForm["E = `1`\nn =`2`\nλ =`3`\nα =`4`\n\
", PaddedForm[N[e], {7, 5}], PaddedForm[n, 2], 
              PaddedForm[λ, {4, 2}], PaddedForm[α, {2, 0}]], 
          ImageSize -> {250, 300}, AspectRatio -> 1.3, 
          PlotStyle -> 
            If[e == Ept[n, λ, α], {Hue[ .0], 
                Thickness[ .01]}, GrayLevel[0]]];

Animation

Do[doplotE[e], {e, EStart, EEnd, EStep}];

Superposition of 4 eigenfunctions

We will construct a superpositition of eigenfunctions with quantum numbers 0,1,2 and 5. The time operator  \(TraditionalForm\`e\^\(\(-i\)\ E\_n\ t\)\) will act onto this wavefunction so the animation shows the time evolution.

Let us assume that the parameters λ and α are integer numbers. Then the eigenvalues can either be an integer number or the half of one. It depends on the parameter λ if it is even or odd as well on the quantum number. All eigenfunctions with  integer eigenvalues have a period of π and all others, such with integer halves as eigenvalues have a period of 4π. So if a superposoition of m eigenfunction is constructed it's period generally is 4π unless the superposition does not include any odd state. After (Subscript[t, 1])=(Subscript[t, 0]) + 2  πthe wavefunction is it's parity transformed wavefunction. For real values of λ and α there is no periodicity anymore.

Set Parameters

λ = 7;
α = 1;
ns = {0, 1, 2, 5};

Hence these are the involved energies:

Ept[ns,λ,α]

coefficients and mean energy

coeff = 
        Table[1/Sqrt[Length[ns]], {i, 
            Length[ns]}];
meanE = 
        Dot[Abs[coeff]^2, Ept[ns, λ, α]];

Range

xLeftRight = 2;
StartTime = 0;
EndTime = 4 Pi;
TimeStep = Pi/64;
lower = meanE*1.02;
upper = meanE* .9;

Wavefunction

Ψ[n_, t_, x_, 
          e_, λ_, α_] := φ[n, x, 
            e, λ, α] (E^(-I Ept[n, λ, α] t));
superp[t_, x_, λ_, α_] := 
      Simplify[Sum[
          coeff[[j]]*Ψ[ns[[j]], t, x, 
              Ept[ns[[j]], λ, α], λ, α], {j, 1, Length [ns]}]]

define plot

doplot[t_, λ_, α_] := 
      Show[{Graphics[
            FilledPlot[{V[x, λ, α], 
                lower}, {x, -xLeftRight, xLeftRight}, 
              Fills -> GrayLevel[0.8], PlotPoints -> 150, 
              DisplayFunction -> Identity]], 
          Graphics[
            QArgColorPlot[
              Evaluate[
                superp[t, x, λ, α]], {x, -xLeftRight, 
                xLeftRight}, PlotPoints -> 120, 
              DisplayFunction -> Identity, QShiftPlot -> meanE, 
              QSquared -> True]]}, PlotRange -> {lower, upper}, 
        Frame -> True, Axes -> {True, None}, 
        PlotLabel -> StringForm["t =`1`", PaddedForm[N[t], {3, 2}]], 
        DisplayFunction -> $DisplayFunction]

Animation

Do[doplot[t, λ, α];, {t, StartTime, EndTime, 
        TimeStep}]