Authors: Amélie Orban (Editor), Inga Schöyen, Alessandro Delmonte, Alexander Rihane.
Based on the project Simulating Single Photon Interference using QuTiP (Quantum Toolbox in Python) carried out in June 2021 at Maastricht University, Maastricht Science Programme.
Last updated: 09/11/21.
import matplotlib.pyplot as plt
import qutip
from IPython.display import Image
from numpy import cos, exp, pi, random, sin, sqrt
from qutip import Qobj, basis
%matplotlib inline
Wave-particle duality is one of the key concepts of quantum mechanics and it illustrates the fact that every quantum particle or object can be described using both particle and wave properties. Neither of these classical concepts can fully describe quantum objects on its own.
In a similar manner to Thomas Young's double-slit experiment, originally carried out in 1801, the wave-particle duality of light can be observed by performing a single-photon interference experiment, or in this case, simulation. The individual particle-like photons interfere with themselves, which is an intrinsically wave-like property, thus exhibiting both particle and wave characteristics at the same time.
This lecture investigates the phenomenon of single-photon interference by creating a simulation of a quantum optical experiment in which single photons go through a polarisation interferometer.
This section presents the setup of the single-photon interference experiment from a theoretical perspective, and details the different optical elements that must be simulated. The mathematical expression of each element will also be presented, in the form of Jones matrices acting as operators on the given quantum state.
The setup and corresponding theory of this single photon interference experiment closely follows the work of Mark Beck in "Quantum Mechanics, Theory and Experiment" (Beck, M. (2012). Quantum mechanics, theory and experiment. Oxford University Press.), more specifically Experiment 6 of Section 3.
A wave plate is used to modify the polarization of a wave, using the fact that its effective index of refraction depends on the polarization of the incident wave.
In this case, a half-wave plate is used, for which the relative phase shift between the fast and slow axes corresponds to a half wavelength shift.
The mathematical expression of a half-wave plate at $45^\circ$ is the following:
\begin{equation} J_{\lambda/2 \hspace{1mm} 45^\circ} = \left[\begin{array}{cc} \cos(2\cdot45) & \sin(2\cdot45) \\ \sin(2\cdot45) & -\cos(2\cdot45) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] \end{equation}A polarization analyzer, also called beam displacing polarizer or polarizing beam splitter, has the property that it can both split a beam into two orthogonal components, and displace each of these components differently.
If the polarization analyzer is positioned perpendicular to the horizontal plane, in such a way that waves propagating horizontally are normally incident, the general waves get split into vertical and horizontal components and it is referred to as a PA$_{HV}$.
If one were to rotate a PA$_{HV}$ by $45^\circ$, a general wave incident on such a polarization analyzer would be analyzed into $|+45\rangle$ and $|-45\rangle$ polarization components, making it a PA$_{45}$.
The setup of the experiment, represented in the figure below, is as follows.
This apparatus acts as an interferometer, and the final super-positioned states should interfere producing an interference pattern and output intensities that depend on the relative phase between them (equivalently, on the path length difference).
Image(filename="images/single-photon-interference-setup.jpg", width=700, embed=True)
Note that a single photon is sent through these optical elements, making the process quantum mechanical in nature. This requires for an important distinction, when it comes to the splitting of the state into orthogonal components.
Code the polarization basis vectors, for the bases $H/V$ and $+/-45$. They can then be used to code any polarization state, since any state can be expressed as a linear combination of basis vectors.
# HV basis
# horizontal polarization
H = basis(2, 0)
# vertical polarization
V = basis(2, 1)
# +45/-45 basis (in terms of HV basis)
# +45 polarization
p45 = 1 / sqrt(2) * (H + V)
# -45 polarization
n45 = 1 / sqrt(2) * (H - V)
It induces a phase shift on the H polarization relative to V polarization, because of the difference in the path lengths of the two arms of the interferometer.
In the case of a single photon, the state doesn't actually get split and there is only a single state throughout the apparatus, so the PA$_{HV}$ only affects the phase of the components.
# Polarization analyzer (HV) n°1
phaseshift1 = pi / 4 # CONSTANT
# should depend on real size of the setup (here: arbitrarily chosen)
PA_HV1 = Qobj([[exp(1j * phaseshift1), 0], [0, 1]])
PA_HV1
In this case, the fast axis makes an angle of $45^\circ$ with the horizontal.
For a single photon, the effect of going through such a half-wave plate is that the $|H\rangle$ and $|V\rangle$ polarization components are switched.
# Half-wave plate
θ = pi / 4 # fast axis orientation
# (!) numpy calculates with rad
halfwave = Qobj([[cos(2 * θ), sin(2 * θ)], [sin(2 * θ), -cos(2 * θ)]])
"""
removes very small elements
(numerical artifacts from the finite precision of the computer)
"""
halfwave.tidyup()
As before, it induces a phase shift on the H polarization relative to V polarization.
# Polarization analyzer (HV) n°2
phaseshift2 = pi / 4 # CHANGE TO CHANGE INTERFERENCE
# should depend on real size of the setup
PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]])
PA_HV2
The PA$_{45}$ acts on the photon such that it has to come out as a $|+45\rangle$ polarized photon or as a $|-45\rangle$ polarized photon.
# Polarization analyzer (45)
# linear Polarizer, transmission axis +45 wrt horizontal
θ = pi / 4
Pp45 = Qobj([[cos(θ) ** 2, cos(θ) * sin(θ)], [cos(θ) * sin(θ), sin(θ) ** 2]])
# linear Polarizer, transmission axis -45 wrt horizontal
θ = -pi / 4
Pn45 = Qobj([[cos(θ) ** 2, cos(θ) * sin(θ)], [cos(θ) * sin(θ), sin(θ) ** 2]])
def PA_45(vector):
p45_comp = Pp45 * vector # retrieve only +45 component
n45_comp = Pn45 * vector # retrieve only -45 component
return p45_comp, n45_comp
Define initial variables:
psi_0 = p45 # define the initial state (+45 vector)
phaseshift2_init = pi / 4 # initial value
phaseshift2_max = 8 * pi
n = 100 # resolution of φ (amount of steps)
step = (
phaseshift2_max - phaseshift2_init
) / n # interval divided by number of small steps we want
N_init = 1000 # number of iterations (range(N) -> 0 to N-1, both included)
# create x and y coords arrays to store the values needed to plot output graph
x_coords = [] # relative phase shift
y1_coords = [] # amount of photons in +45
y2_coords = [] # amount of photons in -45
Create a for loop to run through all of the possible values for the phase shift considered for PA$_{HV2}$:
In this loop, the passage of the photon through the interferometer is simulated using an effective matrix. The PA$_{45}$ is used after that to compute the probability related to each of these polarization states, and eventually to get the results of the experiment.
Inside the first for loop is another loop performing what can be considered as the final measurement. Whether the photon comes out in a $|+45\rangle$ or $|-45\rangle$ polarization state is a random process that depends on the corresponding probabilities.
for i in range(n + 1):
output_p45 = 0
output_n45 = 0
phaseshift2 = phaseshift2_init + i * step
x_coords.append(
(phaseshift2 - phaseshift1) / pi
) # add realtive phase shift to x coords
# create corresponding PA_HV2
PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]])
EffM = PA_HV2 * halfwave * PA_HV1 # define the effective matrix
# apply the effective matrix to the initial state to get the final state
psi_final = EffM * psi_0
psi_p45 = PA_45(psi_final)[0] # retrieve +45 and -45 components
psi_n45 = PA_45(psi_final)[1]
# probab is rounded up to 5 decimals to avoid machine precision artifacts
proba_p45 = round(psi_p45.norm() ** 2, 5)
proba_n45 = round(psi_n45.norm() ** 2, 5)
for j in range(N_init):
"""
generates random number between 1 and 100 (both included),
100 because 100% of the photons need to come out in either
+45 or -45 state
"""
a = random.randint(1, 100)
if a <= proba_p45 * 100:
output_p45 = output_p45 + 1
else:
output_n45 = output_n45 + 1
y1_coords.append(output_p45)
y2_coords.append(output_n45)
Create output plot:
Using the arrays (which are collections of coordinates) that are created as the simulation goes on for different values and setups, a plot can be created to visualize the dependency of the interference on the relative phase shift, or similarly, on the different in path lengths of the two arms of the interferometer.
plt.plot(x_coords, y1_coords, "b.", markersize=9, label="Photons in state +45")
plt.plot(x_coords, y2_coords, "r.", markersize=9, label="Photons in state -45")
legend = plt.legend(loc="upper center", fontsize="x-large")
plt.ylim([-100, N_init + 500])
plt.xlabel("Relative phase shift (multiples of π)")
plt.ylabel("Amount of photons detected")
plt.title("Amount of photons existing in |+/-45> states");
The data shows a clear interference pattern. The measurements oscillate between only measuring photons in the state $|+45\rangle$ or $|-45\rangle$ as the relative phase shift increases.
Alternative cases can be investigated to get a better understanding of what is happening, why this interference pattern appears, and of the quantum nature of this experiment.
In this case, it is shown that, when the superpositon of states collapses (here the superposition refers to the photon taking both paths "at the same time"), the interference pattern disappears. This is done by blocking either the vertical or the horizontal output port of the very first polarization analyzer, thus preventing the photon to be in the previously mentionned superposition of states as it can take only one path. This comes with some modifications, with respect to how the simulation was implemented previously:
# Polarization analyzer (HV) n°1, with V output port BLOCKED
phaseshift1 = pi / 4 # CONSTANT
# should depend on real size of the setup (here: arbitrarily chosen)
PA_HV1vb = Qobj([[exp(1j * phaseshift1), 0], [0, 0]])
psi_0 = p45 # Defining the initial state (+45 vector)
phaseshift2_init = pi / 4 # initial value
phaseshift2_max = 8 * pi
n = 100 # resolution of φ (amount of steps)
step = (
phaseshift2_max - phaseshift2_init
) / n # interval divided by number of small steps we want
# number of iterations (range(N) -> 0 to N-1, both included)
N_init = 1000
x_coords = [] # create x- and y- coords. arrays
# (x = phase shift of 2nd PA_HV,
# y1 = amount of photons in +45,
# y2 = amount of photons in -45)
y1_coords = []
y2_coords = []
for i in range(n + 1):
output_p45 = 0
output_n45 = 0
phaseshift2 = phaseshift2_init + i * step
# add realtive phase shift to x coords
x_coords.append((phaseshift2 - phaseshift1) / pi)
# create corresponding PA_HV2
PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]])
EffM = PA_HV2 * halfwave * PA_HV1vb # Defining the effective matrix
"""
Applying the effective matrix to the initial state to get the final state
"""
psi_final = EffM * psi_0
psi_p45 = PA_45(psi_final)[0] # Determining the probabilities
psi_n45 = PA_45(psi_final)[1]
proba_p45 = round(psi_p45.norm() ** 2, 5)
proba_n45 = round(psi_n45.norm() ** 2, 5)
if (proba_p45 + proba_n45) == 1:
N = N_init # all of the photons get to the end
else:
"""
half of the photons are blocked (should only get N/2 in the ouput)
-> total prob should be 0.5
"""
N = int(N_init / 2)
for j in range(N):
"""
generates random number between 1 and 50 (both included),
50 because 50% of the photons need to come out in either
+45 or -45 state (since other 50% was blocked)
"""
a = random.randint(1, 50)
if a <= proba_p45 * 100:
output_p45 = output_p45 + 1
else:
output_n45 = output_n45 + 1
y1_coords.append(output_p45)
y2_coords.append(output_n45)
plt.plot(x_coords, y1_coords, "b.", markersize=9, label="Photons in state +45")
plt.plot(x_coords, y2_coords, "r.", markersize=9, label="Photons in state -45")
legend = plt.legend(loc="upper center", fontsize="x-large")
plt.ylim([0, 1000])
plt.xlabel("Relative phase shift (multiples of π)")
plt.ylabel("Amount of photons detected")
plt.title("Amount of photons existing in |+/-45> states (V port blocked)");
# Polarization analyzer (HV) n°1, with H output port BLOCKED
PA_HV1hb = Qobj([[0, 0], [0, 1]])
psi_0 = p45 # Defining the initial state (+45 vector)
phaseshift2_init = pi / 4 # initial value
phaseshift2_max = 8 * pi
n = 100 # resolution of φ (amount of steps)
step = (
phaseshift2_max - phaseshift2_init
) / n # interval divided by number of small steps we want
# number of iterations (range(N) -> 0 to N-1, both included)
N_init = 1000
"""
create x- and y- coords. arrays (x = phase shift of 2nd PA_HV,
y1 = amount of photons in +45, y2 = amount of photons in -45)
"""
x_coords = []
y1_coords = []
y2_coords = []
for i in range(n + 1):
output_p45 = 0
output_n45 = 0
phaseshift2 = phaseshift2_init + i * step
# add realtive phase shift to x coords
x_coords.append((phaseshift2 - phaseshift1) / pi)
# create corresponding PA_HV2
PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]])
# Defining the effective matrix
EffM = PA_HV2 * halfwave * PA_HV1vb
"""
Applying the effective matrix to the initial state to get the final state
"""
psi_final = EffM * psi_0
psi_p45 = PA_45(psi_final)[0] # Determining the probabilities
psi_n45 = PA_45(psi_final)[1]
proba_p45 = round(psi_p45.norm() ** 2, 5)
proba_n45 = round(psi_n45.norm() ** 2, 5)
if (proba_p45 + proba_n45) == 1:
N = N_init # all of the photons get to the end
else:
"""
half of the photons are blocked (should only get N/2 in the ouput)
-> total probability should be 0.5
"""
N = int(N_init / 2)
"""
# generates random number between 1 and 50 (both included),
50 because 50% of the photons need to come out in either
+45 or -45 state (since other 50% was blocked)
"""
for j in range(N):
a = random.randint(1, 50)
if a <= proba_p45 * 100:
output_p45 = output_p45 + 1
else:
output_n45 = output_n45 + 1
y1_coords.append(output_p45)
y2_coords.append(output_n45)
plt.plot(x_coords, y1_coords, "b.", markersize=9, label="Photons in state +45")
plt.plot(x_coords, y2_coords, "r.", markersize=9, label="Photons in state -45")
legend = plt.legend(loc="upper center", fontsize="x-large")
plt.ylim([0, 1000])
plt.xlabel("Relative phase shift (multiples of π)")
plt.ylabel("Amount of photons detected")
plt.title("Amount of photons existing in |+/-45> states (H port blocked)");
The results show that, when both pathways are unobstructed, an interference pattern emerges. As the relative phase varies, the detected photons oscillate between the $|+45\rangle$ and $|-45\rangle$ states, clearly displaying interference. Note that, while inside the interferometer, the beams are in orthogonal polarizations, so no interference should occur there, though the components get affected and phase-shifted differently. Interference "happens", or rather becomes visible, at the PA$_{45}$ where both the horizontal and vertical polarisation are projected onto $+45^\circ$ and $-45^\circ$ axes. The degree of interference oscillates with increasing phase difference between the two paths.
When one of the output port from the first polarization analyzer is blocked, the amount of photons detected is constant (apart from some random fluctuations) under a changing relative phase shift and there is no interference pattern. The results are identical for both setups, with the V beam and the H beam blocked. This illustrates that when the superposition collapses, no interference is displayed.
Important takeaways:
qutip.about()
QuTiP: Quantum Toolbox in Python ================================ Copyright (c) QuTiP team 2011 and later. Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li, Jake Lishman, Simon Cross and Asier Galicia. Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng. Original developers: R. J. Johansson & P. D. Nation. Previous lead developers: Chris Granade & A. Grimsmo. Currently developed through wide collaboration. See https://github.com/qutip for details. QuTiP Version: 4.7.6 Numpy Version: 1.22.4 Scipy Version: 1.12.0 Cython Version: 0.29.37 Matplotlib Version: 3.5.2 Python Version: 3.10.4 Number of CPUs: 4 BLAS Info: Generic OPENMP Installed: False INTEL MKL Ext: False Platform Info: Linux (x86_64) Installation path: /usr/share/miniconda3/envs/test-environment/lib/python3.10/site-packages/qutip ================================================================================ Please cite QuTiP in your publication. ================================================================================ For your convenience a bibtex reference can be easily generated using `qutip.cite()`