r"""
Quantum information with polarized light
=================================================
*Author: Nicolás Quesada*
In this tutorial we will explore how to model quantum photonic
experiments involving "internal" degrees of freedom of light, such as
polarization. Note that by default Strawberry Fields deals with
qumodes without an explicit commitment to what degree of freedom those
qumodes correspond to. Using polarization degrees of freedom allows us
to encode quantum information more densely. For example, one can use
the two orthogonal polarizations of a single photon to encode a single qubit.
Polarization is the ability of fields to oscillate in more than one
direction as they propagate. One can think of light as a transverse
field, meaning that it does not oscillate in the direction of
propagation. This implies that in our threedimensional world light has
two directions in which it can oscillate that are perpendicular to its direction of propagation.
Because of this, a photon propagating in a certain direction :math:`k`
can have two orthogonal polarizations. This is schematically shown below:
.. image:: /tutorials/images/polarization.svg
:align: center
:width: 75%
:target: javascript:void(0);

In this tutorial, we will show how polarization "gates" implemented by halfwave plates
and polarizing beamsplitters can be mapped into the usual quantum
operations from Strawberry Fields. Then we will use this knowledge to simulate the generation
of a maximally entangled Bell state encoded in polarization degrees of
freedom.
Manipulating polarization degrees of freedom

It turns out that the polarization degree of freedom of photons can be used
to encode information, and this information is typically manipulated
using socalled halfwave plates and polarizing beamsplitters.
Halfwave plates
~~~~~~~~~~~~~~~~
A halfwave plate (HWP) allows us to rotate the polarization of a photon. For
example,
.. math::
\text{HWP} 1_H 0_V \rangle &= \frac{1}{\sqrt{2}} \left( 1_H 0_V \rangle + 0_H 1_V \rangle \right), \quad \text{(A)}\\
\text{HWP} 0_H 1_V \rangle &= \frac{1}{\sqrt{2}} \left( 1_H 0_V \rangle  0_H 1_V \rangle \right), \quad \text{(B)}
where :math:`\text{HWP}` is a unitary operator representing
a halfwave plate, :math:`1_H \rangle` is a single horizontallypolarized photon
and, similarly, :math:`1_V \rangle` is a single verticallypolarized photon
The states in the right hand sides of the last
equation can be understood as single photons polarized along the diagonal
(:math:`+45^\circ`) and antidiagonal (:math:`45^\circ`) directions, as depicted below:
.. image:: /tutorials/images/HWP.svg
:align: center
:width: 75%
:target: javascript:void(0);

The transformation above is written in the Schrödinger picture for
single photon inputs. It is often more fruitful to write it in the
Heisenberg picture as follows:
.. math::
\text{HWP} ^\dagger a_H^\dagger \text{HWP} &= \frac{1}{\sqrt{2}}\left( a_H^\dagger + a_V^\dagger \right),\\
\text{HWP} ^\dagger a_V^\dagger \text{HWP} &= \frac{1}{\sqrt{2}}\left( a_H^\dagger  a_V^\dagger \right),
where we transformed the creation operators :math:`a_{H/V}^\dagger` of the modes.
Polarizing beamsplitters
~~~~~~~~~~~~~~~~~~~~~~~~
A polarizing beamsplitter (PBS) will transmit light of a certain
polarization while reflecting light from the orthogonal polarization.
This operation acts on two paths each having two polarizations. For two
spatial modes with two polarizations (a total of :math:`4 = 2 \times 2`
creation operators) we can write its Heisenberg action as
.. math::
\text{PBS}^\dagger a_{1,H}^\dagger \text{PBS} &= a_{1,H}^\dagger, \quad \text{(a)}\\
\text{PBS}^\dagger a_{1,V}^\dagger \text{PBS} &= a_{2,V}^\dagger, \quad \text{(b)}\\
\text{PBS}^\dagger a_{2,H}^\dagger \text{PBS} &= a_{2,H}^\dagger, \quad \text{(c)}\\
\text{PBS}^\dagger a_{2,V}^\dagger \text{PBS} &= a_{1,V}^\dagger. \quad \text{(d)}
Note that the horizontal polarizations of both modes are
unaffected (they are "transmitted") and that the vertical polarizations
are swapped (they are "reflected"). The action of the PBS is schematically shown below
for the different path and polarization modes:
.. image:: /tutorials/images/PBS.svg
:align: center
:width: 75%
:target: javascript:void(0);

We can use a PBS to separate two photons with orthogonal polarizations
in the same path into two photons in different paths:
.. math::
\text{PBS}1_H 1_V, 0_H 0_V \rangle = 1_H 0_V, 0_H 1_V \rangle,
where we used the notation :math:`x_H y_V,a_H b_V \rangle`
to indicate :math:`x` (:math:`y`) horizontal (vertical) photons in mode
1, :math:`a` (:math:`b`) horizontal (vertical) photons in mode 2.
Mapping to Strawberry Fields gates

With the notation we developed above we are ready to simulate basic
polarization optical primitives using Strawberry Fields. We simply need
to recall that each optical path with two polarizations can be
represented using 2 qumodes.
For each path :math:`i`, we associate qumodes :math:`2i` and :math:`2i+1` to its horizontal and vertical polarization respectively.
We can the associate Strawberry Fields gates with the polarization primitives discussed
above as follows:
.. math::
\text{HWP}_i &\leftrightarrow \text{BSgate}_{2i,2i+1}\\
\text{PBS}_{i,j} &\leftrightarrow \text{Permute}_{2i+1,2j+1} = \text{Interferometer}\left(\left[\begin{smallmatrix} 0&1 \\1&0 \end{smallmatrix} \right]\right)_{2i+1,2j+1}
To see why the first association is true recall that a
HWP maps the creation operators of two polarizations in the same path to
linear combinations of the same operators, which is precisely what the
``BSgate`` does when applied to modes :math:`2i` and :math:`2i+1`
representing the two orthogonal polarizations of path :math:`i`.
To understand the second expression recall that a
:math:`\text{PBS}` simply swaps the vertical polarization of the two
paths, which correspond precisely to qumodes :math:`2i+1` and
:math:`2j+1`. Finally, note that a swap is readily implemented in Strawberry Fields by
using the ``Interferometer`` gate with argument
:math:`\left[\begin{smallmatrix} 0&1 \\1&0 \end{smallmatrix} \right]`.
Simulating the generation of eventready photon pairs

Zhang et al. [[#zhang]_] propose a method for generating postselected photonic maximally
entangled Bell states by using single photons input into an
interferometer. Their method uses a 4 spatialpath interferometer where
each path takes advantage of the two polarizations (vertical and
horizontal) of a photon. The circuit starts with four horizontal photons
in each of the paths and then proceeds as shown below
.. image:: /tutorials/images/circuit_path.svg
:align: center
:width: 75%
:target: javascript:void(0);

Conditioned in two of the detectors collecting a single photon each and
the other two measuring vacuum, Zhang et al. show that the state of
modes 1 and 4 collapses to a twophoton maximally entangled Bell
state of the form
.. math::
\phi^+ \rangle \propto \left( 1_H 0_V, 1_H 0_V \rangle + 0_H 1_V, 0_H 1_V \rangle \right).
With the rules established in the last section we can translate the
2polarization 4path circuit into a circuit with 8 qumodes as shown below
.. image:: /tutorials/images/circuit.svg
:align: center
:width: 75%
:target: javascript:void(0);

We are now ready to do the simulations!
"""
# Import and preliminaries
import strawberryfields as sf
from strawberryfields.ops import Ket, BSgate, Interferometer
import numpy as np
cutoff_dim = 5 # (1+ total number of photons)
paths = 4
modes = 2 * paths
initial_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
# The ket below corresponds to a single horizontal photon in each of the modes
initial_state[1, 0, 1, 0, 1, 0, 1, 0] = 1
# Permutation matrix
X = np.array([[0, 1], [1, 0]])
# Here is the main program
# We create the input state and then send it through a network of beamsplitters and swaps.
prog = sf.Program(8)
with prog.context as q:
Ket(initial_state)  q # Initial state preparation
for i in range(paths):
BSgate()  (q[2 * i], q[2 * i + 1]) # First layer of beamsplitters
Interferometer(X)  (q[1], q[3])
Interferometer(X)  (q[5], q[7])
BSgate()  (q[2], q[3])
BSgate()  (q[4], q[5])
Interferometer(X)  (q[3], q[5])
BSgate().H  (q[2], q[3])
BSgate().H  (q[4], q[5])
# We run the simulation
eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff_dim})
result = eng.run(prog)
state = result.state
ket = state.ket()
# Check the normalization of the ket.
# This does give the exact answer because of the cutoff we chose.
print("The norm of the ket is ", np.linalg.norm(ket))
######################################################################
# Postselection
# 
#
# Let's consider the case where one horizontallypolarized photon is detected in both paths
# 2 and 3.
#
sub_ket1 = np.round(ket[:, :, 1, 0, 1, 0, :, :], 14) # postselect on correct pattern
p1 = np.round(np.linalg.norm(sub_ket1) ** 2, 14) # Check the probability of this event
print("The probability is ", p1)
print("The expected probability is ", 1 / 32)
# These are the only nonzero components
ind1 = np.array(np.nonzero(np.real_if_close(sub_ket1))).T
print("The indices of the nonzero components are \n ", ind1)
# And these are their coefficients
print("The nonzero components have values ", [sub_ket1[tuple(ind)] for ind in ind1])
######################################################################
# Thus up to normalization the postselected state is indeed :math:`\phi^+ \rangle`.
#
######################################################################
# We can study all the successful postselections. To simplify tensor
# manipulation we will move the modes in which we measure to be the first
# 4 modes of the tensor by using a transposition:
#
# Transpose the ket
ket_t = ket.transpose(2, 3, 4, 5, 0, 1, 6, 7)
# Postselection patterns:
patterns = [
[1, 1, 0, 0],
[1, 0, 1, 0],
[1, 0, 0, 1],
[0, 1, 1, 0],
[0, 1, 0, 1],
[0, 0, 1, 1],
]
######################################################################
# For each pattern we can construct the postselected ket, and find which
# components are nonzero. Note that for each postselection there are only
# two nonzero components, as expected for a Bell state.
#
sub_kets = [np.round(ket_t[tuple(ind)], 15) for ind in patterns]
ps = np.array(list(map(np.linalg.norm, sub_kets))) ** 2
indices = np.array([np.array(np.nonzero(sub_ket)).T for sub_ket in sub_kets])
print(
"The indices of the nonzero components for the six different postselections are \n",
indices,
)
# The successful postselection events occur with the same probability
print("The success probabilities for each pattern are the same \n", ps)
######################################################################
# Conclusion
# 
#
# We have examined how to map the evolution of photonic systems with path
# and polarization degrees of freedom into qumodes evolving under unitary
# operations. The main takeaway is that a system with :math:`N` paths and
# 2 polarization degrees of freedom can be mapped into a system of
# :math:`2N` qumodes. We have also explored in detail how the typical
# optical elements used to couple path and polarization can be mapped to
# qumode simulations. Finally, we used these identifications to
# simulate the generation of eventready (i.e. postselected) Bell states.
#
######################################################################
# References
# 
#
# .. [#zhang]
#
# Q. Zhang, X.H. Bao, C.Y. Lu, X.Q. Zhou, T. Yang, T. Rudolph, and J.W. Pan
# Physical Review A 77, 062316, 2008. doi:10.1103/PhysRevA.77.062316 .