Path integral and transfer matrix

Quantum Mechanics
K. Konishi, G. Paffuti

Compute the path integral defining a transfer matrix between adiacent time slices.

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

Introduction

The (euclidean) path integral can be computed using the technique of trasfer matrix, well known from statistical mechanics.

We start by a discretization of time, in the form (Subscript[t, k])= k Δt. At each time we call (Subscript[q, k]) the coordinate. The path integral, apart from a multiplicative constant, is obtained integrating over intermediate coordinates the exponential of the (euclidean) action, which in discretized form reads

(Subscript[S, E]) = (∑_k) Δt (1/2)((((((Subscript[q, k])- (Subscript[q, k+1]))/(Δt))))^2) +  (∑_k) Δt V[(Subscript[q, k])]

If we consider a Schrodinger picture for the intermediate states it is immediate that the path integral can be written as a matrix element of a power of T defined by

T = 𝒱 𝒦 𝒱

where

𝒱 = Exp[-(1/2)Δt V[Q]];  ⟨q' |𝒦|q⟩ = Exp[-(1/(2 Δt))(q'- q)]

The part depending on the potential energy has been split between the two adiacent times for simmetry. Q is the positionoperator,  with Q|q⟩ = q|q⟩.

The path integral in now written as

Z[q',q] = C⟨q'| (T^N)|q⟩;   T = N Δt.

The numerical implementation of the above procedure can be explicitly done by a discretization in q-space, with a grid of points separated by Δx.

In the following we will focus on a particular case of  (), the periodic case, with q'= q.

In this case as shown in the text, for large T the path integral gives the wave function of the ground state

Z[q,q]→ Exp[-(Subscript[E, 0])T] ((|(Subscript[ψ, 0])(q)|)^2)

The kinetic term

The kinetic term, the matrix 𝒦 above, can be easily written in Mathematica using Outer[] (see Help). An example will be self-explanatory

X = Array[x,3]; Y = Array[y,3];

Outer[Exp[-(((#1 - #2))^2)]&,X,Y]//MatrixForm

Examples

Let us define a grid of points, between -L/2 and L/2. Q are coordinates for these points.

L = 6;
Δx = 0.05;
Q = Range[-L/2,L/2,Δx];

We define a time - step Δt.

Δt = 0.5;

Given a potential V[x] the definition of matrix T is easy.

Harmonic oscillator

Let us consider for example an harmonic oscillator with frequency ω:

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

we have

matV=DiagonalMatrix[Exp[((-1)/2)V[Q]Δt]];
matK = Outer[(((Δx)/(Sqrt[2π Δt]))) Exp[((-1)/(2Δt))(((#1 - #2))^2)]&,Q,Q];
matT = matV.matK.matV;

We choose a number of time steps, here 200 and write the (periodic) path integral, denoted zF below

n = 200;
zF=MatrixPower[matT,n];//Timing

The Feynman theorem allows to extract immediately (ψ^2)for the ground state

ψ2 = Diagonal[zF];

Here is a plot of the function compared with the known analytical result :

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

Show[ListPlot[Transpose[{Q,ψ2/Max[ψ2]}],PlotStyle-> Red],
Plot[psiHarmonic[0,x]^2/((psiHarmonic[0,0])^2),{x,-L/2,L/2}]]

The propagator

Other quantities can be computed in a similar way. Let us consider for instance the propagator ⟨q[t]q[0]⟩.

Normalizing its value to ⟨(q^2)⟩:

matQ = DiagonalMatrix[Q];
norm𝒟 = Tr[matQ.matQ.MatrixPower[matT,n]];
𝒟[0]=1.;

we have

𝒟[t_]:=Tr[matQ.MatrixPower[matT,t].matQ.MatrixPower[matT,n-t]]/norm𝒟

Let us compute its value for a list of times

times = Range[0,10];

valD=Table[𝒟[t],{t,times}]

ListPlot[Transpose[{times,valD}],PlotRange-> All]

We know that the general form of large t behavior is

𝒟[t] ∼ Exp[-((Subscript[E, 1])- (Subscript[E, 0])) t] = Exp[((Subscript[E, 1])- (Subscript[E, 0])) n Δt];

and in particular for the harmonic oscillator the exact result is

𝒟[t] = (1/(2ω))Exp[-ω |t|]

Here is the fit of the exponential slope

data = Log[valD];
line = Fit[data, {1,x},x]

expectedSlope = ω Δt

and this is the result compared with the expectations in usual linear coordinates

Show[
Plot[Exp[-ω n Δt],{n,Min[times],Max[times]},PlotRange-> All],
ListPlot[Transpose[{times,valD}],PlotStyle->{Red,PointSize[0.012]}],
PlotLabel-> "D[t]=<q[t]q[0]> : harmonic oscillator",AxesLabel-> {"t","D[t] / D[0]"}]

Double well

The computation canbe performed for any potential. Let us consider in particular the case of a double well potential

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

We compute again the T matrix.

L = 8;
Δx = 0.05;
Q = Range[-L/2,L/2,Δx];
Δt = 0.8;

matV=DiagonalMatrix[Exp[((-1)/2)V[Q]Δt]];
matK = Outer[(((Δx)/(Sqrt[2π Δt]))) Exp[((-1)/(2Δt))(((#1 - #2))^2)]&,Q,Q];
matT = matV.matK.matV;

We choose a number of time steps, here 200 and write the (periodic) path integral, denoted zF below

n = 350;
zF=MatrixPower[matT,n];

We see the clearl symmetric nature of the state:

ψ2 = Diagonal[zF];
ListPlot[Transpose[{Q,ψ2/Max[ψ2]}],Joined-> True]

The propagator

We repeat the computation of the propagator

matQ = DiagonalMatrix[Q];
norm𝒟 = Tr[matQ.matQ.MatrixPower[matT,n]];
𝒟[0]=1.;

𝒟[t_]:=Tr[matQ.MatrixPower[matT,t].matQ.MatrixPower[matT,n-t]]/norm𝒟
times = Range[0,40];
valD=Table[𝒟[t],{t,times}];

ListPlot[Transpose[{times,valD}],PlotRange-> All]

With this parameter (ω = 1, g = 0.1) we will compute in the chapter on variational methods for ΔE = (Subscript[E, 1])- (Subscript[E, 0])

ΔE=0.056217;
expectedSlope = ΔE Δt

Here is the fit of the exponential slope

data = Transpose[{times,Log[valD]}];
line = Fit[data, {1,x},x]

The reader has to check that T ΔE ≫ 1, in order that the leading contribution to zF  is effectively selected

n Δt ΔE

This is the quality of the fit

Show[ListPlot[data, PlotStyle->Red], Plot[{line}, {x, 0, Max[times]}]]

and this is the result compared with the expectations in usual linear coordinates

Show[
Plot[Exp[-ΔE n Δt],{n,Min[times],Max[times]},PlotRange-> All],
ListPlot[Transpose[{times,valD}],PlotStyle->{Red,PointSize[0.01]}],
PlotLabel-> "D[t]=<q[t]q[0]> : harmonic oscillator",AxesLabel-> {"t","D[t] / D[0]"}]

Higher precision requires more points and smaller Δt.

We leave to the reader as an exercise the study of the usual anharmonic oscillator.

Comment on the procedure

It is quite clear that all the information is contained in the matrix T, the evolution for large times just selects its highest eigenvalue. In practical terms we can consider T as a numerical approximation to the infinitesimal euclidean evolution operator.

To check this let us compute the eigenvalues and eigenstates of T:

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

As the eigenvalues are Exp[- (Subscript[E, n])Δt] the energy gap between the ground state and the first excited state is (compared with the result of the fit)

{Log[((eigH[[2]])/(eigH[[1]]))],line[[2]]/x}

The same is true for the eigenstate, in particular here we show that Feynman theorem select the ground state :

ψ2EigV1 = ((eigV[[1]])^2);

Show[ListPlot[{Transpose[{Q,ψ2EigV1/Max[ψ2EigV1]}],Transpose[{Q,ψ2/Max[ψ2]}]},
Joined-> {False,True},PlotStyle-> {Red,Black}]]