Simulation 103: Monte Carlo Modeling Quantum Mechanics
Handling the probabilistic nature of Quantum Mechanics computationally
Quantum Mechanics is seen as an especially ominous field of physics as it describes reality at a level we don’t perceive on a day to day basis. Not only that, the conclusions Quantum Mechanics gives us are unintuitive and seemingly bizarre. In this article, we will not only demystify some of the underlying facts of Quantum Mechanics but also develop a computational framework to simulate a keystone quantum experiment.
In this article we will:
- Learn a basic overview of Quantum Mechanics
- Go over the theory behind the double slit experiment
- Learn about Monte Carlo methods
- Create a Monte Carlo simulation of the double slit experiment
Quantum Mechanics Overview
Quantum Mechanics is a branch of physics that describes the behavior of matter and energy at the smallest scale, such as atoms and subatomic particles. It provides a framework for understanding the strange and counterintuitive phenomena that occur at such a fundamental level. These phenomena can not be described using classical physics which is what gives Quantum Mechanics its air of mysticism.
For instance, Quantum Mechanics introduces the concept of wave-particle duality, where particles can exhibit both wave-like and particle-like behavior depending on if they are observed. If a particle is not being observed, it does not exist as a discrete object; it exists as a region of probability where it could exist until it is observed. Once observed, it then becomes discrete somewhere in that region. This is an incredibly foreign way to think about existence since we only experience physics at a classical level but there is plenty of experimental evidence proving the quantum understanding of the world is true.
We will go over the most commonly know experiment in Quantum Mechanics and explore the theory behind it.
The Double Slit Experiment
The double slit experiment is by far the most well known experiment in Quantum Mechanics. It experimentally proved wave-particle duality by sending a stream of particles through a pair of slits in a sheet of metal. The particles could go through one slit or the other and land on a secondary sheet (viewing screen) where their final position could be observed.
On the first run of the experiment, a detector was placed to observe the particles before going into the slits. This way, we would know which slit each of the particles went through. The corresponding pattern of particle hits on the secondary sheet of metal lined up with the slits as expected. On the second run of the experiment, the detector was removed. This time, the particle hits did not line up with the slits; a banding pattern was observed. See figure 2 for a visual of the double slit experiment.
So, why did this happen? Why did the particles land where we expected them to when observed but not when unobserved? Here is where wave-particle duality comes into play. When the particles were observed, they were discreate objects so it was like throwing balls through the slits; they would either go through one or the other. When the particles were unobserved they were not discreate objects; they were waves. Now, the particles are like ripples in water going through the slits, seen in figure 3.
When the particles act as waves, they create a banding pattern on the second metal sheet because of wave interference. The particles, now waves, get split into two waves as they pass through the slits. The two waves can then interact with each other and amplify each other when they line up (peak to peak) or cancel each other out when they are misaligned (peak to trough), seen in figure 4.
The peaks of the combined wave become become the regions of highest probability for the particle to exist once it hits the second metal sheet. Conversely, the particle can not exist where the waves have canceled out. This is what creates the banding pattern on the second metal sheet. The bands line up with the peaks in the combined wave and the gaps line up with the areas of cancelation.
We will circle back to some of the theory behind the double slit experiment later on in this article, but we understand enough Quantum Mechanics to know that the behavior of quantum systems is probabilistic. In order to simulate probabilistic systems, we will need Monte Carlo methods; which will be explored in the next section.
Monte Carlo Methods
Monte Carlo methods, named after the Casino de Monte-Carlo in Monaco, are a set of computational methods that use randomness and selectivity to get results. Monte Carlo methods are used when a deterministic answer does not exist or is unfeasible to obtain. These methods are remarkably simple as they only require the rolling of random numbers and then checking the roll against some criteria, but they are incredibly powerful.
We will go over a use case of Monte Carlo methods before using them to simulate a quantum system.
Estimating Pi
We can see the power of Monte Carlo methods with a simple estimation of pi. Here we will run a simulation where we throw darts at a board (place random points). The board will be a circle inside of a square (seen in figure 6) and we will keep track of darts that land inside of the circle as well as the total amount of darts throw.
That’s it, that’s the simulation and everything we need to estimate the value of pi. We can do this by looking at the relationship between the area of our circle and our square on the board seen below in equation 1.
As we can see, the value of pi is equal to four times the ratio between the area of the circle and square. This ratio between areas can be approximated by the ratio of darts that land in the circle vs the number that land in the square. In our set up, the square is the total area of our board so we can look at the number of darts that land in the circle vs the total number of darts. We run the simulation in figure 7 below.
In our simulation, as the number of darts increase, the accuracy of our pi estimation increases in kind. This is just a simple example of how randomness and selectivity can be used to get real results. The code snippet to run this simulation is given below.
import random
def monteCarloPi(numSamples):
dartsInCircle = 0
#Throw random darts
for i in range(numSamples):
x = random.uniform(0, 1)
y = random.uniform(0, 1)
#Check if darts land in circle
if (x**2 + y**2) <= 1:
dartsInCircle += 1
estimate = 4 * dartsInCircle / (numSamples * 1.0)
return estimate
# Example usage:
print(monteCarloPi(10000))Now that we understand what Monte Carlo methods are and how to use them, we will build a simulation of the double slit experiment.
Monte Carlo Quantum Simulation
We will now use the power of Monte Carlo methods to handle the probabilistic nature of Quantum Mechanics. We will take the last two sections of this article and combine them to make a Monte Carlo simulation of the double slit experiment. We first need to define the probability function for a particle existing somewhere on the viewing screen.
Fraunhofer Diffraction Equation
For the sake of brevity, we will not go over the full derivation of the Fraunhofer diffraction equation in this article. We only need to accept it as our probability function in our simulation. In short, the equation takes into account the double slit layout as well as some properties of the particle to find a probability distribution of where our given particle could exist on the second metal sheet. An example of this function is given in figure 8.
Simulation
We first need to write the Fraunhofer Diffraction Equation into code so we can calculate the probability of particle being at a given location.
#Parameters used
#m = Particle mass in kg
#v = Particle velocity in m/s
#lamb = h/(m*v) de broglie wavelength, scaled to be used in this modeled, h = 6.62 × 10^−34 Js.
#a = Slit width in meters
#d = Distance between slits in meters
#l = Distance to screen in meters
def Diffraction(a,d,lamb,l,x):
#Calculating some constants
c = (np.pi*d)/(lamb*l)
k = (np.pi*a)/(lamb*l)
#Probability function
probability = ((np.cos(c))**2)*(((np.sin(k*x))/(k*x))**2)
return probabilityNext we write a function to perform a Monte Carlo roll against our probability function. We first select a random location on the second metal sheet viewing screen. Then we find the probability of the particle existing at that location by evaluating our diffraction equation at that location. Finally, we roll a random number against our probability and if our probability at that location is higher than the random number we accept that location and say a particle exists there.
This Monte Carlo rolling is essentially a sampling from our probability function to estimate the behavior of the probabilistic system. Each roll represents a possible particle going through the slits and manifesting at some location on the viewing screen.
def MonteCarlo(a,d,lamb,l):
#Monte carlo throw for a particle at a given position on our viewing screen
x = np.random.uniform(-.001,.001) #between -1mm and 1 mm
p = Diffraction(a,d,lamb,l,x)
rand = np.random.uniform(0,1)
if (p >= rand):
return xWith our probability function and Monte Carlo roll we can now run our simulation by continuously rolling and recording our results. We are using a feasible set of experimental parameters but they can be changed (may need to change the position roll in the above function if we do).
#Parameters of our system
m= 9.1*10**-31 #Mass of electron
v= 5*10**5 #velocity
lamb = 6.62*10**-34/(m*v) #de broglie wavelength.
a = 150*10**-6 #Slit width
d = 600*10**-6 #Distance between slits
l = 10 #Distance to screen
#The loop that runs the monte carlo to generate the interferance pattern
numParticles = 10**6
particles = []
for i in range(numParticles):
#Rolling a random y value as well to give our bands vertical spread
y = np.random.uniform(0,1)
x = MonteCarlo(a,d,lamb,l)
if x != None:
particles.append([x,y])After running our simulation, we can run the following code snippet to visualize the results seen in figure 9 and compare then to a real interference pattern in figure 10.
#Visualize of the interferance pattern
particles = np.array(particles)
plt.figure(figsize = (20,6))
plt.scatter(particles[:,0],particles[:,1])
plt.title("Double Slit Interferance Pattern")
plt.xlabel("X position")
plt.ylabel("Y position")
plt.show()
plt.figure(figsize = (12,6))
plt.tight_layout()
plt.show()
That looks like a realistic interference pattern! The intensity or clustering of points in our simulation matches what we would observe if we were to actually run the double slit experiment. As another check to our results, we can check the distribution of our points against the theoretical probability. We should expect our simulated distribution to follow the theory since we knew the theoretical distribution beforehand and used it.
It is important to note that in broader cases outside of this article we may not have a probability distribution to work with which is what makes Monte Carlo methods so powerful. Even with an unknown or partially know probability distribution Monte Carlo methods can be used to simulate and sample because they only need a single criterion to roll against. This criterion could come from real data or on the spot calculation when mapping out the entire probability distribution is infeasible.
#Showing that the generated model fits the theoretical model.
plt.figure(figsize = (20,6))
weights = np.ones_like(particles[:,0])/float(len(particles[:,0]))
plt.hist(particles[:,0], bins = 75, weights=weights)#density = True)
x = np.linspace(-.001,.001,500)
plt.plot(x,Diffraction(a,d,lamb,l,x)) #max(p) is there to normalize the function.
plt.title("Theoretical Model Over the Generated Distribution",fontsize = 16)
plt.ylabel("Probability",fontsize = 16)
plt.xlabel("X position",fontsize = 16)
plt.xticks(fontsize = 12)
plt.yticks(fontsize = 12)
plt.show()
Conclusion
In this article, we went over some basic theory of quantum mechanics and the double slit experiment. Additionally, we covered Monte Carlo methods and how to use them to simulate a quantum system. Though our simulation was basic, it shows off all of the tools required to dive deeper into quantum simulation; which we will do in future articles.
Full Code
import numpy as np
import matplotlib.pyplot as plt
#Parameters used
#m = Particle mass in kg
#v = Particle velocity in m/s
#lamb = h/(m*v) "de broglie wavelength", scaled to be used in this modeled, h = 6.62 × 10^−34 Js.
#a = Slit width in meters
#d = Distance between slits in meters
#l = Distance to screen in meters
def Diffraction(a,d,lamb,l,x):
#Calculating some constants
c = (np.pi*d)/(lamb*l)
k = (np.pi*a)/(lamb*l)
#Probability function
probability = ((np.cos(c))**2)*(((np.sin(k*x))/(k*x))**2)
return probability
def MonteCarlo(a,d,lamb,l):
#Monte carlo throw for a particle at a given position on our viewing screen
x = np.random.uniform(-.001,.001)
p = Diffraction(a,d,lamb,l,x)
rand = np.random.uniform(0,1)
if (p >= rand):
return x
#Parameters of our system
m= 9.1*10**-31 #Mass of electron
v= 5*10**5 #velocity
lamb = 6.62*10**-34/(m*v) #de broglie wavelength, scaled to be used in this modeled, h = 1.
a = 150*10**-6 #Slit width
d = 600*10**-6 #Distance between slits
l = 10 #Distance to screen
#The loop that runs the monte carlo to generate the interferance pattern
numParticles = 10**6
particles = []
for i in range(numParticles):
#Rolling a random y value as well to give our bands vertical spread
y = np.random.uniform(0,1)
x = MonteCarlo(a,d,lamb,l)
if x != None:
particles.append([x,y])
#Visualize of the interferance pattern
particles = np.array(particles)
plt.figure(figsize = (20,6))
plt.scatter(particles[:,0],particles[:,1])
plt.title("Double Slit Interferance Pattern", fontsize = 16)
plt.xlabel("X position", fontsize = 16)
plt.ylabel("Y position", fontsize = 16)
plt.xticks(fontsize = 12)
plt.yticks(fontsize = 12)
#plt.show()
plt.figure(figsize = (12,6))
plt.tight_layout()
plt.show()
#Showing that the generated model fits the theoretical model.
plt.figure(figsize = (20,6))
weights = np.ones_like(particles[:,0])/float(len(particles[:,0]))
plt.hist(particles[:,0], bins = 75, weights=weights)#density = True)
x = np.linspace(-.001,.001,500)
plt.plot(x,Diffraction(a,d,lamb,l,x)) #max(p) is there to normalize the function.
plt.title("Theoretical Model Over the Generated Distribution",fontsize = 16)
plt.ylabel("Probability",fontsize = 16)
plt.xlabel("X position",fontsize = 16)
plt.xticks(fontsize = 12)
plt.yticks(fontsize = 12)
plt.show()References
Young’s Double Slit Demonstration Video: https://www.youtube.com/watch?v=nuaHY5lj2AA
Any figures not cited where created by the author
