# Measuring Wigner functions of quantum states of light in the undergraduate
laboratory

Juan Rafael Álvarez<EMAIL_ADDRESS>Laboratorio de Óptica
Cuántica, Universidad de los Andes, A.A. 4976, Bogotá D.C., Colombia.
Université Paris-Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies,
91120, Palaiseau, France.  Clarendon Laboratory, University of Oxford, Parks
Road, Oxford OX1 3PU, United Kingdom  Andrés Martínez Silva  Laboratorio de
Óptica Cuántica, Universidad de los Andes, A.A. 4976, Bogotá D.C., Colombia.
Alejandra Valencia<EMAIL_ADDRESS>Laboratorio de Óptica
Cuántica, Universidad de los Andes, A.A. 4976, Bogotá D.C., Colombia.

###### Abstract

In this work, we present an educational activity aimed at measuring the Wigner
distribution functions of quantum states of light in the undergraduate
laboratory. This project was conceived by students from various courses within
the physics undergraduate curriculum, and its outcomes were used in an
introductory Quantum Optics course at the Universidad de los Andes in Bogotá,
Colombia. The activity entails a two-hour laboratory practice in which
students engage with a pre-aligned experimental setup. They subsequently
employ an open-access, custom-made computational graphical user interface to
reconstruct the Wigner distribution function for various quantum states of
light. Given that the testing phase coincided with the COVID-19 pandemic, we
incorporated the capacity to analyze simulated data into the computational
user interface. The activity is now part of the course syllabus and its
virtual component has proven to be highly valuable for the implementation of
distance learning in quantum optics.

##  I Introduction

Understanding and characterizing the properties of light beams is a key
endeavour in optical physics. Light sources can vary in many ways, having,
among others, different amplitudes, frequencies, phases, and polarizations.
Such variations can be present within the same physical objects and define the
different states that can be assigned to light. In quantum optics, a field of
particular relevance, the measurement of amplitude and phase behaviors in
different states of light plays a critical role. In this context, the Wigner
distribution function (WDF) is a commonly utilized tool for characterizing the
latter behaviors in a phase space [1, 2]. WDFs are typically obtained through
tomographic reconstruction, a technique that retrieves a three-dimensional
image of an object by utilizing partial information from two-dimensional
projections. An established method for performing tomographic reconstructions
is Computed Axial Tomography (CAT), initially introduced in the realm of
medical imaging [3, 4].

In this article, we provide an account of the development and execution of a
laboratory practice designed to reconstruct the WDF for various quantum states
of light, employing a tomographic reconstruction technique. We introduce an
experimental practice that was devised for and by undergraduate physics
students. The development spanned an 18-month period during which both the
experimental and computational facets of the activity were integrated.
Subsequently, the exercise underwent testing by undergraduate and early
graduate students participating in an introductory quantum optics course at
the Universidad de los Andes (Uniandes) in Bogotá, Colombia.

All the computer programming required for this practice was developed in
Python 3, featuring a user-friendly graphical user interface (GUI) that
simplifies its use. The code is freely available and can be accessed at
https://github.com/amartinez1224/quantum-tomography. Given that the testing
phase coincided with the COVID-19 pandemic, the activity was designed to offer
a similar experience for both in-person and virtual experimental sessions,
enhancing accessibility for students participating in virtual classes. Despite
the resumption of in-person classes at the time of writing, this activity
remains valuable, as it continues to serve the needs of distance learning.

We anticipate that this activity can be of valuable use to other educators and
researchers in multiple ways. Firstly, it can serve as pedagogical material to
elucidate the concepts related to the measurement of quadratures, the
understanding of Radon transforms, and their pivotal role in the
reconstruction of Wigner Distribution Functions (WDFs) in quantum optics,
complementing existing references [5, 6, 7, 8]. Secondly, it can function as a
guide for an experimental practice suitable for undergraduate and early
graduate students in quantum optics classes. Lastly, it can act as a versatile
tool for remotely training the quantum workforce through distance learning.

This paper is structured as follows: In Section II, we delineate various tools
for characterizing signals, such as quadratures and Wigner distribution
functions, and elaborate on the methods for their measurement. Section III
offers an overview of diverse quantum states of light, consolidating them in a
comprehensive table. Section IV delves into the details of the experimental
and computational aspects of the educational activity, while Section V reports
on the hands-on and virtual activities developed in collaboration with
undergraduate and early graduate students. Finally, in Section VI, we present
our conclusions and outline potential perspectives.

##  II Tools for state characterization and their measurement

Several conceptual tools are available for describing the properties of light
coming from signals arriving at detectors. This article primarily concentrates
on two of these tools: Quadratures, employed to characterize the behavior of
an electromagnetic wave in terms of conjugate variables, and Wigner
Distribution Functions (WDFs), which can represent signals within a phase
space.

###  II.1 Quadratures

Electric fields associated to electromagnetic waves can be mathematically
described as solutions to the wave equation derived from Maxwell’s equations
[9]. Generally, solutions to the wave equation can be decomposed as the sum of
monochromatic waves with polarization vector $\mathbf{e}$ and frequency
$\omega$. A mathematical approach generally taken to describe one such
monochromatic wave is to use a complex amplitude: the electric field is
considered to be the real part of a vector of complex numbers that evolves in
time, written in the form:

$\mathbf{E}(t)=\operatorname{Re}\left[\tilde{E}_{0}e^{i\omega
t}\right]\mathbf{e},$ (1)

where $\tilde{E}_{0}=E_{0}e^{i\theta}=\left|\tilde{E}_{0}\right|e^{i\theta}$
is a complex number known as the complex amplitude, or phasor, of the field,
and $\theta$ is its initial phase. Another approach to describe the field is
based on the concept of quadratures, where the exponential in Eq. 1 can be
expanded to rewrite $\mathbf{E}(t)$ as:

$\mathbf{E}(t)=\left[X\cos\omega t-Y\sin\omega t\right]\mathbf{e},$ (2)

where $X=\mathrm{Re}\left(\tilde{E}_{0}\right)$ and
$Y=\mathrm{Im}\left(\tilde{E}_{0}\right)$ are two complementary quadratures of
the field.

The quadratures can be thought of as the natural variables for a harmonic
oscillator: In classical mechanics, the Hamiltonian for a harmonic oscillator
with mass $m$ and angular frequency $\omega$ is given by

$H=\frac{{[p(t)]}^{2}}{2m}+\frac{1}{2}m\omega^{2}{[x(t)]}^{2},$ (3)

where $x(t)$ and $p(t)$ denote the position and momentum of the oscillator at
time $t$. The equation of motion associated to this Hamiltonian is given by:

$\displaystyle\ddot{x}$ $\displaystyle=-\omega^{2}x$ (4)

with solutions

$\displaystyle x\left(t\right)$ $\displaystyle=x\left(0\right)\cos\left(\omega
t\right)+\frac{p\left(0\right)}{m\omega}\sin\left(\omega t\right).$ (5)

$x(t)$ has the same form as Eq. 2. Therefore, the quadratures $X$ and $Y$
associated to a monochromatic electric field can be thought of as the position
and momentum of a harmonic oscillator. It is, in this sense that, in the
literature [10], the quadratures of the electromagnetic field are referred to
as the generalized position and momentum of a field, and allow the
representation of optical fields in a diagram of conjugate variables, i.e.,
Fourier transforms of each other, $X$ and $Y$, called a phase space.

The position and momentum variables that have been defined so far correspond
to two conjugate variables that can be used to represent the information
conveyed by $\mathbf{E}(t).$ Nevertheless, said information could also be
represented by other pairs of conjugate variables, called generalized
quadratures $\\{X_{\alpha},X_{\alpha+\pi/2}\\}$, which are defined as

$X_{\alpha}=X\cos\alpha+Y\sin\alpha.$ (6)

Quadratures are measured using different experimental approaches that depend
on the frequency of the electromagnetic field under study: In the range of
frequencies $f=\omega/2\pi$ between kHz and MHz, the use of a lock-in
amplifier for measuring the quadratures of electric fields is a well-
established technique, with commercially available devices that are used
everyday in research laboratories [5, 6, 8].

The working principle of a lock-in amplifier is described as follows:
Following Fig. 1(a), an input signal of interest with a known frequency
$\omega_{S}$ but with unknown amplitude $E_{S}$ and phase $\theta_{S}$,
denoted by
$E_{\text{signal}}\left(t\right)=E_{S}\cos\left(\omega_{S}t+\theta_{S}\right)$,
is mixed with a well-defined local oscillator signal
$E_{\mathrm{LO}}\left(t\right)=E_{R}\cos\left(\omega_{R}t+\theta_{R}\right)$,
in which the frequency $\omega_{R}$, amplitude $E_{R}$ and phase $\theta_{R}$
are well known. The result of this mixing, which mathematically corresponds to
the product of the two signals, can be written as:

$\displaystyle
E_{\text{mix}}(t)=\frac{E_{S}E_{R}}{2}\left[\cos\left(\left(\omega_{S}+\omega_{R}\right)t+\left(\theta_{S}+\theta_{R}\right)\right)\times\right.$
(7)
$\displaystyle\left.\cos\left(\left(\omega_{S}-\omega_{R}\right)t+\left(\theta_{S}-\theta_{R}\right)\right)\right]$

The lock-in amplifier works in the homodyne detection regime, meaning that the
frequencies of the signal and the local oscillator are set to be equal:
$\omega_{S}=\omega_{R}=\omega$. Additionally, the lock-in amplifier uses a
low-pass filter that only keeps the DC component of the mixed signal.
Therefore, after mixing and filtering, Eq. 7 becomes:

$E_{\text{mix}}=\frac{E_{S}E_{R}}{2}\cos\left(\theta_{\mathrm{LO}}\right),$
(8)

where $\theta_{\mathrm{LO}}=\theta_{R}-\theta_{S}$ is the relative phase
between the fields $E_{\mathrm{signal}}$ and $E_{\mathrm{LO}}$.

Since the $X$ quadrature of the signal field is $X_{S}=E_{S}\cos{\theta_{S}}$,
choosing $\theta_{\mathrm{LO}}=0$ and $\theta_{\mathrm{LO}}=\pi/2$ in Eq. 8
enables the retrieval of the values of $X_{S}$ and $Y_{S}$ according to the
following equations:

$\displaystyle E_{\mathrm{mix}}^{(\theta_{\mathrm{LO}}=0)}$
$\displaystyle=\frac{1}{2}E_{S}X_{S},$ (9) $\displaystyle
E_{\mathrm{mix}}^{(\theta_{\mathrm{LO}}=\pi/2)}$
$\displaystyle=\frac{1}{2}E_{S}Y_{S}.$ (10)

This allows to see that, with the appropriate choice of the phase
$\theta_{\mathrm{LO}}$, the lock-in amplifier enables the retrieval of the
quadratures $X_{S}$ and $Y_{S}$.

$\otimes$$E_{\mathrm{\mathrm{LO}}}$$E_{\mathrm{S}}$Low
passfilter($\omega_{S}=\omega_{R}$, set $\theta_{\mathrm{LO}}=0$ as
reference)$X_{S}$$Y_{S}$ (a) Lock-in amplifier

$\mathrm{[S]}$
---
$\mathrm{[LO]}$
---
BSPS $\mathrm{D}_{\mathrm{2}}$
---
$\mathrm{D}_{\mathrm{1}}$
---
(b) Homodyne detector

$X$$Y$$s=0$$s$$\phi$$\mathrm{pr}_{u}(s,\phi)$ (c) WDF and marginal
distributions

Figure 1: (a) Schematic representation of a lock-in amplifier. (b) Scheme for
homodyne detection. (c) Scheme of the WDF and one of its experimentally
accessible marginal distributions, $\mathrm{pr}_{u}(s,\phi)$.

When one is interested in measuring electromagnetic fields in the optical
domain (i.e., with THz frequency), there is a limitation in the speed of the
electronics, which at the time of writing do not exceed the order of a few
GHz. However, there is an optical implementation of homodyne detection that
enables the access to the quadratures of light. This technique has been used
to characterize different quantum states of light [11, 12, 13]. The setup to
perform optical homodyne detection is shown in Fig. 1(b). The local
oscillator, (LO), and the light source to be characterized (S) are combined in
the two input ports of a beam splitter (BS) and then detected using photo-
detectors ($\mathrm{D}_{1}$ and $\mathrm{D}_{2}$). The relative phase between
both fields can be changed by using a phase shifter (PS). The detected
currents are subtracted, eliminating classical fluctuations as they are
correlated in both outputs of the BS. The complex amplitude of the fields at
$\mathrm{D}_{1}$ and $\mathrm{D}_{2}$ can be written as [14]

$\displaystyle\tilde{E}_{\mathrm{D}1/\mathrm{D}2}$
$\displaystyle=\frac{1}{\sqrt{2}}(\tilde{E}_{\mathrm{signal}}\pm
e^{i\theta_{\mathrm{LO}}}\tilde{E}_{\mathrm{LO}}).$ (11)

Since the current at the detectors is proportional to the intensity of light
arriving at them and the detectors have a finite response time compared with
the fast oscillations of the optical fields, the time-averaged substraction of
the currents associated to each detector, $i_{-}$, can be written in terms of
the quadratures of the signal field:

$\displaystyle i_{-}$ $\displaystyle=\left\langle
I_{1}-I_{2}\right\rangle_{t}$ (12)
$\displaystyle\propto\tilde{E}_{\mathrm{LO}}\left[X_{S}\cos\left(\theta_{\mathrm{LO}}\right)+Y_{S}\sin\left(\theta_{\mathrm{LO}}\right)\right].$

Here, $\left\langle\cdot\right\rangle_{t}$ denotes a time averaging.

The result in Eq. (12) shows that the subtracted current allows the
measurement of the quadratures of optical signals. Analogous to the case of
the lock-in amplifier, this is achieved by changing the relative phase between
the fields: For $\theta_{\mathrm{LO}}=0$ and $\theta_{\mathrm{LO}}=\pi/2$, the
values of $X_{S}$ and $Y_{S}$ are retrieved, respectively.

###  II.2 Wigner distribution function (WDF)

One tool that can be used to represent signals in phase space is the Wigner
Distribution Function (WDF). Although this tool was introduced in the early
1930s for analyzing quantum states [15], it can be used to represent any
signal $u(x)$ in the phase space spanned by two conjugate variables $X$ and
$Y$. The WDF of a signal $u(x)$ is defined as:

$W_{u}(X,Y)=\\\ \int
u\left(X-\frac{X^{\prime}}{2}\right)u^{*}\left(X+\frac{X^{\prime}}{2}\right)e^{-2\pi
iX^{\prime}Y}dX^{\prime}.$ (13)

The projection $\mathrm{pr}_{u}(s,\phi)$ of the WDF along the $s$-axis spanned
by the unit vector $\mathbf{e}_{s}=(\cos\phi,\sin\phi)$ is given by:

$\text{pr}_{u}\left(s,\phi\right)=\\\
\int_{-\infty}^{\infty}W_{u}(s\cos\phi-Y\sin\phi,s\cos\phi+Y\sin\phi)dY.$ (14)

The $\mathrm{pr}_{u}(s,\phi)$ are the marginal probability distributions of
obtaining a value for the quadrature $s$. These projections have a physical
meaning: They are the energy distributions of the signal $u(x)$ in the
generalized quadrature $s$. Despite the marginal distributions being positive
and having a physical meaning, the WDF is not necessarily defined as a
positive function. In the context of probability distributions, this makes the
WDF a joint quasiprobability distribution between conjugate variables[15].

Unlike $W_{u}(X,Y)$, the $\mathrm{pr}_{u}(s,\phi)$ are experimentally
accessible. From a set of $\mathrm{pr}_{u}(s,\phi)$, it is possible to use
three-dimensional reconstruction protocols to recover the WDF. Specifically,
by means of a Radon Inverse transform, it is possible to back-project the set
of $\mathrm{pr}_{u}(s,\phi)$ to recover the WDF of a signal[1, 4]. This
corresponds to integrating the different projections weighted by a function
$K(s)$, called the kernel of the Radon transform. Mathematically,

$W_{u}(X,Y)=\\\
\frac{1}{2\pi^{2}}\int_{\phi=0}^{\pi}\int_{s=-\infty}^{+\infty}\text{pr}_{u}(s,\phi)K(X\cos\phi+Y\sin\phi-s)dsd\phi,$
(15)

with $K(s)$ given by

$K\left(s\right)=\frac{1}{2}\int_{-\infty}^{+\infty}|\xi|e^{i\xi
s}\mathrm{~{}d}\xi.$ (16)

At the moment of applying Eq. 15, it is necessary to characterize and limit
the behavior of the kernel $K(s)$ due to the fact that any measurement relies
on the discretization of amplitudes and phases. In order to avoid ill
behaviors [16], the kernel is regularized [1], or filtered, by introducing a
cutoff value, $k_{c}$. When changing the integration limits in Eq. 16 to
$[-k_{c},k_{c}]$ and performing a Taylor expansion for low values of $s$, the
filtered kernel $K_{\text{filter}}$ becomes

$K_{\text{filter}}\left(s\right)=\\\
\begin{cases}\frac{1}{s^{2}}\left(\cos\left(k_{c}s\right)+k_{c}s\sin\left(k_{c}s\right)-1\right)&\left|s\right|>0.1/k_{c},\\\
\frac{k_{c}^{2}}{2}\left(1-\frac{k_{c}^{2}s^{2}}{4}+\frac{k_{c}^{4}s^{4}}{72}-...\right)&\text{otherwise.}\end{cases}$
(17)

By applying Eq. 17 into Eq. 15, the WDF is recovered numerically by having
performed a filtered back-projection algorithm.

##  III Quantum states of light

Understanding certain optical fields, generated experimentally in the last
decades [10], requires a fully quantum-mechanical picture. In this picture, it
is possible to describe the dynamics of a mode of the electromagnetic field
using the Hamiltonian of a quantum harmonic oscillator of frequency $\omega$
as [17]:

$\hat{H}_{\text{opt}}=\hbar\omega\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right),$
(18)

where $\hat{a}^{\dagger}$ and $\hat{a}$ are creation and annihilation
operators, respectively. From these operators, it is possible to define a pair
of conjugate variables known as generalized position and momentum:

$\displaystyle\hat{x}$
$\displaystyle=\frac{1}{2}\left(\hat{a}+\hat{a}^{\dagger}\right),$ (19)
$\displaystyle\hat{p}$
$\displaystyle=-\frac{i}{2}\left(\hat{a}-\hat{a}^{\dagger}\right).$ (20)

The presence of position and momentum observables permits the introduction of
a phase space and the definition of quadratures in the same spirit of those
mentioned in the previous section.

According to quantum mechanics, a physical system can be described by a
density matrix $\hat{\rho}$. The Wigner function of such a quantum state
$\hat{\rho}$ has the same properties as the WDFs introduced above and is given
by

$W_{\hat{\rho}}(x,p)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{ipq}\langle
x-q/2|\hat{\rho}|x+q/2\rangle dq,$ (21)

where $\langle u|\hat{\rho}|v\rangle$ is the matrix element $(u,v)$ when
$\hat{\rho}$ is represented in the continuous position basis.

Statistically, the measurement of the generalized quadrature rotated by an
angle $\theta$: $\hat{s}=\hat{x}\cos\theta+\hat{p}\sin\theta$, on a physical
system, described by $\hat{\rho}$, has a mean given by the trace of the
operator $\hat{s}\hat{\rho}$,

$\left\langle\hat{s}\right\rangle_{\hat{\rho}}=\text{Tr}\left[\hat{s}\hat{\rho}\right],$
(22)

and a standard deviation given by

$\Delta\hat{s}_{\hat{\rho}}=\sqrt{\left\langle\hat{s}^{2}\right\rangle_{\hat{\rho}}-\left\langle\hat{s}\right\rangle_{\hat{\rho}}^{2}}.$
(23)

The conjugate quadratures
$(\hat{s},\hat{s}^{\prime}=-\hat{x}\sin\theta+\hat{p}\cos\theta)$ of a quantum
state satisfy the uncertainty principle:

$\Delta\hat{s}_{\hat{\rho}}\Delta\hat{s}^{\prime}_{\hat{\rho}}\geq\frac{1}{2}.$
(24)

Using the concepts just introduced, this section describes the physical
properties of a gallery of different quantum states of light. In particular,
we describe their density matrices, quadrature distributions, Wigner functions
and phasor diagrams.

###  III.1 Fock states

One respresentation suitable to describe eigenstates of
$\hat{H}_{\mathrm{opt}}$ is the Fock basis, constituted as

${\cal B}=\left\\{\left|n\right\rangle\right\\}_{n=0,1,2,...},$ (25)

where each of the elements is an eigenstate of the number operator
$\hat{n}=\hat{a}^{\dagger}\hat{a}$:
$\hat{n}\left|n\right\rangle=n\left|n\right\rangle$. In this representation, a
state $\left|n\right\rangle$ has a well-defined number of excitations, i.e.,
photons. The density matrix of a quantum state $\hat{\rho}$ in this
representation is

$\rho_{nm}=\left\langle n\right|\hat{\rho}\left|m\right\rangle,$ (26)

which is illustrated in Fig 2 (1a) for the Fock state $|n=1\rangle$.

In a physical system described by the elements of the Fock basis, the
probability distribution for measuring the generalized quadrature $s$ rotated
by an angle $\theta$ is independent of $\theta$ and is given by [18]

$\displaystyle\mathrm{pr}_{\left|n\right\rangle}\left(s,\theta\right)$
$\displaystyle=\left|\langle s|n\rangle\right|^{2}$
$\displaystyle=\frac{1}{2^{n}n!\sqrt{\pi}}e^{-s^{2}}\left(H_{n}\left(s\right)\right)^{2}.$
(27)

Here, $\langle s|n\rangle=\psi_{n}\left(s\right)$, where
$\psi_{n}\left(s\right)$ is the wavefunction associated to the basis state
$|n\rangle$ in the $s$ representation. Additionally, $H_{n}(s)$ is the Hermite
polynomial of $n$th order. This is illustrated in Fig. 2 (2a) for
$|n=1\rangle$. Due to the form of the Hermite polynomial, all the projections
have the shape of a double-hump.

Following Eq. 21 the Wigner function of a Fock state can be written using the
Laguerre polynomials $L_{n}(x)$ as [1]:

$W_{\left|n\right\rangle}(x,p)=\\\
\frac{(-1)^{n}}{\pi}\exp\left(-\left(x^{2}+p^{2}\right)\right)L_{n}\left(2\left(x^{2}+p^{2}\right)\right).$
(28)

This is illustrated in Fig 2 (3a) for $n=1$.

From the Wigner function, it is possible to perform a transverse cut parallel
to the $(x,p)$ plane and draw a corresponding phasor diagram, illustrating the
means and variances of the quadratures: Fock states have a mean value for the
quadrature $s$ of $\left\langle\hat{s}\right\rangle_{\left|n\right\rangle}=0$
and a standard deviation of

$\Delta\hat{s}_{\left|n\right\rangle}=\sqrt{\frac{2n+1}{2}},$ (29)

both of which are independent on the value of $\theta$. For all values of $n$,
the uncertainty principle satisfies

$\Delta\hat{s}\Delta\hat{s}^{\prime}=\frac{2n+1}{2}\geq 1/2,$ (30)

as illustrated in Fig. 2 (4a). For the $|n=1\rangle$ case, the phasor diagram
has a donut shape [10], while its Wigner function attains negative values.

###  III.2 Coherent states

Coherent states correspond to eigenstates of the annihilation operator
$\hat{a}$, and are written as $\left|\alpha\right\rangle$, where $\alpha$ is a
complex number that satisfies

$\displaystyle\hat{a}\left|\alpha\right\rangle=\alpha\left|\alpha\right\rangle.$
(31)

In the Fock basis, coherent states are written as [10]

$\displaystyle\left|\alpha\right\rangle=e^{-\frac{1}{2}\left|\alpha\right|^{2}}\sum_{n=0}^{\infty}\frac{\alpha^{n}}{\sqrt{n!}}\left|n\right\rangle.$
(32) From this expression and using the definition in Eq. (26) the density
matrix in the Fock basis for a coherent state is
$\displaystyle\rho_{nm}=\frac{e^{-\left|\alpha\right|^{2}}\alpha^{n}\alpha^{*m}}{\sqrt{n!m!}},$
(33)

illustrated in Figure 2 (1b) for the case $\alpha=2$.

For a coherent state, the probability distribution for measuing the
generalized quadrature $s$ rotated by an angle $\theta$ is given by

$\text{pr}_{\left|\alpha\right\rangle}\left(s,\theta\right)=\sqrt{\frac{1}{\pi}}e^{-(-s+\text{Re}\left(\alpha\right)\cos\theta+\text{Im}\left(\alpha\right)\sin\theta)^{2}}.$
(34)

An illustration of this is shown in Fig. 2 (2b) for $\alpha=2$.

The Wigner function for a coherent state can be calculated following Eq. 21 to
be a 2D Gaussian distribution:

$W_{\left|\alpha\right\rangle}(x,p)=\frac{1}{\pi}e^{-\left(\left(x-\text{Re}\left(\alpha\right)\right)^{2}+\left(p-\text{Im}\left(\alpha\right)\right)^{2}\right)}.$
(35)

Graphically, this is illustrated in Fig. 2 (3b).

Statistically, it is possible to find the mean value of the generalized
quadrature $\hat{s}$ for the coherent state. This value depends on $\theta$
and is given by

$\left\langle
s\right\rangle_{\left|\alpha\right\rangle}=\text{Re}\left(\alpha\right)\cos\theta+\text{Im}\left(\alpha\right)\sin\theta.$
(36)

The uncertainty in any generalized quadrature is

$\Delta
s^{2}=\left\langle\hat{s}^{2}\right\rangle_{\left|\alpha\right\rangle}-\left\langle\hat{s}\right\rangle_{\left|\alpha\right\rangle}^{2}=\frac{1}{2},$
(37)

for any value of $\theta$. For this reason, the phasor diagram associated to a
coherent state corresponds to a circle displaced from the center by $\alpha$
with a diameter of $1/\sqrt{2}$, as shown in Fig. 2 (4b). Coherent states
correspond to minimum uncertainty states, as they attain the minimum value
allowed by the uncertainty relation: $\Delta s\Delta s^{\prime}=1/2.$

###  III.3 Vacuum

The vacuum state is a particular case of both the coherent states and the Fock
states, where the amplitude is $\alpha=0$ and the photon number is $n=0$.
Therefore, the density matrix is one for the element $\rho_{00}$, as shown in
Fig. 2(1c).

The projections of the vacuum state can be found by setting $n=0$ in Eq. 27 or
$\alpha=0$ in Eq. 34. Therefore,

$\text{pr}_{\left|0\right\rangle}\left(s,\theta\right)=\sqrt{\frac{1}{\pi}}e^{-s^{2}}.$
(38)

This is shown in Fig. 2(2c). According to Eq. 21 and setting $n=0$ in Eq. 28
and $\alpha=0$ in Eq. 35, the Wigner function of the vacuum state is given by:

$W_{\left|0\right\rangle}(x,p)=\frac{1}{\pi}e^{-\left(x^{2}+p^{2}\right)}.$
(39)

Here, the displacement in the phasor diagram is zero for both quadratures, as
shown in Fig. 2 (3c). Just as before, the uncertainties in the quadratures are
given by Eq. 37, corresponding to minimum-uncertainty states, represented in
the phasor diagram as a circle centered at the origin with diameter
$1/\sqrt{2}$.

###  III.4 Thermal states

Thermal states characterize photons at a frequency $\omega$ radiated from a
blackbody at a temperature $T$ and thermodynamic beta $\beta=1/k_{B}T$. The
density matrix in the Fock basis for a thermal state is given by [19]:

$\hat{\rho}_{\text{th}}=(1-e^{-\beta\hbar\omega})\sum_{n=0}^{\infty}e^{-n\beta\hbar\omega}|n\rangle\langle
n|.$ (40)

Since the thermal distribution has a mean number of photons of the form

$\langle n\rangle=\frac{1}{e^{\beta\hbar\omega}-1},$ (41)

the density matrix can be rewritten as

$\hat{\rho}_{\text{th}}=\sum_{n=0}^{\infty}\frac{\langle
n\rangle^{n}}{\left(\langle n\rangle+1\right)^{n+1}}|n\rangle\langle n|.$ (42)

As revealed in Fig. 2 (1d), the density matrix of a thermal state is diagonal
and therefore represents a statistical mixture of Fock states.

Using the linearity of Eq. 21, it is possible to calculate the Wigner function
of a thermal state by writing the weighted sum of the Wigner functions of all
the different Fock states:

$W_{\text{th}}\left(x,p\right)=\sum_{n=0}^{\infty}\frac{\langle
n\rangle^{n}}{\left(\langle
n\rangle+1\right)^{n+1}}W_{\left|n\right\rangle}\left(x,p\right)\text{ }.$
(43)

This equation can be rewritten taking into account that the generating
function of the Laguerre polynomials is given by
$\sum_{n=0}^{\infty}t^{n}L_{n}(u)=\frac{1}{1-t}e^{-tu/(1-t)}$ and setting
$t=-\frac{\langle n\rangle}{\langle n\rangle+1}$ and
$u=2\left(x^{2}+p^{2}\right)$, the Wigner function for the thermal state
becomes:

$W_{\text{th}}\left(x,p\right)=\frac{1}{\sqrt{2\pi\sigma_{x}^{2}}}e^{-x^{2}/2\sigma_{x}^{2}}\frac{1}{\sqrt{2\pi\sigma_{p}^{2}}}e^{-p^{2}/2\sigma_{p}^{2}},$
(44)

where $\sigma_{x}=\sigma_{p}=\sigma=\sqrt{\frac{2\langle n\rangle+1}{2}}$. Eq.
44 reveals that the Wigner function for a thermal state is the product of two
Gaussian functions in the quadratures $x$ and $p$ centered at the origin. The
independence between $x$ and $p$ reveals a lack of correlations betweeen the
quadratures. This lack of correlations implies that the projection in any
generalized quadrature will have the same mean, centered at zero, and
variance, given by $\sigma$. The projections are shown in Fig. 2(2d) and the
Wigner function is shown in Fig. 2(3d). In the phasor diagram (Fig. 2(4d)), a
thermal state corresponds to a circle of non-minimum uncertainty centered
around zero.

###  III.5 Squeezed states

Squeezed states have different uncertainties in each of their quadratures. Any
general state given by the density matrix $\hat{\rho}$ can be transformed
using the squeezing operator

$\hat{S}(\zeta)=\exp\left[\frac{1}{2}\left(\zeta^{*}\hat{a}^{2}-\zeta\hat{a}^{\dagger
2}\right)\right],$ (45)

where $\zeta$ is a complex squeezing parameter. The squeezing operator
modifies the density matrix $\hat{\rho}$ of any state as
$\hat{\rho}_{\mathrm{sq}}=\hat{S}\hat{\rho}\hat{S}^{\dagger}$, rescaling the
Wigner function of the original state $\hat{\rho}$ [19]:

$\begin{aligned} &W_{\hat{S}(\hat{\rho})}(x,p)\\\
&=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{ipq}\left\langle
x-\frac{q}{2}\left|\hat{S}\hat{\rho}\hat{S}^{\dagger}\right|x+\frac{q}{2}\right\rangle
dq\\\
&=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{ipq}\mathrm{e}^{\zeta}\left\langle\mathrm{e}^{\zeta}\left(x-\frac{q}{2}\right)|\hat{\rho}|\mathrm{e}^{\zeta}\left(x+\frac{q}{2}\right)\right\rangle
dq\\\
&=W_{\hat{\rho}}\left(\mathrm{e}^{\zeta}x,\mathrm{e}^{-\zeta}p\right).\end{aligned}$
(46)

The rescaling is different in each of the quadratures, compressing one and
expanding the other one, as indicated by the different sign in the
exponentials in the last line of Eq. 46.

One example of a squeezed state is the squeezed vacuum. Since vacuum is a pure
state, the effect of the squeezing operator can be written as
$\hat{S}\left(\zeta\right)\left|0\right\rangle$. Writing
$\zeta=re^{2i\delta}$, $r$ can be seen as the squeezing amplitude, while
$\delta$ corresponds to the angle along which squeezing is performed.

In the Fock basis, the squeezed vacuum can be written as [20]

$\hat{S}\left(re^{2i\delta}\right)\left|0\right\rangle=\frac{1}{\sqrt{\cosh(r)}}\\\
\sum_{m=0}^{\infty}(-1)^{m}\frac{\sqrt{(2m)!}}{2^{m}m!}e^{2im\delta}\tanh(r)^{m}\left|2m\right\rangle.$
(47)

From this expression, the density matrix for squeezed vacuum is calculated to
be

$(\rho_{\text{sq}})_{nm}=\frac{1}{\cosh(r)}(-1)^{m+n}\\\
\frac{\sqrt{(2m)!(2n)!}}{2^{m}m!2^{n}n!}e^{2i(m-n)\delta}\tanh(r)^{m+n}$ (48)

for even values of $n$ and $m$. This implies that the density matrix in the
Fock representation for the squeezed vacuum only has values different from
zero in even values of $n$ and $m$, as shown in Fig. 2 (1e).

Following Eq. 46, the Wigner function of squeezed vacuum corresponds to a
Gaussian function with different widths in the quadratures, rotated by an
angle $\delta$. Just like vacuum states, squeezed vacuum has a zero mean
around any generalized quadrature. However, the variances along the
quadratures $x$ and $p$ are given by:

$\displaystyle\Delta x^{2}$
$\displaystyle=\frac{1}{4}\left[\cosh\left(2r\right)-2\sinh r\cosh r\cos
2\delta\right],$ (49) $\displaystyle\Delta p^{2}$
$\displaystyle=\frac{1}{4}\left[\cosh\left(2r\right)+2\sinh r\cosh r\cos
2\delta\right].$ (50)

For example, if we perform squeezing along the quadrature $p$, i.e.,
$\delta=\pi/2$, the variances will be reduced to

$\Delta x^{2}=\frac{1}{4}e^{+2r},\quad\text{ and }\quad\Delta
p^{2}=\frac{1}{4}e^{-2r}.$ (51)

Although the squeezing is not equally distributed in both quadratures, $\Delta
x\Delta p$ still attains minimum uncertainty, while interestingly, one of the
uncertainties on the quadratures can be made smaller than that of the vacuum
state. The marginal distributions of this state are shown in Fig. 2 (2e),
while the Wigner function is shown in Fig. 2 (3e). As a consequence, the
phasor diagram in Fig. 2 (4e) is an ellipse rotated by an angle $\delta$ whose
major and minor axes have been stretched.

Figure 2: Column (a) shows, for a Fock state with $n=1$, the density matrix
(1), marginal distributions (2), Wigner function (3) and phasor diagram
representations (4). Likewise, Column (b) shows, for a coherent state with
$\alpha=2$, the same representations. Column (c) does the same for a vacuum
state ($\alpha=n=0$), and Column (d) shows the same representations for a
thermal state with $\langle n\rangle=10$, and finally, Column (e) shows, for
squeezed vacuum with squeezing factor $\zeta=2e^{i\pi/2}$, the different
representations.

##  IV Experimental Implementation by students

The following section describes the experimental and computational procedures
constituting an educational activity for recovering the Wigner function of
different quantum states of light. These procedures were implemented by
undergraduate students interested in learning the topic, and was documented in
a bachelor monography, mandatory to obtain the degree of physicist at
Uniandes. Additionally, the created experimental setup remains in the
laboratory and is the one used by students in the quantum optics class.

Laser
---
633nm
---
Dove prism
---
$\mathrm{BS}_{2}$$\mathrm{BS}_{1}$PZT  $\frac{\lambda}{2}$
---
$\frac{\lambda}{2}$
---
[Signal]
---
[LO]
---
$L$
---
$L$
---
$\mathrm{PD}_{\mathrm{2}}$
---
$\mathrm{PD}_{\mathrm{1}}$
---

Figure 3: Experimental setup for optical homodyne detection.

###  IV.1 Experimental activity

To recover the experimentally accessible quadratures, an optical set-up for
homodyne detection, such as the one depicted in Fig. 3, was implemented. A
stable $633\mathrm{nm}$ He-Ne laser with a power close to 0.5 mW and linear
polarization provides the light source for both S and LO. The laser is split
using a beam splitter, labeled $\mathrm{BS}_{1}$. Both beams are then
recombined using another beam splitter, $\mathrm{BS}_{2}$. The half-wave
plates before $\mathrm{BS}_{2}$ ensure the proper polarization matching
between the two beams. Finally, two lenses with a 30mm focal length are used
to focus the light into two silicon photodiodes (Thorlabs FDS100). To change
the LO phase, a dove prism is mounted into a piezo-electric (PZT) platform,
changing the optical path difference between both inputs of the homodyne
detector. The subtraction of detector outputs can be electronically
implemented either by using a subtracting circuit or an oscilloscope. For the
present implementation, the latter is used: An oscilloscope (Tektronix DPO
4054) records each photodiode output and the subtracted data, taking ten
thousand data points during an acquisition time of 2 seconds, and the process
is repeated for a new position of the piezo-electric, changed in steps of
$\Delta x=0.01\ \mathrm{\mu m}$ in a range of $0.6\ \mathrm{\mu m}$ for a
stepwise change of the local oscillator phase of
$\Delta\phi_{\mathrm{LO}}=\pi/12$.

A pipeline for the tomographic reconstruction of the Wigner function is shown
in Fig. 4: First, different values of $i_{-}$ are sampled (Fig. 4(a)), second,
they are tallied as different histograms (Fig. 4(b)), third, histograms are
recorded for different values of $\phi_{\mathrm{LO}}$ (Fig. 4(c)), and
finally, using a filtered back-projection algorithm, the Wigner function is
recovered (Fig. 4(d)).

(a) Samples for fixed $\phi_{\mathrm{LO}}$ (b) Histogram @
$\phi_{\mathrm{LO}}=0$ (c) Set of $\mathrm{pr}_{u}(s,\phi_{\mathrm{LO}})$ (d)
Recovery of Wigner

(e) Graphical User Interface (GUI)

Figure 4: Pipeline for Wigner function reconstruction. (a) Raw data retrieval
for a fixed phase $\phi_{\mathrm{LO}}$. (b) The experimental outcomes are
tallied to generate a histogram of $i_{-}$ values for a fixed
$\phi_{\mathrm{LO}}$. (c) Set of measured marginal distributions
$\mathrm{pr}_{u}(s,\phi_{\mathrm{LO}})$ for different values of
$\phi_{\mathrm{LO}}$. (d) Reconstructed Wigner function for a coherent state
with $\alpha=2$. (e) Graphical User Interface (GUI) for tomographic
reconstruction of the Wigner function.

###  IV.2 Tomographic reconstruction and simulation

The recovery of the Wigner function from experimentally measured data was
developed in Python by implementing a Graphical User Interface (GUI), shown in
Fig. 4(e). The GUI receives the set of marginal distributions as a function of
the LO phase, and computes the reconstruction of the Wigner function
performing a filtered back-projection algorithm. This program is openly
accessible, and can be found at https://github.com/amartinez1224/quantum-
tomography.

The software created for the tomographic reconstruction of the Wigner function
can also be used to analyze simulated data. This feature was introduced due to
COVID restrictions that forbade student access to the laboratory. These
simulations can be used to emulate the measurement of more exotic states of
light which are not experimentally feasible to achieve in an educational
laboratory with a HeNe laser, and constitutes an important tool for distance
learning.

The GUI uses the information from marginal distributions, simulated or
measured. The file to be introduced must be formatted appropriately,
containing three sets of data embedded in a JSON file as indicated in the
aforementioned website. Two arrays, s and phi, are row vectors containing all
the possible values of $i_{-}$ and the phase $\theta$, respectively. A third
component, pr, is formatted as a two-dimensional array of dimensions
size(s)*size(phi) containing the frequencies of the histograms,
$\mathrm{pr}_{u}(s,\theta)$. The GUI requires the user to define the ranges of
$X$ and $Y$ to display the tomography. This is done by using the inputs Xmin,
Xmax, Ymin, Ymax and Density. The filter of the Radon kernel, $k_{c}$, can be
tuned by changing the parameter Kc. Finally, the Wigner function can be used
to obtain the density matrix of a quantum state of light in the photon number
basis.

Additionally, the interface enables color map and perspective settings, which
can be changed by using the computer mouse, and enables the user to save the
figures and the reconstructed tomographic data.

##  V Hands-on activity

The different tools introduced in the previous section were condensed in a
laboratory practice developed by students and for students. It was
subsequently used in the “Quantum optics: theory and practice” course at
Uniandes. The laboratory activities consisted in attending the laboratory in
groups of maximum three students, who visited the lab in different pre-booked
two-hour shifts. The optical setup was already implemented, and students
controlled the piezoelectric to change the phase of the local oscillator and
perform the different quadrature measurements. The data was taken home to
process and hand in a report.

As mentioned before, since this practice was originally implemented during one
of the lockdowns of the COVID-19 pandemic, some students could not attend the
laboratory but could generate their own data using the part of the program
used to simulate data. In the in-person laboratory, students recovered the
Wigner function for a vacuum state (in which one of the ports for the homodyne
detector was blocked) and for a coherent state (a laser with a power of
approximately 0.5 mW). The students working remotely also reconstructed more
exotic states, such as single-photon states and squeezed vacuum. Some examples
of the Wigner functions recovered by in-person students are shown in Fig. 5
(a) and (b), and one example of a Wigner function recovered by a student
following the course remotely is shown in Fig. 5 (c).

(a)  (b)  (c)

Figure 5: (a) and (b) show experimental results obtained by students for the
vacuum and the coherent states. (c) Shows the Wigner function obtained from
data simulated by a student to emulate the quadrature measurements of a
squeezed state.

##  VI Conclusions and perspectives

This paper has presented the development of an undergraduate educational
activity introducing the reconstruction of Wigner distribution functions using
optical homodyne tomography. Interestingly, this activity was developed by
undergraduate students in three semesters in different mandatory courses of
the physics undergraduate curriculum at Uniandes. This activity was later used
in the “Quantum Optics: theory and practice” course, which is open for
advanced undegraduate students, as well as early graduate students.

Since the use of the educational activity coincided with the COVID-19
lockdowns, the option of processing real experimental and simulated data was
included. This activity is now running in-person, and the virtual part opens
the possibility of implementing distance learning in a quantum optics course.
Future iterations of this activity could improve the experimental results
obtained by the students: Implementing a Maximum Likelihood method [21] would
enable the recovery of smoother Wigner functions, and a thorough
characterization of detector performance, in particular detector clearance and
common mode rejection ratio, would enable students to interpret and correct
their experimental data adequately.

###### Acknowledgements.

J.R.A. acknowledges funding by the European Union Horizon 2020 (MSCA
765075-LIMQUET, FET 899544-PHOQUSING), and from the Plan France 2030 through
the project ANR-22-PETQ-0006. A.V. ackowledges financial support from the
Facultad de Ciencias at the Universidad de los Andes.

## References

  * Leonhardt and Paul [1995] U. Leonhardt and H. Paul, Progress in Quantum Electronics 19, 89 (1995).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, Reviews of Modern Physics 81, 299 (2009).
  * Cormack [1973] A. M. Cormack, Physics in Medicine and Biology 18, 195 (1973).
  * Maier  _et al._ [2018] A. Maier, S. Steidl, V. Christlein, and J. Hornegger, eds.,  _Medical Imaging Systems_ (Springer International Publishing, 2018).
  * Wolfson [1991] R. Wolfson, American Journal of Physics 59, 569 (1991).
  * Libbrecht  _et al._ [2003] K. G. Libbrecht, E. D. Black, and C. M. Hirata, American Journal of Physics 71, 1208 (2003).
  * Case [2008] W. B. Case, American Journal of Physics 76, 937 (2008).
  * DeVore  _et al._ [2016] S. DeVore, A. Gauthier, J. Levy, and C. Singh, American Journal of Physics 84, 52 (2016).
  * Griffiths [2017] D. J. Griffiths,  _Introduction to Electrodynamics_ (Cambridge University Press, 2017).
  * Furusawa [2015] A. Furusawa, in  _SpringerBriefs in Mathematical Physics_ (Springer Japan, 2015) pp. 1–67.
  * Beck  _et al._ [1993] M. Beck, D. T. Smithey, and M. G. Raymer, Physical Review A 48, R890 (1993).
  * Breitenbach  _et al._ [1997] G. Breitenbach, S. Schiller, and J. Mlynek, Nature 387, 471 (1997).
  * Lvovsky  _et al._ [2001] A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, Physical Review Letters 87, 10.1103/physrevlett.87.050402 (2001).
  * Bachor and Ralph [2019] H.-A. Bachor and T. C. Ralph,  _A Guide to Experiments in Quantum Optics_ (Wiley, 2019).
  * Wigner [1932] E. Wigner, Physical Review 40, 749 (1932).
  * Roach [1991] G. F. Roach, ed.,  _Inverse Problems and Imaging_ , Pitman Research Notes in Mathematics Series No. 245 (Longman Scientific & Technical ; Wiley, Harlow, Essex, England : New York, 1991).
  * Cohen-Tannoudji  _et al._ [1997] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms_ (Wiley, 1997).
  * Beck [2013] M. Beck, Contemporary Physics 54, 77 (2013).
  * Leonhardt [2009] U. Leonhardt,  _Essential Quantum Optics_ (Cambridge University Press, 2009).
  * Gerry and Knight [2004] C. Gerry and P. Knight,  _Introductory Quantum Optics_ (Cambridge University Press, 2004).
  * Lvovsky [2004] A. I. Lvovsky, Journal of Optics B: Quantum and Semiclassical Optics 6, S556 (2004).