In [ ]:
circuit_custom_unitary.draw(output='mpl')
Out[ ]:
Links The Jupyter notebook can be found at: https://github.com/robhendrik/HOMsimulation The blog post related to this notebook can be found at http://armchairquantumphysicist.com/2023/05/22/simulating-the-hom-effect-on-quantum-computer/
In this notebook we simulate the HOM effect on a quantum computer. The HOM effect is an optical quantum effect where two photons enter a beamsplitter and somehow interact such that they end up together in one of the output ports. Classical physics would predict that there is an equal likelihood of finding the photons together in one output, or finding them separately in different outputs. The quantum effect is what is actually observed by experiment (see references for details).
In this notebook we build a quantum algorithm that correctly predicts the experimentally observationn and run this algorithm on a quantum computer. There is no real quantum advantage (the effect can easily be simulated by a regular computer). The exercise is for learning and for fun.
Contents
References
# Imported libraries
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.lines
import numpy as np
import math
from IPython.display import Image # for displaying picture
from qutip import * # for quantum simulations
import pandas as pd # used for easily printing tables
# Settings for this notebook
length_of_fock_state = 3 # max number of photons in a channel is 2, so possible states are |0>, |1> or |2>
no_of_channels = 2 # the beamsplitter has two input and two output channels, so we work with two channels that are''mixed' by the beamsplitter
precision =2 # used for printing output)
pd.set_option("display.precision", precision)
In this notebook we use a very simple optical component: A beamsplitter. This beamsplitter has two input ports and two output ports. For single photons we find the expected, and rather boring outcome that the photon either passes through the beamsplitter, or is reflected by the beamsplitter with a probability depending on the reflectance or transmittance of the beamsplitter. For a 50%/50% beamsplitter you would find an input photon with equal likelihood in any of the output ports of the beamsplitter.
However, when using two photons it gets a lot more interesting. We see that when the photons enter into different ports and are fundamentally 'indistinguishable' an effect occurs leading to the two photons always being together in any of the outputs. This 'bunching' is known at the HOM effect, named after Hong, Ou and Mandel (see their paper "Measurement of subpicosecond time intervals between two photons by interference" which mentioned in the references for this notebook)
To model the beamsplitter we define the matrix below. Here the transmission coefficient is $t$ and the reflection coefficient is $r$. These coefficients act on the ampitudes of the electric field. The probability for photons to be transmitted or reflected is the square of the absolute value of these amplitudes. So transmission probability is $|t|^{2}$ and reflection probability is $|r|^{2}$. If for instance $t = r = \frac{1}{\sqrt{2}}$ the probability for a single photon to be reflected is $r^{2} = \frac{1}{2}$ and the probability to be transmitted is $t^{2} = \frac{1}{2}$. So this would be a 50%/50% beamsplitter equally splitting the light between the output ports. The matrix describing the general beamsplitter is
$\begin{pmatrix} r & t \\ t & -r \\ \end{pmatrix}$
For the 50%/50% beamsplitter the matrix looks like
$\begin{pmatrix} \frac{1}{2} \sqrt{2} & \frac{1}{2} \sqrt{2} \\ \frac{1}{2} \sqrt{2} & -\frac{1}{2} \sqrt{2} \\ \end{pmatrix} = \frac{1}{2} \sqrt{2} \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ \end{pmatrix}$
For more information on how to describe a beamsplitter with a matrix check out https://en.wikipedia.org/wiki/Beam_splitter or https://www.pas.rochester.edu/~howell/mysite2/Tutorials/Beamsplitter2.pdf.
To describe the light we use Fock states (see https://en.wikipedia.org/wiki/Fock_state). In a Fock state there is a well defined amount or photons in each channel. So if we look at the two input ports of the beamsplitter we could have state $|00>$ with zero photons in either input, or we could have state $|21>$ with 2 photons in one port and one photon in the other port. The HOM effects is such that if we have input state $|11>$ (one photon in each input) we measure output states $|02>$ and $|20>$ with 50% probability each. We never measure $|11>$ at the output.
For the beamsplitter we have two input channels and two output channels. When using Fock states the beamsplitter 'transforms' input state $|nm>$ into output state $|uv>$. Here $n$ and $m$ are the number of photons in input channels 1 and 2, and $u$ and $v$ are the number of photons in output channels 1 and 2. We can make a transition matrix describing how this transformation happens. This matrix will have for instance a coeffient at indices '10' and '01' which indicates the amplitude for the transition from state $|10>$ at the input to state $|01>$ at the output. If this is amplitude is for instance $\frac{1}{\sqrt{2}}$ the probability to end up in state $|01>$ would be $(\frac{1}{\sqrt{2}})^{2}$ or $\frac{1}{2}$.
As we deal with lossless beamsplitter the total amount of photons in the input and output have to be the (this is the law of conservation of energy, we cannot create photons out of nothing and without loss they can also not disappear). So an input state with one photon ($|10>$ or $|01>$) can only transform to an output state with one photon ($|10>$ or $|01>$). The total probability should be 100%. So if we know that the probability to go from input $|10>$ to output $|01>$ is $\frac{1}{2}$ we automatically know that the probability to from input $|10>$ to output $|10>$ is also $\frac{1}{2}$. We also know that the probability to go from input $|10>$ to for instance output $|10>$, $|11>$ or $|20>$ is zero as the total number of photons is not kept the same.
The transition matrix contains amplitudes and the transition probability is the square of the absolute value of the amplitude. The probability is a real number ranging from 0 to 1. The amplitude can be complex number. We need to calculate with amplitudes to be able to include quantum interference where a probability is reduced due to two amplitudes cancelling eachother for the different ways to get to an component. As explained in https://armchairquantumphysicist.com/2023/04/19/the-hom-effect-explained/ this is the basis for the HOM effect. If we would ignore the amplitude and only calculate with probabilities we would not be able to describe the interesting quantum effects.
Below we define the Python function to calcute the Fock state transition matrix. Here $r$ and $t$ are reflection and transmission amplitudes. For the 50%/50% beamsplitter we take $t = r = \frac{1}{\sqrt{2}}$. The name for this function in our Python code is:
create_Fock_coefficients(r, t, make_unitary)
If the paramater make_unitary is set to False we use the transition coefficients defined by the physical system. This leads to a practical issue when modelling. To keep the size of the model finite we have to set a maximum to the number of photons per channel. We do this by setting the parameter length_of_fock_state. If set this to for instance 3 the maximum photon count per channel is 2 (we have 3 options per channel, either 0,1 or 2 photons).
If we limit the photon count per channel to for instance 2, we can have input state $|22>$. This means the output state could be $|31>$ or $|40>$, with photon count per channel that is higher than the earlier set value of 2. But if we extend the basis to include this and allow photon count of 4 per channel the input could be $|44>$ and the output should include $|80>$. We can of course have a different basis as input and output, but that is not practical. Also the system would not be 'revesible' which poses an issue later if we want to model on a quantum computer. So for practical purposes we can set the option make_unitary to True. In that case the coefficients will be articially set if the total photon count at either input or output is above the maximum defined. So for a maximum photon count per channel set to 2 and make_unitary set to True the transition from $|21>$ to $|21>$ will be artificially set to 1 and the transition from $|21>$ to $|03>$ will be zero (so we can leave state $|03>$ out of consideration). This is not the real physics. A real beamsplitter will have a probability to go from $|21>$ to $|03>$. Purely from practical purposes we exclude this. For total photon counts up to 2 the matrix will describe the real system, for higher counts the matrix is artificial and does not describe the real system with make_unitary set to True.
Below the two transition matrices for the physical system with make_unitary = False. In some rows and columnesn(for $|21>$, $|12>$ and $|22>$ ) the squares of the amplitudes do not add up to 1.
r = t = 1/math.sqrt(2)
make_unitary = False
Fock_coefficients = create_Fock_coefficients(r, t, make_unitary)
00 01 02 10 11 12 20 21 22
00 1.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.0
01 0.0 0.71 0.00 0.71 0.00 0.00 0.00 0.00 0.0
02 0.0 0.00 0.50 0.00 0.71 0.00 0.50 0.00 0.0
10 0.0 -0.71 0.00 0.71 0.00 0.00 0.00 0.00 0.0
11 0.0 0.00 -0.71 0.00 0.00 0.00 0.71 0.00 0.0
12 0.0 0.00 0.00 0.00 0.00 -0.35 0.00 0.35 0.0
20 0.0 0.00 0.50 0.00 -0.71 0.00 0.50 0.00 0.0
21 0.0 0.00 0.00 0.00 0.00 -0.35 0.00 -0.35 0.0
22 0.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.5
When we model with with make_unitary = True we artificially set the transition for states with in total more than two photons to 1, and the it becomes a 'nice' transition matrix.
r = t = 1/math.sqrt(2)
make_unitary = True
Fock_coefficients = create_Fock_coefficients(r, t, make_unitary)
00 01 02 10 11 12 20 21 22
00 1.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0 0.0
01 0.0 0.71 0.00 0.71 0.00 0.0 0.00 0.0 0.0
02 0.0 0.00 0.50 0.00 0.71 0.0 0.50 0.0 0.0
10 0.0 -0.71 0.00 0.71 0.00 0.0 0.00 0.0 0.0
11 0.0 0.00 -0.71 0.00 0.00 0.0 0.71 0.0 0.0
12 0.0 0.00 0.00 0.00 0.00 1.0 0.00 0.0 0.0
20 0.0 0.00 0.50 0.00 -0.71 0.0 0.50 0.0 0.0
21 0.0 0.00 0.00 0.00 0.00 0.0 0.00 1.0 0.0
22 0.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0 1.0
An alternative approach would be to limit the set of input and output states only to a specific set (i.e., only look at states with up to 2 photons and exclude the rest). However, when modelling with a register of (qu)bits it is convenient to have a number of states that is a power of 2. So we can easily deal with a basis of 2, 4, 8 or 16 states, but have to take specific care when dealing with for instance a basis consisting of 6 states. Right now we have nine states, so the next step is to condense to 8 states by removing one of the states with total number of photons larger than 2 (for those states we anyway have set the transition amplityde articificially to 1, so there is no physical meaning to it).
def create_Fock_coefficients(r, t, make_unitary):
Fock_coefficients = dict([])
for n in range(length_of_fock_state):
for m in range(length_of_fock_state):
# n and m are the photon counts for the input state |nm>.
# The for loops run through all possible input states |nm> with photon count n and m less than 'length_of_fock_state'
input_state = str(n) + str(m)
coeff = [[0 for v in range(2*length_of_fock_state)] for u in range(2*length_of_fock_state)] # initialize at zero
for i in range(0,n+1):
for j in range(0,m+1):
# u and v are the photon counts in the output state |uv>.
# The loops are such that we run through all possible values of u and v such that n+m = u+v
# (so total photon count at input same as at output)
# Note that the same value for u and v occur multiple times in the for loop, more ways to creat one output state!
v = (n+m) - (i+j)
u = (i+j)
# the factor from the binomial expansion
data = math.comb(n,i)*math.comb(m,j)*(r**(n-i))*(t**(m-j))*(t**(i))*((-r)**(j))
# factor in sqrt n for boson raising operation c|n> = sqrt(n+1)|n+1> and l|n> = sqrt(n)|n-1>. <n|a*a|n> = n
data = data * math.sqrt(math.factorial(u)*math.factorial(v)/(math.factorial(n)*math.factorial(m)))
# total number of photons cannot exceed (length-1) for each of the output channels
coeff[u][v] = coeff[u][v] + data
# if make_unitary is set to True we have to set all coefficients for output states with
if make_unitary:
if n + m >= length_of_fock_state or n + m >= length_of_fock_state:
if (n == u) and (m == v):
coeff[u][v] = 1.0
else:
coeff[u][v] = 0.0
Fock_coefficients[input_state] = {str(u) + str(v): coeff[u][v] for v in range(length_of_fock_state) for u in range(length_of_fock_state)}
return Fock_coefficients
r = t = 1/math.sqrt(2)
make_unitary = True
Fock_coefficients = create_Fock_coefficients(r, t, make_unitary)
print("Transition matrix between Fock states in 2 channels, with total photon count of maximum 2")
print(pd.DataFrame([[Fock_coefficients[input_state][output_state] for input_state in Fock_coefficients.keys()] for output_state in Fock_coefficients.keys()], columns = list(Fock_coefficients.keys()), index=list(Fock_coefficients.keys())))
Transition matrix between Fock states in 2 channels, with total photon count of maximum 2
00 01 02 10 11 12 20 21 22
00 1.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0 0.0
01 0.0 0.71 0.00 0.71 0.00 0.0 0.00 0.0 0.0
02 0.0 0.00 0.50 0.00 0.71 0.0 0.50 0.0 0.0
10 0.0 -0.71 0.00 0.71 0.00 0.0 0.00 0.0 0.0
11 0.0 0.00 -0.71 0.00 0.00 0.0 0.71 0.0 0.0
12 0.0 0.00 0.00 0.00 0.00 1.0 0.00 0.0 0.0
20 0.0 0.00 0.50 0.00 -0.71 0.0 0.50 0.0 0.0
21 0.0 0.00 0.00 0.00 0.00 0.0 0.00 1.0 0.0
22 0.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0 1.0
In this notebook we want to model the HOM effect on ultimately a real quantum computer. This means we need to map the Fock states discussed above on a (qu)bit state. Above transition matrix has a basis of 9 states. So we would need 4 bits, with $2^{4} = 16$ states. As we know that some states in the Fock state basis are artificially set anyway we can also limit the basis to 8 states and use a register of 3 (qu)bits.
Anyway we will need a mapping table where we map a Fock state on a qubit register. This is a choice we make. There is no physics prescribing whether we map for state $|01>$ on qubit state '001' or '100'.
Let's take the simplest approach and follow the ordering of basis states implemented in Python QuTiP anyway. We first show how the Fock states are mapped on the QuTiP bases state via create_lookup_Fock_states(length_of_fock_state, no_of_channels). We then can do the same and map the QuTiP 'QubitCircuit' states on the same basis via create_lookup_QubitCircuit_states(number_of_qubits)
Lookup table for a system with two channels and maximum two photons in each channel (so either 0, 1 or 2 photons per channel)
Fock state <==> Qutip basis state
|00> <==> 100000000
|10> <==> 010000000
|20> <==> 001000000
|01> <==> 000100000
|11> <==> 000010000
|21> <==> 000001000
|02> <==> 000000100
|12> <==> 000000010
|22> <==> 000000001
If we ignore the 9th state in the Fock state basis (which we do not intend to use anyway, and already did set to an artificial transition coefficient to make the matrix unitary) we can keep the order of the basis and come to the mapping defined in merge_lookups(map, first_lookup, second_lookup)
Lookup table from Fock states to Qubits
Fock states <==> QubitCircuit qubits
|00> <==> 000
|01> <==> 001
|02> <==> 010
|10> <==> 011
|11> <==> 100
|12> <==> 101
|20> <==> 110
|21> <==> 111
With this mapping we can define the transition matrix between qubit states to model the beamsplitter (and model the HOM effect)
Transition matrix between registers of 3 qubits.
This matrix models the transition between Fock states in 2 channels, with total photon count of maximum 2
000 001 010 011 100 101 110 111
000 1.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0
001 0.0 0.71 0.00 -0.71 0.00 0.0 0.00 0.0
010 0.0 0.00 0.50 0.00 -0.71 0.0 0.50 0.0
011 0.0 0.71 0.00 0.71 0.00 0.0 0.00 0.0
100 0.0 0.00 0.71 0.00 0.00 0.0 -0.71 0.0
101 0.0 0.00 0.00 0.00 0.00 1.0 0.00 0.0
110 0.0 0.00 0.50 0.00 0.71 0.0 0.50 0.0
111 0.0 0.00 0.00 0.00 0.00 0.0 0.00 1.0
def create_lookup_Fock_states(length_of_fock_state, no_of_channels):
digits = [0]*(length_of_fock_state**no_of_channels)
states = [[]]
for c in range(no_of_channels):
states_new = []
for state in states:
for n in range(length_of_fock_state):
states_new.append(state + [n])
states = states_new[:]
lookup = dict([])
for state in states:
word = ''.join([str(n) for n in state])
bit_state = basis(length_of_fock_state,state[0])
for c in range(1,no_of_channels):
bit_state = tensor(basis(length_of_fock_state, state[c]), bit_state)
data = bit_state.full(order="C")
bit_string = ''.join((['1' if b == 1 else '0' for b in data]))
lookup[bit_string] = word
return lookup
print('Lookup table for a system with two channels and maximum two photons in each channel (so either 0, 1 or 2 photons per channel)')
print('Fock state <==> Qutip basis state')
for (key, value) in sorted(create_lookup_Fock_states(length_of_fock_state, no_of_channels).items(), key=lambda x: x[0], reverse=True):
print("|{}> <==> {}".format(value,key))
Lookup table for a system with two channels and maximum two photons in each channel (so either 0, 1 or 2 photons per channel) Fock state <==> Qutip basis state |00> <==> 100000000 |10> <==> 010000000 |20> <==> 001000000 |01> <==> 000100000 |11> <==> 000010000 |21> <==> 000001000 |02> <==> 000000100 |12> <==> 000000010 |22> <==> 000000001
from qutip.qip.circuit import QubitCircuit, Gate
number_of_qubits = 3
def create_lookup_QubitCircuit_states(number_of_qubits):
digits = [0]*(2**number_of_qubits)
states = [[]]
for c in range(number_of_qubits):
states_new = []
for state in states:
for n in range(2):
states_new.append(state + [n])
states = states_new[:]
lookup = dict([])
for state in states:
word = ''.join([str(n) for n in state])
bit_state = basis(2,state[0])
for c in range(1,number_of_qubits):
bit_state = tensor(bit_state, basis(2, state[c]))
data = bit_state.full(order="C")
bit_string = ''.join((['1' if b == 1 else '0' for b in data]))
lookup[bit_string] = word
return lookup
print('Lookup table for a system with 3 qubits')
print('QubitCircuit <==> Qutip basis state')
for (key, value) in sorted(create_lookup_QubitCircuit_states(number_of_qubits).items(), key=lambda x: x[0], reverse=True):
print("|{}> <==> {}".format(value,key))
Lookup table for a system with 3 qubits QubitCircuit <==> Qutip basis state |000> <==> 10000000 |001> <==> 01000000 |010> <==> 00100000 |011> <==> 00010000 |100> <==> 00001000 |101> <==> 00000100 |110> <==> 00000010 |111> <==> 00000001
def merge_lookups(map, first_lookup, second_lookup):
# merge two lookups. If the same value is present in both dictionaries the new lookup has item {key in first lookup: key in second lookup}
lookup = dict([])
first_keys = list(first_lookup.keys())
second_keys = list(second_lookup.keys())
for index in map:
lookup[first_lookup[first_keys[index]]] = second_lookup[second_keys[map[index]]]
return lookup
first_lookup = create_lookup_Fock_states(length_of_fock_state, no_of_channels)
second_lookup = create_lookup_QubitCircuit_states(number_of_qubits)
map = [index for index in range(min(len(first_lookup), len(second_lookup)))]
print('Lookup table from Fock states to Qubits')
print('Fock states <==> QubitCircuit qubits')
for (key, value) in merge_lookups(map, first_lookup, second_lookup).items():
print("|{}> <==> {}".format(key, value))
Lookup table from Fock states to Qubits Fock states <==> QubitCircuit qubits |00> <==> 000 |01> <==> 001 |02> <==> 010 |10> <==> 011 |11> <==> 100 |12> <==> 101 |20> <==> 110 |21> <==> 111 Fock states <==> QubitCircuit qubits |00> <==> 000 |01> <==> 001 |02> <==> 010 |10> <==> 011 |11> <==> 100 |12> <==> 101 |20> <==> 110 |21> <==> 111
# Calculate the transitions coefficients for the Fock states
r = t = 1/math.sqrt(2)
make_unitary = True
Fock_coefficients = create_Fock_coefficients(r, t, make_unitary)
# Define the matrix describing the HOM effect as a transition matrix between the qubit states
first_lookup = create_lookup_Fock_states(length_of_fock_state, no_of_channels)
second_lookup = create_lookup_QubitCircuit_states(number_of_qubits)
map = [index for index in range(min(len(first_lookup), len(second_lookup)))]
lookup_table = merge_lookups(map, first_lookup, second_lookup)
reverse_lookup_table = {value:key for key, value in lookup_table.items()}
qubit_transition_matrix = np.array(
[[Fock_coefficients[Fock_ver][Fock_hor] for Fock_hor in lookup_table.keys()] for Fock_ver in lookup_table.keys()]
)
print("Transition matrix between registers of 3 qubits.")
print("This matrix models the transition between Fock states in 2 channels, with total photon count of maximum 2")
print(pd.DataFrame(qubit_transition_matrix, columns = lookup_table.values(), index=lookup_table.values()))
Transition matrix between registers of 3 qubits.
This matrix models the transition between Fock states in 2 channels, with total photon count of maximum 2
000 001 010 011 100 101 110 111
000 1.0 0.00 0.00 0.00 0.00 0.0 0.00 0.0
001 0.0 0.71 0.00 -0.71 0.00 0.0 0.00 0.0
010 0.0 0.00 0.50 0.00 -0.71 0.0 0.50 0.0
011 0.0 0.71 0.00 0.71 0.00 0.0 0.00 0.0
100 0.0 0.00 0.71 0.00 0.00 0.0 -0.71 0.0
101 0.0 0.00 0.00 0.00 0.00 1.0 0.00 0.0
110 0.0 0.00 0.50 0.00 0.71 0.0 0.50 0.0
111 0.0 0.00 0.00 0.00 0.00 0.0 0.00 1.0
We use Qiskit to create a quantum circuit that simulates the beamsplitter. The circuit will have three qubit and three classical bits (to store the measurment results). We have a helper function to cycle through all input states in our lookup table and determine the probability of the various outcomes (this function is get_all_probabilities(circuit, lookup_table) ). We also define a helper function to plot the unitary matrix of the quantum circuit (print_circuit_unitary(circuit, lookup_table)).
First we create a custom gate which exactly contains the transition matrix defined before as qubit_transition_matrix.:
circuit.unitary(Operator(qubit_transition_matrix), circuit.qubits, label='HOM simulation')
We get the expected outcome when running the circuit:
Input Fock state: |00> (qubit |000>)
---- Outcome Fock state: |00> (qubit |000>) with probability 100%
Input Fock state: |01> (qubit |001>)
---- Outcome Fock state: |10> (qubit |011>) with probability 49%
---- Outcome Fock state: |01> (qubit |001>) with probability 51%
Input Fock state: |02> (qubit |010>)
---- Outcome Fock state: |02> (qubit |010>) with probability 25%
---- Outcome Fock state: |11> (qubit |100>) with probability 50%
---- Outcome Fock state: |20> (qubit |110>) with probability 24%
Input Fock state: |10> (qubit |011>)
---- Outcome Fock state: |10> (qubit |011>) with probability 50%
---- Outcome Fock state: |01> (qubit |001>) with probability 50%
Input Fock state: |11> (qubit |100>)
---- Outcome Fock state: |02> (qubit |010>) with probability 50%
---- Outcome Fock state: |20> (qubit |110>) with probability 50%
Input Fock state: |12> (qubit |101>)
---- Outcome Fock state: |12> (qubit |101>) with probability 100%
Input Fock state: |20> (qubit |110>)
---- Outcome Fock state: |11> (qubit |100>) with probability 50%
---- Outcome Fock state: |02> (qubit |010>) with probability 26%
---- Outcome Fock state: |20> (qubit |110>) with probability 23%
Input Fock state: |21> (qubit |111>)
---- Outcome Fock state: |21> (qubit |111>) with probability 100%
If we then check the unitary matrix for the total circuit we see indeed the correct results
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 0.71 . -0.71 . . . .
010 . . 0.5 . -0.71 . 0.5 .
011 . 0.71 . 0.71 . . . .
100 . . 0.71 . . . -0.71 .
101 . . . . . 1.0 . .
110 . . 0.5 . 0.71 . 0.5 .
111 . . . . . . . 1.0
So this is not very exciting. We have simply created a custom gate by using the transition coefficients between the different Fock states and by using a lookup table mapping the Fock states on the qubit states. However, with this infrastructure setup we can now try to build a quantum circuit from more standard gates (so not the custom gate) to come to the same result.
# Create the quantum circuit using Qiskit
from qiskit import *
#from qiskit import QuantumCircuit, execute
from qiskit.quantum_info.operators import Operator
#from qiskit import Aer
from qiskit_aer import AerSimulator
# Define a function to loop through all qubit states in the lookup_table and return the probabilities of the results for a given quantum circuit
def get_all_probabilities(circuit, lookup_table):
outcomes = dict([])
for qubit_state in lookup_table.values():
number_of_shots = 1000
simulator = AerSimulator()
temp_q = QuantumRegister(circuit.num_qubits,'q')
temp_c = ClassicalRegister(circuit.num_clbits,'c')
if len(qubit_state) != circuit.num_qubits:
print('Length of qubit state in lookup table to matching the quantum circuit')
temp_circuit = QuantumCircuit(temp_q,temp_c)
temp_circuit.initialize(qubit_state, temp_circuit.qubits)
temp_circuit = temp_circuit.compose(circuit)
temp_circuit.measure([i for i in range(circuit.num_qubits)], [i for i in range(circuit.num_clbits)])
job = simulator.run(temp_circuit, shots=number_of_shots)
result = job.result()
counts = result.get_counts(temp_circuit)
for word in counts.keys():
counts[word] = counts[word]/number_of_shots
outcomes.update({qubit_state : counts})
return outcomes
# Define a function to get an print the overall unitary for a given circuit
def print_circuit_unitary(circuit, lookup_table):
precision = 2
backend = Aer.get_backend('unitary_simulator')
job = execute(circuit, backend)
result = job.result()
matrix = np.asarray(result.get_unitary(circuit,precision)).round(precision)
table = [[0 for j in range(len(lookup_table.values()))] for i in range(len(lookup_table.values()))]
for index_h, bit_h in enumerate(lookup_table.values()):
for index_v, bit_v in enumerate(lookup_table.values()):
number = matrix[int(bit_h,2)][int(bit_v,2)].real
if number != 0:
table[index_h][index_v] = matrix[int(bit_h,2)][int(bit_v,2)].real
else:
table[index_h][index_v] = '.'
df = pd.DataFrame(table, columns = lookup_table.values(), index=lookup_table.values())
print(df)
return
# Define a quantum circuit with a custom unitary build from the known transition matrix
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_custom_unitary = QuantumCircuit(q,c)
circuit_custom_unitary.unitary(Operator(qubit_transition_matrix), circuit_custom_unitary.qubits, label='HOM simulation')
outcomes = get_all_probabilities(circuit_custom_unitary, lookup_table)
for input_qubit_state, outcome in outcomes.items():
print("\nInput Fock state: |{}> (qubit |{}>)".format(reverse_lookup_table[input_qubit_state], input_qubit_state))
for output_qubit_state, probability in outcome.items():
print(
"---- Outcome Fock state: |{}> (qubit |{}>) with probability {:.0%}".format(
reverse_lookup_table[output_qubit_state],
output_qubit_state,
probability)
)
Input Fock state: |00> (qubit |000>) ---- Outcome Fock state: |00> (qubit |000>) with probability 100% Input Fock state: |01> (qubit |001>) ---- Outcome Fock state: |10> (qubit |011>) with probability 51% ---- Outcome Fock state: |01> (qubit |001>) with probability 49% Input Fock state: |02> (qubit |010>) ---- Outcome Fock state: |11> (qubit |100>) with probability 49% ---- Outcome Fock state: |02> (qubit |010>) with probability 26% ---- Outcome Fock state: |20> (qubit |110>) with probability 25% Input Fock state: |10> (qubit |011>) ---- Outcome Fock state: |10> (qubit |011>) with probability 51% ---- Outcome Fock state: |01> (qubit |001>) with probability 49% Input Fock state: |11> (qubit |100>) ---- Outcome Fock state: |02> (qubit |010>) with probability 50% ---- Outcome Fock state: |20> (qubit |110>) with probability 50% Input Fock state: |12> (qubit |101>) ---- Outcome Fock state: |12> (qubit |101>) with probability 100% Input Fock state: |20> (qubit |110>) ---- Outcome Fock state: |11> (qubit |100>) with probability 52% ---- Outcome Fock state: |20> (qubit |110>) with probability 23% ---- Outcome Fock state: |02> (qubit |010>) with probability 25% Input Fock state: |21> (qubit |111>) ---- Outcome Fock state: |21> (qubit |111>) with probability 100%
print_circuit_unitary(circuit_custom_unitary, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 0.71 . -0.71 . . . . 010 . . 0.5 . -0.71 . 0.5 . 011 . 0.71 . 0.71 . . . . 100 . . 0.71 . . . -0.71 . 101 . . . . . 1.0 . . 110 . . 0.5 . 0.71 . 0.5 . 111 . . . . . . . 1.0
circuit_custom_unitary.draw(output='mpl')
We want to create a quantum circuit to simulate the beamsplitter and the HOM effect using standard gates. The goal is to create the same unitary matrix as we had from the custom gate defined in previous section. Obvsiouly there are automated functions available to do this (like Qiskit transpile), but is insightfull to try by hand. An easy and step-by-step way to resolve this manually is to start from the custom gate and then extend with standard gates to turn the matrix into the identity matrix. So we want to get to $I = M \cdot Gate_1 \cdot Gate_2 \cdot Gate_3 \ldots$ (where $I$ is the identity matrix, $M$ the custom unitary and $Gate_i$ a standard gate). So we add gates to the custom gate to create a circuit which does nothing. We then know that the added gates are the reverse of the custom gate, and that reversing the added standard gates should replicate the custom gate. So
$\begin{array}{lcl} I & = & M \cdot M^{-1} \\[5pts] I & = & M \cdot Gate_1 \cdot Gate_2 \cdot Gate_3 \ldots Gate_{n-1} \cdot Gate_n \\[5pts] M^{-1} & = & Gate_1 \cdot Gate_2 \cdot Gate_3 \ldots Gate_{n-1} \cdot Gate_n \\[5pts] M & = & (Gate_1 \cdot Gate_2 \cdot Gate_3 \ldots Gate_{n-1} \cdot Gate_n)^{-1} \\[5pts] M & = & Gate_n^{-1} \cdot Gate_{n-1}^{-1} \ldots Gate_3^{-1} \cdot Gate_2^{-1} \cdot Gate_1^{-1} \end{array}$
We start by the full unitary matrix which comes from our 'custom'
circuit.unitary(Operator(qubit_transition_matrix), circuit.qubits, label='HOM simulation')
This is the orginal 'custom gate'
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 0.71 . -0.71 . . . .
010 . . 0.5 . -0.71 . 0.5 .
011 . 0.71 . 0.71 . . . .
100 . . 0.71 . . . -0.71 .
101 . . . . . 1.0 . .
110 . . 0.5 . 0.71 . 0.5 .
111 . . . . . . . 1.0
The idea is now to add gates after this custom gate to ulimately get to the identity transformation. We can start by removing the mixing between states '001' and '011'.
circuit.x(2)
circuit.unitary(CCH, [q[1], q[0], q[2]], label='CCH')
circuit.x(2)
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 1.0 . . . . . .
010 . . 0.5 . -0.71 . 0.5 .
011 . . . -1.0 . . . .
100 . . 0.71 . . . -0.71 .
101 . . . . . 1.0 . .
110 . . 0.5 . 0.71 . 0.5 .
111 . . . . . . . 1.0
We see that for states '001' and '011' the matrix has simplified and these states only transition to themselves. The idea is to continue untill we have all 1's on the diagonal. Let's remove the factor '0.5' and '-0.5'
circuit.x(0)
circuit.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit.x(0)
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 1.0 . . . . . .
010 . . 0.71 . . . 0.71 .
011 . . . -1.0 . . . .
100 . . 0.71 . . . -0.71 .
101 . . . . . 1.0 . .
110 . . . . -1.0 . . .
111 . . . . . . . 1.0
A few more steps to get to a fully diagonal matrix
circuit.x(0)
circuit.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit.toffoli(2,0,1)
circuit.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit.x(0)
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 1.0 . . . . . .
010 . . 1.0 . . . . .
011 . . . -1.0 . . . .
100 . . . . -1.0 . . .
101 . . . . . 1.0 . .
110 . . . . . . 1.0 .
111 . . . . . . . 1.0
Basically we are there, but to get to the full identity matrix let's also remove the minus sign for states '011' and '100'
circuit.cnot(control_qubit=2, target_qubit=0)
circuit.cnot(control_qubit=2, target_qubit=1)
circuit.cp(math.pi, control_qubit=0, target_qubit=1)
circuit.cnot(control_qubit=2, target_qubit=0)
circuit.cnot(control_qubit=2, target_qubit=1)
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 1.0 . . . . . .
010 . . 1.0 . . . . .
011 . . . 1.0 . . . .
100 . . . . 1.0 . . .
101 . . . . . 1.0 . .
110 . . . . . . 1.0 .
111 . . . . . . . 1.0
We now have added standard gates to our custom gate to get to the identity matrix. By inverting this sequence of standard gates we can build our custom gate.
# Define the double controlled Hadamard gate as a 'helper'. later we see how to build this from more elementary gates.
s = math.sqrt(2)/2
CCH = Operator([
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, s, s],
[0, 0, 0, 0, 0, 0, s, -s]
])
# Define a quantum circuit with a custom unitary build from the known transition matrix
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit = QuantumCircuit(q,c)
circuit.unitary(Operator(qubit_transition_matrix), circuit.qubits, label='HOM simulation')
circuit.barrier()
print("This is the orginal 'custom gate'")
print_circuit_unitary(circuit, lookup_table)
This is the orginal 'custom gate'
000 001 010 011 100 101 110 111
000 1.0 . . . . . . .
001 . 0.71 . -0.71 . . . .
010 . . 0.5 . -0.71 . 0.5 .
011 . 0.71 . 0.71 . . . .
100 . . 0.71 . . . -0.71 .
101 . . . . . 1.0 . .
110 . . 0.5 . 0.71 . 0.5 .
111 . . . . . . . 1.0
circuit.x(2)
circuit.unitary(CCH, [q[1], q[0], q[2]], label='CCH')
circuit.x(2)
circuit.barrier()
print_circuit_unitary(circuit, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 1.0 . . . . . . 010 . . 0.5 . -0.71 . 0.5 . 011 . . . -1.0 . . . . 100 . . 0.71 . . . -0.71 . 101 . . . . . 1.0 . . 110 . . 0.5 . 0.71 . 0.5 . 111 . . . . . . . 1.0
circuit.x(0)
circuit.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit.x(0)
circuit.barrier()
circuit.x(0)
circuit.toffoli(2,0,1)
circuit.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit.x(0)
circuit.barrier()
print_circuit_unitary(circuit, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 1.0 . . . . . . 010 . . 1.0 . . . . . 011 . . . -1.0 . . . . 100 . . . . -1.0 . . . 101 . . . . . 1.0 . . 110 . . . . . . 1.0 . 111 . . . . . . . 1.0
circuit.cnot(control_qubit=2, target_qubit=0)
circuit.cnot(control_qubit=2, target_qubit=1)
circuit.cp(math.pi, control_qubit=0, target_qubit=1)
circuit.cnot(control_qubit=2, target_qubit=0)
circuit.cnot(control_qubit=2, target_qubit=1)
circuit.barrier()
print_circuit_unitary(circuit, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 1.0 . . . . . . 010 . . 1.0 . . . . . 011 . . . 1.0 . . . . 100 . . . . 1.0 . . . 101 . . . . . 1.0 . . 110 . . . . . . 1.0 . 111 . . . . . . . 1.0
We can draw the circuit that starts with the custom gate and ultimately result in the identity gate. If we only take the gates behind the custom gate and invert their order we should get the custom gate back.
circuit.draw(output='mpl')
# Define a quantum circuit with a custom unitary build from the known transition matrix
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_manual = QuantumCircuit(q,c)
circuit_manual.cnot(control_qubit=2, target_qubit=0)
circuit_manual.cnot(control_qubit=2, target_qubit=1)
circuit_manual.cp(math.pi, control_qubit=0, target_qubit=1)
circuit_manual.cnot(control_qubit=2, target_qubit=0)
circuit_manual.cnot(control_qubit=2, target_qubit=1)
circuit_manual.x(0)
circuit_manual.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit_manual.toffoli(2,0,1)
circuit_manual.unitary(CCH, [q[2], q[0], q[1]], label='CCH')
circuit_manual.x(0)
circuit_manual.x(2)
circuit_manual.unitary(CCH, [q[1], q[0], q[2]], label='CCH')
circuit_manual.x(2)
print_circuit_unitary(circuit_manual, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 0.71 . -0.71 . . . . 010 . . 0.5 . -0.71 . 0.5 . 011 . 0.71 . 0.71 . . . . 100 . . 0.71 . . . -0.71 . 101 . . . . . 1.0 . . 110 . . 0.5 . 0.71 . 0.5 . 111 . . . . . . . 1.0
Indeed this circuit of standard gates exacly replicates the custom gate. The circuit is drawn below.
circuit_manual.draw(output='mpl')
The disadvantage of our circuit is that it contains double controlled gates. We can easily replace them by single controlled gates using the right decomposition. First decompose the double controlled Hadamard gate.
from qiskit.circuit.library import HGate
gate = HGate().power(1/2).control(1)
from qiskit import QuantumCircuit
circuit_CCH = QuantumCircuit(3)
circuit_CCH.append(gate, [1,0])
circuit_CCH.cnot(target_qubit=1, control_qubit=2)
circuit_CCH.append(gate.inverse(), [1,0])
circuit_CCH.cnot(target_qubit=1, control_qubit=2)
circuit_CCH.append(gate, [2,0])
circuit_CCH.draw('mpl')
print_circuit_unitary(circuit_CCH, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 1.0 . . . . . . 010 . . 1.0 . . . . . 011 . . . 1.0 . . . . 100 . . . . 1.0 . . . 101 . . . . . 1.0 . . 110 . . . . . . 0.71 0.71 111 . . . . . . 0.71 -0.71
Then we decompose the Toffoli gate, or double controlled NOT gate into only 2-qubit gates:
from qiskit.circuit.library import XGate
gate = XGate().power(1/2).control(1)
from qiskit import QuantumCircuit
circuit_CCX = QuantumCircuit(3)
circuit_CCX.append(gate, [1,0])
circuit_CCX.cnot(target_qubit=1, control_qubit=2)
circuit_CCX.append(gate.inverse(), [1,0])
circuit_CCX.cnot(target_qubit=1, control_qubit=2)
circuit_CCX.append(gate, [2,0])
circuit_CCX.draw('mpl')
print_circuit_unitary(circuit_CCX, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 1.0 . . . . . . 010 . . 1.0 . . . . . 011 . . . 1.0 . . . . 100 . . . . 1.0 . . . 101 . . . . . 1.0 . . 110 . . . . . . . 1.0 111 . . . . . . 1.0 .
# Build our circuit only including 2 - qubit gates by decomposing the 3-qubit gates used before
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_manual_2 = QuantumCircuit(q,c)
circuit_manual_2.cnot(control_qubit=2, target_qubit=0)
circuit_manual_2.cnot(control_qubit=2, target_qubit=1)
circuit_manual_2.cp(math.pi, control_qubit=0, target_qubit=1)
circuit_manual_2.cnot(control_qubit=2, target_qubit=0)
circuit_manual_2.cnot(control_qubit=2, target_qubit=1)
circuit_manual_2.x(0)
circuit_manual_2.compose(circuit_CCH, [2,0,1], inplace=True)
circuit_manual_2.compose(circuit_CCX, [1,0,2], inplace=True)
circuit_manual_2.compose(circuit_CCH, [2,0,1], inplace=True)
circuit_manual_2.x(0)
circuit_manual_2.x(2)
circuit_manual_2.compose(circuit_CCH, [1,0,2], inplace=True)
circuit_manual_2.x(2)
print_circuit_unitary(circuit_manual_2, lookup_table)
000 001 010 011 100 101 110 111 000 1.0 . . . . . . . 001 . 0.71 . -0.71 . . . . 010 . . 0.5 . -0.71 . 0.5 . 011 . 0.71 . 0.71 . . . . 100 . . 0.71 . . . -0.71 . 101 . . . . . 1.0 . . 110 . . 0.5 . 0.71 . 0.5 . 111 . . . . . . . 1.0
circuit_manual_2.draw(output='mpl')
In this section we started of with a 3 qubit quantum circuit and a single 'custom' gate. The circuit is called
circuit_custom_unitaryThe custom gate was defined by the unitary matrix describing the transition between Fock states for our beamsplitter. With the propose mapping in thelookup_table. We then manually decomposed the circuit by using three qubit gates (the double controlled Hadamard gate) and the double controlled NOT gate (aka Toffoli gate). This circuit is calledcircuit_manualWe then manually further decomposed the circuit by using only 2 qubit gates. The result is calledcircuit_manual_2
We can continue and for instance include the allowed coupling between qubits in a specific execution of the quantum computer (for instance, by adding swap gates we can make sure only neighbouring qubits control eachother, and re-write gates where qubit 0 controls qubit 2). However, before execution on a quantum device the circuit is transpiled anyway and those algorithms are more specific and more efficient than we can be manually. So we stop now with the circuit that simulates the beamsplitter using only standard, single and two qubit gates.
Next we want to run our circuits on a quantum computer. The code below starts by re-defining the three circuits we want to use. In this way we decouple from any changes made in the notebook above and also make it easier to copy into the "IBM Quantum Lab" environment. But before we can define the circuits we re-create some helper functions like the operator for the double controlled Hadamard gate, the decomponisitions of the double control gates, our target transition matrix and the lookup table between Fock states and qubit states. All this we used above, but is redefined to make a 'clean start' which can easily be copied into a different environment.
# First get all the libraries sorted
# Create the quantum circuit using Qiskit
from qiskit import *
#from qiskit import QuantumCircuit, execute
from qiskit.quantum_info.operators import Operator
#from qiskit import Aer
from qiskit_aer import AerSimulator
# Imported libraries
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.lines
import numpy as np
import math
from IPython.display import Image # for displaying picture
import pandas as pd # used for easily printing tables
Below we define the operator for the double controlled Hadamard gate. If the two control qubits are set to 1 the Hadamard gate is applied to the third qubit.
# Define the double controlled Hadamard gate as a 'helper'.
s = math.sqrt(2)/2
CCH_operator = Operator([
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, s, s],
[0, 0, 0, 0, 0, 0, s, -s]
])
And we define replacement circuit for the double controlled Hadamard gate and for the Toffoli gate (double controlled NOT gate). This decomposition we use to make a circuit that avoid 3 qubit gates and only relies on gates connecting two qubits (next to single qubit gates)
# Define the circuit to replace the double controlled NOT gate (aka Toffoli gate) by a circuit using only 2 bit gates.
from qiskit.circuit.library import XGate
gate = XGate().power(1/2).control(1)
from qiskit import QuantumCircuit
circuit_CCX = QuantumCircuit(3)
circuit_CCX.append(gate, [1,0])
circuit_CCX.cnot(target_qubit=1, control_qubit=2)
circuit_CCX.append(gate.inverse(), [1,0])
circuit_CCX.cnot(target_qubit=1, control_qubit=2)
circuit_CCX.append(gate, [2,0])
circuit_CCX.draw('mpl')
# Define the circuit to replace the double controlled Hadamard gate by a circuit using only 2 bit gates.
from qiskit.circuit.library import HGate
gate = HGate().power(1/2).control(1)
from qiskit import QuantumCircuit
circuit_CCH = QuantumCircuit(3)
circuit_CCH.append(gate, [1,0])
circuit_CCH.cnot(target_qubit=1, control_qubit=2)
circuit_CCH.append(gate.inverse(), [1,0])
circuit_CCH.cnot(target_qubit=1, control_qubit=2)
circuit_CCH.append(gate, [2,0])
circuit_CCH.draw('mpl')
And we define the overall transition matrix for our beamsplitter. This matrix has been calculated above. Below we define it directly based on values rather than calculate the values. This is to easier copy to code to different environments, like "IBM Quantum Lab", without having to run the full calculation.
qubit_transition_matrix = [
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.7071067811865475, 0.0, -0.7071067811865475, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.4999999999999999, 0.0, -0.7071067811865475, 0.0, 0.4999999999999999, 0.0],
[0.0, 0.7071067811865475, 0.0, 0.7071067811865475, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.7071067811865475, 0.0, 0.0, 0.0, -0.7071067811865475, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.4999999999999999, 0.0, 0.7071067811865475, 0.0, 0.4999999999999999, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
]
And finally we define the look-up table between the Fock states and the qubit states.
# Define the matrix describing the HOM effect as a transition matrix between the qubit states
lookup_table = {'00': '000', '01': '001', '02': '010', '10': '011', '11': '100', '12': '101', '20': '110', '21': '111'}
reverse_lookup_table = {value:key for key, value in lookup_table.items()}
Now we are ready to define the three circuits we want to run on a real quantum computer. These circuits are:
The circuit called
circuit_custom_unitary. This is based on a single gate which is the 'custom gate' defined by the unitary matrix describing the transition between Fock states for our beamsplitter.
The circuit made by manually decomposing the custom gate using three qubit gates (the double controlled Hadamard gate and the double controlled NOT gate). This circuit is called
circuit_manual
The third circuit results from further decomposition until we use only 2 qubit gates. This is called
circuit_manual_2
# Redefine the circuit based on unitary matrix
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_custom_unitary = QuantumCircuit(q,c)
circuit_custom_unitary.unitary(Operator(qubit_transition_matrix), circuit_custom_unitary.qubits, label='HOM simulation')
# Redefine the circuit with "manual transpile"
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_manual = QuantumCircuit(q,c)
circuit_manual.cnot(control_qubit=2, target_qubit=0)
circuit_manual.cnot(control_qubit=2, target_qubit=1)
circuit_manual.cp(math.pi, control_qubit=0, target_qubit=1)
circuit_manual.cnot(control_qubit=2, target_qubit=0)
circuit_manual.cnot(control_qubit=2, target_qubit=1)
circuit_manual.x(0)
circuit_manual.unitary(CCH_operator, [q[2], q[0], q[1]], label='CCH')
circuit_manual.toffoli(2,0,1)
circuit_manual.unitary(CCH_operator, [q[2], q[0], q[1]], label='CCH')
circuit_manual.x(0)
circuit_manual.x(2)
circuit_manual.unitary(CCH_operator, [q[1], q[0], q[2]], label='CCH')
circuit_manual.x(2)
# Redefine the circuit with "manual transpile" and onlu use 2-qubit gates
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
circuit_manual_2 = QuantumCircuit(q,c)
circuit_manual_2.cnot(control_qubit=2, target_qubit=0)
circuit_manual_2.cnot(control_qubit=2, target_qubit=1)
circuit_manual_2.cp(math.pi, control_qubit=0, target_qubit=1)
circuit_manual_2.cnot(control_qubit=2, target_qubit=0)
circuit_manual_2.cnot(control_qubit=2, target_qubit=1)
circuit_manual_2.x(0)
circuit_manual_2.compose(circuit_CCH, [2,0,1], inplace=True)
circuit_manual_2.compose(circuit_CCX, [1,0,2], inplace=True)
circuit_manual_2.compose(circuit_CCH, [2,0,1], inplace=True)
circuit_manual_2.x(0)
circuit_manual_2.x(2)
circuit_manual_2.compose(circuit_CCH, [1,0,2], inplace=True)
circuit_manual_2.x(2)
<qiskit.circuit.instructionset.InstructionSet at 0x1f2fd6ae500>
We run these circuits first on a simulated backend and compare the results to an ideal, noiseless, situation. For the noiseless reference we usen Qiskit's build in 'AerSimulator' as backend. We compare results to the 'FakeLimaV2' backend (see https://qiskit.org/documentation/apidoc/providers_fake_provider.html for documentation).
Two different visualizations are used to compare the results. First we use a two-dimensional plot where the gray scale indicates the transition probability. Black is a 100% probability and white is 0% probability. We can compare how this picture looks for the noiseless simulation and for the simulations on the fake backend. We clearly see the white spots in the noiseless simulation disappear for the simulations on the fake backend. This illustrates the noise we have to deal with on a real quantum computer. Note that we used a non-linear scaling of the color bar to make the low noise values visible. We do also see that the result is best when using the custom gate as input. It gets worse if we use our 'manually transpiled circuits'. The reason is that the circuit is transpiled anyway for optimal use on the backend and the algorithms in Qiskit are doing a better job at this than we did. Giving these algorthims the orginal circuit allows for more efficient optimization.
#Define the function to run a circuit on a given back-end. The function returns a matrix with observed transition probabilities.
def run_circuit_on_backend(circuit_to_run, backend):
number_of_shots = 1000
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
initial_states = lookup_table.values()
list_of_circuits = []
for state in initial_states:
circuit_in_list = QuantumCircuit(q,c)
circuit_in_list.initialize(state, q)
circuit_in_list.compose(circuit_to_run, [0, 1, 2], inplace=True)
circuit_in_list.measure_all()
transpiled_circuit = transpile(circuit_in_list, backend)
list_of_circuits.append(transpiled_circuit)
job = backend.run(list_of_circuits, shots=number_of_shots)
matrix = []
for index in range(len(initial_states)):
res = job.result().get_counts(index)
matrix_line = []
for output_state in lookup_table.values():
word = output_state + ' 000'
if word in res.keys():
matrix_line.append(1.0*res[word]/number_of_shots)
else:
matrix_line.append(0)
matrix.append(matrix_line)
return matrix
from qiskit.providers.fake_provider import FakeLimaV2
result_circuit_noiseless = run_circuit_on_backend(circuit_custom_unitary, AerSimulator())
result_circuit_custom_unitary = run_circuit_on_backend(circuit_custom_unitary, FakeLimaV2())
result_circuit_manual = run_circuit_on_backend(circuit_manual, FakeLimaV2())
result_circuit_manual_2 = run_circuit_on_backend(circuit_manual_2, FakeLimaV2())
plt.rcParams['figure.figsize'] = [15, 6]
fig, ax = plt.subplots(1,4)
colormap = 'binary'
fig.suptitle(' Results from running on \'fake\' quantum computer ', fontsize=30)
ax[0].set_title('Noiseless result \n using \'AerSimulator()\'\n')
ax[1].set_title('\'FakeLimaV2()\' \n using custom unitary \n no manual transpilation\n')
ax[2].set_title('\'FakeLimaV2()\'\n using manual transpilation \n (3 qubit gates)\n')
ax[3].set_title('\'FakeLimaV2()\'\n using manual transpilation \n (2 qubit gates)\n')
for i in range(4):
ax[i].set_xticks([i for i in range(len(lookup_table.values()))])
ax[i].set_yticks([i for i in range(len(lookup_table.values()))])
ax[i].set_xticklabels(lookup_table.values())
ax[i].set_yticklabels(lookup_table.values())
# We use a powernorm to scale the color coding. A linear scale would hide the noise too much and differences would not be visible.
im = ax[0].matshow(result_circuit_noiseless, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[1].matshow(result_circuit_custom_unitary, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[2].matshow(result_circuit_manual, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[3].matshow(result_circuit_manual_2, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
cbar_ax = fig.add_axes([0.92, 0.25, 0.01, 0.55])
fig.colorbar(im, cax=cbar_ax)
<matplotlib.colorbar.Colorbar at 0x1f2fd822f10>
width = 0.2
state_11, state_02, state_20 = list(lookup_table.keys()).index('11'), list(lookup_table.keys()).index('02'), list(lookup_table.keys()).index('20')
colors = ['blue', 'red', 'green', 'purple']
for i, result in enumerate([result_circuit_noiseless, result_circuit_custom_unitary, result_circuit_manual, result_circuit_manual_2]):
plt.bar([x+(i-2)*width for x in range(3)],
[
result[state_11][state_02],
result[state_11][state_11]+0.002, # For visual effect lift the value slightly above zero
result[state_11][state_20]
],
width,
color = colors[i]
)
plt.xticks([x for x in range(3)], [
'Output state |02>',
'Output state |11>',
'Output state |20>'
])
plt.ylabel("Probability")
custom_lines = [matplotlib.lines.Line2D([0], [0], color='blue', lw=4),
matplotlib.lines.Line2D([0], [0], color='red', lw=4),
matplotlib.lines.Line2D([0], [0], color='green', lw=4),
matplotlib.lines.Line2D([0], [0], color='purple', lw=4)
]
plt.legend(custom_lines, ['Noiseless result \n using \'AerSimulator()\'',
'\'FakeLimaV2()\' \n using custom unitary \n no manual transpilation',
'\'FakeLimaV2()\' \n using manual transpilation \n (3 qubit gates)',
'\'FakeLimaV2()\' \n using manual transpilation \n (2 qubit gates)'
])
text = 'Probabilities of finding the output states for input Fock state |11> \n'
text += 'For ideal simulation of the HOM effect we would see 50%/50% in Fock states |02> and |20> \n'
text += 'The observed probability for Fock state |11> is indicative of the noise level in the simulation'
plt.title(text)
plt.show()
This is the code we use to run the circuits on the IMB quantum backend. We use the 'ibmq_lima'device (see https://quantum-computing.ibm.com/services/resources?tab=systems) and use 'runtime' for execution. This part of the code we execute on the IBM server via IBM quantum Lab. This avoids complexity with SSL licenses and credentials.
We run the three circuits sequentially. As the waiting time can be up to a few days we fetch the results from the IMB dashboard. Once the circuits are defined in IBM quantum lab we use below code to run the circuits.
circuit_to_run = circuit_custom_unitary
#circuit_to_run = circuit_manual_2
#circuit_to_run = circuit_manual_2
backend = service.get_backend('ibmq_lima')
q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
initial_states = lookup_table.values()
list_of_circuits = []
for state in initial_states:
circuit_in_list = QuantumCircuit(q,c)
circuit_in_list.initialize(state, q)
circuit_in_list.compose(circuit_to_run, [0, 1, 2], inplace=True)
circuit_in_list.measure_all()
transpiled_circuit = transpile(circuit_in_list, backend)
list_of_circuits.append(transpiled_circuit)
from qiskit_ibm_runtime import (
QiskitRuntimeService,
Session,
Sampler,
Estimator,
Options,
)
service = QiskitRuntimeService(channel="ibm_quantum")
options = Options(optimization_level=3, resilience_level=1)
sampler = Sampler(session=backend)
job = sampler.run(list_of_circuits)
print(job.job_id())
The results are fetched (after a few days) from the IBM site. For processing we copy the raw 'quasi distribution' below. This quasi distribution is effectively the observed probability for transitioning from input state to output state.
# These are the results fetched from the IBM Q backend.
quasi_dist_circuit_custom_unitary = [
{
"100000": 0.046623767723592974,
"101000": 0.01695406132793082,
"110000": 0.05433316345302092,
"111000": 0.021928979879237007,
"000000": 0.709766005255246,
"001000": 0.07422933951805274,
"010000": 0.03890223489003641,
"011000": 0.037262447952883145
},
{
"100000": 0.04594482059609681,
"101000": 0.054433876903609665,
"110000": 0.030507078335165985,
"111000": 0.025431099348189573,
"000000": 0.04775156979721547,
"001000": 0.4119664100787966,
"010000": 0.05662204287414959,
"011000": 0.32734310206677614
},
{
"100000": 0.4399491825582809,
"101000": 0.04200221417018358,
"110000": 0.22551600940887095,
"111000": 0.027574929845079164,
"000000": 0.0262353457114064,
"001000": 0.03214812688287603,
"010000": 0.16427487970292212,
"011000": 0.04229931172038081
},
{
"100000": 0.022648463882807897,
"101000": 0.035437904483281206,
"110000": 0.02269237385834693,
"111000": 0.05958430643868795,
"000000": 0.035983399180089866,
"001000": 0.39785693181127785,
"010000": 0.042192469917452916,
"011000": 0.3836041504280552
},
{
"100000": 0.03933878578193142,
"101000": 0.04757490981889651,
"110000": 0.35161318962671034,
"111000": 0.03425113895757403,
"000000": 0.03430253744599728,
"001000": 0.031710149266330014,
"010000": 0.41816712805466677,
"011000": 0.04304216104789359
},
{
"100000": 0.051092380700449885,
"101000": 0.7270709623238296,
"110000": 0.04865507988864886,
"111000": 0.03188108249796108,
"000000": 0.024504249423010258,
"001000": 0.03986782983164607,
"010000": 0.02426303955900279,
"011000": 0.05266537577545146
},
{
"100000": 0.31527996870053965,
"101000": 0.044449606105537255,
"110000": 0.27856762065744406,
"111000": 0.03518078052565727,
"000000": 0.05782993523156828,
"001000": 0.013475061539461649,
"010000": 0.23620554583180453,
"011000": 0.01901148140798729
},
{
"100000": 0.02465127998176131,
"101000": 0.02987416650923599,
"110000": 0.054016069456806454,
"111000": 0.7369467534130348,
"000000": 0.0240689596554278,
"001000": 0.06796398508189393,
"010000": 0.020670779068689107,
"011000": 0.04180800683315045
}
]
quasi_dist_circuit_manual = [
{
"100000": 0.10321835254492363,
"101000": 0.08497262294731132,
"110000": 0.08763947875232597,
"111000": 0.06414240472433555,
"000000": 0.2974638099838497,
"001000": 0.11384291855230176,
"010000": 0.14111411603990376,
"011000": 0.10760629645504817
},
{
"100000": 0.08891459734538953,
"101000": 0.09908440188329849,
"110000": 0.05463783761804659,
"111000": 0.07416629481070046,
"000000": 0.1174325092527995,
"001000": 0.3084171661167175,
"010000": 0.0826519480354107,
"011000": 0.1746952449376373
},
{
"100000": 0.17386521343141711,
"101000": 0.10844393853071635,
"110000": 0.1573194324325035,
"111000": 0.13256350493953134,
"000000": 0.11545007777801793,
"001000": 0.07987324889232715,
"010000": 0.1473453786162634,
"011000": 0.0851392053792232
},
{
"100000": 0.05489932941939224,
"101000": 0.08872741382294479,
"110000": 0.08886309718968884,
"111000": 0.0912999946754674,
"000000": 0.05607678736705861,
"001000": 0.16641192406343605,
"010000": 0.1054514123553473,
"011000": 0.3482700411066647
},
{
"100000": 0.12068748498320793,
"101000": 0.10575186720460086,
"110000": 0.16160339247601047,
"111000": 0.12126600090078103,
"000000": 0.09568648256494078,
"001000": 0.07843728700709345,
"010000": 0.24050406306915142,
"011000": 0.07606342179421408
},
{
"100000": 0.1068225458077494,
"101000": 0.25294977216842407,
"110000": 0.08237760952895498,
"111000": 0.15169908105258387,
"000000": 0.11444007589952813,
"001000": 0.10121642775803791,
"010000": 0.10444941512081898,
"011000": 0.08604507266390257
},
{
"100000": 0.2078484031525109,
"101000": 0.10534972627049384,
"110000": 0.20038661352963885,
"111000": 0.09989742577338971,
"000000": 0.10830125923083957,
"001000": 0.07206205989551906,
"010000": 0.1336264675366393,
"011000": 0.0725280446109688
},
{
"100000": 0.08331075064223938,
"101000": 0.14204448391538696,
"110000": 0.14381954743387385,
"111000": 0.3086926360840581,
"000000": 0.04817783952790376,
"001000": 0.08995501614846173,
"010000": 0.08990419729159436,
"011000": 0.09409552895648177
}
]
quasi_dist_circuit_manual_2 = [
{
"100000": 0.0903492397786909,
"101000": 0.07237151669164554,
"110000": 0.11162422881253728,
"111000": 0.06342475974200189,
"000000": 0.37198014850658045,
"001000": 0.1118341301004222,
"010000": 0.06637827963427219,
"011000": 0.11203769673384951
},
{
"100000": 0.04734387119624626,
"101000": 0.07833319853800583,
"110000": 0.05815418145162998,
"111000": 0.056313922252337056,
"000000": 0.1842697305063958,
"001000": 0.19763133586454038,
"010000": 0.14887512530655114,
"011000": 0.2290786348842935
},
{
"100000": 0.19494303951603584,
"101000": 0.08579843971446288,
"110000": 0.20733509208443945,
"111000": 0.12345796582081167,
"000000": 0.06847937588367675,
"001000": 0.11980019128934068,
"010000": 0.13138982085594586,
"011000": 0.06879607483528684
},
{
"100000": 0.05495280988159646,
"101000": 0.05113106204727031,
"110000": 0.06904746315947934,
"111000": 0.07878753040561266,
"000000": 0.08634789602806647,
"001000": 0.24541866734016754,
"010000": 0.12450718052707514,
"011000": 0.2898073906107321
},
{
"100000": 0.058744129550529926,
"101000": 0.08870217969878873,
"110000": 0.23094620891703896,
"111000": 0.10472158702632826,
"000000": 0.050022379427530926,
"001000": 0.06910758175669694,
"010000": 0.3189735723593165,
"011000": 0.07878236126376975
},
{
"100000": 0.09110930563857574,
"101000": 0.5127014786853866,
"110000": 0.0639195307145584,
"111000": 0.07207940866439336,
"000000": 0.05550377905738124,
"001000": 0.08443051311261295,
"010000": 0.04429206344634052,
"011000": 0.07596392068075124
},
{
"100000": 0.3744961840032315,
"101000": 0.08653387901063407,
"110000": 0.1565386200850637,
"111000": 0.096419412139166,
"000000": 0.058755893070804734,
"001000": 0.05359217689160483,
"010000": 0.12304792594968196,
"011000": 0.05061590884981317
},
{
"100000": 0.0700629332931067,
"101000": 0.08576318079514962,
"110000": 0.09732946749700226,
"111000": 0.37503807068908407,
"000000": 0.06871478397233899,
"001000": 0.06297604825515009,
"010000": 0.1481192331227688,
"011000": 0.09199628237539945
}
]
We have the date return from IBM quantum. Now we can visualize in the same way we visualized the data on the 'fake' backend. We make a heatmap with the full result where we start with the noiseless result as reference, and the plot the results from the three circuits.
The two-dimensional grey scale shows the noise level. The left-most picture is the noisless simulation included for reference. Clearly we see the noise added by the quantum computer. We also see, just as predicted from our runs on the fake quantum computer, that the circuit based on the custom gate performs best.
If we look at the quantum HOM effect specifically we can see how much the result on the real quantum computer shows that the outcome $ |11>$ is suppressed. For all runs on the real quantum computer we see a non-zero probability to go from $ |11>$ input to $ |11>$ output, but especially for the circuit based on our custom gate there is a nice contrast between the high probability expected in the $ |02>$ or $ |20>$ and the low probability in the $ |11>$.
# Create a function to process the quasi probability distribution in to a matrix, similar to the output we used from the simulators.
def data_from_quasi_dist(quasi_dist, lookup_table):
matrix = [[0 for val in lookup_table.values()] for distribution in quasi_dist]
for dist_index, distribution in enumerate(quasi_dist):
for key_word,value in distribution.items():
key = int(key_word,2)
word = format(int(key/8), '03b')
key_index = list(lookup_table.values()).index(word)
matrix[dist_index][key_index] = value
return matrix
data_circuit_custom_unitary = data_from_quasi_dist(quasi_dist_circuit_custom_unitary, lookup_table)
data_circuit_manual = data_from_quasi_dist(quasi_dist_circuit_manual, lookup_table)
data_circuit_manual_2 = data_from_quasi_dist(quasi_dist_circuit_manual_2, lookup_table)
plt.rcParams['figure.figsize'] = [15, 6]
fig, ax = plt.subplots(1,4)
fig.suptitle(' Results from simulation on real quantum computer ', fontsize=30)
ax[0].set_title('Noiseless result (from simulator)')
ax[1].set_title('Real quantum computer \n using custom unitary')
ax[2].set_title('Real quantum computer \n using manual transpilation \n (3 qubit gates)')
ax[3].set_title('Real quantum computer \n using manual transpilation \n (2 qubit gates)')
for i in range(4):
ax[i].set_xticks([i for i in range(len(lookup_table.values()))])
ax[i].set_yticks([i for i in range(len(lookup_table.values()))])
ax[i].set_xticklabels(lookup_table.values())
ax[i].set_yticklabels(lookup_table.values())
# We use a powernorm to scale the color coding. A linear scale would hide the noise too much and differences would not be visible.
im = ax[0].matshow(result_circuit_noiseless, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[1].matshow(data_circuit_custom_unitary, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[2].matshow(data_circuit_manual, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
im = ax[3].matshow(data_circuit_manual_2, cmap = colormap, norm=matplotlib.colors.PowerNorm(gamma = 0.5, vmin=0, vmax=1))
cbar_ax = fig.add_axes([0.92, 0.25, 0.01, 0.55])
fig.colorbar(im, cax=cbar_ax)
<matplotlib.colorbar.Colorbar at 0x1f2fef57650>
width = 0.2
state_11, state_02, state_20 = list(lookup_table.keys()).index('11'), list(lookup_table.keys()).index('02'), list(lookup_table.keys()).index('20')
colors = ['blue', 'red', 'green', 'purple']
for i, result in enumerate([result_circuit_noiseless, data_circuit_custom_unitary, data_circuit_manual, data_circuit_manual_2]):
plt.bar([x+(i-2)*width for x in range(3)],
[
result[state_11][state_02],
result[state_11][state_11]+0.002, # For visual effect lift the value slightly above zero
result[state_11][state_20]
],
width,
color = colors[i]
)
plt.xticks([x for x in range(3)], [
'Output state |02>',
'Output state |11>',
'Output state |20>'
])
plt.ylabel("Probability")
custom_lines = [matplotlib.lines.Line2D([0], [0], color='blue', lw=4),
matplotlib.lines.Line2D([0], [0], color='red', lw=4),
matplotlib.lines.Line2D([0], [0], color='green', lw=4),
matplotlib.lines.Line2D([0], [0], color='purple', lw=4)
]
plt.legend(custom_lines, ['Noiseless result (from simulator)',
'Real quantum computer \n using custom unitary',
'Real quantum computer \n using manual transpilation \n (3 qubit gates)',
'Real quantum computer \n using manual transpilation \n (2 qubit gates)'
])
text = 'Probabilities of finding the output states for input Fock state |11> \n'
text += 'For ideal simulation of the HOM effect we would see 50%/50% in Fock states |02> and |20> \n'
text += 'The observed probability for Fock state |11> is indicative of the noise level on the quantum computer'
plt.title(text)
plt.show()
