Boson Bunching
A prominent example of boson bunching is the Hong-Ou-Mandel effect, where the bunching of two photons arises from a destructive quantum interference between the trajectories where they both either cross a beam splitter or are reflected. This effect takes its roots in the indistinguishability of identical photons. You can read more about it [1]. More generally one can send in \(n\) photons into an optical interferometer via \(m\) modes and one then measures in which modes the photons are ending up. For the example discussed here there are 7 photons sent into 7 modes and we are interested in the probability that all 7 photons end up in the first two modes. The interferometer is described by an unitary matrix and the incoming photons can be fully indistinguishable, fully distinguishable or partially distinguishable.
Introduction
The paper [2] shows that, unlike in the Hong-Ou-Mandel effect, bunching gets maximized with partially distinguishable photons. This disproofs the following conjecture:
Conjecture 1
(Generalized Bunching). Consider any input state of classically correlated photons. For any linear interferometer \(\hat{U}\) and any nontrivial subset \(\mathcal{K}\) of output modes, the probability that all photons are found in \(\mathcal{K}\) is maximal if the photons are (perfectly) indistinguishable.
Evidence for this conjecture was given [3].
The probability for such a case is given by:
where \(n\) is the number of input photons, \(S\) the distinguishability matrix and \(H\) a matrix which depends on \(\hat{U}\) and \(\mathcal{K}\). The symbol \(\odot\) is the elementwise multiplication.
The distinguishability matrix \(S\) is defined by:
where \(<i|\) is the photon in the \(i^{th}\) input mode \(S_{ij} = 1\) if \(\text{photon}_i\) and \(\text{photon}_j\) are indistinguishable and \(0\) if \(\text{photon}_i\) and \(\text{photon}_j\) are fully distinguishable. So for fully indistinguishable photons, \(S\) is the matrix with ones in all entries.
The matrix \(H\) is defined by:
with \(i,j \in \{1, ... , n\}\) where \(n\) is the number of input photons. If we now compare the probability for fully indistinguishable photons to partially distinguishable photons we can rewrite conjecture 1 like this:
This version of the conjecture has been proven wrong with a 7 by 7 counterexample by [4].
In this notebook we are constructing an optical circuit based on this counter example and input the partially distinguishable photons to show that this indeed amplifies the bunching probability.
Setting up the required packages.
[1]:
# Perceval
import perceval as pcvl
# Other
import numpy as np
import math
import matplotlib.pyplot as plt
Constructing the unitary matrix
The counter example to the conjecture is \(H = A\) and \(S^T = A^T\) with
and \(\omega = exp(\frac{2*\pi *i}{5})\), the \(5^{th}\) root of unitiy.
Thus we have to construct the unitary circuit fom the definition of \(H\). There are infinitly many such unitaries and we therefore generate a random one which returns \(H=A\).
[2]:
# number of input photons and modes (works for n ≥ 4) See [2] for the generalisation of the above to n photons.
n = 7
# definition of w as the (n-2)th root of unity
w = np.exp((2*math.pi*1j)/(n-2))
def create_vectors(n: int, w: complex) -> tuple:
# first two orthonormal vectors from which the Unitary matrix for the circuit is built
v1 = np.zeros(n, dtype='complex_')
v1[0] = 1
v1[1] = 0
for i in range(2, n):
v1[i] = 1/math.sqrt(2)
v2 = np.zeros(n, dtype='complex_')
v2[0] = 0
v2[1] = 1
for i in range(2, n):
v2[i] = np.conj(w**(i-2))/math.sqrt(2)
return v1, v2 # these are the complex conjugate of the first two rows of the matrix M for n=7
def normalize(vector: np.ndarray) -> np.ndarray:
return vector / np.linalg.norm(vector)
def make_unitary(n: int, w: complex) -> np.matrix:
vector1, vector2 = create_vectors(n, w)
orthonormal_basis = [normalize(vector1), normalize(vector2)]
for _ in range(n-2):
# Create a random n-dimensional complex vector
random_vector = np.random.rand(n) + 1j*np.random.rand(n)
# Subtract projections onto existing basis vectors to make it orthogonal
for basis_vector in orthonormal_basis:
random_vector -= np.vdot(basis_vector,
random_vector) * basis_vector
# Normalize the orthogonal vector to make it orthonormal
orthonormal_basis.append(normalize(random_vector))
# Define the unitary matrix from the calculated vectors
mat = []
for i in range(n):
mat.append(orthonormal_basis[i])
unitary_matrix = np.matrix(mat)
return unitary_matrix
unitary_matrix = make_unitary(n, w)
Building the Circuit
Now that we have the unitary matrix we can use perceval to create a circuit from it. We do this by decomposing the unitary into Mach-Zender Interferometers. This gives us a 7 modes interferometer.
[3]:
# building Mach-Zender Interferometer block of the circuit
mzi = pcvl.catalog["mzi phase last"].build_circuit(
phi_a=pcvl.Parameter("φ_a"), phi_b=pcvl.Parameter("φ_b"))
# convert Unitary matrix into perceval language
unitary = pcvl.Matrix(unitary_matrix)
# create circuit
circuit_rand = pcvl.Circuit.decomposition(unitary, mzi,
phase_shifter_fn=pcvl.PS,
shape=pcvl.InterferometerShape.TRIANGLE)
pcvl.pdisplay(circuit_rand, recursive=True)
[3]:
Inputting photons of different degrees of distinguishability
Now that we have the circuit which represents the correct unitary matrix and therefore the correct \(H\) from the conjecture, we need to feed in the right photons following the distinguishably matrix \(S=A\). In the [2] the input photons are described via polarization. In this simulation it is not specified how the photons are distinguishable. The first two photons are fully distinguishable. The other five photons are each a linear combination of the first two photons and form an equally spaced star shape when drawn on the Bloch sphere.
Here one can see the two fully distinguishable photons with red arrows and the five linear combinations in the equator plane.
[4]:
def create_inputs(n: int, w: complex) -> tuple:
# make a list of fully distinguishable photons of the form |{a:i}> that means we give the photon an attribute a with value i
states_dist = []
for i in range(1, n+1):
x = "|{{a:{}}}>".format(i)
states_dist.append(pcvl.StateVector(x))
# make a list of all input photons which are superpositions of |{a:1}> and |{a:2}>
# start with the pure states |{a:1}> and |{a:2}>
states_par_dist = [states_dist[0], states_dist[1]]
for i in range(2, n):
x = states_par_dist[0] + states_par_dist[1] * \
w**(i-2) # add the states |{a:1}> + w^(i-2)|{a:2}>
states_par_dist.append(x)
# Input states
# initialise the variables
indistinguishable_photons = []
partially_distinguishable_photons = 1
distinguishable_photons = 1
# fill the states
for i in range(n):
indistinguishable_photons.append(1) # gives the state |1,1, ... , 1>
partially_distinguishable_photons = partially_distinguishable_photons * \
states_par_dist[i]
distinguishable_photons = distinguishable_photons * \
states_dist[i] # gives the state |{a:1},{a:2}, ... ,{a:n}>
return pcvl.BasicState(indistinguishable_photons), partially_distinguishable_photons, distinguishable_photons
# indistinguishable photons, partially distinguishable photons, fully distinguishable photons respectively
indistinguishable_photons, partially_distinguishable_photons, distinguishable_photons = create_inputs(
n, w)
To see if photon bunching gets amplified for partially distiguishable photons we need to compare it to the case of indistinguishable photons and for interest we are also comparing it to fully distinguishable photons. We are expecting that the probability for partially distinguishable photons will be highest followed by the indistinguishable photons and fully distinguishable photons having the lowest probability.
[5]:
# simulating boson sampling
processor = pcvl.Processor("SLOS")
processor.set_circuit(circuit_rand)
simulator = pcvl.SimulatorFactory().build(processor)
[6]:
def calc_prob(distinguishable_type, n: int, simulator) -> tuple:
def reverse(lst):
new_lst = lst[::-1]
return new_lst
# initialize the probabilities
probability_photons = 0
probability_distribution_photons = []
# summing over all cases where all photons end up in only two modes
# range over half the distribution because of symmetry
for i in range(math.ceil((n+1)/2)):
probability = simulator.probability(
distinguishable_type, pcvl.BasicState([i, n-i]+[0] * (n-2)))
probability_distribution_photons.append(probability)
probability_photons += probability
X = []
for i in range(math.ceil((n+1)/2), n+1):
X.append(probability_distribution_photons[i-math.ceil((n+1)/2)])
probability_photons += probability_distribution_photons[i-math.ceil(
(n+1)/2)]
probability_distribution_photons = probability_distribution_photons + \
reverse(X)
for i in range(n+1):
if probability_photons != 0:
probability_distribution_photons[i] = probability_distribution_photons[i] / \
probability_photons
return probability_distribution_photons, probability_photons
# Save the Probabilities so we don't have to run the CalcProb function multiple times
probability_distribution_indistinguishable, probability_indistinguishable = calc_prob(
indistinguishable_photons, n, simulator)
probability_distribution_distinguishable, probability_distinguishable = calc_prob(
distinguishable_photons, n, simulator) # high runtime
probability_distribution_partially_distinguishable, probability_partially_distinguishable = calc_prob(
partially_distinguishable_photons, n, simulator) # high runtime
print(
f"The probability for all {n} indistinguishable photons ending up in the first two modes is: {probability_indistinguishable*100:.3f} %")
print(
f"The probability for all {n} distinguishable photons ending up in the first two modes is: {probability_distinguishable*100:.3f} %")
print(
f"The probability for all {n} partially distinguishable photons ending up in the first two modes is: {probability_partially_distinguishable*100:.3f} %")
The probability for all 7 indistinguishable photons ending up in the first two modes is: 0.699 %
The probability for all 7 distinguishable photons ending up in the first two modes is: 0.016 %
The probability for all 7 partially distinguishable photons ending up in the first two modes is: 0.751 %
From the results above we can clearly see that the implemented partial distinguishablity amplified the bunching probability from 0.699% to 0,751%. This however does not proof that our set up gives the maximum bunching probability. As it only disproofs the conjecture with a counter example.
Plotting the results
The first plot shows two-mode bunching probabilities for the three different scenarios. Here we can see the results match our expectations.
The second plot shows the photon-number probability distribution for the first two output modes. The distributions are normalized, i.e., the probabilities are conditioned on events where all 7 photons end up in the first two modes (two-mode bunching events). This plot is interesting because it shows that the shape of the distribution of indistinguishable photons is very different to the shape of the distributions of partially and fully distinguishable photons.
[7]:
def plot(n: int, probability_distribution_indistinguishable, probability_indistinguishable, probability_distribution_partially_distinguishable, probability_partially_distinguishable, probability_distribution_distinguishable, probability_distinguishable):
X = []
for i in range(n+1):
# creates a list of all possible tuple such that a+b=n for tuple (a,b)
X.append("({},{})".format(i, n-i))
output_modes = X
colours = ['#58c4e1', '#946cba', '#383e48'] # light blue, purple and gray
def add_labels(x, y):
for i in range(len(x)):
if y[i] > 0.0005:
plt.text(i, y[i]/2, round(y[i], 5), ha='center', color='white')
else:
plt.text(i, 1.25*y[i], round(y[i], 5), ha='center')
cases = ["indistinguishable",
"partially distinguishable", "distinguishable"]
probabilities = [probability_indistinguishable,
probability_partially_distinguishable, probability_distinguishable]
plt.bar(cases, probabilities, color=colours)
plt.ylabel('probability')
plt.title("Comparing the total bunching probabilities for {} photons".format(n))
add_labels(cases, probabilities)
plt.show()
plt.scatter(output_modes, probability_distribution_indistinguishable,
label='indistinguishable', color=colours[0])
plt.scatter(output_modes, probability_distribution_partially_distinguishable,
label='partially distinguishable', color=colours[1])
plt.scatter(output_modes, probability_distribution_distinguishable,
label='distinguishable', color=colours[2])
plt.ylabel('probability')
plt.xlabel('number of photons in first two modes')
plt.title("Normalized bunching distributions")
plt.legend()
plt.show()
plot(n, probability_distribution_indistinguishable, probability_indistinguishable, probability_distribution_partially_distinguishable,
probability_partially_distinguishable, probability_distribution_distinguishable, probability_distinguishable)
To see the that this set up only works for at least 7 photons we can compare the cases of 4,5,6 and 7 photons. You can see this if you run the above for n = 4,5 and 6. Running it for 8 or higher values for n will take a long time.
The probability for all 8 indistinguishable photon ending up in the first two modes is: 0.240 %
The probability for all 8 distinguishable photon ending up in the first two modes is: 0.002 %
The probability for all 8 partially distinguishable photon ending up in the first two modes is: 0.294 %
References
[1] Agata M. Branczyk, Hong-Ou-Mandel Interference (2017) in arXiv:1711.00080 (http://arxiv.org/abs/1711.00080)
[2] Benoit Seron, Leonardo Novo and Nicolas J. Cerf, Boson bunching is not maximized by indistinguishable particles (2022) in arXiv:2203.01306 (https://arxiv.org/abs/2203.01306)
[3] V. S. Shchesnovich, Universality of Generalized Bunching and Efficient Assessment of Boson Sampling (2016) in arXiv:1509.01561 (https://arxiv.org/abs/1509.01561)
[4] Stephen W. Drury, A counterexample to a question of Bapat and Sunder (2016) in Electronic Journal of Linear Algebra Vol 31 69-70 (https://journals.uwyo.edu/index.php/ela/article/view/1631/1631)