Feynman’s Path Integral Formulation Actually Explained (Part 1)

Motivation

When writing a book, one must choose what to emphasize, what to deemphasize, what to cut, and what to mention in passing. Feynman’s path integral formulation didn’t seem to fit into my book in all its glory, as it is a conceptually complicated method to calculate things.

In my book, I have two chapters on the Schrödinger equation, first As an energy balance and then As a wave equation. I managed to briefly mention Feynman’s formulation in the second of those chapters. I was able to do this because Feynman’s method gives the exact same answers as Schrödinger’s equation. They both reflect the same fundamental law of the universe; both describe the time evolution of wave functions.

This may surprise some of you because it’s common for physicists to claim that “particles” go through an infinitude of paths on the way to a location. Usually, images of various chaotic paths are shown to illustrate the idea. Those same authors rarely mention that the stories they are telling refer to an intermediate calculation that is part of a larger calculation that gives you the wave function at some future time.

What I want to do is give you a fairly good idea of Feynman’s path integral formulation, but that won’t happen with a single blog post. It will take several. With this blog post, I will zoom out and give you the context of where path integrals come into the calculation of a wave function. This means that I will show you the highest level equation in Feynman’s path integral formulation (the whole method), but I won’t show you how to calculate a path integral (part of Feynman’s method). Please note the difference between those two terms.

With a bit of foresight, I predict that future blog posts will be most educational if I show you how to code the complicated path integrals. So, in this post, to prepare you for that, I will show you how I coded in MATLAB the time evolution of a one-dimensional, free electron, which appears appears in my book Physicists at Fault: Why you don’t understand quantum mechanics, yet. Thus, I assume you can read some basic code. If you cannot, MATLAB has free tutorials that can give you the basics in a few hours.

The Top-Level Equation

If I were writing this as a chapter in my book, then I would make Feynman’s path integral version of the wave function as simple as possible, at least at first. It would probably look like this the first time you saw it:
$$ \Psi = \int^\infty_{-\infty}K\cdot\Psi_0\cdot dx’$$
This is the equation I will unpack and code into MATLAB today. The \(\Psi\) is the wave function that we will find.

The Starting Wave Function

Let’s now turn to the integral on the right-hand side of the equation above. The \(\Psi_0\) has a subscript of zero, indicating that it is the starting wave function that exists at the time \(t=0\). We have the freedom to say what shape we will start with, so I chose an easy one. I pick a Gaussian function traveling to the right with some momentum.

First, I will show you a Gaussian with no momentum. It is given by the following equation and looks like a simple bump:
$$\left(\frac{1}{2\pi\cdot \sigma^2}\right)^{1/2}\cdot\exp{\frac{-(x’)^2}{4\cdot \sigma^2}}$$
The little mark beside the x is pronounced prime, so we say “eks-prime.” I include this symbol instead of a regular \(x\) that indicates the position along the \(x\)-axis because, in the calculation, it is the integration variable and will disappear in the end. The Greek lowercase sigma (\(\sigma\)) in the exponential sets the width of the Gaussian bump. The sigma in the fraction out front of the exponential helps set the height of the Gaussian.

Tacking on a complex exponential function gives the starting wave function some momentum \(p\), which shows up as some characteristic oscillations in the wave function:
$$\Psi_0=\left(\frac{1}{2\pi\cdot \sigma^2}\right)^{1/4}\cdot\exp{\frac{-(x’)^2}{4\cdot \sigma^2}}\cdot\exp{(i\cdot p\cdot x’\,/\,\hbar)}$$
This final equation is our starting wave function, and it looks like the following image.

This function is coded as follows. Note that I treat Plank’s constant (\(\hbar\)) as equal to one, so it doesn’t appear anywhere in the script.

%initial wave function
Psi0=@(xprime) exp( -xprime.^2/4/(sigma.^2) ) / ( 2*pi*sigma^2 )^(1/4) .* exp(1i*p*xprime);

Even if you don’t know how to code in MATLAB, for the most part, it is easy to see which parts of this line of code correspond to the equation. The part at the beginning, Psi0=@(xprime), makes this a named function I can call later with a variable ($ latex x’$) I must supply. I will show you soon how I pre-declare the rest of the variables in lines of code higher up in the script.

The Path Integral

The function \(K\) is commonly called a kernel or a propagator. Because of Feynman, it is also called a path integral. For each situation, one must solve a complicated integral to get a form we can work with. As I said, I plan to show you how to calculate this in later posts. Today, I give you the result of the calculation for a free electron:
$$K\;=\; \left(\frac{m}{2\pi\cdot i \cdot\hbar\cdot t}\right)^{1/2}\cdot \exp{\frac{i\cdot (x-x’)^2}{2\cdot\hbar\cdot t}}$$
Since you are reading a physics blog, I think most of you will recognize the electron’s mass \(m\). The equation gets coded into MATLAB as

%Path integral/Kernel
K=@(xprime,xx,t) sqrt(1/2/pi/1i./t) .* exp( 1i*( (xx-xprime).^2/2./t ) );

Since the last two functions combine to form the integrand of the integral, I combine them into one function:

%integrand
KPsi0=@(xprime,xx,t) K(xprime,xx,t) .* Psi0(xprime);

Declaring the Values of Variables and Whatnot

With that, the essential physics is coded in. If you want to run the code, I need to show you some code that comes before this that sets up some preliminary variables. Reading this should be straightforward, even for the uninitiated.

%%%%%%%%%%%%%%
% One-dimensional, free-particle evolution
% via Feynman Path Integral Formulation
% code written by Dr. Nathan Armstrong
%%%%%%%%%%%%%%

%Makes a number of plots equal rows x cols
plotrows=3;
plotcols=2;
numplots=plotrows*plotcols;

tl1 = tiledlayout( plotrows, plotcols, 'padding', 'tight');

%limits of plots
endza=-20;
endzb=40;
dx=0.1;
xxx=endza:dx:endzb;

%set initial time
%because of division by t in the kernel, the integral function has warnings
%for values much smaller than this
tinit=0.1;

%final time
tfinal=16;

dt=(tfinal - tinit) / (numplots-1);

%width of Gaussian peak
sigma=2;

%momentum of wavefunction
p=1;

%counter for time step display
n=1;

For-loops to Make the Plots at Various Snapshots in Time

In a moment, I will show you a for-loop that contains 23 lines of code. The purpose of the loop is to make one snapshot of a wave function for each loop. You don’t need to understand all the code inside the loop, and I won’t explain it all. Most of it is used to make the plots. I include all 23 lines so that you can run the code if you want. You can run it on MATLAB Online by creating an account. If you don’t want to paste it together from what I show you in this blog post, you can find my copy here.

The most important part of the for-loop is on line 51. This line repeatedly calculates the first equation of this post to generate a complex number for each point in space. This line is:

Psi=integral( @(xprime) KPsi0(xprime,xxx,t), -endzb*3,endzb*3,'ArrayValued',true );

For the reader who is a novice coder, this line brings up an interesting point. In that line, I call the function named integral. It requires the upper and lower bounds of integration, which in the first equation on this page are plus and minus infinity. (From my previous post, you may know I have an opinion on infinities.) We can’t numerically integrate to infinity because infinity is a concept and not a number. We must choose finite values. Because I know how spread out the function will be in the end, I know I can make it only three times the size of the plot. Inside those two bounds are the significant contributions to the integral. Should you take my code and modify the final time, Gaussian width, or momentum, the bounds of integration may need to be made larger to maintain accuracy.

The full for-loop that makes the plots is:

for t=tinit:dt:tfinal
    nexttile

    Psi=integral( @(xprime) KPsi0(xprime,xxx,t), -endzb*3,endzb*3,'ArrayValued',true );

    plot(xxx,real(Psi),'-','linewidth',2.5)
    hold on
    plot(xxx,imag(Psi),'-','linewidth',2.5)
    hold off

    if t==tinit
        ybound=1.2*max([real(Psi),imag(Psi)]);
        lgd=legend('Real \Psi','Imaginary \Psi','location','southeast','edgecolor','none');
        fontsize(lgd,12,'points')
    end

    ylim([-ybound ybound])

    set(gca,'fontsize', 19)
    set(gca,'linewidth',2)

    str = {['Time step: ' num2str(n)]};
    text(endzb-27,0.85*ybound,str,'fontsize',15)

    n=n+1;
end

You now have seen the full script that creates plots similar to what appears in my book. Running this script, you will see the following series of snapshots of a one-dimensional, free electron.

On the Physics

Now, I haven’t explained the physics of the wave functions you see in the images. Nor have I mentioned the advantages and disadvantages of Feynman’s method over that of the Schrödinger equation. You will find information on both of those in my book. You can sign up to be notified when it goes on sale by clicking the button in the side panel. Additionally, I will soon be making videos explaining the wave functions, so stay tuned for that.