Eigenvalue problems

Main idea and formulation in the linear algebra

The word "eigenvalue" stems from the German word "Eigenwert" that can be translated into English as "Its own value" or "Inherent value". This is a value of a parameter in the equation or system of equations for which this equation has a nontrivial (nonzero) solution. Mathematically, the simplest formulation of the eigenvalue problem is in the linear algebra. For a given square matrix A one has to find such values of λ, for which the equation (actually the system of linear equations)

A.X==λX

has a nontrivial solution for a vector (column) X. Moving the right part to the left, one obtains the equation

(A-λI).X==0,

where I is the identity matrix having all diagonal elements one and nondiagonal elements zero. This matrix equation has nontrivial solutions only if its determinant is zero,

Det[A-λI]==0.

This is equivalent to a Nth order algebraic equation for λ, where N is the rank of the mathrix A. Thus there are N different eigenvalues (Subscript[λ, n]) (that can be complex), for which one can find the corresponding eigenvectors (Subscript[X, n]). Eigenvectors are defined up to an arbitrary numerical factor, so that usually they are normalized by requiring

(Subscript[X, n]^(T* )).(Subscript[X, n])==1,

where (X^(T* )) is the row transposed and complex conjugate to the column X.

It can be proven that eigenvectors that belong to different eigenvalues are orthogonal, so that, more generally than above, one has

(Subscript[X, m]^(T* )).(Subscript[X, n])==(Subscript[δ, mn]).

Here (Subscript[δ, mn]) is the Kronecker symbol,

(Subscript[δ, mn])=(Piecewise[{{1,, m==n}, {0,, m≠n.}}])

An important class of square matrices are Hermitean matrices that satisfy

(A^(T* ))==A.

Eigenvalues of Hermitean matrices are real. A real Hermitean matrix is just a symmetric matrix, (A^T)==A. For such matrices eigenvectors can be chosen real.

Matrix eigenvalue problem in Mathematica

Mathematica offers a solver for the matrix eivenvalue problem. If one is interested in eigenvalues only, one can use the command Eigenvalues[...]. Eigenvectors are computed by Eigenvectors[...], while both eigenvalues and eigenvectors are computed by the command Eigensystem[...]. Let us illustrate how it works for a real symmetric matrix

A=(({{a, b}, {b, -a}}));

Its eigenvalues are given by

Eigenvalues[A]

Its eigenvectors are given by

Eigenvectors[A]

that is,

(Subscript[X, i_]):=Eigenvectors[A][[i]]

These two eigenvectors are orthogonal to each other

(Subscript[X, 1]).(Subscript[X, 2]) //Simplify

However, they are not normalized

(Subscript[X, 1]).(Subscript[X, 1]) //Expand//Factor

To see which eigenvector corresponds to each eigenvalue, one has to use the command Eigensystem

ESys=Eigensystem[A]

The first part of this List are eigenvalues and the second part are eigenvectors. One can better see the correspondence in the form

TableForm[Transpose[ESys]]

Mathematica also solves matrix eigenvalue problems numerically, that is the only way to go for big matrices. For instance,

ESys=Eigensystem[A/.{a->1.,b->2.}]

The numerical eigenvectors

(Subscript[X, i_]):=ESys[[2]][[i]]

are orthonormal

(Subscript[X, 1]).(Subscript[X, 1]) 
(Subscript[X, 1]).(Subscript[X, 2])

For a complex Hermitean matrix eigenvalues are indeed real, although eigenvectors are complex

TableForm[Transpose[Eigensystem[(({{2, I}, {I, -2}}))]]]

({{-(Sqrt[3]), ({{I (-2+(Sqrt[3]))}, {1}})}, {(Sqrt[3]), ({{-I (2+(Sqrt[3]))}, {1}})}})

Eigenvalue problem for systems of linear ODEs on time

The importance of the eigenvalue problem in physics (as well as in engineering and other areas) is that it arises on the way of solution of systems of linear ordinary differential equations with constant coefficients. We have already obtained the solution for the harmonic oscillator on this way in the chapter on differential equations.

Every linear ODE or a system of ODEs can be represented in the basic matrix form with a constant matrix A

X'[t]+A.X[t]==0,

X being a vector. (We drop the inhomogeneous term.) Searching for the solution in the form

X[t]=(Subscript[X, 0])(e^(-λt)),

one arrives at the eigenvalue problem Eq. () with X⇒(Subscript[X, 0]). After finding eigenvalues (Subscript[λ, n]) and normalized eivenvectors (Subscript[X, 0n]) by linear algebra, one can write down the general solution of the equation as a linear superposition of all these solutions,

X[t]=(∑_(n=1)^N)(Subscript[C, n])(Subscript[X, 0n])(e^(-(Subscript[λ, n])t)),

where (Subscript[C, n]) are integration constants that can be found from the initial conditions.

Eigenvalue problems for PDE

In physical problems described by partial differential equations, eigenvalue problems usually arise due to boundary conditions.

Standing waves in a pipe

Consider, as an example, the wave equation for the pressure change (see Waves) in the 1d region 0≤x≤L,

(∂_t^2)δP-(c^2)(∂_x^2)δP.

If we consider a pipe with both ends open to the atmosphere, the boundary conditions are

δP[0,t]==0 ,             δP[L,t]==0

because the pressure at the open end (practically) merges with the constant atmospheric pressure. Searching for δP in the form

δP[x,t]=ψ[x]Cos[ωt+(Subscript[ϕ, 0])]

one obtains the stationary wave equation

(∂_x^2)ψ+(k^2)ψ==0,             k=ω/c,

k being the wave vector. This is an eigenvalue problem because this equation has nontrivial solutions that satisfy the boundary conditions only for some values of k and thus of ω. The general solution of the ODE above is

ψ[x]=(Subscript[C, 1])Sin[kx]+(Subscript[C, 2])Cos[kx].

Since Cos[kx] does not satisfy the BC at x=0, the solution simplifies to

ψ[x]=C Sin[kx]

that describes a standing wave. Next, the BC at x=L requires Sin[kL]=0 from which one obtains the eigenvalues of the wave vector

k=(Subscript[k, n])=((πn)/L),      n=1,2,3,⃛

In terms of the wave length of the standing wave λ=2π/k one has

(Subscript[λ, n])=((2L)/n),      n=1,2,3,⃛

Standing wave with n=1 is called fundamental wave, whereas those with n=2,3,⃛ are called overtones or harmonics. For the frequencies of these waves f=ω/(2π)=c/λ one has

(Subscript[f, n])=(cn/(2L)),      n=1,2,3,⃛

The general solution of the wave equation for a pipe with both ends open is a linear superposition of all these solutions,

δP[x,t]=(∑_(n=1)^∞)(Subscript[C, n])Sin[(Subscript[k, n])x]Cos[(Subscript[ω, n])t+(Subscript[ϕ, n])],      (Subscript[k, n])=((Subscript[ω, n])/c)=((πn)/L),      n=1,2,3,⃛

The coefficients (Subscript[C, n]) and phases (Subscript[ϕ, n]) in this solution are arbitrary. Similar results can be obtained for a pipe with both ends closed (such as flute). If one end is closed and one is open (clarinet), the solution is somewhat different and only odd overtones exist. (Excersize). While the phases (Subscript[ϕ, n]) are irrelevant for our ears (Ohm's phychoacoustical law), the amplitudes (Subscript[C, n]) define the quality of the sound via the relative weight of the overtones. This depends on how the music instrument is constructed and how it is played.

Eigenfunctions corresponding to different eigenvalues are orthogonal,

(∫_0^L)Sin[(Subscript[k, m])x]Sin[(Subscript[k, n])x]dx=(L/2)(Subscript[δ, mn]).

One can define the normalized eigenfunctions

(Subscript[ψ, n])[x]=(Sqrt[(2/L)])Sin[(Subscript[k, n])x]

that satisfy

(∫_0^L)(Subscript[ψ, m])[x](Subscript[ψ, n])[x]dx=(Subscript[δ, mn]).

Several lowest eigenfunctions are plotted below.

L=1;          ψnx[n_,x_]=(Sqrt[(2/L)])Sin[((π n)/L)x];
Plot[Table[ψnx[n,x],{n,1,4}],{x,0,L}]

General formulation of the eigenvalue problem for PDE

In general, the eigenvalue problem for PDE can be formulated in the form

(OverHat[L])ψ[r]==λψ[r],

where (OverHat[L]) is a differential operator. The best example is the stationary Schrödinger equation for a quantum particle

(OverHat[H])ψ==Eψ,                      (OverHat[H])=-(((ℏ^2)Δ)/(2m))+U[r],

where (OverHat[H]) is the Hamilton operator. For a particle bound to a potential well, this equation has solutions satisfying boundary conditions only for select discrete values of the energy, E=(Subscript[E, n]),  n=0,1,2,⃛ These energy eigenvalues are labeled by the index n in the increasing order, and n=0 corresponds to the minimal energy, the so-called ground-state energy.

When a quantum particle undergoes a transition from the energy level (energy eigenstate) m to a lower energy level n, it emits a quantum of electromagnetic radiation (a photon), to conserve the total energy. The energy of the photon ε is related to its frequency ω as ε=ℏω (Max Planck), and from the energy conservation follows ℏω = (Subscript[E, m])-(Subscript[E, n]). Spectroscopic investigations of the light emitted (and absorbed) by simple quantum systems such as hydrogen atoms have shown that only particular discrete frequencies obeying certain simple laws are present in the spectrum. This was at a  contradiction with classical mechanics and electrodynamics that predicted a continuous energy spectrum. Physicists started to look for a mathematical tool that describes discreteness and made use of the eigenvalue problem. Werner Heisenberg proposed his matrix quantum mechanics operating with eigenvalues of matrices describing quantum systems. Erwin Schrödinger obtained his famous Schrödinger equation in which discrete energy levels of a quantum system arise as eigenvalues of a differential operator. The approaches of Heisenberg and Schrödinger are equivalent.

Solving the eigenvalue problem for a general differential operator, especially in 2d and 3d, is not an easy task, and there is no Mathematica solver for this problem at the moment. If the differential operator and boundary conditions possess special properties such as symmetry (e.g., spherical or cylindrical symmetry), the solution can be searched for in the form of a product of functions of different variables, so that the solution simplifies. This case is called "separation of variables".

Similarly to eigenvectors of matrices, eigenfunctions of differential operators can be complex. It can be shown that eigenfunctions corresponding to different eigenvalues are orthogonal. Also eigenfunctions corresponding to discrete eigenvalues can be normalized. Allowing for complex eigenfunctions, the orthonormality conditions takes the form

∫(Subscript[ψ, m]^( * ))[r](Subscript[ψ, n])[r]dV=(Subscript[δ, mn]),

where integration is carried out over the volume. An important mathematical theorem states that the whole set of eigenfunctions (Subscript[ψ, n])[r] of a differential operator satisfying boundary conditions forms a complete basis, so that any function f[r] satisfying these BC can be expanded over the eigenfunctions as

f[r]=(∑_(n=1)^∞)(Subscript[c, n])(Subscript[ψ, n])[r].

The expansion coefficients can be found using the orthonormality condition that yields

∫(Subscript[ψ, m]^( * ))[r]f[r]dV=∫(∑_(n=1)^∞)(Subscript[c, n])(Subscript[ψ, m]^( * ))[r](Subscript[ψ, n])[r]dV=(∑_(n=1)^∞)(Subscript[c, n])(Subscript[δ, mn])=(Subscript[c, m]).

Quantum particle in a one-dimensional potential box

The stationary Schrödinger equation () for a quantum particle in a one-dimensional potential box 0≤x≤L has the form

(∂_x^2)ψ+(k^2)ψ==0,             k=(Sqrt[2mE])/ℏ.

Everywhere in the box U=0 while outside the box U=∞, so that the wave function outside the box is zero and by continuity it has to be zero at the box walls,

ψ[0,t]==0 ,             ψ[L,t]==0.

Mathematically this eigenvalue problem is equivalent to the problem of standing waves in a pipe considered above, except for another definition of the wave vector k. It has the same solution for the normalized eigenfunctions

(Subscript[ψ, n])[x]=(Sqrt[(2/L)])Sin[(Subscript[k, n])x],              (Subscript[k, n])=((πn)/L),      n=1,2,3,⃛

The corresponding energy eigenvalues are given by

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

These energy levels are discrete, and a transition m→n should be accompanied by emission of a photon of frequency

(Subscript[ω, mn])=(((Subscript[E, m])-(Subscript[E, n]))/ℏ)=(((π^2)ℏ )/(2(mL^2)))((m^2)-(n^2)).

Measuring all these frequencies, one can figure out the energy levels.

Quantum particle in a three-dimensional potential box and density of states

Energy levels

Let us consider now a quantum particle in 3d potential box 0≤x≤(Subscript[L, x]), 0≤y≤(Subscript[L, y]), 0≤z≤(Subscript[L, z]). One can easily check that the solution with separated variables

ψ[x,y,z]∝Sin[(Subscript[k, x])x]Sin[(Subscript[k, y])y]Sin[(Subscript[k, z])z]

satisfies the stationary Schrödinger equation

Δψ+(k^2)ψ==0

since (Subscript[k, x]^2)+(Subscript[k, y]^2)+(Subscript[k, z]^2)=(k^2). Also this solution satisfies the BC at x,y,z=0, where ψ vanishes. Taking into account three BC at the other sides of the box finally yields

(Subscript[ψ, (Subscript[n, x])(Subscript[n, y])(Subscript[n, z])])[x,y,z]=((2^(3/2))/(Sqrt[(Subscript[L, x])(Subscript[L, y])(Subscript[L, z])]))Sin[(Subscript[k, (Subscript[xn, x])])x]Sin[(Subscript[k, (Subscript[yn, y])])y]Sin[(Subscript[k, (Subscript[zn, z])])z],

where the quantized wave vectors are given by

(Subscript[k, (Subscript[αn, α])])=((π )/(Subscript[L, α]))(Subscript[n, α]),              (Subscript[n, α])=1,2,3,⃛,              α=x,y,z.

The discrete energy eigenvalues are parametrized by three integers,

(Subscript[E, (Subscript[n, x])(Subscript[n, y])(Subscript[n, z])])=(((π^2)(ℏ^2))/(2m))(((Subscript[n, x]^2)/(Subscript[L, x]^2))+((Subscript[n, y]^2)/(Subscript[L, y]^2))+((Subscript[n, z]^2)/(Subscript[L, z]^2))).

Another important case of separation of variables is the hydrogen atom which possesses a spherical symmetry. Rewriting Schrödinger equation in spherical coordinates allows to separate variables r, θ, and ϕ and find analytical solutions for the energy levels and corresponding wave functions.

Density of states

For potential boxes of sufficiently large sizes the quantum energy levels are very close to each other, so that one can consider the density of states ρ[E] as the number of energy levels d(Subscript[N, E]) in the energy interval dE,

d(Subscript[N, E])=ρ[E]dE.

For the 3d box the density of states is given by

ρ[E]=(V/(((2π))^2))(((((2m)/(ℏ^2))))^(3/2))(Sqrt[E]),               V=(Subscript[L, x])(Subscript[L, y])(Subscript[L, z]).

It is instructive to plot a hystogram of the discrete energy levels given by Eq. () and to compare it with the formula above.

Lx=91;                    Ly=100;                   Lz=110;
V=Lx Ly Lz;
nxMax=Lx;             nyMax=Ly;             nzMax=Lz;  
nLevels=nxMax nyMax nzMax 
ℏ=1;     M=1.;
ρ[EE_]=(V/(((2π))^2))(((((2M)/(ℏ^2))))^(3/2))(Sqrt[EE]);
Enxnynz[nx_,ny_,nz_]=(((π^2)(ℏ^2))/(2M))(((nx^2)/(Lx^2))+((ny^2)/(Ly^2))+((nz^2)/(Lz^2)));
EMin=0;
EMax=(1/3)Enxnynz[nxMax,nyMax,nzMax];
NE=100;   (* Number of energy bins \[DifferentialD]E *)
dE=(EMax-EMin)/NE;  (* Width of energy bin \[DifferentialD]E *)
EList=Flatten[Table[Enxnynz[nx,ny,nz],{nx,0,nxMax},{ny,0,nyMax},{nz,0,nzMax}]];  (* List of energy eigenvalues *)
EBinCounts=(1/dE)BinCounts[EList,{EMin,EMax,dE}];  (* Number of eigenvalues in energy bins \[DifferentialD]E divided by \[DifferentialD]E *)
DOSList=Table[{nE dE,EBinCounts[[nE]]},{nE,0,NE}];   (* Same with the energy values for the energy bins *)
Show[
ListPlot[DOSList,PlotStyle->{Orange,PointSize[0.02]},AxesLabel->{"E","ρ[E]"}],
Plot[ρ[EE],{EE,EMin,EMax},PlotStyle->{Black,Thick},AxesLabel->{"E","ρ[E]"}]
]

Quantum harmonic oscillator

Consider a one-dimensional quantum harmonic oscillator that is described by the Hamilton function

H=((p^2)/(2m))+U[x],      U[x]=(((Subscript[mω, 0]^2)(x^2))/2)

where (Subscript[ω, 0]) is the oscillator's frequency. The stationary Schrödinger equation Eψ==(OverHat[H])ψ for the oscillator reads

(((ℏ^2)/(2m))((d^2)/(d(x^2)))+E-U[x])ψ==0.

This is a second-order linear ODE with a variable coefficient. Boundary conditions require that ψ decreases to zero at x=±∞. Solution of such differential equations is a more complicated task than for those with constant coefficients, and the analytical solution is usually expressed via special functions, if they are known. In this case an involved analytical solution yields the energy levels

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

and normalized eigenfunctions

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

where (Subscript[H, n])[x] are Hermite polynomials defined by

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

In Mathematica one can calculate Hermite polynomials straightforwardly and fast enough as

f[x_]=(E^(-(x^2)));
MyHermiteH[n_,x_]:=(((-1))^n)(E^(x^2))Derivative[n][f][x]//Simplify

MyHermiteH[0,x]
MyHermiteH[1,x]
MyHermiteH[2,x]
MyHermiteH[3,x]

etc., although there is a built-in function HermiteH[n, x]

HermiteH[0,x]
HermiteH[1,x]
HermiteH[2,x]
HermiteH[3,x]

One can see that eigenfunctions of the quantum harmonic oscillator are, in fact, elementary functions. Still, Hermite and many other polynomials are put into the class of special functions by tradition. Let us plot the energy eigenvalues and eighefunctions together with the potential.

Parameters={m->1,ℏ->1,(Subscript[ω, 0])->1};
  U[x_]=((m (Subscript[ω, 0]^2)(x^2))/2)/.Parameters
En[n_]=ℏ (Subscript[ω, 0])(n+1/2)/.Parameters
ψnx[n_,x_]=(1/(Sqrt[(2^n)n!]))(((((m (Subscript[ω, 0]))/(π ℏ))))^(1/4))Exp[-((m (Subscript[ω, 0])(x^2))/(2ℏ))]HermiteH[n,(Sqrt[((m (Subscript[ω, 0]))/ℏ)])x]/.Parameters

xMax=6;   nMax=16;
Show[
Plot[U[x],{x,-xMax,xMax},PlotStyle->{Black,Thick}], (* Potential energy *)
Plot[Table[En[n],{n,0,nMax}],{x,-xMax,xMax},PlotStyle->Dashed], (* Energy eigenvalues *)
Plot[Evaluate[Table[0.7ψnx[n,x]+En[n],{n,0,nMax}]],{x,-xMax,xMax},PlotStyle->Black] (* Eigenfunctions at the height of their eigenvalues *)
]

It is convenient to put eigenfunctions at the height of their eigenvalues for better assosiation and separation. The numeric factor at eigenfunctions has been introduced to avoid overlap. Command Evaluate makes plotting faster by first evaluating eigenfunctions and then plotting them. Otherwise they are evaluated a every plot point.

One can see that even and odd eigenfunctions alternate. Similarly to plane quantum waves discussed in the chapter on Schrödinger equation, wave functions oscillate in the classically allowed region E-U[x]>0 and exponentially decrease in the classically prohibited region E-U[x]<0.

As we know, a general solution of a second-order differential equation is a superposition of two independent functions. As in the case of a plane wave, in the classically prohibited region one of these functions is exponentially increasing while the other is exponentially decreasing. In this and similar problems, eigenfunctions decrease on both sides of the potential well. Other functions increase in the prohibited region and have to be rejected. It turns out that the solutions of the stationary Schrödiger equation that decrease on both sides exist only if E=(Subscript[E, n]). For other energy values, both solutions of the DE increase at least on one side of the well and thus violate the boundary conditions.

This fact can be used to numerically find eigenvalues and eigenfunctions by the shooting method. One can start at some (Subscript[x, left]) left from the well setting any boundary conditions, say ψ[(Subscript[x, left])]=1 and  ψ'[(Subscript[x, left])]=0 and solve the DE for any E to the right until some (Subscript[x, right]). If ψ[(Subscript[x, right])] is exponentially small, E is one of the energy eigenvalues, and if ψ[(Subscript[x, right])] is exponentially large, E is not an eigenvalue.

Parameters={m->1,ℏ->1,(Subscript[ω, 0])->1};
U[x_]=((m (Subscript[ω, 0]^2)(x^2))/2)/.Parameters ;
En[n_]=ℏ (Subscript[ω, 0])(n+1/2)/.Parameters 

xRight=10;    xLeft=-xRight;    

IConds={ψ[xLeft]==1,ψ'[xLeft]==0};

Sol[EE_]:=NDSolve[Join[{((ℏ^2)/(2m))ψ''[x]+(EE-U[x])ψ[x]==0 /.Parameters },IConds],ψ,{x,xLeft,xRight}];
ψx[EE_?NumericQ,x_]:=ψ[x]/.Sol[EE][[1]](* Without ?NumericQ FindMinimum below crashes *)

Plotting shows that for our parameter choice the energy eigenvalues indeed are 0.5, 1.5, 2.5, etc., that is, (Subscript[E, n])=n+1/2.

Plot[Log[Abs[ψx[EE,xRight]]],{EE,0,4}]

Now one can plot eigenfunctions for the found eigenvalues

xMax=5;  xMin=-xMax;
Plot[Evaluate[ψx[0.5,x]],{x,xMin,xMax},PlotRange->All]
Plot[Evaluate[ψx[1.5,x]],{x,xMin,xMax},PlotRange->All]
Plot[Evaluate[ψx[2.5,x]],{x,xMin,xMax},PlotRange->All]

These are essentially the same eigenfunctions as the analytical result considered above, only they are not normalized.

One could find the exact eigenvalues by minimizing Log[Abs[ψx[EE, xRight]]] starting from different values of EE.

FindMinimum[Log[Abs[ψx[EE,xRight]]],{EE,0}]

Finding exact eigenvalues by numerical minimization is important because this shooting method can be applied to potentials that do not allow for an analytical solution. For instance, if one replaces (x^2)⇒(x^4) in U[x], one obtains, for the same parameters,

Plot[Log[Abs[ψx[EE,xRight]]],{EE,0,5}]

These nonequidistant energy eigenvalues are not described by any analytical formula, so that finding them numerically is important. Note that power potentials U[x]∝((((x/(Subscript[x, 0]))))^n) approach the potential box -(Subscript[x, 0])<x<(Subscript[x, 0]) in the limit n->∞, where the analytical results have been found above.

Tunneling of a quantum particle

Let us now consider the potential U[x]=-(ax^2)+(bx^4) with a, b>0 that has two symmetric minima at x=±(Subscript[x, 0])=±(Sqrt[a/(2b)]). It is convenient to represent this potential in the form

U[x]=((m (Subscript[ω, 0]^2))/8)(((((x^2)-(Subscript[x, 0]^2)))^2)/(Subscript[x, 0]^2)).

Its second derivative at the minima is given by

U[x_]=((m (Subscript[ω, 0]^2))/8)(((((x^2)-(x0^2)))^2)/(x0^2));
(∂_(x,x))U[x]/.x->x0

that is the same as for a harmonic oscillator with frequency (Subscript[ω, 0]). A classical particle of mass m would perform harmonic oscillations with this frequency near the minima of this potential.

Now the solution

Parameters={m->1,ℏ->1,(Subscript[ω, 0])->1,(Subscript[x, 0])->3.5};
U[x_]=((m (Subscript[ω, 0]^2))/8)(((((x^2)-(Subscript[x, 0]^2)))^2)/(Subscript[x, 0]^2))/.Parameters ;
xRight=7;    xLeft=-xRight;    

IConds={ψ[xLeft]==1,ψ'[xLeft]==0};

Sol[EE_]:=NDSolve[Join[{((ℏ^2)/(2m))ψ''[x]+(EE-U[x])ψ[x]==0 /.Parameters },IConds],ψ,{x,xLeft,xRight}];
ψx[EE_?NumericQ,x_]:=ψ[x]/.Sol[EE][[1]]; (* Without ?NumericQ FindMinimum below crashes *)

xMax=6;  xMin=-xMax;    EMax=2;
Plot[U[x],{x,xMin,xMax}];
Plot[Log[Abs[ψx[EE,xRight]]],{EE,0,EMax},PlotRange->All]

reveals that low-lying energy eigenvalues form doublets. To see the lowest doublet better one can zoom in:

Plot[Log[Abs[ψx[EE,xRight]]],{EE,0.475,0.479},PlotRange->All]

One can find both eigenvalues

E0=EE/.FindMinimum[Log[Abs[ψx[EE,xRight]]],{EE,0}][[2]]
E1=EE/.FindMinimum[Log[Abs[ψx[EE,xRight]]],{EE,0.5}][[2]]

And their spitting

Δ=E1-E0

To normalize the eigenfunctions, one at first has to calculate the norms of both eigenfunctions

N0=N[(∫_xLeft^xRight)((ψx[E0,x])^2)\[DifferentialD]x];         N1=N[(∫_xLeft^xRight)((ψx[E1,x])^2)\[DifferentialD]x];
ψx0[x_]=(1/(Sqrt[N0]))ψx[E0,x];              ψx1[x_]=(1/(Sqrt[N1]))ψx[E1,x];

Below the two normalized eigenfunctions are plotted together with the scaled potential energy and a energy levels.

xMax=7;  xMin=-xMax;
Show[
Plot[Evaluate[{ψx0[x],ψx1[x]}],{x,xMin,xMax},PlotRange->All,PlotStyle->{Blue,Red}],Plot[0.1U[x],{x,-xMax,xMax},PlotStyle->{Black,Thick},PlotRange->{0,1}], (* Potential energy *)
Plot[0.1E0,{x,-xMax,xMax},PlotStyle->{Black,Dashed}],
Plot[0.1E1,{x,-xMax,xMax},PlotStyle->{Black,Dashed}]
]

One can see that the eigenfunctions of the ground-state doublet are even and odd. According to a quantum-mechanical theorem, the eigenfunction of the ground state (Subscript[ψ, 0])[x] does not change its sign on x. They are localized in the wells with little probalility to be under the barrier, x~0. These two eigenfunctions can be interpreted as even and odd linear combinations of the states in the left and right wells that are hybridized via tunneling, an essentially quantum phenomenon. With increasing the width or/and the height of the barrier, the separation between energy eigenvalues in tunneling doublets (tunneling splitting) decreases very fast, and it becomes difficult to calculate the splitting numerically. Fortunately, in this case tunneling splittting can be calculated analytically with the help of the quasiclassical approximation.

Matrix quantum mechanics

Eigenfunctions (Subscript[χ, n])[r] of a Hamiltonian satisfying boundary conditions form a complete basis that can be used to expand the solution of Schrödinger equation. The Hamiltonian to which the eigenfunctions belong can be the Hamiltonian of the problem or any other Hamiltonian. Let us call it (Subscript[(OverHat[H]), 0]). The eigenvalue expansion has the form

Ψ[r,t]=(∑_(n=0)^∞)(Subscript[a, n])[t](Subscript[χ, n])[r],

similarly to the pressure expansion for the standing waves in a pipe, Eq. (). This expansion transforms the time-dependent Schrödinger equation to a system of ODEs for the coefficients (Subscript[a, n])[t]. Substituting Eq. () into the Schrödinger equation, one obtains

(∑_(n=0)^∞)iℏ(Subscript[(OverDot[a]), n])[t](Subscript[χ, n])[r]==(∑_(n=0)^∞)(Subscript[a, n])[t](OverHat[H])(Subscript[χ, n])[r],

Multiplicating this equation (Subscript[χ, m]^( * ))[r]and integrating it over the space taking into account orthonormality of the eigenfunctions, one obtains

iℏ(Subscript[(OverDot[a]), m])[t]==(∑_(n=0)^∞)(Subscript[H, mn])(Subscript[a, n])[t],           (Subscript[H, mn])≡∫(Subscript[χ, m]^( * ))[r](OverHat[H])(Subscript[χ, n])[r]dV.

This is a system of a (generally) infinite number of linear homogeneous ODEs that are coupled to each other via the matrix elements (Subscript[H, mn]). If  (OverHat[H])=(Subscript[(OverHat[H]), 0]), one has (Subscript[H, mn])=(Subscript[ϵ, m])(Subscript[δ, mn]), where (Subscript[ϵ, m]) are eigenvalues corresponding to the eigenfunctions (Subscript[χ, m]). In this case equations decouple, iℏ (Subscript[(OverDot[a]), m]) [t] ==(Subscript[ϵ, m])(Subscript[a, m])[t], and the solution is (Subscript[a, m])[t]=(Subscript[a, m])[0]Exp[-((Subscript[iϵ, m])/ℏ)t]. However, this is a trivial case, and usually the Hamiltonian of the system (OverHat[H]) is more complicated than (Subscript[(OverHat[H]), 0]), so that the matrix (Subscript[H, mn]) is non-diagonal.

Eq. () can be written in the matrix form

iℏ(OverDot[A])[t]==H.A[t] ,

where A[t] is the vector composed of (Subscript[a, m])[t]. If matrix H is time-independent, its solution is the matrix exponential

A[t]=Exp[-((iH)/ℏ)t].A[0]

describing time evolution of the initial state A[0]. This formula can be used to obtain the solution with the help of Mathematica. It should be noted that Heisenberg proposed his matrix quantum mechanics independently of Schrödinger equation.

Stationary states are described by the time dependence

A[t]=Exp[-((iE)/ℏ)t]C.

Substituting this into Eq. (), one arrives at the matrix eigenvalue problem

EC==H.C

or, explicitly,

(Subscript[Ec, m])==(∑_(n=0)^∞)(Subscript[H, mn])(Subscript[c, n])

for the energy eigenvalues E that is equivalent to Eq. (). If matrix H is infinite, there is an infinite number of eigenvalues. Practically, the functional basis used to build matrices is cut, so that only N lowest basis eigenfunctions of (Subscript[(OverHat[H]), 0]) are used and the matrices are of rank N+1. Some quantum-mechanical problems such that spin problems are described by finite matrices from the very beginning. An advantage of the matrix formalism is that it allows to obtain all eigenvalues and eigenfunctions as a result of a single computation.

Particle in a double - well potential

Let us apply the matrix formalism to the problem of the quartic double-well potential, Eq. (), and choose the basic formed by the eigenstates of the Hamiltonian (Subscript[(OverHat[H]), 0]) of the harmonic oscillator, Eq. (). The Hamiltonian of the problem can be written as

(OverHat[H])=(Subscript[(OverHat[H]), 0])+(OverHat[V]),                (OverHat[V])=((m (Subscript[ω, 0]^2))/2)((((((x^2)-(Subscript[x, 0]^2)))^2)/(4(Subscript[x, 0]^2)))-(x^2)).

Such a separation of the Hamiltonian into two parts is typical, and in the case (OverHat[V]) is small, there are perturbative analytical methods based on expansion in powers of (OverHat[V]).  However, in our case (OverHat[V]) is not small and the problem will be solved numerically. With the help of (Subscript[(OverHat[H]), 0mn])=(Subscript[ϵ, m])(Subscript[δ, mn]) Eq. () becomes

(E-(Subscript[ϵ, m]))(Subscript[c, m])==(∑_(n=0)^∞)(Subscript[V, mn])(Subscript[c, n]),

where (Subscript[ϵ, m])=(Subscript[ℏω, 0])(m+1/2) are the energy eigenvalues of the harmonic oscillator and

(Subscript[V, mn])≡(∫_(-∞)^∞)(Subscript[χ, m])[x](OverHat[V])(Subscript[χ, n])[x]dx,

where (Subscript[χ, n])[x] are given by Eq. (). There is an analytical expression for these matrix elements for general m and n, however, Mathematica 7 cannot find it. Still, Mathematica can analytically calculate (Subscript[V, mn]) for any particular values of m and n that we need to build the Hamiltonian matrix. We will be labeling the eigenvalues of the problem with the index ν, so that (Subscript[E, ν]) are energy eigenvalues, (Subscript[c, νn]) are the basis expansion coefficients, and

(Subscript[ψ, ν])[x]≡(∑_(n=0)^∞)(Subscript[c, νn])(Subscript[χ, n])[x]

are eigenfunctions of the problem. Below is Mathematica code that solves the matrix eigenvalue problem for a quartic potential using the harmonic oscillator basis.

(* Definitions *)
NBasis=20;  (* Number of the lowest oscillator eigenfunctions in the basis *)
M=1;          (* Mass of the particle *)
ℏ=1;          (* Planck's constant *)
(Subscript[ω, 0])=1;       (* Frequency of oscillation near the minima of U[x] *)
(Subscript[x, 0])=3.5;   (* Positions of the minima of U[x]  *)
U[x_]=((M (Subscript[ω, 0]^2))/8)(((((x^2)-(Subscript[x, 0]^2)))^2)/(Subscript[x, 0]^2));    (* Potential energy of the particle *)
U0[x_]=((M (Subscript[ω, 0]^2)(x^2))/2);     (* Potential energy of the harmonic oscillator *)
ϵn[n_]=ℏ (Subscript[ω, 0])(n+1/2);   (* Energy levels of the harmonic oscillator *)
χnx[n_,x_]=(1/(Sqrt[(2^n)n!]))(((((M (Subscript[ω, 0]))/(π ℏ))))^(1/4))Exp[-((M (Subscript[ω, 0])(x^2))/(2ℏ))]HermiteH[n,(Sqrt[((M (Subscript[ω, 0]))/ℏ)])x];  (* Eigenfunctions of the harmonic oscilator *)
V[x]=U[x]-U0[x];    (* Perturbation Hamiltonian *)
Vmn[m_,n_]:=If[Abs[m-n]==0||Abs[m-n]==2||Abs[m-n]==4,(∫_(-∞)^∞)χnx[m,x]V[x]χnx[n,x]\[DifferentialD]x,0];    (* Most of matrix elements of V[x] are zero, so do not calculate them *)
Hmn[m_,n_]:=ϵn[n]KroneckerDelta[m,n]+Vmn[m,n];  (* Elements of Hamiltonian matrix *)

(* Creating Hamiltonian matrix and solving matrix eigenvalue problem *)
Timing[ HMatrix=Table[Hmn[m,n],{m,0,NBasis},{n,0,NBasis}]; ]
TableForm[Take[HMatrix,9,9]]   (* Show part of the Hamiltonian matrix *)
Timing[ ES=Eigensystem[HMatrix]; ]
ES=Transpose[Sort[Transpose[ES]]]; (* Eigensystem sorted in accending order *)
EVal=First[ES]    (* Eigenvalues *)
EVecT=Last[ES];  (* This is a matrix of eigenvectors lying horizontally *)
EVec=Transpose[EVecT]; (* This is a matrix of eigenvectors standing vertically *) 
TableForm[Chop[Take[EVecT.EVec,9,9]]]  (* Check that eigenvectors are orthonormal *)

(* Finalizing the solution *)
Eν[ν_]:=EVal[[ν+1]];      (* Energy eigenvalues of the problem *)
Δ=Eν[1]-Eν[0];         (* Ground-state splitting *)
c[ν_,n_]:=EVecT[[ν+1,n+1]]    (* Expansion coefficients of the eigenfunctions of the problem over the harmonic oscillator basis *)
ψνx[ν_,x_]:=(∑_(n=0)^NBasis)c[ν,n]χnx[n,x]     (* Eigenfunctions of the problem *)

The splitting of the ground-state doublet

Δ=Eν[1]-Eν[0]

is in an excellent accordance with the value obtained above by the shooting method. The eigenfunctions of the ground-state doublet, plotted together with the scaled potential energy and a group of low-lying energy levels

xMax=7;  xMin=-xMax;
Show[
Plot[Evaluate[{ψνx[0,x],ψνx[1,x]}],{x,xMin,xMax},PlotRange->All,PlotStyle->{Blue,Red}],
Plot[0.1U[x],{x,-xMax,xMax},PlotStyle->{Black,Thick},PlotRange->{0,1}], (* Potential energy *)
Plot[Table[0.1Eν[ν],{ν,0,12}],{x,-xMax,xMax},PlotStyle->{Black,Dashed}]
]

are the same as above, except for the irrelevant signs. The splitting of the ground-state doublet is not seen in this plot, whereas the spliting of the second doublet is seen. The states above the barrier are not grouped into doublets, and their separation increases with energy.

Sparce matrices

We have seen above that the main time of the calculation goes to fill the Hamiltonian matrix with elements, whereas finding its eigensystem is very fast. Filling a matrix of rank N requires (N^2) operations, even if many its elements are zero. In our Hamiltonian matrix, only the elements with |m-n|=0,2,4 are nonzero. Such matrices in which most of the elements are zero are called sparce matrices. In Mathematica a special fast method of filling sparce matrices is employed. First, a matrix with all elements equal to zero is created by one step. Then nonzero elements are inserted at their places. Actually, it is only nonzero elements and the matrix dimensions that are stored in the memory. Thus, both memory usage and speed are dramatically improved. In the case of nearly-diagonal matrices, there are ~N nonzero elements, so that the time to fill such a sparce matrix if proportional to N rather than to (N^2).

Mathematica command to create a sparce matrix object is SparceArray. Let us create a 5×5 matrix with diagonal elements equal to 1 and subdiagonal elements equal to (Sqrt[2]).

MySparceMatrix=SparseArray[{{i_,j_}/;Abs[i-j]==0->1,{i_,j_}/;Abs[i-j]==1->(Sqrt[2])},{5,5}]

Here /; means "under the condition". The output object is named but its content is hidden because it is not a standard list. It can be visualized with the command MatrixForm:

MatrixForm[MySparceMatrix]

One can apply all standard commands such as Eigensystem to spare matrix objects. Let us now compare the time needed to create the Hamiltonian matrix above as a standard (dense) matrix and sparse matrix. (The preceding code must be run before so that Mathematica knows definitions.) For NBasis = 100 one obtains

NBasis=40;

Timing[ HMatrix=Table[Hmn[m,n],{m,0,NBasis},{n,0,NBasis}]; ]
MatrixForm[Take[HMatrix,5,5]] 

Timing[ HMatrix=SparseArray[{{m_,n_}/;Abs[m-n]==0||Abs[m-n]==2||Abs[m-n]==4->Hmn[m-1,n-1]},{NBasis+1,NBasis+1}]; ]
MatrixForm[Take[HMatrix,5,5]]

For higher NBasis the gain in using sparce matrices will theoretically increase. This example is not the best, however, because the difference in speed between the two methods is marginal. The reason is that the standard method also uses the optimized definition of the matrix elements, in which their lengthy calculation via the integral is done only in the case they are nonzero, otherwise zeros are being put into the Hamiltonian matrix. A real gain in using sparce matrix objects can be achieved only for really big matrices in which calculation of nonzero elements is fast.

Quantum rotator and pendulum

Let us apply matrix quantum mechanics to the  problem of a quantum pendulum that is described by the Hamiltonian

(OverHat[H])=(((ℏ^2)((OverHat[l])^2))/(2I))+U[ϕ],       (OverHat[l])≡-i(∂/(∂ϕ)),        U[ϕ]=Mga(1-Cos[ϕ]),

where ℏ(OverHat[l]) is the quantum-mechanical angular momentum operator for a rotation around a fixed axis by the angle ϕ, I is the moment of inertia, M is the mass, and a is the distance between the center of mass and the pivot point. Classical pendulum performs oscillations for E<2Mga and rotations in positive or negative direction for E>2Mga. The frequency of small oscillations near the minimum of the potential energy is

(Subscript[ω, 0])=(Sqrt[(Mga/I)]).

In the case U[ϕ]=0 pendulum simplifies to rotator.

Since adding 2πn to the angle does not change the physical position of the pendulum, the wave function satisfies the periodicity condition

Ψ[ϕ,t]=Ψ[ϕ+2πn,t],       n=±1,±2,⃛

Thus we consider the wave functions in the interval -π≤ϕ≤ϕ and require Ψ[-π, t] = Ψ[π, t]. A suitable complete set of basis functions satisfying this boundary conditions consists of the eigenfunctions of the angular momentum

(OverHat[l])(Subscript[χ, m])=(Subscript[mχ, m])

that are given by

(Subscript[χ, m])=(1/(Sqrt[2π])) (e^(imϕ)),          m=0,±1,±2,⃛

One can check that these functions are orthonormal. They are eigenfunctions of the stationary Schrödinger equation Eψ==(OverHat[H])ψ for the rotator with the rotational energies

(Subscript[E, m])=(((ℏ^2)(m^2))/(2I)).

For the pendulum, these functions are not eigenfunctions but can be used as a basis to formulate the matrix Schrödinger equation (). The matrix elements of the Hamiltonian are given by

(Subscript[H, mn])=(∫_(-π)^π)(Subscript[χ, m]^( * ))[ϕ](OverHat[H])(Subscript[χ, n])[ϕ]dϕ=(∫_(-π)^π)(e^(-imϕ))(-((ℏ^2)/(2I))((∂^2)/(∂(ϕ^2)))+Mga(1-Cos[ϕ]))(e^(inϕ))dϕ =((((ℏ^2)(m^2))/(2I))+Mga)(Subscript[δ, mn])-(Mga/(2π))(∫_(-π)^π)(e^(-imϕ))Cos[ϕ](e^(inϕ))dϕ.

Using Cos[ϕ]=((e^(iϕ))+(e^(-iϕ)))/2, for the last term one obtains Mga((Subscript[δ, m,n+1])+(Subscript[δ, m,n-1]))/2, so that the matrix element becomes

(Subscript[H, mn])=((((ℏ^2)(m^2))/(2I))+Mga)(Subscript[δ, mn])-(Mga/2)((Subscript[δ, m,n+1])+(Subscript[δ, m,n-1])).

Below is a Mathematica code for solving the eigenvalue problem for the quantum pendulum. Among the parameters, at least one number has to be real (such as M=1.), so that Mathematica does not try to find eigenvalues analytically.

(* Definitions *)
NBasis=300;  (* Number of the lowest rotator eigenfunctions in the basis *)
M=1.;            (* Mass of the pendulum *)
II=300;          (* Moment of inertia of the pendulum *)
ℏ=1;           (* Planck's constant *)
g=1;           (* Gravity acceleration *)
a=1;            (* Lever arm of the gravity force *)
(Subscript[ω, 0])=(Sqrt[((M g a)/II)]);       (* Frequency of small oscillations *)
χmϕ[m_,ϕ_]=(1/(Sqrt[2π])) (E^(I m ϕ));     (* Eigenfunctions of the rotator *)
ϵn[m_]=((((ℏ^2)(m^2))/(2II))+M g a);
Vmn[m_,n_]=-((M g a)/2)(KroneckerDelta[m,n+1]+KroneckerDelta[m,n-1]);
Hmn[m_,n_]:=ϵn[m]KroneckerDelta[m,n]+Vmn[m,n];  (* Elements of Hamiltonian matrix *)

(* Creating Hamiltonian matrix and solving matrix eigenvalue problem *)
Timing[ HMatrix=Table[Hmn[m,n],{m,-NBasis/2,NBasis/2},{n,-NBasis/2,NBasis/2}]; ]
Print["Part of the Hamiltonian matrix:"]
TableForm[Take[HMatrix,9,9]]   (* Show part of the Hamiltonian matrix *)
Timing[ ES=Eigensystem[HMatrix]; ]
ES=Transpose[Sort[Transpose[ES]]]; (* Eigensystem sorted in accending order *)
EVal=First[ES] ;   (* Eigenvalues *)
EVecT=Last[ES];  (* This is a matrix of eigenvectors lying horizontally *)
EVec=Transpose[EVecT]; (* This is a matrix of eigenvectors standing vertically *) 
TableForm[Chop[Take[EVecT.EVec,9,9]]] ; (* Check that eigenvectors are orthonormal *)

(* Finalizing the solution *)
Eμ[μ_]:=EVal[[μ+1]];      (* Energy eigenvalues of the problem *)
c[μ_,m_]:=EVecT[[μ+1,m+1+NBasis/2]]    (* Expansion coefficients of the eigenfunctions of the problem over the rotator basis *)
ψμϕ[μ_,ϕ_]:=(∑_(m=-NBasis/2)^(NBasis/2))c[μ,m]χmϕ[m,ϕ]     (* Eigenfunctions of the problem *)

Print["Energy levels:"]
LevelsToShow=70;              LevelsToShow=Min[ LevelsToShow,NBasis];
Take[EVal,LevelsToShow]
Show[
ListPlot[Take[EVal,LevelsToShow],AxesLabel->{"μ","(Subscript[E, μ])"}],
Plot[2M g a,{x,0,LevelsToShow}]
]

One can see that the low-lying energy levels are equidistant, as the levels of the harmonic oscillator. Levels with (Subscript[E, μ])<2Mga are oscillatory. Levels above this energy correspond to rotations over the barrier. These levels form doublets because of the two posisble directions of the rotation. However, these doublets are split because of the essentially quantum-mechanical effect, the overbarrier reflection. For instance,

Eμ[50]-Eμ[49]
Eμ[52]-Eμ[51]
Eμ[54]-Eμ[53]

Splitting of level doublets, as the overbarrier reflection, decreases very fast with the energy. At high energies one approximately has

(Subscript[E, μ])=((((ℏ^2)(μ/2))^2)/(2I)),

the energy levels of the rotator. The eigenfunctions of the problem are similar to those of the harmonic oscillator at low energies. However, at high energies the probability density ((|(Subscript[ψ, μ])|)^2) is oscillating because eigenfunctions describe superpositions of rotations in different directions.

Plot[Evaluate[{((Abs[ψμϕ[0,π x]])^2),((Abs[ψμϕ[20,π x]])^2),((Abs[ψμϕ[60,π x]])^2)}],{x,-1,1},PlotRange->All,PlotPoints->30,PlotStyle->{{Black,Thick},{Black,Thin},{Red,Thin}},AxesLabel->{"ϕ/π","((\[LeftBracketingBar](Subscript[ψ, μ])[ϕ]\[RightBracketingBar])^2)"}]

To better see mixing of rotations in different directions, one can plot the coefficients (Subscript[c, μm]) for one the doublets

Eμ[54]
Eμ[53]
Eμ[54]-Eμ[53]
Plot[{c[53,Round[m]],c[54,Round[m]]},{m,-LevelsToShow,LevelsToShow},PlotRange->All,PlotStyle->{{Blue,Dashed},{Black}},AxesLabel->{"m","(Subscript[c, μm])"}]

One can see that for one of the eigenfunctions the coefficients (Subscript[c, μm]) are even in m (solid line) while for the other they are odd in m (dashed line), that is, both of them are mixtures of rotations in different directions.