Monte Carlo path integral

Quantum Mechanics
K. Konishi, G. Paffuti

Write a Monte Carlo simulation for one dimensional quantum systems. Apply the method to the computation of the propagator for harmonic and anharmonic oscillator.

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

Introduction

Euclidean path integral

The euclidean version of the path integral

Z = ∫𝒟 q Exp[- (Subscript[S, E])],

is identical to a partition function of a statistical system with β H = (Subscript[S, E]). To compute () we will use periodic boundary conditions and take a discretization of time, with a step Δt.  Our "statistical variables" will be (Subscript[q, 1])... (Subscript[q, n]), the coordinate variables at the times (Subscript[t, k])= k Δt.

The time integration is traded for a sum over discrete times:

∫dt F[t]→ (∑_k) Δt F[(Subscript[t, k])].

The discrete version of the action is (we use m = ℏ = 1)

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

In writing this expression we used periodicity.

The Monte Carlo method

The naive method

In principle we could consider (Subscript[q, i])as stochastic variables, and Exp[-S] their probability density, or, equivalently, we can compute  () by a "standard" Monte Carlo technique, the same used in the Mathematica routines: just take a random value for the (Subscript[q, i])and accept them with probability Exp[-S]. This is the same as computing the area subtended by a curve drawing randomly points in the plane and accepting only those fitting the area.

Let us make an example. Suppose we want to compute

(∫_(-1)^(+1))(1-(x^2))d x = (4/3) ≃1.33.

We draw a plot of the function and frame it with a rectangle, of area 2.

Plot[1-(x^2),{x,-1,1},Filling-> Bottom,Frame-> True]

Let us generate random points in the rectangle.

numRandomPoints = 500;
SeedRandom[1234];
xR =RandomReal[{-1,1},numRandomPoints];
yR = RandomReal[{0,1},numRandomPoints];

The points belonging to the area we are interested in are

selectedPoints=Pick[Range[numRandomPoints],Negative[ yR - (1-(xR^2))]];

Plot[1-(x^2),{x,-1,1},Filling-> Bottom,Frame-> True,Epilog-> Point/@Transpose[{xR[[selectedPoints]],yR[[selectedPoints]]}]]

The approximate value of the area is given by the fraction (selected Points)/(total Points) * (Area Rectangle)

areaRectangle=2;
area =N@ Length[selectedPoints]/numRandomPoints*areaRectangle

This method can be generalized easily to any dimension. However in this elementary form this method cannot be used to compute () : by the very nature of a system at statistical equilibrium the probability distribution of the configurations (i.e. the coordinates in our case) is sharply peaked around the equilibrium, if we just try to extract the interesting configurations starting form a uniform random distribution we have no chance, an enormous fraction of our points is just rejected in the selection mechanism.

The way to overcome this problem is to construct a Markov process which, in fact, simulate our intuitive approach of "local thermalization" for a system.

Monte Carlo procedure

We give a brief description of the method, referring to the literature for more details (see for example [1]).

Let us suppose to have a configuration C (in our case a set of values for the coordinates) and suppose we can construct a transition probability process, i.e. a function which gives the probability to find a new configuration C'once C is given: T(C,C'). If the process T satisfies the so-called  detailed balance condition (and an ergodicity condition, not relevant in our case, for which we refer to [1])

T[C,C'] Exp[-S[C]] =T[C',C] Exp[-S[C']]

then the procedure T[(Subscript[C, n]),(Subscript[C, n+1])]generate a sequence distributed with probability Exp[-S], as n→∞, i.e. for any observable A[C]

(Underscript[Lim, N→∞]) (1/N)(∑_(n=1)^N)A[(Subscript[C, n])] = ∫Exp[-C] A[C].

The main point is that we can change a single coordinate in the procedure T, in this way the action as a small change and we have a reasonable accepting rate.

More precisely the (simplest) procedure we can adopt is the following. Fix the attention on a single coordinate in S, let's say (Subscript[q, i]). Let us choose a random new value  for this coordinate, (Subscript[q, i]^′). In this process the action has changed

ΔS = S[(Subscript[q, i]^′)] - S[(Subscript[q, i])]

We accept this change with a conditional probability Exp[-S] in this way: choose a random number 0 < r < 1

If  r < Exp[-ΔS] then  (Subscript[q, new])= q'  ;   If  r > Exp[-ΔS] then  (Subscript[q, new])= q  ;

We see in particular that is S has lowered, then ΔS < 0 and Exp[-ΔS] > 1 the change is surely accepted.

It is possible to show that the algorithm () satisfy detailed balance condition.

Heat bath method

Even with this improvement the above method would have some difficulty to be used in a reasonable computer-time in Mathematica . We can slightly modify it by changing the starting point. These kind of  methods usually are known as "heat bath methods".

A given coordinate (Subscript[q, i])appears in the action in the form

(1/(Δt))(Subscript[q, i]^2) +Δt V[(Subscript[q, i])] - (Subscript[q, i])((Subscript[q, i+1])+(Subscript[q, i-1]))

The first part can be handled in an exact way, instead of generating a random q to do the acceptance test we generate a q with a probability distribution

P[q] = C Exp[-((1/(Δt))(q^2) +Δt V[q])]

As an effective action we have just to test the "force" term

- (Subscript[q, i])((Subscript[q, i+1])+(Subscript[q, i-1]))

This term is sometimes called the hopping term, for obvious reasons. Below we will show how to generate in the general case points with the distribution ().

Example

Let us make an example which will be relevant for the harmonic oscillator path integral. Let us suppose we want to compute the mean values of x and (x^2)for an action x + (x^2)/2 (it is a toy model with time reduced to a single point).

int[k_]:=Integrate[Exp[-x^2/2]Exp[-x]x^k,{x,-∞,∞}];
normalization = int[0];

Table[int[k],{k,1,2}]/normalization

The heat bath method works first by generating a number with distribution Exp[-(x^2)/2], this is provided directly by Mathematica with the function NormalDistribution[μ, σ] which generate numbers with probability density

PDF[NormalDistribution[μ,σ],x]

In our case μ = 0, σ = 1. The function updateX below translate in Mathematica the procedure described above, as input we have the "old" value of x, as output the new value, passed through the acceptance test

SeedRandom[1234];

σ = 1;
β = 1.;
updateX[x_]:= Module[{q0,qnew,deltaS,rn,force},
q0=x;
qnew=RandomReal[ NormalDistribution[0,σ]];
force = 1.;
deltaS = (qnew-q0)*force;
rn = RandomReal[{0,1}];
q0=If[Negative[rn- Exp[-β deltaS]],qnew,q0];
q0]

We make a number of "measurement", i.e. generation of the configuration, separated by a number numTherm of steps, to avoid possible statistical correlations.

numTherm=10; numMeas = 1000;

This is the starting value of q

q= RandomReal[ NormalDistribution[0,σ]];

Here we produce a set of configurations :

valQ=Table[Do[q=updateX[q],{numTherm}];q,{numMeas}];//Timing

The statistical averages are

Mean/@{valQ,(valQ^2)}

in agreement with the exact values. Clearly the agreement improves for longer updating.

Harmonic oscillator

As a first example we will compute the propagator for an harmonic oscillator. We know the exact result

D[t] = ⟨q[t]q[0]⟩= (1/(2ω))Exp[-ω |t|]

The computation is almost identical to the previous example for the heat bath method. The coordinate (Subscript[q, i]) appears in the action in the form

(Subscript[q, i]^2)((1/2)(1/(Δt))2 + (1/2)(ω^2)Δt) -  (1/(Δt))(Subscript[q, i])((Subscript[q, i+1])+(Subscript[q, i]))

The first term is generated by the function NormalDistribution, with mean μ=0 and σ  given by

(1/(σ^2)) = (1/(Δt))2 + (ω^2)Δt .

The second term is subjected to acceptance test.

It is apparent that for too small Δt the free-particle term dominate, in this case a very long and accurate simulation has to be done to distinguish the small contribution due to ω (the potential energy term). For this reason below we will take Δt relatively large.

A crucial observation : even points are subjected only to the force of odd points and viceversa, this means that we can update at once all even points and after that all odd points, this parallelism greatly reduces computational work.

The notations

The first and the last point are called below i1 and i2. As we need forward and backward shifts to compute the force we can use the function Mod to limit ourselves to numbers between 1 and npt (the total number of time points) or just add a point to the left and one to the right in a periodic way. We choose this alternative. In this way our vectors q have npt+2 components. We update the points from 2 to npt+1  (our i1 and i2) while the points q[[1]] and q[[npt+2]] are fixed by periodicity.

Parameters

We choose 200 points.

npt=200;
ω = 1.;
Δt = 1./(2 ω);
T = npt *Δt

We fix the random generator for reproducibility of the results.

SeedRandom[1234];

σ = 1/(Sqrt[(2./(Δt)) + (ω^2)Δt]) ;
β = 1./Δt;

This is our convention for points, we explicitly select even and odd points

i1=2; i2=npt+1;iF = npt+2;
even = Range[i1+1,i2,2];
odd = Range[i1,i2-1,2];

The initial configuration and the periodicity conditions are now given :

q = RandomReal[NormalDistribution[0,σ],npt+2];
{q[[1]],q[[iF]]}={ q[[i2]],q[[i1]]};

The "force" term is written as a function, this is not very efficient but it is done to allow the reader an easier analysis of the program

forceEven[q_]:=   - q[[even]]*(q[[even+1]]+ q[[even-1]]);
forceOdd[q_]:=  - q[[odd]]*(q[[odd+1]]+ q[[odd-1]]);

This is the updating procedure. We first generate even coordinates, then the odd ones. The output is the new configuration. The procedure follow the lines of the example above, the only difference being the necessity of writing the periodicity condition.

update[qInput_]:= Module[{q0,qnew,deltaS,rn},
q0 = qInput;
qnew=q0;
qnew[[even]] = RandomReal[NormalDistribution[0,σ],npt/2];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceEven[qnew]-forceEven[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[even,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
qnew=q0;
qnew[[odd]] = RandomReal[NormalDistribution[0,σ],npt/2];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceOdd[qnew]-forceOdd[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[odd,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
q0];

Given a set of configurations we compute some statistical average: q, (q^2), and q[t]q[0]. The last quantity is computed by averaging over the whole time axis points differing by t in time, in this way we use the periodicity properties to depress statistical errors.

makeCorrelations[Q_]:= Module[{num,numdist,meanQ,meanQ2,corr},
{num,numdist} = Dimensions[Q];
meanQ = Table[Mean[Q[[All,i]]],{i,1,numdist}];
meanQ2= Table[Mean[(((Q[[All,x]]-meanQ[[x]]))^2)],{x,1,numdist}];
corr=Table[
Mean[Table[Mean[(Q[[All,i]]-meanQ[[i]])*(Q[[All,i+x]]-meanQ[[i+x]])],
{i,1,npt/2}]],{x,0,npt/2-4}];
{meanQ,meanQ2,corr}]

The propagator

numTherm=10; numMeas = 2000;

This is the updating procedure

valQ=Table[Do[q=update[q],{numTherm}];q,{numMeas}];//Timing

and now we compute correlations and mean values.

cd=makeCorrelations[valQ];

We expect ⟨q⟩=0 and ⟨(q^2)⟩ = 1/2ω for an harmonic oscillator:

{Mean[cd[[1]]],2 ω Mean[cd[[2]]]}

This is a plot of the correlation (i.e. the propagator). We see an exponential damped behavior

ListPlot[cd[[3]][[1;;10]]/cd[[3]][[1]],PlotStyle-> PointSize[0.015],PlotRange-> All]

It is easy to extract the slope

data = Log[cd[[3]][[1;;10]]/cd[[3]][[1]]];
line = Fit[data, {1,x},x]

The expected value is ω Δt  ( in the exponent ω t = ω Δt n)

ω Δt

The quality of the fit is this one

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

and here is the same graph in usual linear coordinates

normψ = cd[[3]][[1;;12]][[1]];
Show[
Plot[Exp[-ω n Δt],{n,0,11},PlotRange-> All],
ListPlot[Transpose[{Range[0,11],((cd[[3]][[1;;12]])/(normψ))}],PlotStyle->{Red,PointSize[0.012]}],
PlotLabel-> "D[t]=<q[t]q[0]> : harmonic oscillator",AxesLabel-> {"t","D[t] / D[0]"}]

The anharmonic oscillator

The model, in natural units is :

H = ((p^2)/(2m)) + (1/2)(ω^2)(q^2) + (1/2)g (q^4).

The procedure is in principle equal to the previous one, but we need a substitution for NormalDistribution. In effect the coordinate (Subscript[q, i])now appears as

(Subscript[q, i]^2)((1/(Δt)) + (1/2)(ω^2)Δt) +(1/2)g Δt (Subscript[q, i]^4) -  (1/(Δt))(Subscript[q, i])((Subscript[q, i+1])+(Subscript[q, i]))

i.e. in a non gaussian form. Let us explain how an arbitrary distribution can be generated.

Generation of the distribution

In principle the procedure is the following : do a change of variables

dξ = P[x] dx

If ξ is generated with a uniform distribution then  x[ξ] is generated with density P[x], as dξ/dx = P[x]. To do this in practice we have to compute

ξ = (∫_(-∞)^x)P[τ]dτ

and invert this relation. This cannot be done analytically for a quartic polynomial in the exponent but there are not particular problems in Mathematica if we use interpolating functions.

First let us choose the parameters for the problem

npt=200;
ω = 1.;
g=1.;
Δt = 1./(2 ω);
T = npt *Δt;

Then we define the quartic polynomial and the integral to be done for the change of variables:

α2 = (1/2)(1/(Δt))2 + (1/2)(ω^2)Δt; α4=(1/2)g Δt;
norm=NIntegrate[Exp[-α2 (x^2)- α4 (x^4)],{x,-∞,+∞}];
gFun[z_]:= NIntegrate[Exp[-α2 (x^2)- α4 (x^4)],{x,-∞,z}]/norm;

We compute g[z] for a bunch of values of z. To interpolate the inverse function is just sufficient to switch x and y in the call to Interpolation

points=Table[{gFun[z],z},{z,-2,2,0.01}];
points = Union[points];

ξFun = Interpolation[points]

Warning : Interpolation needs a unique value for the result of the interpolation (otherwise we cannot define a function). As our function gFun saturates rapidly both for large and small values the interpolation interval must not be to large, otherwise numerical approximation give the same value for different z, which is a problem for the inverse function. The situation is clearer if we plot gFun[z]

Plot[gFun[z],{z,-4,4}]

Computation

From now on the computation is identical to the previous one, except the substitution of NormalDistribution with ξFun in the update procedure.

SeedRandom[1234];

β = 1./Δt;
i1=2; i2=npt+1;iF = npt+2;
q = ξFun[RandomReal[{0,1},npt+2]];
{q[[1]],q[[iF]]}={ q[[i2]],q[[i1]]};

even = Range[i1+1,i2,2];
odd = Range[i1,i2-1,2];

forceEven[q_]:=   - q[[even]]*(q[[even+1]]+ q[[even-1]]);
forceOdd[q_]:=  - q[[odd]]*(q[[odd+1]]+ q[[odd-1]]);

update[qInput_]:= Module[{q0,qnew,deltaS,rn},
q0 = qInput;
qnew=q0;
qnew[[even]] = ξFun[RandomReal[{0,1},npt/2]];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceEven[qnew]-forceEven[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[even,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
qnew=q0;
qnew[[odd]] = ξFun[RandomReal[{0,1},npt/2]];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceOdd[qnew]-forceOdd[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[odd,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
q0];

Here is the actual computation

numTherm=10; numMeas = 1000;

valQ=Table[Do[q=update[q],{numTherm}];q,{numMeas}];//Timing

cd=makeCorrelations[valQ];

We check that ⟨q⟩ ≃ 0.

{Mean[cd[[1]]], Mean[cd[[2]]]}

This is the propagator.

ListPlot[cd[[3]][[1;;10]]/cd[[3]][[1]],PlotStyle-> PointSize[0.015],PlotRange-> All]

data = Log[cd[[3]][[1;;5]]/cd[[3]][[1]]];
line = Fit[data, {1,x},x]

On general grounds the time constant for the decay must be (for large t, see Problems and text):

((Subscript[E, 1])- (Subscript[E, 0])) t = ((Subscript[E, 1])- (Subscript[E, 0])) n Δt

With our parameters (ω = 1, g = 1) we can compute easily the first few Hamiltonian eigenvalues (see chapter on variational methods):

ene={0.6961758207650811`,2.3244063521066423`,4.327524978879773`};

Then we expect for the time constant

δene =(ene[[2]]-ene[[1]])Δt

This is the quality of the fit

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

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

normψ = cd[[3]][[1;;12]][[1]];
Show[
Plot[Exp[-δene n],{n,0,11},PlotRange-> All],
ListPlot[Transpose[{Range[0,11],((cd[[3]][[1;;12]])/(normψ))}],PlotStyle->{Red,PointSize[0.012]}],
PlotLabel-> "D[t]=<q[t]q[0]> : anharmonic oscillator",AxesLabel-> {"t","D[t] / D[0]"}]

Double well

The same computation can be performed for a double well oscillator, i.e. for (ω^2) negative in the distribution. We now that for small g the first two states are almost degenerate, in this case a larger Δt simplifies the computation.

npt=200;
ω = 1.;
g=0.1;
Δt = 3./( ω);
T = npt *Δt;

Here we generate the distribution

α2 = (1/2)(1/(Δt))2 - (1/2)(ω^2)Δt; α4=(1/2)g Δt;
norm=NIntegrate[Exp[-α2 (x^2)- α4 (x^4)],{x,-∞,+∞}];
gFun[z_]:= NIntegrate[Exp[-α2 (x^2)- α4 (x^4)],{x,-∞,z}]/norm;

The shape is somewhat larger, as expected, then in the interpolation procedure we will allow a larger interval for z:

Plot[gFun[z],{z,-4,4}]

points=Table[{gFun[z],z},{z,-5,5,0.01}];
points = Union[points];

ξFun = Interpolation[points]

The computation

SeedRandom[1234];

β = 1./Δt;
i1=2; i2=npt+1;iF = npt+2;
q = ξFun[RandomReal[{0,1},npt+2]];
{q[[1]],q[[iF]]}={ q[[i2]],q[[i1]]};

even = Range[i1+1,i2,2];
odd = Range[i1,i2-1,2];

forceEven[q_]:=   - q[[even]]*(q[[even+1]]+ q[[even-1]]);
forceOdd[q_]:=  - q[[odd]]*(q[[odd+1]]+ q[[odd-1]]);

update3[qInput_]:= Module[{q0,qnew,deltaS,rn},
q0 = qInput;
qnew=q0;
qnew[[even]] = ξFun[RandomReal[{0,1},npt/2]];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceEven[qnew]-forceEven[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[even,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
qnew=q0;
qnew[[odd]] = ξFun[RandomReal[{0,1},npt/2]];
{qnew[[1]],qnew[[iF]]}={ qnew[[i2]],qnew[[i1]]};
deltaS = forceOdd[qnew]-forceOdd[q0];
rn = RandomReal[{0,1},npt/2];
changedSites = Pick[odd,Negative[rn- Exp[-β deltaS]]];
q0[[changedSites]]= qnew[[changedSites]];
{q0[[1]],q0[[iF]]}={ q0[[i2]],q0[[i1]]};
q0];

numTherm=20; numMeas = 1000;

valQ=Table[Do[q=update[q],{numTherm}];q,{numMeas}];//Timing

cd=makeCorrelations[valQ];

{Mean[cd[[1]]],2 ω Mean[cd[[2]]]}

ListPlot[cd[[3]][[1;;10]]/cd[[3]][[1]],PlotStyle-> PointSize[0.015],PlotRange-> All]

data = Log[cd[[3]][[1;;8]]/cd[[3]][[1]]];
line = Fit[data, {1,x},x]

For (ω^2)→ -1 and g = 0.1 the variational method gives for ΔE = (Subscript[E, 1])- (Subscript[E, 0])

ΔE=0.056217;

then we expect a time constant

δene =ΔE Δt

Here is the comparison with the Monte Carlo results :

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

Here is the comparison between M results and expected behavior

normψ = cd[[3]][[1;;12]][[1]];
Show[
Plot[Exp[-δene n],{n,0,11},PlotRange-> All],
ListPlot[Transpose[{Range[0,11],((cd[[3]][[1;;12]])/(normψ))}],PlotStyle->{Red,PointSize[0.012]}],
PlotLabel-> "D[t]=<q[t]q[0]> : double well",AxesLabel-> {"t","D[t] / D[0]"}]

References

1) G. Parisi, Statistical Field Theory, Addison-Wesley Publishing Company, 1988