Cubic phase states using a circuit with three modes

Author: Nicolás Quesada

In this tutorial we study how to prepare a cubic phase state \(| \psi \rangle \sim | 0 \rangle +i a \sqrt{\frac{3}{2}} | 1 \rangle + i a | 3 \rangle\) by heralding two modes of a three mode circuit using photon-number-resolving detectors. This idea was proposed in “Production of photonic universal quantum gates enhanced by machine learning” Phys. Rev. A 100, 012326 (2019) by Krishna Kumar Sabapathy, Haoyu Qi, Josh Izaac, and Christian Weedbrook.

import numpy as np
from qutip import wigner, Qobj, wigner_cmap

import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import cm

import strawberryfields as sf
from strawberryfields.ops import *
from thewalrus.quantum import state_vector, density_matrix

Ideal preparation

Here we setup some basic parameters, like the value of the photon-number-resolving detectors we will use to herald and the amount of squeezing and displacement to use. These parameters are taken from the reference cited above.

# the Fock state measurement of mode 0 to be post-selected
m1 = 1
# the Fock state measurement of mode 1 to be post-selected
m2 = 2
sq_r = [0.71, 0.67, -0.42]
# squeezing phase
sq_phi = [-2.07, 0.06, -3.79]
# displacement magnitudes
d_r = [-0.02, 0.34, 0.02]
# beamsplitter theta
bs_theta1, bs_theta2, bs_theta3 = [-1.57, 0.68, 2.5]
# beamsplitter phi
bs_phi1, bs_phi2, bs_phi3 = [0.53, -4.51, 0.72]

Now we setup a 3-mode quantum circuit in Strawberry Fields and obtain the covariance matrix and vector of means of the Gaussian state.

nmodes = 3
prog = sf.Program(nmodes)
eng = sf.Engine("gaussian")

with prog.context as q:
    for k in range(3):
        Sgate(sq_r[k], sq_phi[k]) | q[k]
        Dgate(d_r[k]) | q[k]

    BSgate(bs_theta1, bs_phi1) | (q[0], q[1])
    BSgate(bs_theta2, bs_phi2) | (q[1], q[2])
    BSgate(bs_theta3, bs_phi3) | (q[0], q[1])

state = eng.run(prog).state
mu = state.means()
cov = state.cov()

# Here we use the sf circuit drawer and standard linux utilities
# to generate an svg representing the circuit
file, _ = prog.draw_circuit()
filepdf = file[0:-3]+"pdf"
filepdf = filepdf.replace("circuit_tex/","")
filecrop = filepdf.replace(".pdf","-crop.pdf")
name = "cubic_circuit.svg"
!pdflatex  $file > /dev/null 2>&1
!pdfcrop $filepdf > /dev/null 2>&1
!pdf2svg $filecrop $name

Here is a graphical representation of the circuit. It is always assumed that the input is vacuum in all the modes.

We can now inspect the covariance matrix and vector of means. Note that the vector of means is non-zero since we did use displacement gates in the circuit above.

print(np.round(mu,10))
print(np.round(cov,10))

[-0.50047867  0.37373598  0.01421683  0.26999427  0.04450994  0.01903583]
[[ 1.57884241  0.81035494  1.03468307  1.14908791  0.09179507 -0.11893174]
 [ 0.81035494  1.06942863  0.89359234  0.20145142  0.16202296  0.4578259 ]
 [ 1.03468307  0.89359234  1.87560498  0.16915661  1.0836528  -0.09405278]
 [ 1.14908791  0.20145142  0.16915661  2.37765137 -0.93543385 -0.6544286 ]
 [ 0.09179507  0.16202296  1.0836528  -0.93543385  2.78903152 -0.76519088]
 [-0.11893174  0.4578259  -0.09405278 -0.6544286  -0.76519088  1.51724222]]

We now use The Walrus to obtain the Fock representation of the heralded Gaussian state when mode 1 is heralded in the value \(n=1\) in the variable psi. We also calculate the probability of success in heralding in the variable p_psi.

cutoff = 10
psi = state_vector(mu, cov, post_select={0: m1, 1: m2}, normalize=False, cutoff=cutoff)
p_psi = np.linalg.norm(psi)
psi = psi/p_psi
print("The probability of successful heralding is ", np.round(p_psi**2,5))

The probability of successful heralding is  0.02016

We now plot the photon-number distribution of the heralded state. Note that the state has zero support on the Fock state with \(i=2\) photons.

plt.bar(np.arange(cutoff),np.abs(psi)**2)
plt.xlabel("$i$")
plt.ylabel(r"$p_i$")
plt.show()

We can now plot the Wigner function of the heralded state,

grid = 100
xvec = np.linspace(-5,5,grid)
Wp = wigner(Qobj(psi), xvec, xvec)
wmap = wigner_cmap(Wp)
sc1 = np.max(Wp)
nrm = mpl.colors.Normalize(-sc1, sc1)
fig, axes = plt.subplots(1, 1, figsize=(5, 4))
plt1 = axes.contourf(xvec, xvec, Wp, 60,  cmap=cm.RdBu, norm=nrm)
axes.contour(xvec, xvec, Wp, 60,  cmap=cm.RdBu, norm=nrm)
axes.set_title("Wigner function of the heralded state");
cb1 = fig.colorbar(plt1, ax=axes)
fig.tight_layout()
plt.show()

and a cut of the Wigner function along \(x=0\).

plt.plot(xvec, Wp[:,grid//2])
plt.title(r"$W(0,p)$")
plt.xlabel(r"p")
plt.show()

Adding Loss

We can now study what happens when loss in the heralding arm is increased. We will add loss in the heralding arm by varying the efficiency of the detector which is parametrized by a transmission \(\eta\) ranging from \(\eta = 50\%\) to \(\eta = 100\%\) (ideal operation).

eta_vals = np.arange(1.,0.45,-0.05)
fidelities = np.zeros_like(eta_vals)
success_p = np.zeros_like(eta_vals)
nmodes = 3
for i,eta in enumerate(eta_vals):
    prog = sf.Program(nmodes)
    eng = sf.Engine("gaussian")
    with prog.context as q:
        for k in range(3):
            Sgate(sq_r[k], sq_phi[k]) | q[k]
            Dgate(d_r[k]) | q[k]

        BSgate(bs_theta1, bs_phi1) | (q[0], q[1])
        BSgate(bs_theta2, bs_phi2) | (q[1], q[2])
        BSgate(bs_theta3, bs_phi3) | (q[0], q[1])
        LossChannel(eta)|q[0]
        LossChannel(eta)|q[1]
    state = eng.run(prog).state
    mu = state.means()
    cov = state.cov()
    rho = density_matrix(mu, cov, post_select={0: m1, 1: m2}, normalize=False, cutoff=cutoff)
    success_p[i] = np.real_if_close(np.trace(rho))
    fidelities[i] = np.real_if_close(psi.conj() @ rho @ psi/success_p[i])

We now plot the probability of success of the heralding scheme as a function of the transmission,

plt.plot(eta_vals, success_p)
plt.xlabel(r"Transmission $\eta$")
plt.ylabel(r"Success probability")
plt.show()

and similarly study the fidelity of the heralded state with respect to the ideal (lossless) state.

plt.plot(eta_vals, fidelities)
plt.xlabel(r"Transmission $\eta$")
plt.ylabel(r"Fidelity with ideal state")
plt.show()

We can also look at the photon number distribution of the nonideal state. Notice that now there is a two photon component.

plt.bar(np.arange(cutoff),np.real_if_close(np.diag(rho/np.trace(rho))))
plt.xlabel(r"$i$")
plt.ylabel(r"$p_i$")
plt.show()

Now we plot the Wigner function of the heralded state for \(\eta = 50\%\).

sc1 = np.max(Wp)
W = wigner(Qobj(rho/np.trace(rho)), xvec, xvec)
nrm = mpl.colors.Normalize(-sc1, sc1)
fig, axes = plt.subplots(1, 2, figsize=(7.5, 3),sharey=True)
axes[0].contourf(xvec, xvec, Wp, 60,cmap=cm.RdBu, norm=nrm)
axes[0].contour(xvec, xvec, Wp, 60,  cmap=cm.RdBu, norm=nrm)
axes[0].set_title(r"$W(x,p)$ for $\eta = 1.0$");
cb1 = fig.colorbar(plt1, ax=axes[0])
plt2 = axes[1].contourf(xvec, xvec, W ,60, cmap=cm.RdBu, norm=nrm)
plt2 = axes[1].contour(xvec, xvec, W ,60, cmap=cm.RdBu, norm=nrm)
axes[1].set_title(r"$W(x,p)$ for $\eta = 0.5$");
cb2 = fig.colorbar(plt1, ax=axes[1])
fig.tight_layout()
plt.show()

and also a cut of the Wigner function along \(x=0\)

plt.plot(xvec, W[:,grid//2], label=r"$\eta=0.5$")
plt.plot(xvec, Wp[:,grid//2], label=r"$\eta=1.0$")
plt.title(r"$W(0,p)$")
plt.xlabel(r"$p$")
plt.legend()
plt.show()

%reload_ext version_information
%version_information qutip, strawberryfields, thewalrus

Software	Version
Python	3.7.5 64bit [GCC 7.3.0]
IPython	7.10.1
OS	Linux 4.15.0 72 generic x86_64 with debian stretch sid
qutip	4.4.1
strawberryfields	0.13.0-dev
thewalrus	0.11.0-dev
Mon Dec 30 10:14:55 2019 EST

Note

Click here to download this gallery page as an interactive Jupyter notebook.

gallery/cubic

 

Download Python script

 

Download Notebook

 

View on GitHub