# Energy of the quantum harmonic oscillator in the Monte-Carlo path integral approach

## The problem

The partition function for the quantum harmonic oscillator can be written in the path integral formulation as $$Z\propto\int Dx(\tau)\exp\left(-\frac{S_E}{\hbar}\right)=\int Dx(\tau)\exp\left(-\frac{1}{\hbar}\int_0^{\beta \hbar}d\tau\left(\frac{1}{2}m\dot x^2(\tau)+\frac{1}{2}m\omega^2 x^2(\tau)\right)\right) .$$

To attack this problem numerically, we can make this quantity discrete by choosing a sequence of $$N$$ points $$x_1$$,...,$$x_N$$ on the path $$x(\tau)$$ with distance $$a$$ between each other; rescaling these variables by an appropriate factor, we can write $$\frac{S_E^{discrete}}{\hbar}=\sum_{j=1}^N\left(\frac{1}{2\eta}(y_{j+1}-y_j)^2+\frac{\eta}{2}y_j^2\right)$$ where $$\eta= a\omega$$. The continuum is then recovered in the limit of small $$\eta$$.

In this discrete form, the internal energy of the system transforms as follows: $$\frac{U^{cont}}{\hbar\omega}=\partial_\beta\log Z=\hbar\omega\left(\frac{1}{2}+\frac{1}{e^{N\eta}-1}\right) \ \Rightarrow \ \frac{U^{discrete}}{\hbar\omega}=\frac{1}{2\eta}-\frac{1}{2\eta^2}\langle \Delta y^2\rangle+\frac{1}{2}\langle y^2\rangle$$ where the averages are first taken over a single path, and then over all paths generated in the Monte Carlo simulation according to the distribution function $$p(y_1,...,y_N)=e^{-S_E^{discrete}/\hbar}.$$

My goal is to compute the internal energy and compare it to the theoretical value for various values of $$N$$, in function of the inverse temperature $$\beta=N\eta$$.

## The algorithm

In order to sample paths that follow this distribution, we can employ a MCMC method based on the Metropolis algorithm:

1. Initialize the discrete path in a cold configuration (all zeros) or hot configuration (random numbers)

2. For each Monte Carlo sweep, loop on all $$j=1,...,N$$ and propose to modify $$y_j$$ with a uniformly sampled $$y_p$$ in the interval $$[y_j-δ,y_j+δ]$$ for some parameter $$δ$$.

3. The update is accepted with probability $$r=e^{-\Delta S_{E}^{discrete}}$$ where the difference in the action (between the original and the modified path) can be easily computed as
\begin{align*} \Delta S_{E}^{discrete}&=(y^p)^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)-\frac{y^p}{\eta}(y_{j_0+1}+y_{j_0-1})-y_{j_0}^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)\\&-\frac{1}{2}y_{j_0}(y_{j_0+1}+y_{j_0-1}). \end{align*}

## The code

I report here the code that I wrote of a MCMC simulation of the quantum harmonic oscillator in the path integral formulation (the first cell contains the needed functions, the second does the plotting):

import numpy as np
import matplotlib.pyplot as plt
from numpy.random import default_rng, Generator
import time as tm

def metropolis_ho(
path: np.ndarray,
rand: Generator,
eta: float,
delta: float,
ntimes: int,
equilibrium_start: int = 10_000,
) -> tuple[
np.ndarray,  # observables 1
np.ndarray,  # observables 2
]:
"""Monte Carlo Simulation"""
n = len(path)

# Initialize arrays of observables
obs1 = np.empty(ntimes)
obs2 = np.empty(ntimes)

# Useful constants
c1 = 1/eta
c2 = c1 + eta/2

# Iterate loop on all sites
for i in range(ntimes):
for j in range(n):
for _ in range(3):
# Set y as j-th point on path
y = path[j]

# Propose modification
y_p = rand.uniform(y - delta, y + delta)

# Calculate accept probability
force = path[(j + 1) % n] + path[j - 1]
p_ratio = c1*force*(y_p - y) + c2*(y*y - y_p*y_p)

# Accept-reject
if rand.random() < np.exp(p_ratio):
path[j] = y_p

# Average of y^2 on the path
obs1[i] = np.dot(path, path)/n

# Average of Delta y^2 on the path
diff = path - np.roll(path, -1)
obs2[i] = np.dot(diff, diff)/n

# Get rid of non-equilibrium states and decorrelate
n_term = equilibrium_start
obs1 = obs1[n_term:]
obs2 = obs2[n_term:]

return obs1, obs2

def U(obs1,obs2,eta):
"""Computes internal energy"""

c1=1/(2*eta)
c2=1/(2*eta**2)

av1=np.average(obs1)
av2=np.average(obs2)

return c1-c2*av2+av1/2

def exact_U(n,eta):
""""Theoretical expectation for internal energy"""

return 0.5+1./(np.exp(n*eta)-1)

begin=tm.time()

#Initialize arrays to plot U(n)
n_array=np.asarray([5,10,15,20,25,30,40,50])
U_array=[]

for n in n_array:

#Initialize path
hot = True
rand: Generator = default_rng(seed=0)

if hot:
start = rand.uniform(low=-1, high=1, size=n)
else:
start = np.zeros(n)

#Run MCMC
obs1, obs2 = metropolis_ho(
path=start,
rand=rand,
eta=0.01,
delta=1,
ntimes=80000,
equilibrium_start=1000,
)

#Compute energy
U_array.append(U(obs1,obs2,eta=0.01))

plt.plot(np.dot(n_array,0.01),U_array,linestyle='None',marker='s',color='k')

#Set arrays for plotting exact U
n_array=np.arange(0.01,300,0.1)
exact_U_array=[]

for n in n_array:

#Compute theoretical value
exact_U_array.append(exact_U(n,eta))

plt.plot(np.dot(eta,n_array),exact_U_array,'b')
plt.ylim(0,12)
plt.xlim(0,3)
plt.show()

end=tm.time()
print(f"Total runtime: {end-begin}.")


## Results

Even without error bars, it is pretty evident that my results are just not that accurate for the small value of $$\eta$$ that I have chosen, even though they roughly follow the theoretical curve. This makes me think my code is largely correct, but something is amiss. Can the code be improved or fixed in order to obtain a closer fit?

• i think it can be hard to judge if you don't plot the error bars. Error bar shows how confident we are about the range where true value is most likely to be. On the one hand, your method can be principally correct, and the result just needs to average on more Monte Carlo threads. On the other hand, your code or theory may contain a minor error which does not affect results much when temperature is high. It requires error bars to tell which has more credibility. Commented Oct 19, 2022 at 5:08
• @AndyChen Maybe this requires a separate question, but I'm still in the process of figuring out how to properly calculate the error bars, since $U$ is the sum of two averages obtained with Monte Carlo simulations and therefore there is a high degree of autocorrelation. I'm not sure how to deal with this in Python... Commented Oct 19, 2022 at 7:50
• Do you know how to get error bars of Monte Carlo sampling with any other model like the Ising model or integrating a function $\int{f(x)p(x)d(x)}$ where $p(x)$ is a probability distribution? The methods to get error bars for different MC simulations are conceptually the same, and you can easily extend the doing to this path integral case even if it has more complicated notations. Commented Oct 20, 2022 at 6:23
• xkcd, Label the axes! More seriously: what is your convergence criterion? Gelman-Rubin? Commented Mar 14, 2023 at 8:46

There are two problems I see with your algorithm in the $$\eta\rightarrow0$$ limit that do not make it wrong but are probably the cause of large statistical errors:
1. Your third equation counts on cancellation of two terms that go to infinity, which is a rule of thumb'' sign that your procedure might not be numerically stable (more careful analysis would reveal statistical error that grows to infinity with $$\eta\rightarrow0$$). The best way to solve the problem is by using centroid virial estimator published here: https://doi.org/10.1063/1.2013257