"""
Studying realistic bosonic qubits using the bosonic backend
===========================================================
.. note::
This tutorial is third in a series on the ``bosonic`` backend. For an introduction
to the backend, see :doc:`part 1 `. To learn about
how the backend can be used to sample non-Gaussian states, see
:doc:`part 2 `. These tutorials accompany
our research paper [[#bourassa2021]_].
In a previous :doc:`tutorial `, we briefly introduced
Gottesman-Kitaev-Preskill (GKP) states [[#gottesman2001]_] as an example of
non-Gaussian states which can be simulated using the ``bosonic``
backend. Here, we use the ``bosonic`` backend to take a deep dive into
the advantages of using GKP states as a means for encoding a qubit in a
photonic mode. GKP qubits have a key property that we will
explore in this tutorial: a universal set of qubit gates and measurements can be
applied using the already-familiar Gaussian gates and measurements from
previous :doc:`tutorials `.
After providing more intuition for what ideal and practical GKP states
look like in phase space, we will show how to apply qubit gates and
measurements, and how such operations perform under realistic, noisy
conditions. While some of these simulations have been performed before,
many of these, especially the two-qubit Clifford gates and the T-gate
teleportation circuits have eluded simulation because of their computational
cost. The ``bosonic`` backend brings these simulations to within our reach.
We assume that readers of this tutorial have some familiarity with phase
space methods like the :doc:`Wigner function ` and
`basic CV gates `_,
as well as `basic qubit gates `_.
"""
######################################################################
# GKP states as qubits
# --------------------
#
######################################################################
# In an ideal world
# ~~~~~~~~~~~~~~~~~
#
######################################################################
# For the GKP encoding, the Wigner function of the ideal qubit
# :math:`|0\rangle_{\rm gkp}` state consists of a linear combination of
# Dirac delta functions centred at half-integer multiples of
# :math:`\sqrt{\pi\hbar}` in phase space [[#gottesman2001]_]:
#
# .. math:: W^0_{\rm gkp}(q,p) = \sum_{s,t = -\infty}^{\infty} (-1)^{st}\delta\left(p-\frac{s\sqrt{\pi\hbar}}{2}\right)\delta\left(q-t\sqrt{\pi\hbar}\right).
#
# The GKP qubit :math:`|1\rangle_{\rm gkp}` state can be obtained by
# shifting the :math:`|0\rangle_{\rm gkp}` state by
# :math:`\sqrt{\pi\hbar}` in the position quadrature [[#gottesman2001]_].
#
# Below is a visualization of the ideal :math:`|0\rangle_{\rm gkp}` and :math:`|1\rangle_{\rm gkp}`
# states in phase space using the ``bosonic`` backend [[#footnote1]_]. Each blue (red) dot
# represents a positive (negative) delta function.
#
# Import relevant libraries
import strawberryfields as sf
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors, colorbar
# Set the scale for phase space
sf.hbar = 1
scale = np.sqrt(sf.hbar * np.pi)
# Create a GKP |0> state
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP() | q
eng = sf.Engine("bosonic")
gkp_0 = eng.run(prog).state
# Create a GKP |1> state
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP() | q
sf.ops.Xgate(np.sqrt(np.pi * sf.hbar)) | q
eng = sf.Engine("bosonic")
gkp_1 = eng.run(prog).state
# Get the phase space coordinates of the delta functions for the two states
q_coords_0 = gkp_0.means().real[:, 0]
p_coords_0 = gkp_0.means().real[:, 1]
q_coords_1 = gkp_1.means().real[:, 0]
p_coords_1 = gkp_1.means().real[:, 1]
# Determine whether the delta functions are positively or negatively weighted
delta_sign_0 = np.sign(gkp_0.weights().real)
delta_sign_1 = np.sign(gkp_1.weights().real)
# Plot the locations and signs of the deltas
fig,ax = plt.subplots(1, 2)
ax[0].scatter(q_coords_0 / scale,
p_coords_0 / scale,
c=delta_sign_0,
cmap=plt.cm.RdBu, vmin=-1.5, vmax=1.5)
ax[1].scatter(q_coords_1 / scale,
p_coords_1 / scale,
c=delta_sign_0,
cmap=plt.cm.RdBu, vmin=-1.5, vmax=1.5)
for i in range(2):
ax[i].set_xlim(-3.5, 3.5)
ax[i].set_ylim(-3.5, 3.5)
ax[i].set_xlabel(r'$q$ (units of $\sqrt{\pi\hbar}$ )', fontsize=15)
ax[i].set_aspect("equal")
ax[0].set_title(r'$|0\rangle_{\rm gkp}$ Wigner function', fontsize=15)
ax[1].set_title(r'$|1\rangle_{\rm gkp}$ Wigner function', fontsize=15)
ax[0].set_ylabel(r'$p$ (units of $\sqrt{\pi\hbar}$ )', fontsize=15)
fig.tight_layout()
plt.show()
######################################################################
# Like standard qubit
# :math:`|0\rangle` and :math:`|1\rangle` states, these two states are orthogonal:
# wherever the :math:`|0\rangle_{\rm gkp}` state has a column of
# delta functions with alternating signs, the :math:`|1\rangle_{\rm gkp}` state
# has a column of all positive delta functions, and vice versa. This
# results in a zero overlap between the :math:`|0\rangle_{\rm gkp}` and
# :math:`|1\rangle_{\rm gkp}` Wigner functions, meaning they are orthogonal.
#
# The Wigner function for a general qubit state
#
# .. math:: |\psi\rangle = \cos{\tfrac{\theta}{2}} \vert 0 \rangle_{\rm gkp} + e^{-i \phi} \sin{\tfrac{\theta}{2}} \vert 1 \rangle_{\rm gkp}
#
# also consists of a weighted linear combination of delta functions
# at half-integer multiples of :math:`\sqrt{\pi\hbar}` in phase space,
# with the weights depending :math:`\theta` and :math:`\phi` [[#alvarez2020]_, [#bourassa2021]_].
#
# One might have already noticed a problem:
# since their Wigner functions consist of a linear combination of
# deltas, ideal GKP states have infinite energy, cannot be normalized, and are
# therefore unphysical!
#
######################################################################
# Let's get physical
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
######################################################################
# Phyiscal, normalizable and finite-energy GKP states can be constructed
# by applying a Fock damping operator to the ideal states:
#
# .. math:: |k^\epsilon\rangle_{\rm gkp} = e^{-\epsilon \hat{n}} |k\rangle_{\rm gkp}, k=0,1,
#
# where :math:`\hat{n}` is the photon number operator. For
# any realistic application calling for GKP states, finite-energy GKP
# states will need to be employed.
#
# Through a mathematical derivation best enjoyed with coffee â and a snack đȘ
# [[#bourassa2021]_], we find the effect of the finite-energy operator on
# the GKP Wigner functions to be as follows:
#
# 1. Each delta function is replaced with a Gaussian function of non-zero variance.
#
# 2. The location of each of these new Gaussian peaks is contracted towards the origin.
#
# 3. The weights of the Gaussians are damped the further the peak is located from the origin.
#
# As we saw in the :doc:`introductory tutorial `, Wigner functions that can be
# expressed as a linear combination of Gaussian functions in phase space
# are the bread and butter of the ``bosonic`` backend. Now we have shown
# exactly how GKP states fall neatly into that form.
#
# For these new sets of Gaussian peaks, the variance, the location shift,
# and the damping of the weights all depend on the finite-energy parameter :math:`\epsilon`. Let's
# visualize how the Wigner function changes as we vary :math:`\epsilon`:
#
# Choose some values of epsilon
epsilons = [0.05, 0.1, 0.2, 5]
# Pick a region of phase space to plot
quad = np.linspace(-3.5, 3.5, 200) * scale
wigners = []
for epsilon in epsilons:
# Create a GKP |0> state
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(epsilon=epsilon) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
# Calculate the Wigner function
wigner = gkp.wigner(mode=0, xvec=quad, pvec=quad)
wigners.append(wigner)
# Plot the results
fig, axs = plt.subplots(1, 5, figsize=(16, 4), gridspec_kw={"width_ratios": [1, 1, 1, 1, 0.05]})
cmap = plt.cm.RdBu
cmax = np.real_if_close(np.amax(np.array(wigners)))
norm = colors.Normalize(vmin=-cmax, vmax=cmax)
cb1 = colorbar.ColorbarBase(axs[4], cmap=cmap, norm=norm, orientation="vertical")
for i in range(4):
axs[i].contourf(
quad / scale,
quad / scale,
wigners[i],
levels=60,
cmap=plt.cm.RdBu,
vmin=-cmax,
vmax=cmax,
)
axs[i].set_title(r'$\epsilon =$'+str(epsilons[i]), fontsize=15)
axs[i].set_xlabel(r'$q$ (units of $\sqrt{\pi\hbar}$ )', fontsize=15)
axs[0].set_ylabel(r'$p$ (units of $\sqrt{\pi\hbar}$ )', fontsize=15)
cb1.set_label("Wigner function", fontsize=15)
fig.tight_layout()
plt.show()
######################################################################
# We see that as :math:`\epsilon` increases, the variance of each peak
# grows, the peaks get closer and closer to the origin, and the weights
# (which govern how large the Wigner function gets) drop exponentially
# away from the origin. In the limit as :math:`\epsilon\rightarrow\infty`,
# the Fock damping gets so strong that all we are left with is the vacuum
# state, as can be seen in the last plot.
#
# Since the delta spikes are now Gaussians, the
# :math:`|0^\epsilon\rangle_{\rm gkp}` and
# :math:`|1^\epsilon\rangle_{\rm gkp}` states are no longer exactly
# orthogonal, which will lead to a small but inherent readout error
# when distinguishing the computational-basis states.
# As compared to ideal GKP qubits, finite-energy qubits are noisier;
# consequently, we can refer to this type of imperfection as "finite-energy noise".
#
# Now that we've established that realistic GKP states fall under the
# simulation capabilities of the ``bosonic`` backend, let's go through the
# set of qubit measurements and gates for a GKP state. We'll see that the
# only resources one needs in addition to GKP states are Gaussian
# transformations and measurements [[#baragiola2019]_], which we also know to be easily
# simulated using the ``bosonic`` backend.
#
######################################################################
# Qubit Pauli measurements: CV homodyne measurements
# ---------------------------------------------------
#
######################################################################
# The first example of correspondence between GKP qubit operations and
# Gaussian resources are Pauli measurements. For GKP qubits, they can be
# implemented by performing homodyne measurements and processing the
# outcomes. The association is:
#
# ===== == ============
# Pauli Homodyne
# ===== == ============
# Z :math:`q`
# X :math:`p`
# Y :math:`q-p`
# ===== == ============
#
# The outcomes of these homodyne measurements are rounded or "binned" to
# the nearest :math:`n\sqrt{\pi\hbar}`, :math:`n\in\mathbb{Z}`. Then, the
# parity of :math:`n` is taken to determine the value of the Pauli
# measurement. Note that the Pauli Y measurements can be achieved by
# performing a homodyne measurement along :math:`(q - p)/\sqrt{2}` and rescaling by :math:`\sqrt{2}`.
#
# As an example, we look at the marginal distributions of homodyne outcomes
# for GKP :math:`|+^\epsilon\rangle_{\rm gkp}` to confirm Pauli measurements can be performed with homodyne.
# Alternatively, we could have directly simulated samples like in the :doc:`last tutorial `.
#
# The angles :math:`\theta,\phi` that specify the qubit state can be passed to the
# ``sf.ops.GKP()`` operation via the ``state`` parameter.
#
# The homodyne outcomes that fall into blue (red) bins can be interpreted as a Pauli operator
# eigenvalue measurement with outcome :math:`+1` (:math:`-1`).
# The bins can be used to provide the expectation value of each Pauli operator:
# we integrate the marginal distributions, each multiplied by an appropriate weighting function
# that changes the sign of the distribution depending on the Pauli operator eigenvalue associated
# with the bin.
#
# Create a GKP |+> state
plus = [np.pi / 2, 0]
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.08) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
def calc_and_plot_marginals(state, mode):
'''Calculates and plots the p, q-p, and q quadrature marginal
distributions for a given circuit mode. These can be used to determine
the Pauli X, Y, and Z outcomes for a GKP qubit.
Args:
state (object): a strawberryfields ``BaseBosonicState`` object
mode (int): index for the circuit mode
'''
# Calculate the marginal distributions.
# The rotation angle in phase space is specified by phi
marginals = []
phis = [np.pi / 2, -np.pi / 4, 0]
quad = np.linspace(-5, 5, 400) * scale
for phi in phis:
marginals.append(gkp.marginal(mode, quad, phi=phi))
# Plot the results
paulis = ["X", "Y", "Z"]
homodynes = ["p", "q-p", "q"]
expectations = np.zeros(3)
fig, axs = plt.subplots(1, 3, figsize=(12, 4))
for i in range(3):
if i == 1:
# Rescale the outcomes for Pauli Y
y_scale = np.sqrt(2 * sf.hbar) / scale
axs[i].plot(quad * y_scale, marginals[i] / y_scale, 'k-')
axs[i].set_xlim(quad[0] * y_scale, quad[-1] * y_scale)
# Calculate Pauli expectation value
# Blue bins are weighted +1, red bins are weighted -1
bin_weights = 2 * (((quad * y_scale - 0.5) // 1) % 2) - 1
integrand = (marginals[i] / y_scale) * bin_weights
expectations[i] = np.trapz(integrand, quad * y_scale)
else:
axs[i].plot(quad / scale, marginals[i] * scale, 'k-')
axs[i].set_xlim(quad[0] / scale, quad[-1] / scale)
# Calculate Pauli expectation value
# Blue bins are weighted +1, red bins are weighted -1
bin_weights = 2 * (((quad / scale - 0.5) // 1) % 2) - 1
integrand = (marginals[i] * scale) * bin_weights
expectations[i] = np.trapz(integrand, quad / scale)
# Color the qubit bins blue and red
for j in range(-10, 10):
axs[i].axvspan((2 * j - 0.5), (2 * j + 0.5), alpha=0.2, facecolor='b')
axs[i].axvspan((2 * j + 0.5), (2 * j + 1.5), alpha=0.2, facecolor='r')
axs[i].set_title("Homodyne data for Pauli " + paulis[i] +
"\n" + r'$\langle$'+paulis[i]+r'$\rangle$='+
str(np.around(expectations[i],2)))
axs[i].set_xlabel(homodynes[i] + r' (units of $\sqrt{\pi\hbar}$ )', fontsize=15)
axs[0].set_ylabel("Marginal distribution", fontsize=15)
fig.tight_layout()
plt.show()
calc_and_plot_marginals(gkp, 0)
######################################################################
# We can see from the :math:`p`-homodyne measurement outcomes that the
# results are mostly within the blue bins, so there will be a high chance
# of reading out a :math:`+1` eigenvalue from a Pauli X measurement, as would be
# expected for a standard qubit :math:`|+\rangle` state.
#
# For the other two homodyne measurements, the peaks of the marginal
# distributions are effectively evenly distributed between :math:`+1` and :math:`-1` bins,
# just like Pauli Y and Z measurements on a standard qubit
# :math:`|+\rangle` state.
#
# One of the strengths of the bosonic backend is the ease of applying
# Gaussian transformations to states. For qubits encoded in CV systems, a
# major source of noise is loss or dissipation. Let's see how the
# distributions of homodyne outcomes change after we apply the
# :class:`~strawberryfields.ops.LossChannel` operation:
#
# Create a GKP |+> state and apply loss
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.08) | q
sf.ops.LossChannel(0.85) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
#Calculate and plot marginals
calc_and_plot_marginals(gkp, 0)
######################################################################
# We see that the peaks in the homodyne distribution get broadened and
# shifted towards the origin, resulting in outcomes falling outside of the
# correct bins. This corresponds to qubit readout errors.
#
# From this, we see the ``bosonic`` backend can be used to study noise and develop
# strategies to deal with realistic effects on bosonic qubits.
#
######################################################################
# Single-qubit Clifford gates
# ---------------------------
#
######################################################################
# With Pauli measurements under our belt, we can move on to the
# single-qubit Clifford gates. As a summary, the association between
# GKP single-qubit Clifford gates and qubits is [[#footnote2]_]:
#
# ============ == ==================================================
# Qubit gates CV gates
# ============ == ==================================================
# Pauli X :math:`q`-displacement by :math:`\sqrt{\pi\hbar}`
# Pauli Z :math:`p`-displacement by :math:`\sqrt{\pi\hbar}`
# Hadamard Rotation by :math:`\pi/2`
# Phase Quadratic phase gate of strength 1
# ============ == ==================================================
#
# We'll demonstrate this mapping in the next few subsections.
#
######################################################################
# Qubit Pauli X and Z gates: CV :math:`q` and :math:`p` displacements
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# The simplest examples are the qubit `Pauli X `_
# and `Z `_ gates,
# which for the GKP encoding correspond to `displacements `_
# in phase space by :math:`\sqrt{\pi\hbar}` along :math:`q` and :math:`p`, respectively.
# These can be implemented using the :class:`~strawberryfields.ops.Xgate` and :class:`~strawberryfields.ops.Zgate` operations.
#
# Below, we take a GKP qubit with :math:`\theta=\phi=np.pi/4` and plot its marginal distributions along
# :math:`p`, :math:`q-p`, and :math:`q` before and after applying a
# bit-phase flip. Recall that, once binned, this data provides the outcome
# of the Pauli measurements: blue (red) bins correspond to Pauli :math:`+1` (:math:`-1`) eigenvalues.
#
# Create a GKP state
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=[np.pi / 4, np.pi / 4], epsilon=0.08) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
#Calculate and plot marginals
calc_and_plot_marginals(gkp, 0)
######################################################################
# Run it again but apply a bit-phase flip this time
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=[np.pi / 4, np.pi / 4], epsilon=0.08) | q
sf.ops.Xgate(np.sqrt(np.pi * sf.hbar)) | q
sf.ops.Zgate(np.sqrt(np.pi * sf.hbar)) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
#Calculate and plot marginals
calc_and_plot_marginals(gkp, 0)
######################################################################
# For :math:`p` and :math:`q` (associated with Pauli X and
# Z), we see that the regions of the distributions originally in
# the :math:`+1` (blue) bins get shifted to :math:`-1` (red) bins and vice versa. The
# distribution of :math:`q-p` remains unchanged. This confirms the gates
# act as expected, since the probabilities of :math:`+1` and :math:`-1` outcomes from Pauli
# X and Z measurements get swapped by the bit-phase flip, while the Pauli
# Y is invariant!
#
######################################################################
# Qubit Hadamard gate: CV rotation
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# The GKP `Hadamard `_
# gate can be implemented by a :math:`\pi/2`
# `rotation `_
# in phase space. We've already indirectly seen this gate in action when
# toggling between different homodyne measurements. In the previous
# example, if a :math:`\pi/2` rotation were applied to the state, the
# distribution of outcomes in :math:`q` and :math:`p` are interchanged,
# which means the the statistics of the Pauli X and Z
# outcomes are swapped as well, just as expected from a qubit Hadamard
# gate!
#
######################################################################
# Qubit phase gate: CV quadratic phase gate
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# The GKP `qubit phase `_
# gate maps to the `CV quadratic phase `_
# gate, which can be implemented in ``Strawberry Fields`` with the existing
# :class:`~strawberryfields.ops.Pgate`. This is a Gaussian operation, but unlike
# rotations, the CV phase gate requires a `squeezing operation `_.
#
# Let's see the GKP phase gate in action on a qubit
# :math:`|+^\epsilon\rangle_{\rm gkp}` state:
#
# Create a GKP |+> state and apply a phase gate
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.08) | q
sf.ops.Pgate(1) | q
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
calc_and_plot_marginals(gkp, 0)
######################################################################
# As expected, the homodyne outcomes associated with the Pauli Y
# measurement now fall entirely in the :math:`+1` (blue) bins, reflecting that a
# standard qubit phase gate converts
# :math:`\hat{S}|+\rangle = |+i\rangle`.
#
# As mentioned, the CV phase gate requires inline squeezing. In photonics,
# a feasible technique to do this is to apply measurement-based squeezing
# [[#filip2005]_], which is essentially a method to teleport a squeezing
# operation onto a state. Briefly, the steps to do this are:
#
# 1. Interfere an ancillary squeezed vacuum state at a beamsplitter with the target state.
#
# 2. Perform a homodyne measurement on the ancillary mode.
#
# 3. Apply a feedforward displacement to the target mode based on the homodyne outcome.
#
# The quality of measurement-based squeezing depends on the level of
# squeezing in the ancillary mode, and the efficiency of the homodyne
# detection.
#
# While all backends of ``Strawberry Fields`` have implemented the ideal
# inline squeezing operation through :class:`~strawberryfields.ops.Sgate`, to be able to study the
# effects of realistic noise on bosonic qubits, for the ``bosonic``
# backend, we have also implemented measurement-based squeezing as
# :class:`~strawberryfields.ops.MSgate`.
#
# :class:`~strawberryfields.ops.MSgate` performs the teleportation circuit described above. The user
# can specify whether they would like the map corresponding to the average
# application of the teleportation circuit, or whether they want to
# simulate a single-shot run and keep the homodyne sample on the ancillary
# mode.
#
# Below we examine how the quality of the GKP qubit phase gate changes
# once we replace the inline squeezing in the CV phase gate with
# measurement-based squeezing. We look first at the average map of the
# gate.
#
# First, we define a phase gate that uses measurement-based squeezing.
# ``r_anc`` and ``eta_anc`` are the level of squeezing and detection efficiency for
# the ancilla. ``avg`` specifies whether to apply the average or single-shot map.
#
def MPgate(q, s, r_anc, eta_anc, avg):
'''Applies a CV phase gate to a given mode using measurement-based squeezing.
Args:
q (object): label for the mode in the strawberryfields program
s (float): the ideal phase gate parameter
r_anc (float): level of squeezing in the ancillary mode
eta_anc (0 <= float <= 1): detector efficiency for the ancillary mode
avg (bool): True (False) applies the average (single-shot) mapping
'''
r = np.arccosh(np.sqrt(1 + s ** 2 / 4))
theta = np.arctan(s / 2)
phi = - np.pi / 2 * np.sign(s) - theta
sf.ops.MSgate(r, phi, r_anc, eta_anc, avg) | q
sf.ops.Rgate(theta) | q
######################################################################
# Next, we create a GKP :math:`|+^\epsilon\rangle_{\rm gkp}` state and
# apply a realistic phase gate (average map):
#
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.08) | q
MPgate(q, 1, r_anc=1.2, eta_anc=0.95, avg=True)
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
calc_and_plot_marginals(gkp, 0)
######################################################################
# What happened? The peaks stayed essentially in the same places as in the
# previous simulation, but now they've become broader! Since more of each
# peak falls outside the correct bins, this will have an impact on
# likelihood of obtaining the correct Pauli measurement outcomes.
#
# Let's take a look at the single-shot implementation of the
# measurement-based squeezing operation in the phase gate by setting ``avg`` to ``False``:
#
prog = sf.Program(1)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.08) | q
MPgate(q, 1, r_anc=1.2, eta_anc=0.95, avg=False)
eng = sf.Engine("bosonic")
results = eng.run(prog)
gkp = results.state
calc_and_plot_marginals(gkp, 0)
######################################################################
# On the face of it, these marginal distributions look very different!
# Luckily, the Pauli outcomes only depend on how much of the marginal
# distributions fall inside the correct bins, and less on the shapes and
# heights of the peaks within those bins. We notice that the homodyne
# distribution used for the Pauli Y measurement still falls primarily in
# the :math:`+1` (blue) bins, although the peaks have broadened as compared to the
# ideal CV phase gate.
#
# As promised, the homodyne outcome used for the feedforward operation can
# be retrieved from the ``Results`` object output by the circuit:
#
print(results.ancilla_samples)
######################################################################
# The dictionary key ``0`` denotes the mode where measurement-based
# squeezing was applied, and the list entry contains the result of the
# homodyne measurement on the ancillary state used in the feedforward
# displacement of the measurement-based squeezing gate.
#
# We've again seen through this example that the ``bosonic`` backend is
# useful for studying realistic noise effects present in photonic circuits
# applied to bosonic qubits. Next we move to two-qubit Clifford gates.
#
######################################################################
# Two-qubit Clifford gates
# -----------------------------------
#
######################################################################
# Now that we understand single-qubit Clifford gates, we can move to
# two-qubit entangling gates. Another boon to the GKP encoding is the fact
# that certain two-qubit entangling gates such as the qubit
# `CNOT `_
# and `CZ `_ gates
# correspond to deterministic Gaussian operations, specifically the CV
# `CX `_
# and `CZ `_ gates:
#
# ============ == =======================
# Qubit gates CV gates
# ============ == =======================
# CNOT CX gate of strength 1
# CZ CZ gate of strength 1
# ============ == =======================
#
# As a simple demonstration of a two-qubit entangling gate, we apply a GKP
# CNOT gate on :math:`|+^\epsilon\rangle_{\rm gkp}` and :math:`|0^\epsilon\rangle_{\rm gkp}`.
# For standard qubits, this yields the entangled state
# :math:`\frac{1}{\sqrt{2}}(|0,0\rangle+|1,1\rangle)`.
#
# As a sanity check that this is indeed the state produced from the CNOT
# gate, we project the first qubit into the :math:`|0\rangle` state by
# measuring :math:`q`-homodyne (corresponding to a qubit Pauli Z
# measurment) and post-select on :math:`q_1=0`, which forces the outcome
# of the first mode into a :math:`+1` (blue) bin. When we plot the marginal
# distribution of the second qubit after the measurement on the first qubit, we see that the
# distribution of :math:`q_2` outcomes fall mainly into the :math:`+1` (blue) bins,
# meaning the Pauli Z outcome on the second qubit will also be :math:`+1`,
# as expected.
#
# Below, we create :math:`|+^\epsilon\rangle_{\rm gkp}\otimes|0^\epsilon\rangle_{\rm gkp}` and apply a CNOT gate.
# Then, we measure :math:`q`-homodyne on the first mode, post-selecting on 0
# via the ``select`` parameter:
#
prog = sf.Program(2)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.1) | q[0]
sf.ops.GKP(epsilon=0.1) | q[1]
sf.ops.CXgate(1) | (q[0],q[1])
sf.ops.MeasureHomodyne(0, select=0) | q[0]
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
calc_and_plot_marginals(gkp, 1)
######################################################################
# We can do the same simulation, but this time postselect
# :math:`q_1=\sqrt{\pi\hbar}`, forcing the
# Pauli Z outcome of the first qubit to be :math:`-1`, since that value falls in a red bin.
# When we plot the marginal
# distributions of the second qubit, we see that the distribution of
# outcomes in :math:`q_2` now fall in the :math:`-1` (red) bins, indicating the
# correct qubit state was created by the two-qubit gate.
#
prog = sf.Program(2)
with prog.context as q:
sf.ops.GKP(state=plus, epsilon=0.1) | q[0]
sf.ops.GKP(epsilon=0.1) | q[1]
sf.ops.CXgate(1) | (q[0],q[1])
sf.ops.MeasureHomodyne(0, select=np.sqrt(np.pi * sf.hbar)) | q[0]
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
calc_and_plot_marginals(gkp, 1)
######################################################################
# The CV `CX `_
# gate used to implement the GKP qubit CNOT gate requires
# beamsplitters, rotations and inline squeezing.
# While we did not examine it here, the same realistic
# noise effects we examined for the other operations -- be it loss or
# measurement-based squeezing -- could be simulated for two qubit gates
# using the ``bosonic`` backend.
#
######################################################################
# Non-Clifford T gate
# -------------------
#
######################################################################
# To complete the set of universal gates for GKP qubits, we need one
# non-Clifford gate. As it turns out, the GKP
# `qubit T `_
# gate can be applied through gate teleportation by using an ancillary GKP qubit magic state
# (a non-Pauli eigenstate) and Gaussian operations [[#gottesman2001]_].
#
# The steps to apply the T gate build on all the previous examples we've
# seen:
#
# 1. Entangle the ancillary magic state and the data state with a qubit CNOT gate, i.e. a CV CX gate.
#
# 2. Perform a Pauli Z measurement on the ancillary mode via a binned CV :math:`q`-homodyne measurement.
#
# 3. Based on the outcome of the Pauli measurement, apply a feedforward qubit phase gate (a CV phase gate) to the data mode.
#
# Ideally, the GKP qubit magic state one would use is given by
# :math:`|M\rangle_{\rm gkp}=\frac{1}{\sqrt{2}}(e^{-i\pi/8}|0\rangle_{\rm gkp} + e^{i\pi/8}|1\rangle_{\rm gkp})`.
# Next, we will see the effect of using a finite-energy GKP
# :math:`|M^{\epsilon}\rangle_{\rm gkp}` for gate teleportation.
#
# We apply the qubit T gate to the state
# :math:`|M^{\epsilon}\rangle_{\rm gkp}` itself, knowing that for
# standard qubits :math:`\hat{T}|M\rangle=|+i\rangle`. For this
# simulation, we will postselect on a homodyne outcome for the ancillary
# mode that will cause the feedforward phase gate to be applied.
#
# Define a binning function to get the Pauli Z outcome for the feedforward
def gkp_binning(x):
'''Bins a homodyne outcome to a binary 0 or 1 according to the GKP bin structure.
Pauli +1 (-1) corresponds to 0 (1).
Args:
x (float): homodyne outcome value
'''
term_1 = (x + np.sqrt(np.pi) / 2) / np.sqrt(4 * np.pi)
term_2 = (x + np.sqrt(np.pi) / 2) // np.sqrt(4 * np.pi)
return (term_1 - term_2) // 0.5
# Create GKP magic states and perform the T gate teleportation circuit
magic = [np.pi / 2, - np.pi / 4]
prog = sf.Program(2)
with prog.context as q:
sf.ops.GKP(state=magic, epsilon = 0.1) | q[0]
sf.ops.GKP(state=magic, epsilon = 0.1) | q[1]
sf.ops.CXgate(1) | (q[0], q[1])
sf.ops.MeasureHomodyne(0, select=np.sqrt(np.pi * sf.hbar)) | q[1]
sf.ops.Pgate(gkp_binning(q[1].par)) | q[0]
eng = sf.Engine("bosonic")
gkp = eng.run(prog).state
calc_and_plot_marginals(gkp, 0)
######################################################################
# We find that the output homodyne distributions match what we expect
# for the :math:`|+i\rangle` qubit state: the
# distribution along :math:`q-p` falls mainly within :math:`+1` (blue) bins, while
# the other quadratures fall evenly between red and blue bins.
#
# An interesting effect to note is that the peaks are broader in some
# quadratures than others. This demonstrates that through the gate
# teleportation process, the data qubit picked up some of the
# finite-energy noise from the ancillary magic state. These types of realistic
# effects are easy to study using the ``bosonic`` backend!
#
######################################################################
# Conclusion
# ----------
#
######################################################################
# In this tutorial, we showed the translation between the universal set of
# GKP qubit gates and their CV Gaussian phase space implementations. We
# performed simulations using the ``bosonic`` backend, a tool specialized to
# study how bosonic qubits transform under the types of Gaussian
# transformations and measurements necessary for the GKP encoding.
#
######################################################################
# References and Footnotes
# ----------
#
# .. [#bourassa2021]
#
# J. Eli Bourassa, NicolĂĄs Quesada, Ilan Tzitrin, Antal SzĂĄva, Theodor Isacsson,
# Josh Izaac, Krishna Kumar Sabapathy, Guillaume Dauphinais, and Ish Dhand.
# Fast simulation of bosonic qubits via Gaussian functions in phase space.
# `arXiv:2103.05530 `_, 2021.
#
# .. [#gottesman2001]
#
# Daniel Gottesman, Alexei Kitaev, and John Preskill. Encoding a qubit in an oscillator.
# doi:10.1103/PhysRevA.64.012310.
#
# .. [#footnote1]
#
# For the advanced reader: when we call ``sf.ops.GKP()``, we are still only creating a
# finite-energy GKP state. However, the locations of peaks in that state are close enough
# to the ideal state, and the signs of the weights match the ideal state, so we can use it
# for a pedagogical illustration of the ideal GKP states.
#
# .. [#alvarez2020]
#
# L. GarcĂa-Ălvarez, A. Ferraro, and G. Ferrini. From the Bloch Sphere to Phase-Space
# Representations with the GottesmanâKitaevâPreskill Encoding. doi:10.1007/978-981-15-5191-8_9.
#
# .. [#baragiola2019]
#
# Ben Q. Baragiola, Giacomo Pantaleoni, Rafael N. Alexander, Angela Karanjai, and Nicolas C. Menicucci.
# All-Gaussian Universality and Fault Tolerance with the Gottesman-Kitaev-Preskill Code.
# doi:10.1103/PhysRevLett.123.200502.
#
# .. [#footnote2]
#
# A note for the keen reader: An extra advantage of the GKP encoding is that there is in
# fact a one-to-many mapping of qubit to CV gates. Taking Pauli gates as an exmaple, a
# displacement by any odd multiple of :math:`\sqrt{\pi\hbar}` along :math:`q` (:math:`p`) will effect
# a Pauli X (Z) gate. The direction of displacements can be varied in the circuit compilation to
# avoid too much energy being injected into the state, which would happen if the displacements
# were always in the same direction.
#
# .. [#filip2005]
#
# Radim Filip, Petr Marek, and Ulrik L. Andersen. Measurement-induced continuous-variable
# quantum interactions. doi:10.1103/PhysRevA.71.042308.