Fourier series

Quantum Mechanics
 K. Konishi, G. Paffuti

Use a Fourier basis to study variational problems.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Introduction

General remarks

One of the most intuitive approximation for the computation of energies of bound states is to put the system in a box of length L. If L is sufficiently large the exponentially decreasing wave functions can be taken to be zero at the boundaries without appreciable error.

For wave functions confined in a box a natural set of basis functions is given by the free particle eigenstates, which are normalizable in a finite box:

(Subscript[f, n])(x) = (Sqrt[(2/L)])Sin[((n π x)/L)];    n = 1...dim

These functions have an immediate connections to Fourier series and for these series various powerful mathematical tools have been invented, in particular techniques like the Fast Fourier Transform ( FFT ) give efficient and accurate tools to manipulate such objects. In the following we will see that exactly these mathematical structures are useful in solving the Schrödinger equation.

The role of the parameter L is similar to a variational parameter, albeit of a special nature: varying L in effect change the domain of our problem. In any case our guiding principle will be always the same: look for a stability region in L for the eigenvalues.

In the box our wave functions admit the expansion

ψ = (∑_n) (Subscript[c, n]) (Subscript[f, n])[x]

and (spectral) coefficients (Subscript[c, n]) for Hamiltonian eigenstates must satisfy (we will use natural units):

(Subscript[H, nk]) (Subscript[c, k]) = E (Subscript[c, n]) ;           H = ((p^2)/2) + V[x]

This equation will be truncated to a certain order, dim. A crucial point is the possibility of computing efficiently matrix elements, to this point is devoted the next subsection.

Matrix elements

Matrix elements for the kinetic term are trivial, those for the potential can be written in term of integrals with Cos factors :

(n|(p^2)|m) = (∫_0^L)(Subscript[f', n])(Subscript[f', m])dx = (((((π n)/L)))^2)(Subscript[δ, nm]) ;
(n|V|m)=(∫_0^L)(Subscript[f, n])V(x) (Subscript[f, m])dx =  (1/L)(∫_0^L)V(x)(Cos[((π(n-m)x)/L)] - Cos[((π(n+m)x)/L)])dx .

It is convenient to introduce the notation

(Subscript[F, s])= (1/L)(∫_0^L)V(x)Cos[((π s x)/L)] ⟹  (n|V|m) = ((Subscript[F, n-m]) - (Subscript[F, n+m]))

To compute these integrals let us use the trapezoidal rule. We know from our study of the Gauss formulas that this rule is indeed for trigonometric polynomials, up to a certain degree depending on the dimension of the grid used. Our potential will not be periodic on the interval L.  The formulas are by no means exact, but in any case they are an approximation to the integral.

In general, see [*], if h = L/n is the grid spacing, n being the number of intervals, and with (Subscript[x, 0]) and (Subscript[x, n]) the extrema of the interval:

(∫_(Subscript[x, 0])^(Subscript[x, n]))f[x]𝕕x = h(((Subscript[f, 0])/2)+ (Subscript[f, 1])+ ... + (Subscript[f, n-1])+ ((Subscript[f, n])/2)) - n ((h^3)/12)f''[ξ]

ξ being an internal point to the interval. For fixed L the error decrease at least as 1/n^2, clearly if we take L too large with respect to n we go into troubles.

Let us translate the trapezoidal rule in the notation that will be used in the notebook. Take N = dim+1 intervals and a function f[x] defined at internal points plus boundary points, a total of N+1 = dim+2 points:

L = N Δx ;   (Subscript[x, r])= (r-1) Δx ;  r = 1 ... N, N+1.  ;              (Subscript[f, r]) ≡ f((Subscript[x, r]))

The trapezoidal rule, neglecting the error, is

(∫_0^L)f[x] 𝕕x ≃Δx{ (1/2)((Subscript[f, 1])+(Subscript[f, 2]))+(1/2)((Subscript[f, 2])+(Subscript[f, 3]))+... + (1/2)((Subscript[f, N])+(Subscript[f, N+1]))}= Δx{(1/2)((Subscript[f, 1])+(Subscript[f, N+1]))+(∑_(r=2)^N) (Subscript[f, r])}

Applying () to our integrals (), with

L = N Δx ;   (Subscript[x, r])= (r-1) Δx ;  r = 1 ... N, N+1.  ;              (Subscript[V, r]) ≡ V((Subscript[x, r]))

we have :

(Subscript[F, s])=((Δx)/(N Δx)) ((1/2)V(0) +  (1/2)V(L)(((-1))^s) )+(1/(N Δx))(∑_(r=2)^N)(Subscript[V, r])Cos[((π s (Subscript[x, r]))/(N Δx))]Δx=
(1/N) ((1/2)(Subscript[V, 1]) +  (1/2)(Subscript[V, N+1])(((-1))^s) )+(1/N)(∑_(r=2)^N)(Subscript[V, r])Cos[((π s (r-1))/N)]= (1/(Sqrt[2 N]))FourierDCT[V,1][[s+1]]

The last expression is just the discrete Fourier Cosine Transform (one of the version of this). In Mathematica language is the function FourierDCT[ V, 1]. There exist 4 version of this function, changing for normalization and indexing, we will use the first one. There exist an analogous  discrete Fourier Sine Transform, its definition is reported below. Let us note that we have written our formulas for n+1 points in the case of DCT and for n-1 points in the case of DST, as this is the form that we will need.

(Subscript[c, s])=(Sqrt[(2/n)])(((Subscript[a, 1])/2)+(1/2)(((-1))^(s-1))(Subscript[a, n+1])+(∑_(j=2)^n)cos( (((j-1)(s-1)π)/n))(Subscript[a, j ]))                                 (for n+1 points),   FourierDCT[ a, 1]

(Subscript[c, s])=(Sqrt[(2/n)])(∑_(j=1)^(n-1))sin(((jsπ)/n))(Subscript[a, j]).                                                                                                                 (for n-1 points)    FourierDST[ a, 1]

The DCT transform in our notation has only N + 1 components. To compute matrix elements with higher indexes in () we can use the relation, which follows from trigonometric identities:

(Subscript[F, s])= (Subscript[F, 2N-s]).

As  DCT vectors :

(Subscript[F, s]) = (Subscript[F, 2N-s])  ⇒  DCT[V][[s+1]] = DCT[V][[2N-s+1]]  ⇒  DCT[V][[k]]= DCT[V][[2N-k+2]]

The last combination can be easily written by dropping the last element of the vector ( the component n+1) and reverting the order with Reverse. An example can help. Let us take n=7; then the components are n+1 = 8. The component number 9 is given, from () by 2 n - 9 + 2 = 7. This is the result obtained with the Reverse trick:

Join[Range[8],Reverse[Drop[Range[8],-1]]]

We see that component 9 is indeed equal to component 7, the rest follows the natural descending order.

N.B.

To make some order in the indices let us note that matrix V has indices that run from 1 to n-1 = dim and describe the modes of the basis. Their relation with DCT comes from (). In particular diagonal elements are, apart from a multiplicative factor 1/(Sqrt[2n]):

V[[i,i]] = (Subscript[F, 0])- (Subscript[F, 2i])= DCT[[1]] - DCT[[2i+1]]

while non-diagonal factors are

V[[i,j]] = (Subscript[F, i-j])- (Subscript[F, i+j])= DCT[[i-j+1]] - DCT[[i+j+1]]

Indices higher than n + 1 in DCT are obtained with the Reverse trick. Indices i and j in these relations run from 1 to n-1, the resulting V matrix being  a dim×dim matrix.

All this means that to compute matrix elements we need just to call the Mathematica function FourierDCT. This is done in the function getMatV below

Functions

The function getMatV needs as input the values of the potential on the extended grid, i.e. with the extrema, it gives as output the matrix ⟨n|V|k⟩

getMatV[Vx_]:=  Module[{mat,fV,nV,fV1,dim},
dim = Length[Vx]-2;
mat = ConstantArray[0.,{dim,dim}];
nV = Length[Vx]-1;
fV = FourierDCT[Vx,1];
fV1 = Drop[fV,-1];
fV = Join[fV,Reverse[fV1]];
Table[mat[[i,j]] = fV[[i-j+1]] - fV[[i+j+1]],{i,1,dim},{j,1,i-1}];
mat = mat+Transpose[mat];
Table[mat[[i,i]] = fV[[1]] - fV[[2 i+1]],{i,1,dim}];
(1/(Sqrt[2 nV]))mat]

Hermite polynomial basis used for check

We copy here the functions of notebook [*] that will be used to perform a check on the results for anharmonic oscillator.

hermiteMatrices[dim_,parity_Integer:-1]:=
Module[{x2,p2,x4,x1}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[parity == 0 || parity == 1,Return[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

getscale[n_,g_,μ_]:= x/.FindRoot[(2n+1)(x-(μ/(x^3)))-((3g)/(x^5))(2(n^2)+2n+1),{x,(((3g))^(1/3))}]

An example and a test

Harmonic oscillator : eigenvalues

As a test let us take the usual harmonic oscillator

V[x_]:= (1/2)(x^2)

Here we fix the parameters of the grid,  Vx is the list of values of the potential on the extended grid, i.e. with the extrema included.

L = 20; dim = 200; (* All internal points *)
Npt = dim+1;
dx = L/Npt;
xptTot = -L/2+L/(Npt)Range[0,Npt]; (* N+1 points *)
Vx = V[xptTot];

The kinetic part of the Hamiltonian is trivial and we have immediately :

p2 =N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
matV = getMatV[Vx];
ham =  (1/2)p2 + matV;

We can compute eigenvalues and eigenvectors

{eigH,eigV}= Eigensystem[ham];
order=Ordering[eigH];
eigH= eigH[[order]];
eigV = eigV[[order]];

ene1=eigH[[Range[50]]]

A reasonable agreement with theoretical results is seen, given in this case by 0.5 + Range[0, N].

Eigenfunctions

The expected theoretical eigenfunctions are

psiTh[n_,x_]:= (1/(π^(1/4)))(1./(Sqrt[(2^n)(n)!]))HermiteH[n,x]Exp[-x^2/2];

To compute the values of ψ on our grid of points we have to sum the components of the eigenvector, (Subscript[c, i])with the corresponding eigenfunctions. Using Npt = dim+1, L = Npt Δx we have

ψ[(Subscript[x, k])]=(∑_(j=1)^dim)(Subscript[c, j]) (Subscript[f, j])[(Subscript[x, k])]= (Sqrt[(2/(Npt Δx))])(∑_(j=1)^(Npt-1))(Subscript[c, j]) Sin[j (π/(Npt Δx)) k Δx]

This is just the definition of FourierDST apart an overall factor 1/(Sqrt[Δx]) (see eq.() ):

ψ=FourierDST[eigV[[1]],1]/(Sqrt[dx]);

Here is the comparison between the theoretical expectation and numerical result :

Show[Plot[psiTh[0,x],{x,-5,5},PerformanceGoal-> "Speed",PlotStyle-> Red],
ListPlot[Transpose[{xptTot[[2;;dim+1]],-ψ}]],PlotRange-> {{-5,5},{0,0.8}}]

We stress that also the absolute normalization is reproduced correctly.

Test of the DCT for the potential

For a quadratic potential we can compute exactly the matrix elements in the box. Let us see what the error induced by taking FourierDCT, i.e. the trapezoidal rule, is.

For s ≠ 0 :

Integrate[(1/(2ℒ))(-ℒ/2+ξ)^2 Cos[π s (ξ/ℒ)],{ξ,0,ℒ},Assumptions-> s∈ Integers]

while for s = 0

Integrate[(1/(2ℒ))(-ℒ/2+ξ)^2 ,{ξ,0,ℒ},GenerateConditions-> False]

Defining the function

fF[0]:= 1/24;
fF[s_Integer]:=((1+(((-1))^s))/(2(π^2) (s^2)))

We can write exact matrix elements (in the box) as (see eq. () ):

matVInt=Table[(L^2)N[fF[i-j]- fF[i+j]],{i,1,dim},{j,1,dim}];

Here is the difference between some of the matrix elements. As expected the errors grow with the mode number.  This is due mathematically to

the error form of trapezoidal formula: it involves a second derivative and the second derivative of Cos[k x] grows like k^2.

Max[Abs[Flatten[N@matVInt[[1;;#1,1;;#1]]-matV[[1;;#1,1;;#1]]]]]&/@{5,10,20,50}

These are eigenvalues computed with exact matrix elements

ene2=Sort[Eigenvalues[(p2/2)+N@matVInt]][[Range[50]]]

The difference between these and the previous ones is :

Max[Abs[ene1[[Range[#1]]]-ene2[[Range[#1]]]]]&/@Range[10,50,10]

It is seen that the lower modes are practically good, while the error grows with mode number. The situation clearly goes better with a finer grid.

Anharmonic oscillator

Anharmonic oscillator

Let us study within this approach our standard target model : the anharmonic oscillator. Here are the parameters for the potential

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

and the choice of the dimension for the basis

dim = 200; (* All internal points *)
Npt = dim+1;

Matrices for L = 1

The potential is a sum of homogeneous terms,  therefore it is sufficient to compute matrices for L = 1 and multiply for the needed power of L

L=1;
dx = L/Npt;
xptTot = -L/2+L/(Npt)Range[0,Npt];
p2L1 = SparseArray@ N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
V2L1 = 1/2 (xptTot^2);
matV2L1 = SparseArray@ getMatV[1./2 (xptTot^2)];
matV4L1 =  SparseArray@ getMatV[1./2 g xptTot^4];

We first give a quick look at the L - dependence of eigenvalues

LRange= Range[4,40,2];
numEig=Floor[dim/2];

eneL=Table[
p2 =1/L^2 p2L1;
matV = (L^2) matV2L1 + (L^4)matV4L1;
ham =  (1/2)p2 + matV;
Sort[Eigenvalues[ham]][[Range[numEig]]],{L,LRange}];//Timing

ListPlot[Table[Transpose[{LRange,eneL[[All,j]]}],{j,1,numEig}],Joined-> True,
PlotRange-> {{-1,40},{0,300}},AxesOrigin-> {0,0}]

A very good stability region is found around  L ∼20.

Let us note that at large L  the eigenvalues collapse in pairs, this peculiar effect will be discussed at the end of this section.

The energies for L = 20 can be extracted from the previous data with

L20=First[Flatten[Position[LRange,22]]]

eneFFT=eneL[[L20]][[1;;50]];

Let us compare these results with the ones obtained in a more usual Hermite polynomials basis. We use functions introduced in notebook [*] to generate the data:

{p2,x2,x4,x1}=N/@hermiteMatrices[dim];

μ=1;
α = getscale[dim,g,μ];

eneVar=Sort[Eigenvalues[Normal[ (α^2)(p2/2)+ (μ/(2 (α^2)))x2+(g/(2 (α^4)))x4]]];

ListPlot[eneFFT - eneVar[[1;;50]]]

The agreement is excellent.

Comments on scale dependence

First of all we notice that the apparent degeneracy for large L in the previous calculations is surely due to discretization : for finite L the eigenvalue's problem is a standard Sturm-Liouville problem on a finite interval which do not have degenerate eigenvalues.

Check

Let us check this point with the harmonic oscillator by using exact values for matrix elements using the function fF defined above.

dim = 20; L=1;
Npt = dim+1;
dx = L/Npt;
xptTot = -L/2+L/(Npt)Range[0,Npt];
p2L1 = SparseArray@ N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
V2L1 = 1/2 (xptTot^2);

This is the FFT version for matrix elements

matV2L1 = SparseArray@getMatV[1./2 (xptTot^2)];

this the exact version

matVInt=Table[(L^2)N[fF[i-j]- fF[i+j]],{i,1,dim},{j,1,dim}];

Let us compute separately the two set of eigenvalues at different L

LRange= Range[4,50,2];
numEig=Floor[dim/2];

eneLFFT=Table[
p2 =1/L^2 p2L1;
ham =  (1/(2(L^2)))p2L1 +(L^2) matV2L1;
Sort[Eigenvalues[ham]],{L,LRange}];

eneLExact=Table[
p2 =1/L^2 p2L1;
ham =  (1/(2(L^2)))p2L1 +(L^2) matVInt;
Sort[Eigenvalues[ham]],{L,LRange}];

Here are the two respective variations with L of the eigenvalues. No collapse is present in the second case.

GraphicsRow[{ListPlot[Table[Transpose[{LRange,eneLFFT[[All,j]]}],{j,1,numEig}],Joined-> True],
ListPlot[Table[Transpose[{LRange,eneLExact[[All,j]]}],{j,1,numEig}],Joined-> True]},PlotRange-> {{-1,40},{0,50}},AxesOrigin-> {0,0},ImageSize-> 8 72]

L dependence

Let us consider the number of grid points, N, fixed. As L goes large the kinetic term, proportional for the k-th mode to k^2/(L^2)do not contribute to the eigenvalues, or better can be considered as a small perturbation with respect the potential term, so let us consider the matrix ⟨i|V|j⟩  in the case of FFT. The discrete Fourier Transform is a unitary operator, then the eigenvalues are the same as those computed in the original x - representation. In this representation V is diagonal with eigenvalues evidently given by V[(Subscript[x, k])]. For instance for the last values computed in the previous check for harmonic oscillator we see a clear correspondence:

((Max[LRange])^2)Sort[N@(V2L1[[2;;dim+1]])]

Last[eneLFFT]

The diagonal matrix (Subscript[V, k]) is clearly doubly degenerate for symmetric potential as (Subscript[x, k])and - (Subscript[x, k]) both belong to the grid.  The exception being grids with an odd number of internal points, where the single eigenvalue V[0] is non degenerate. This is why couples of eigenvalues collapse for large L. The reader can check that, for instance, with a potential (x^2) + 0.1 (1+ Sin[x]) that  this do not happen.

This considerations raise an interesting point; it is reasonable to assume that the point where approximation breaks down is determined by perturbation theory where kinetic term acts as a perturbation with respect to V.  The "energy splitting" in this approach are of the order (Subscript[V, k+1])- (Subscript[V, k]), at least for monotonic potentials, and for the k-th mode we can reasonably assume that the critical L is given by

((π^2)/2)((k^2)/(L^2))∼  (Subscript[δV, k]);

This order of magnitude can be useful to have an idea priori on critical L.

Double well

We apply the same technique to the double well version of the anharmonic oscillator. The main point is to show that quasi - degeneracy between ground state and first-excited level appear only for sufficiently large L, otherwise the levels are quite well separated. This is nothing but a manifestation of the Heisenberg principle, for small L the contribution of kinetic energy dominates and there is no effective barrier hence tunnel effect.

μ  = -1 ; g = 0.1;
V[x_]:= (1/2)μ x^2 + (1/2)g x^4
 Plot[V[x],{x,-1.5(Sqrt[((Abs[μ])/g)]),1.5(Sqrt[((Abs[μ])/g)])},AxesLabel-> {"x","V[x]"}]

dim = 100; (* All internal points *)
Npt = dim+1;

Matrices for L = 1

L=1;
dx = L/Npt;
xptTot = -L/2+L/(Npt)Range[0,Npt];
p2L1 =N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
matV2L1 = getMatV[μ/2 (xptTot^2)];
matV4L1 =  getMatV[1/2 g xptTot^4];

LRange= Range[3,20,0.25];
numEig=10;

eneL=Table[
p2 =1/L^2 p2L1;
matV = (L^2) matV2L1 + (L^4)matV4L1;
ham =  (1/2)p2 + matV;
Sort[Eigenvalues[ham]][[Range[numEig]]],{L,LRange}];

ListPlot[Table[Transpose[{LRange,eneL[[All,j]]}],{j,1,numEig}],Joined-> True,
PlotRange-> {{-1,14},{-1,2}},AxesOrigin-> {0,0},AxesLabel-> {"L","(Subscript[E, n])"}]

Hydrogen atom

The variables

We take this example to show where problems can arise. In Coulomb systems the point r=0 is a singular one. The direct method used above cannot be applied as V[r] is singular at the boundary. One can consider intervals [ϵ,L+ϵ] in the r variable, but this from introduces another free parameter, and does not cure the main problem: with a uniform grid the singular region is poorly sampled, and poor results are expected. We leave as an exercise to explore this kind of approach.

A better possibility is to make a change of variables in such a way that the singular point is sent to ∞. There are several possibilities, as in the analogous case of the Chebyshev approximation;  here we choose the one which works for the semiclassical version of the problem (and is alternative to the approach presented in the notebook on Chebyshev approximation).

The Schrödinger equation for the reduced radial function is, in natural units

-(1/2)((d^2)/(dr^2))φ + (1/2)((ℓ(ℓ+1))/(r^2))φ- (1/r) φ = E φ

With the change of variables

r = Exp[x];  φ[r] = Exp[x/2] u[x]

the equation is transformed into

-(1/2)((d^2)/(dx^2))u +( (1/2)((ℓ(ℓ+1))/(x^2))+(1/8))u- (e^x) u = E  (e^(2x))u

The singular point has been moved  to x = -∞. We will apply FFT technique to this equation. For simplicity we will limit to s-states, ℓ = 0.

The computation

Let us look at the scale dependence of the eigenvalues

dim = 100; (* All internal points *)
Npt = dim+1;

LRange= Range[15,55,5];
numEig=10;

eneL=Table[
xptTot = -L/2+L/(Npt)Range[0,Npt];
p2 =N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
matVE1 = getMatV[Exp[xptTot]];
matVE2 = getMatV[Exp[2 xptTot]];
ham =  (1/2)p2 + (1/8)IdentityMatrix[dim]-matVE1;
Sort[Eigenvalues[{ham,matVE2}]][[Range[numEig]]],{L,LRange}];

We plot  j E[[j]] + j to clearly distinguish eigenvalues and their stability range:

ListPlot[Table[Transpose[{LRange,j*eneL[[All,j]]+(j)}],{j,1,numEig}],Joined-> True,
PlotRange-> {{0,50},{0,10}},AxesOrigin-> {0,0},AxesLabel-> {"L","(Subscript[E, n])"}]

We see that L = 25 is a good value and we can perform a more accurate computation for this value of L.

L=25;
dim = 250;

Npt = dim+1;
xptTot = -L/2+L/(Npt)Range[0,Npt];
p2 =N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
matVE1 = getMatV[Exp[xptTot]];
matVE2 = getMatV[Exp[2 xptTot]];
ham =  (1/2)p2 + (1/8)IdentityMatrix[dim]-matVE1;
eneH250=Sort[Eigenvalues[{ham,matVE2}]][[Range[numEig]]];

The approximation can be extracted from the following list of values  for n^2(Subscript[E, n]):

NumberForm[eneH250[[1;;numEig]]*Range[numEig]^2,8]