# Toolbox for non-classical state calculations

Filippus S. Roux<EMAIL_ADDRESS>National Metrology Institute of South
Africa, Meiring Naudé Road, Brummeria 0040, Pretoria, South Africa

###### Abstract

Computational challenges associated with the use of Wigner functions to
identify non-classical properties of states are addressed with the aid of
generating functions. It allows the computation of the Wigner functions of
photon-subtracted states for an arbitrary number of subtracted photons. Both
the formal definition of photon-subtracted states in terms of ladder operators
and the experimental implementation with heralded photon detections are
analyzed. These techniques are demonstrated by considering photon subtraction
from squeezed thermal states as well as squeezed Fock states. Generating
functions are also used for the photon statistics of these states. These
techniques reveal various aspects of the parameter dependences of these
states.

Keywords: quantum optics, photon subtraction, Wigner function, generating
function, squeezed thermal state

##  1 Introduction

The advantage of quantum information systems over their classical equivalents
can be related to the properties of the quantum states that serve as resources
for such systems [1, 2, 3]. Therefore, the “quantumness” of these states is an
important albeit elusive property, which can be defined or quantified in
different ways.

Such quantum states are often distinguished from classical states based on
some non-classical properties. An example of such a non-classical property is
the presence of regions where the state’s Wigner function in terms of the
particle-number degrees of freedom111A Wigner function computed for the
spatiotemporal degrees of freedom could have negative regions, but that would
not indicate non-classical properties. is negative. States with such
negativity are useful resource states for quantum information processing [4,
5, 6]. The experimental techniques to prepare such non-classical states are
reviewed in [7].

Most photonic states that are produced with existing technology have Gaussian
Wigner functions, which do not contain negative regions on phase space. As a
result, readily producible photonic states need additional processing to
prepare non-classical photonic states. Negativity is generally obtained by
applying some form of post-processing or post-selection on Gaussian states
[7].

One such process is photon subtraction, where one or more photons are removed
from a given input state [8]. Such photon-subtracted states can have regions
on phase space where the Wigner function is negative. However, the input state
needs to be a special state to produce such negativity [9]. Photon-subtracted
squeezed vacuum states have Wigner functions with negative regions. Their
Wigner functions resemble those of Schrödinger cat states. For example, the
three-dimensional view of the Wigner function of a squeezed vacuum state with
five photons subtracted, shown in figure 1, demonstrates its resemblence with
the Wigner function of a Schrödinger cat state. For a Schrödinger cat state
produced as the superposition of two coherent states, its Wigner function
consists of two Gaussian functions displaced to geometrically opposite
locations with respect to the origin and a Gaussian function at the origin
modulated by oscillations along a direction perpendicular to the
displacements. The Wigner function in figure 1 has the same qualitative
features. Photon subtraction from squeezed vacuum states has indeed been
proposed as a way to produce Schrödinger cat states [10].

Figure 1: Wigner function of a squeezed vacuum states with five photons
subtracted (not to scale).

Another way to produce Wigner negativity, closely related to photon
subtraction, is the opposite process: photon addition [11]. The formal
definition implies an interchange of creation and annihilation operators, and
the experimental implementation replaces the beam splitter with an optical
parametric amplifier, which can be seen as replacing the SU(2) process by a
SU(1,1) process. Here, we focus on photon subtraction since it suffices to
demonstrate the use of generating functions. However, similar techniques can
be used to aid the analytical calculations associated with photon addition.

While a large set of photon-subtracted states exists, depending on the number
of photons that are subtracted and the type of input state that is being used,
many of these cases are challenging to analyze, especially when multiple
photons are subtracted or when the input state is a mixed state. Such
challenges are often addressed with the aid of numerical simulations [12, 9].
To address these challenges without resorting to numerical simulations, we
propose a generating function approach for such analyses.

Formally, photon subtraction involves the application of ladder operators to
the state $\hat{\rho}\rightarrow\hat{a}\hat{\rho}\hat{a}^{{\dagger}}$. In
practice, the photon subtraction process is performed by detecting a specific
number of photons in a portion of the state, which heralds the existence of
the photon-subtracted state in the remainder [7].

Here, we investigate the formal process, as well as its experimental
implementation. We compute the Wigner functions after the subtraction process
is applied to an input state to see if it contains regions where the Wigner
function is negative. While many of these states are qualitatively similar,
their quantitative properties are different. For coherent states, neither the
formal subtraction process, nor any experimental implementation thereof can
produce a Wigner function with negative regions. Therefore, the input state
needs to be squeezed [9].

The formal photon-subtracted squeezed vacuum state has been compared to a
squeezed Fock state [7]. Here, we use a generating function for the Wigner
functions of Fock states to investigation this comparison.

We consider only the particle-number degrees of freedom in these analyses.
However, the approach presented here can also be applied in cases where the
spatiotemporal degrees of freedom are included. Thus, one can generalize it to
continuous variables in terms of the symplectic formalism [13, 14] for a
finite number of discrete modes. Alternatively, one can use it in a Wigner
functional formalism [15, 16, 17], which includes all the spatiotemporal
degrees of freedom without any restrictions or truncations.

##  2 Generating functions

The operations that are performed to produce photon-subtracted states can
either be modeled formally in terms of ladder operators or practically in
terms of projection operators. Both these approaches are generalized to
represent the subtraction of an arbitrary number of photons. Here, we
alleviate the burden of such calculations for arbitrary numbers of photon
subtractions by introducing generating functions.

To determine whether the Wigner functions of photon-subtracted states have
negative regions on phase space, the calculations can be converted into the
Moyal formalism [18, 19, 20] where all the operators are represented in terms
of their Wigner functions. Products of operators then become star products of
their Wigner functions.

Instead of using the Wigner functions of the ladder operators and projection
operators in these phase space calculations, we use the Wigner functions of
the generating functions of these operators. It alleviates the phase space
calculations by producing generating functions for the Wigner functions of
states with arbitrary numbers of subtracted photons.

###  2.1 The concept of a generating function

A generating function is a form of mathematical tagging, where multiple
mathematical entities are combined into one mathematical entity in which the
different entities are individually tagged in some way. The tagging makes it
possible to extract the individual entities from the combined entity. In the
case of a generating function, the tagging is done with the aid of a
generating parameter, raised to different powers and then multiplied with the
different entities (functions) before being summed. For example,

$\mathcal{G}(J)=\sum_{n=0}^{\infty}J^{n}f_{n},$ (1)

where $J$ is the generating parameter and $f_{n}$ represents the different
entities or functions in this case. The individual entities in a generating
function are extracted with the aid of derivatives in the same way that one
would calculate the terms in a Taylor series expansion. For the above example,

$f_{n}=\frac{1}{n!}\left.\partial_{J}^{n}\mathcal{G}\right|_{J=0}.$ (2)

The benefit of such combined entities is that one can apply a linear operation
to all the entities by applying the linear operation once to the combined
entity. For instance, applying a linear operation represented by
$L\\{\cdot\\}$ to the above example, we get

$L\\{\mathcal{G}(J)\\}=\sum_{n=0}^{\infty}J^{n}L\\{f_{n}\\}.$ (3)

The results of the linear operation on individual entities can be extracted
from this result in the same way as the individual entities are extracted from
the generating function.

###  2.2 Generating functions for ladder operators

The generating functions for multiple ladder operators are given by

$\hat{G}=\exp(\eta^{*}\hat{a})~{}~{}~{}{\rm
and}~{}~{}~{}\hat{G}^{{\dagger}}=\exp(\hat{a}^{{\dagger}}\eta),$ (4)

where $\hat{a}$ and $\hat{a}^{{\dagger}}$ are ladder operators (annihilation
and creation operators, respectively), and $\eta=J+\rmi K$ is a complex
generating parameter. The derivatives with respect to these parameters can be
defined in terms of derivatives with respect to their real and imaginary
parts:

$\partial_{\eta}\equiv\frac{1}{2}\left(\frac{\partial}{\partial
J}-\rmi\frac{\partial}{\partial K}\right)~{}~{}~{}{\rm
and}~{}~{}~{}\partial_{\eta^{*}}\equiv\frac{1}{2}\left(\frac{\partial}{\partial
J}+\rmi\frac{\partial}{\partial K}\right).$ (5)

It then follows that

$\partial_{\eta}\eta^{*}=\partial_{\eta^{*}}\eta=0.$ (6)

As a result, the complex generating parameter and its complex conjugate act as
independent generating parameters. Multiple ladder operators are produced from
these generating functions by

$\left(\hat{a}\right)^{n}=\left.\partial_{\eta^{*}}^{n}\hat{G}\right|_{\eta^{*}=0}~{}~{}~{}{\rm
or}~{}~{}~{}\left(\hat{a}^{{\dagger}}\right)^{n}=\left.\partial_{\eta}^{n}\hat{G}^{{\dagger}}\right|_{\eta=0}.$
(7)

The Wigner functions of the generating functions in equation (4) are given by

$W_{\hat{G}}=\exp(\eta^{*}\alpha)~{}~{}~{}{\rm
and}~{}~{}~{}W_{\hat{G}^{{\dagger}}}=\exp(\alpha^{*}\eta),$ (8)

where

$\alpha\equiv\frac{1}{\sqrt{2}}(q+\rmi p).$ (9)

The Wigner functions for specific powers of the ladder operators are obtained
by applying equation (7) to the Wigner functions in equation (8). For example,
by replacing $\hat{G}\rightarrow W_{\hat{G}}$ in equation (7), we obtain
$\alpha^{n}$.

###  2.3 Generating function for projection operators

The photon-number-resolving detection process is represented by projection
operators for fixed numbers of photons

$\hat{P}_{n}=\left|{n}\right\rangle\left\langle{n}\right|.$ (10)

The probability to detect $n$ photons is given by

$P(n)=\tr\\{\hat{P}_{n}\hat{\rho}\\}=\left\langle{n}\right|\hat{\rho}\left|{n}\right\rangle.$
(11)

However, it is more convenient to perform the calculations with the aid of a
generating function for these projection operators. It is defined by

$\mathcal{P}=\sum_{n=0}^{\infty}J^{n}\hat{P}_{n}=\sum_{n=0}^{\infty}\left|{n}\right\rangle
J^{n}\left\langle{n}\right|,$ (12)

where $J$ is the generating parameter, so that

$\hat{P}_{n}=\left.\frac{1}{n!}\partial_{J}^{n}\mathcal{P}\right|_{J=0}.$ (13)

When the generating parameter is set equal to 1, the generating function
reproduces the identity

$\left.\mathcal{P}\right|_{J=1}=\sum_{n=0}^{\infty}\left|{n}\right\rangle\left\langle{n}\right|=\mathds{1},$
(14)

where $\mathds{1}$ is the identity operator.

A generating function for the probability to detect $n$ photons from a state
$\hat{\rho}$ is obtained from the trace

$\mathcal{F}=\tr\\{\mathcal{P}\hat{\rho}\\}=\sum_{n=0}^{\infty}J^{n}\tr\\{\hat{P}_{n}\hat{\rho}\\}=\sum_{n=0}^{\infty}J^{n}P(n).$
(15)

This generating function does not only provide the means to produce the
individual probabilities. It can also be used to compute the expectation value
for the number of photons in the state. To calculate this expectation value,
we perform

$\left.\partial_{J}\mathcal{F}\right|_{J=1}=\sum_{n=0}^{\infty}nP(n)=\langle
n\rangle.$ (16)

The variance in the number of photons in the state, which is given by

$\sigma^{2}=\langle n^{2}\rangle-\langle n\rangle^{2},$ (17)

can be calculated with the aid of the second moment, which is directly
obtained from the generating function as

$\left.\partial_{J}\left(J\partial_{J}\mathcal{F}\right)\right|_{J=1}=\sum_{n=0}^{\infty}n^{2}P(n)=\langle
n^{2}\rangle.$ (18)

Since we intend to perform these calculations with the aid of Wigner
functions, we need the Wigner function of the generating function for these
projection operators. It is given by

$\mathcal{W}_{\mathcal{P}}=\frac{2}{1+J}\exp\left(-2|\alpha|^{2}\frac{1-J}{1+J}\right).$
(19)

Since the projection operator is also the density operator for a Fock state,
the Wigner function in equation (19) also serves as a generating function for
the Wigner functions of Fock states. This Wigner function shares some
properties with the original generating function. By setting $J=1$, we obtain
the sum over the Wigner functions of all the projection operators. It gives

$\mathcal{W}_{\mathcal{P}}(J=1)=1,$ (20)

which is the Wigner function for the identity. When we apply a derivative to
the Wigner function in equation (19) before setting $J=1$, the result is the
Wigner function for the number operator

$\left.\partial_{J}\mathcal{W}_{\mathcal{P}}\right|_{J=1}=|\alpha|^{2}-{\textstyle\frac{1}{2}}\equiv
W_{\hat{n}},$ (21)

which produces the expectation value for the number of photons in a state when
it is traced with the Wigner function of that state.

##  3 Mixed squeezed vacuum states

As an example, we’ll consider the case where the input state is a mixed
squeezed vacuum state. Such a state is obtained when a homogeneous loss (as
found when passing the state through a neutral density filter) is applied to a
pure squeezed vacuum state. It can also be produced as a squeezed thermal
state.

The Wigner function for a pure squeezed vacuum state is

$W_{{\rm sq}}=2\exp\left(-2|\alpha|^{2}A-\alpha^{2}B^{*}-\alpha^{*2}B\right),$
(22)

where $A^{2}-|B|^{2}=1$. The constants are parameterized in terms of a
complex-valued squeezing parameter $\xi=|\xi|\exp(\rmi\varphi)$, by

$A=\cosh(2|\xi|),~{}~{}~{}~{}~{}B=\exp(\rmi\varphi)\sinh(2|\xi|).$ (23)

They are related to the average number of photons in the squeezed vacuum state
$\langle{n}\rangle$ by

$A=1+2\langle{n}\rangle,~{}~{}~{}~{}~{}B=2\exp(\rmi\varphi)\sqrt{\langle{n}\rangle}\sqrt{\langle{n}\rangle+1}.$
(24)

Henceforth, we’ll assume that $\varphi=0$, which means the Wigner function of
the squeezed state is elongated along the $p$-direction and squeezed along the
$q$-direction.

When a pure squeezed vacuum state passes through a neutral density filter with
transmission coefficient $t<1$, the Wigner function becomes

$W_{{\rm
msq}}=\frac{2}{\sqrt{1+2ta(A-1)}}\exp\left[\frac{-2|\alpha|^{2}(a+tA)-t\alpha^{2}B^{*}-t\alpha^{*2}B}{1+2ta(A-1)}\right],$
(25)

where $a=1-t$ is the absorption coefficient. The resulting state is mixed,
having a purity of

${\rm purity}=\frac{1}{\sqrt{1+2ta(A-1)}}.$ (26)

Another way to produce a mixed squeezed vacuum state is to apply squeezing to
a thermal state. The Wigner function for a thermal state is

$W_{{\rm th}}=2T\exp\left(-2T|\alpha|^{2}\right),$ (27)

where $T<1$ is the purity of the state. For $T=1$, we obtain the Wigner
function of a vacuum state.

The squeezing process is done with the aid of a Bogoliubov transformation. It
is a unitary process that transforms the arguments of the Wigner function of
the state on which it is applied. Hence, the Bogoliubov transformation
produces

$W_{\hat{\rho}}(\alpha)\rightarrow W_{\hat{\rho}}(U\alpha+V\alpha^{*}),$ (28)

where $U$ and $V$ are constants expressed in terms of the complex-valued
squeezing parameter $\xi=|\xi|\exp(\rmi\varphi)$ by

$U=\cosh(|\xi|),~{}~{}~{}~{}~{}V=\exp(\rmi\varphi)\sinh(|\xi|),$ (29)

so that $U^{2}-|V|^{2}=1$.

After applying a Bogoliubov transformation to the thermal state, we obtain a
squeezed thermal state with a Wigner function expressed as

$W_{{\rm
sts}}=2T\exp\left(-2T|\alpha|^{2}A_{0}-T\alpha^{2}B_{0}^{*}-T\alpha^{*2}B_{0}\right).$
(30)

The constants $A_{0}$, $B_{0}$, and $B_{0}^{*}$ are produced by the Bogoliubov
transformation and are given in equation (23) in terms of the squeezing
parameter.

The expressions in equations (25) and (30) are formally equivalent. It can be
shown by setting

$\eqalign{T=\frac{1}{\sqrt{1+2ta(A-1)}},\cr
A_{0}=\frac{a+tA}{\sqrt{1+2ta(A-1)}},\cr
B_{0}=\frac{tB}{\sqrt{1+2ta(A-1)}},\cr
B_{0}^{*}=\frac{tB^{*}}{\sqrt{1+2ta(A-1)}},}$ (31)

in equation (30) to obtain equation (25). Since the Bogoliubov transformation
is a unitary process, the purity of the state is not affected by the
transformation.

##  4 Formal photon subtraction

The formal representation of the process for the subtraction of $n$ photons
from a state $\hat{\rho}$ is defined by

$\hat{\rho}_{{\rm
fs}}=\frac{\left(\hat{a}\right)^{n}\hat{\rho}\left(\hat{a}^{{\dagger}}\right)^{n}}{\tr\left\\{\left(\hat{a}\right)^{n}\hat{\rho}\left(\hat{a}^{{\dagger}}\right)^{n}\right\\}}.$
(32)

The denominator serves to normalize the state after the subtraction process,
which does not preserve the normalization of the state.

The Wigner function of the state obtained from the formal photon subtraction
process is obtained with the aid of a triple star product. It involves the
Wigner functions of the input state and those of the generating functions for
the ladder operators given in equation (8). By computing this triple star
product, we obtain a generating function for the unnormalized Wigner function
of the photon-subtracted states in terms of the Wigner function of the input
state. The expression reads

$\displaystyle\mathcal{W}_{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}}$
$\displaystyle=\int\exp[(\alpha^{*}-\alpha_{1}^{*})\alpha_{2}-\alpha_{2}^{*}(\alpha-\alpha_{1})]W_{\hat{G}}[{\textstyle\frac{1}{2}}(\alpha_{1}+\alpha+\alpha_{2})]$
(33) $\displaystyle~{}~{}~{}\times
W_{\hat{\rho}}(\alpha_{1})W_{\hat{G}^{{\dagger}}}[{\textstyle\frac{1}{2}}(\alpha_{1}+\alpha-\alpha_{2})]\
\rmd^{2}\alpha_{1}\ \rmd^{2}\alpha_{2}$
$\displaystyle=\exp\left(\eta^{*}\alpha+\alpha^{*}\eta+{\textstyle\frac{1}{2}}|\eta|^{2}\right)W_{\hat{\rho}}\left(\alpha+{\textstyle\frac{1}{2}}\eta\right),$

where

$\rmd^{2}\alpha_{n}=\frac{1}{2\pi}\ \rmd q_{n}\ \rmd p_{n}.$ (34)

For single-photon subtraction, we compute

$\left.\partial_{\eta}\partial_{\eta^{*}}\mathcal{W}_{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}}\right|_{\eta^{*}=\eta=0}={\textstyle\frac{1}{2}}\left({\textstyle\frac{1}{2}}\partial_{\alpha}\partial_{\alpha^{*}}+\alpha\partial_{\alpha}+\alpha^{*}\partial_{\alpha^{*}}+1+2|\alpha|^{2}\right)W_{\hat{\rho}},$
(35)

which gives a differential operation that produces the (unnormalized) Wigner
function for the single-photon-subtracted state when it is applied to the
Wigner function of an arbitrary input state. In a similar way, one can compute
the differential operations required for larger numbers of photon
subtractions.

The trace of the generating function

$\mathcal{W}_{\tr\\{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}\\}}=\int\mathcal{W}_{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}}\
\rmd^{2}\alpha,$ (36)

gives a generating function for the traces of the photon-subtracted states.
These traces serve as the inverses of the required normalization constants.
The Wigner functions for the normalized $n$-photon-subtracted states
(according to the formal process) are then given by

$W_{\hat{\rho}-n}=\left.\frac{\partial_{\eta^{*}}^{n}\partial_{\eta}^{n}\mathcal{W}_{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}}}{\partial_{\eta^{*}}^{n}\partial_{\eta}^{n}W_{\tr\\{\hat{G}\hat{\rho}\hat{G}^{{\dagger}}\\}}}\right|_{\eta=\eta^{*}=0}.$
(37)

###  4.1 Formal photon-subtracted squeezed thermal states

Applying equation (33) to the Wigner function in equation (30), we obtain a
generating function for photon subtraction from a squeezed thermal state,
according to the formal definition of the process. It is given by

$\displaystyle\mathcal{W}_{{\rm fs-sts}}=$ $\displaystyle
2T\exp\left[-2T|\alpha+{\textstyle\frac{1}{2}}\eta|^{2}A-T(\alpha+{\textstyle\frac{1}{2}}\eta)^{2}B^{*}-T(\alpha^{*}+{\textstyle\frac{1}{2}}\eta^{*})^{2}B\right.$
(38)
$\displaystyle\left.+\eta^{*}\alpha+\alpha^{*}\eta+{\textstyle\frac{1}{2}}|\eta|^{2}\right],$

where we dropped the subscript 0 of the constants. The generating function for
the traces of these states, obtained by computing the trace of the generating
function in equation (38), is given by

$\mathcal{W}_{{\rm tr\\{fs-
sts\\}}}=\exp\left(\frac{2|\eta|^{2}(A-T)-\eta^{2}B^{*}-\eta^{*2}B}{4T}\right).$
(39)

A state with $n$ subtracted photons is obtained from this generating function
with the aid of equation (37). For one subtracted photon, the normalized
Wigner function is

$\displaystyle W_{{\rm sts}-1}=$
$\displaystyle\left[\left(1-T^{2}\right)|\alpha|^{2}+(TA-1)\left(2T|\alpha|^{2}A+T\alpha^{2}B^{*}+T\alpha^{*2}B-{\textstyle\frac{1}{2}}\right)\right]$
(40)
$\displaystyle\times\frac{4T^{2}}{A-T}\exp\left(-2T|\alpha|^{2}A-T\alpha^{2}B^{*}-T\alpha^{*2}B\right).$

The region at the center, where $\alpha=0$ is negative, provided that $TA>1$.

Figure 2: Wigner functions of formal photon-subtracted squeezed vacuum states
($T=1$) with zero to five subtracted photons, plotted along the squeezed
direction. The input squeezed vacuum state has $\langle{n}\rangle=1$, so that
$A=3$ and $B=B^{*}=2.83$.

The Wigner functions of formal photon-subtracted squeezed vacuum states
($T=1$) and photon-subtracted squeezed thermal states with purity $T=0.9$, for
zero to five subtracted photons, are shown in figures 2 and 3, respectively.
In both figures, the parameter values for $A$, $B$, and $B^{*}$ are the same,
given by equation (24) for a squeezed vacuum state with $\langle{n}\rangle=1$.
All the curves are plotted along the line where $p=0$ as a function of $q$,
which is the squeezed direction. The curves in figures 2 and 3 show the
positive and negative regions produced by the photon subtractions, each having
a width that is smaller than the width of a coherent state. The number of
negative regions is equal to the number of subtracted photons. For an odd
number of subtracted photons, there is always a negative region located at the
origin. In figure 2, the amplitudes of the curves start from the same maximum
or minimum value, but in figure 3 the amplitudes decay for larger numbers of
subtracted photons. This trend indicates that mixing has a detrimental effect
on the amplitudes of the oscillations produced by the photon subtractions. The
extent of the sensitivity of these oscillations to a reduction in the purity,
is made apparent by the fact that, even for the slightly reduced value of
$T=0.9$, the decay in the amplitudes of the oscillations is already
significant.

Figure 3: Wigner functions of formal photon-subtracted squeezed thermal states
with zero to five subtracted photons, plotted along the squeezed direction for
$T=0.9$, $A=3$ and $B=B^{*}=2.83$.

##  5 Heralded photon subtraction

###  5.1 Experimental process

The formal definition of the photon subtraction process does not correspond
directly to a physical process; the physical implementation of the photon
subtraction process is different from its formal description. It uses
projective measurements in conjunction with a beam splitter to simulate the
action of the ladder operators. Nevertheless, one can model the physically
implemented photon subtraction process to investigate it analytically and
compare it with the formal description.

The experimental implementation of a photon subtraction process involves a
heralded process whereby the detection of a specific number of photons in a
portion of the state signals the existence of the required photon-subtracted
state in the remainder. For this purpose, a portion of the state is separated
from the remainder of the state with the aid of an unbalanced beam splitter.
The portion that is separated off, is then subjected to a measurement to
detect a specific number of photons. Ideally, a photon-number-resolving
detector is needed for this purpose, but the generating function approach
allows one to investigate the case where the detector is not photon-number
resolving.

The two input ports of the beam splitter respectively receive the input state
and a vacuum state. The unitary process that represents the beam splitter
transforms the arguments of the product of the Wigner functions of these two
states, leading to

$W_{\hat{\rho}}(\alpha)W_{{\rm vac}}(\beta)\rightarrow
W_{\hat{\rho}}(\sqrt{1-\zeta}\alpha+\rmi\sqrt{\zeta}\beta)W_{{\rm
vac}}(\sqrt{1-\zeta}\beta+\rmi\sqrt{\zeta}\alpha),$ (41)

where $\zeta$ represents the reflectivity of the beam splitter. The
measurement is done on the portion of the state represented by $\beta$.

The projective measurement is performed with the aid of the generating
function for the projection operators in equation (19), with
$\alpha\rightarrow\beta$. After multiplying it with the Wigner function after
the beam splitter, we compute the trace by integrating over $\beta$:

$\displaystyle\mathcal{W}_{{\rm hs}}=$ $\displaystyle\int
W_{\hat{\rho}}(\sqrt{1-\zeta}\alpha+\rmi\sqrt{\zeta}\beta)W_{{\rm
vac}}(\sqrt{1-\zeta}\beta+\rmi\sqrt{\zeta}\alpha)$ (42)
$\displaystyle\times\frac{2}{1+J}\exp\left(-2|\beta|^{2}\frac{1-J}{1+J}\right)\
\rmd\beta.$

The result is a generating function for the Wigner functions of the states
that are heralded by the detection of certain numbers of photons. However,
these Wigner functions are not normalized, because the heralding process does
not preserve the trace. The trace of the generating function in equation (42),
which implies an integration over $\alpha$, provides a generating function for
the traces of states,

$\mathcal{W}_{\tr\\{{\rm hs}\\}}=\int\mathcal{W}_{{\rm hs}}\ \rmd^{2}\alpha,$
(43)

as it was done for the formal process in equation (36). The individual
normalized Wigner functions are then generated by

$W_{\hat{\rho}-n}=\left.\frac{\partial_{J}^{n}\mathcal{W}_{{\rm
hs}}}{\partial_{J}^{n}\mathcal{W}_{\tr\\{{\rm hs}\\}}}\right|_{J=0}.$ (44)

like the way it is done for the formal process in equation (37).

For the evaluation of the integral in equation (42), we need to specify the
Wigner function of the input state $W_{\hat{\rho}}$. However, it is often
assumed that the beam splitter separates off only a small portion of the
state, so that $\zeta\ll 1$. In that case, we can expand equation (42) up to
first order in $\zeta$. The resulting expression then allows us to evaluate
the integral over $\beta$, without having to specify the input state. The
result reads

$\displaystyle\mathcal{W}_{{\rm hs}}\approx$
$\displaystyle{\textstyle\frac{1}{2}}\zeta\left[{\textstyle\frac{1}{2}}(1+J)\partial_{\alpha}\partial_{\alpha^{*}}+J\alpha\partial_{\alpha}+J\alpha^{*}\partial_{\alpha^{*}}\right]W_{\hat{\rho}}$
(45)
$\displaystyle+\left[1+{\textstyle\frac{1}{2}}(1+J)\zeta-(1-J)\zeta|\alpha|^{2}\right]W_{\hat{\rho}},$

It only allows single-photon detection at this level of the expansion, which
leads to

$\left.\partial_{J}\mathcal{W}_{{\rm
hs}}\right|_{J=0}\approx{\textstyle\frac{1}{2}}\zeta\left({\textstyle\frac{1}{2}}\partial_{\alpha}\partial_{\alpha^{*}}+\alpha\partial_{\alpha}+\alpha^{*}\partial_{\alpha^{*}}+1+2|\alpha|^{2}\right)W_{\hat{\rho}}.$
(46)

A comparison with the differential operation in equation (35) produced by the
formal process shows that, in the limit $\zeta\ll 1$, the two processes
produce the same state for single photon subtraction.

The detection of $n$ photons assumes that the detector is photon-number
resolving. However, it is also possible to use a non-photon-number-resolving
detector. The ensemble averaged state that would be produced by heralding the
detection of an arbitrary number of photons is obtained by removing the
probability to observe the vacuum. In terms of the generating function, it is
given by

$W_{{\rm npnr}}=\mathcal{N}\left(\left.\mathcal{W}_{{\rm
hs}}\right|_{J=1}-\left.\mathcal{W}_{{\rm hs}}\right|_{J=0}\right),$ (47)

where $\mathcal{N}$ is a required normalization factor. Under the assumption
that $\zeta\ll 1$, the state that is obtained from a non-photon-number-
resolving detection is the same as the one obtained in equation (46).

###  5.2 Heralded photon-subtracted squeezed thermal states

The generating function for heralded photon subtraction from a squeezed
thermal state is obtained by substituting the Wigner function for the squeezed
thermal state, given in equation (30), into equation (42) and evaluating the
integral over $\beta$. The resulting generating function reads

$\displaystyle\mathcal{W}_{{\rm hs-sts}}=$
$\displaystyle\frac{2T}{\sqrt{\mathcal{K}}}\exp\left[-\frac{(1-\zeta)T}{\mathcal{K}}\left(2|\alpha|^{2}A+\alpha^{2}B^{*}+\alpha^{*2}B\right)-\frac{\zeta}{\mathcal{K}}(1-J+T^{2}+T^{2}J)|\alpha|^{2}\right.$
(48)
$\displaystyle\left.-\frac{\left(1-J^{2}\right)\zeta^{2}}{2\mathcal{K}}|\alpha|^{2}(2TA-1-T^{2})\right].$

where

$\mathcal{K}\equiv
1+(1+J)(TA-1)\zeta-{\textstyle\frac{1}{4}}(1+J)^{2}(2TA-1-T^{2})\zeta^{2}.$
(49)

The generating function for the associated traces reads

$\mathcal{W}_{{\rm tr\\{hs-
sts\\}}}=\left[1+\frac{(1-J)(A-T)\zeta}{T}-\frac{(1-J)^{2}(2TA-1-T^{2})\zeta^{2}}{4T^{2}}\right]^{-1/2}.$
(50)

The expressions of the normalized Wigner functions of the photon-subtracted
states are more complicated for the heralded process than for the formal
process. In the limit of a small reflectivity $\zeta\rightarrow 0$, the
expression for one subtracted photon becomes the same as in equation (40). For
larger numbers of subtracted photons, the Wigner functions are qualitatively
similar, but slight quantitative differences appear, even in the limit of a
small reflectivity $\zeta\rightarrow 0$.

Hence, the experimental process produces similar results as the formal
definition of the process. Moreover, even for finite $\zeta$, the Wigner
functions of the heralded states have negative regions for arbitrary amounts
of squeezing.

For a quantitative comparison with the formal process, we plot the curves for
the heralded photon-subtracted squeezed vacuum states and the heralded photon-
subtracted squeezed thermal states with zero to five subtracted photons as
functions of $q$ along the line where $p=0$ in figures 4 and 5, respectively.
The assumed reflectivity of the beam splitter is $\zeta=0.2$. The other
parameters are the same as those for the curves in figures 2 and 3,
respectively. These curves have the same qualitative appearance as those for
the formal process, but they are quantitatively different. The effect of a
nonzero reflectivity is that it broadens the curves, reducing the effect of
the squeezing. By comparing figures 5 and 3, we also notice that a nonzero
reflectivity slightly reduces the sensitivity with respect to reduced purity
in the amplitudes of the oscillations produced by the photon subtractions.

Figure 4: Wigner functions of heralded photon-subtracted squeezed vacuum
states ($T=1$) with zero to five subtracted photons, plotted along the
squeezed direction, with $\zeta=0.2$. The input squeezed vacuum state has
$\langle{n}\rangle=1$, so that $A=3$ and $B=B^{*}=2.83$.

Figure 5: Wigner functions of heralded photon-subtracted squeezed thermal
states with zero to five subtracted photons, plotted along the squeezed
direction for $\zeta=0.2$, $T=0.9$, $A=3$ and $B=B^{*}=2.83$.

##  6 Squeezed Fock states

Under certain conditions, such as for small squeezing, the Wigner functions of
photon-subtracted squeezed vacuum states resemble those of squeezed Fock
states. We can use generating functions to investigate such squeezed Fock
states.

In this case, the generating function is not associated with a photon
subtraction process, but instead with the different Fock states. A generating
function for the normalized Wigner functions of the Fock states is given by

$\mathcal{W}_{{\rm
Fock}}=\frac{2}{1+\nu}\exp\left(-2\frac{1-\nu}{1+\nu}|\alpha|^{2}\right),$
(51)

where $\nu$ is the generating parameter. The Wigner functions of the Fock
states are produced by

$W_{\left|{n}\right\rangle}=\left.\frac{1}{n!}\partial_{\nu}^{n}\mathcal{W}_{{\rm
Fock}}\right|_{\nu=0}.$ (52)

The generating function for the Wigner functions of squeezing Fock states is
obtained by applying a Bogoliubov transformation, as discussed in section 3,
to the generating function for the Wigner functions of the Fock states in
equation (51). The result reads

$\mathcal{W}_{{\rm sq-
Fock}}=\frac{2}{1+\nu}\exp\left[-\frac{1-\nu}{1+\nu}\left(2|\alpha|^{2}A+\alpha^{2}B^{*}+\alpha^{*2}B\right)\right],$
(53)

where $A$ and $B$ are defined in equation (23).

The Wigner function of the squeezed single-photon Fock state is given by the
same expression as the single-photon-subtracted squeezed vacuum state in
equation (40) with $T=1$. For larger numbers of photons, the expressions are
different.

Figure 6: Wigner functions of squeezed Fock states with the same squeezing as
in figure 2, for zero to five initial photons (prior to squeezing), plotted
along $q$ for $p=0$.

In figure 6, we show the curves for the Wigner functions of squeezed Fock
states along the line where $p=0$ for the first six Fock states, including the
vacuum state. The amount of squeezing is such that the squeezed vacuum state
has $\langle{n}\rangle=1$. The curves in figure 6 show that the Wigner
functions of these squeezed Fock state produce qualitatively similar positive
and negative regions as found for the photon-subtracted squeezed vacuum states
shown in figure 2, but apart from the vacuum and single photon Fock states,
they are quantitatively different.

##  7 Photon statistics

###  7.1 Generating function

The properties of a quantum state can be revealed by measuring the photon
statistics of the state. The photon statistics is given by the probability
distribution of the photon number in the state. The projection operators,
their generating function and its Wigner function are provided in section 2.3.

To obtain a generating function for the photon statistics of a state, we
compute the trace of the state, multiplied by the generating function for the
photon number projection operators. In terms of Wigner functions, the
generating function for the photon statistics of the photon-subtracted states
is obtained by multiplying the generating function for the photon-subtracted
states by the generating function for the projection operators, given in
equation (19), but replacing $J\rightarrow K$ to have different generating
parameters for the photon subtractions and the photon statistics,
respectively. Then we integrate over $\alpha$ to obtain a double generating
function for both the number of subtracted photons and the photon statistics
of such a state.

###  7.2 Formal photon-subtracted squeezed thermal states

In the case of the formal photon subtraction process applied to the squeezed
thermal state, the generating function for the photon statistics is obtain by
multiplying equation (38) by equation (19) and integrate over $\alpha$. The
result is

$\displaystyle\mathcal{R}_{{\rm fs}}=$
$\displaystyle\frac{2T}{\sqrt{(1+K)^{2}T^{2}+2(1-K^{2})TA+(1-K)^{2}}}$ (54)
$\displaystyle\times\exp\left[-\frac{-2TK|\eta|^{2}A+T\eta^{2}B^{*}+T\eta^{*2}B-(1-K-T^{2}-T^{2}K)|\eta|^{2}}{(1+K)^{2}T^{2}+2(1-K^{2})TA+(1-K)^{2}}\right].$

It still needs to be normalized with the aid of equation (39).

The expression in equation (54) is a double generating function for the photon
statistics of a squeezed thermal state with an arbitrary number of subtracted
photons. The number of subtracted photons is determined by the number with
derivatives with respect to $\eta$ and $\eta^{*}$. The probability for
observing a given number of photons in the resulting state is determined by
the number with derivatives with respect to $K$.

This double generating function for the photon statistics can be used to
obtain a generating function for the average number of photons in a photon-
subtracted squeezed thermal state. It is given by

$\displaystyle\mathcal{M}_{{\rm fs}}$
$\displaystyle\equiv\left.\partial_{K}\mathcal{R}_{{\rm fs}}\right|_{K=1}$
(55)
$\displaystyle=\left[\frac{A-T}{2T}+\frac{A^{2}-1+(A-T)^{2}}{4T^{2}}|\eta|^{2}-\frac{A-T}{4T^{2}}\eta^{2}B^{*}-\frac{A-T}{4T^{2}}\eta^{*2}B\right]$
$\displaystyle~{}~{}~{}\times\exp\left(\frac{2|\eta|^{2}(A-T)-\eta^{2}B^{*}-\eta^{*2}B}{4T}\right).$

For the case with no subtracted photons, the average number of photons is

$\langle{n}\rangle=\frac{A-T}{2T}.$ (56)

For $A=1$, we obtain the relationship between the purity of a thermal state
and its average number of photons, and for $T=1$ we obtain a similar
relationship for the squeezed vacuum state.

###  7.3 Heralded photon-subtracted squeezed thermal states

Applying the calculation of the generating function for the photon statistics
to the generating function for the heralded photon-subtracted squeezed thermal
states, we obtain a double generating function given by

$\displaystyle\mathcal{R}_{{\rm hs}}$
$\displaystyle=2T\left\\{\left[1-(K-\zeta K+\zeta
J)^{2}\right]\left(2TA-1-T^{2}\right)\right.$ (57)
$\displaystyle~{}~{}~{}~{}\left.-2(J-K)\left(1-T^{2}\right)\zeta+2\left(1-K+T^{2}+T^{2}K\right)\right\\}^{-1/2}.$

The generating function for the average number of photons in a heralded
photon-subtracted squeezed thermal state is

$\displaystyle\mathcal{M}_{{\rm hs}}$
$\displaystyle\equiv\left.\partial_{K}\mathcal{R}_{{\rm hs}}\right|_{K=1}$
(58) $\displaystyle=\frac{2T(1-\zeta)\left[(1-\zeta+\zeta
J)^{2}\left(2TA-1-T^{2}\right)+1-T^{2}\right]}{\left\\{\left[1-(1-\zeta+\zeta
J)^{2}\right]\left(2TA-1-T^{2}\right)+2(1-J)\left(1-T^{2}\right)\zeta+4T^{2}\right\\}^{3/2}}.$

In the limit of small reflectivity, the average number of photons for no
subtracted photons is the same as given in equation (56).

##  8 Conclusions

Once the development of quantum information technology moves beyond the proof-
of-principle stage and the detailed design of quantum information systems
becomes more important, the analyses of such systems will need to address all
aspects of physical implementations. Here, we provide analytical tools to aid
such detailed quantitative analyses. Although, we consider only the particle-
number degrees of freedom here, generating functions can also be used in
analyses that incorporate the spatiotemporal degrees of freedom. Such analyses
are generally much more complex and the use of generating functions should
alleviate the complexity of the calculation in such comprehensive analyses.

The use of generating functions makes it possible to compute the Wigner
functions for an arbitrary number of photon subtractions. It allows different
models to be investigated, including the formal photon subtractions, squeezed
Fock states and the experimentally implemented heralding process to produce
such states through photon detection with or without photon-number-resolving
detectors.

The generating functions also allow one to compute the photon statistics of
such states as generating functions of generating functions. They facilitate
the direct calculations of the average number of photons in the states or the
variance and standard deviation in the photon number.

Generating functions for the marginal probability distributions can be
directly obtained via the integration of the generating functions for the
Wigner functions. The complex parameter $\alpha$ is expressed in terms of $q$
and $p$, as defined in equation (9). The resulting expression of the
generating function is then integrated over either $q$ or $p$ to obtain a
generating function for the marginal probability distributions as a function
of the remaining quadrature variable.

## Acknowledgement

This work was supported in part by funding from the National Research
Foundation of South Africa (Grant Numbers: 118532).

## References

## References

  * [1] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto.  Efficient classical simulation of continuous variable quantum information processes.  Phys. Rev. Lett., 88:097904, 2002.
  * [2] A. Mari and J. Eisert.  Positive Wigner functions render classical simulation of quantum computation efficient.  Phys. Rev. Lett., 109:230503, 2012.
  * [3] S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves.  Sufficient conditions for efficient classical simulation of quantum optics.  Phys. Rev. X, 6:021039, 2016.
  * [4] E. Chitambar and G. Gour.  Quantum resource theories.  Rev. Mod. Phys., 91:025001, 2019.
  * [5] S. D. Bartlett and B. C. Sanders.  Universal continuous variable quantum computation: requirement of optical nonlinearity for photon counting.  Phys. Rev. A, 65:042304, 2002.
  * [6] J. Eisert, S. Scheel, and M. B. Plenio.  Distilling Gaussian states with Gaussian operations is impossible.  Phys. Rev. Lett., 89:137903, 2002.
  * [7] A. I. Lvovsky, P. Grangier, A. Ourjoumtsev, V. Parigi, M. Sasaki, and R. Tualle-Brouri.  Production and applications of non-Gaussian quantum states of light.  arXiv preprint arXiv:2006.16985, 2020.
  * [8] A. Biswas and G. S. Agarwal.  Nonclassicality and decoherence of photon-subtracted squeezed states.  Phys. Rev. A, 75:032104, 2007.
  * [9] M. Walschaers, V. Parigi, and N. Treps.  Practical framework for conditional non-Gaussian quantum state preparation.  PRX Quantum, 1:020305, 2020.
  * [10] M. Dakna, T. Anhut, T. Opatrnỳ, L. Knöll, and D.-G. Welsch.  Generating Schrödinger-cat-like states by means of conditional measurements on a beam splitter.  Phys. Rev. A, 55:3184, 1997.
  * [11] A. Zavatta, S. Viciani, and M. Bellini.  Quantum-to-classical transition with single-photon-added coherent states of light.  Science, 306:660–662, 2004.
  * [12] P. R. Sharapova, G. Frascella, M. Riabinin, A. M. Pérez, O. V. Tikhonova, S. Lemieux, R. W. Boyd, G. Leuchs, and M. V. Chekhova.  Properties of bright squeezed vacuum at increasing brightness.  Phys. Rev. Research, 2:013371, 2020.
  * [13] S. L. Braunstein and P. Van Loock.  Quantum information with continuous variables.  Rev. Mod. Phys., 77:513–577, 2005.
  * [14] G. Adesso, S. Ragy, and A. R. Lee.  Continuous variable quantum information: Gaussian states and beyond.  Open Syst. Inf. Dyn., 21:1440001, 2014.
  * [15] F. S. Roux.  Evolution equation for multi-photon states in turbulence.  J. Phys. A: Math. Theor., 52:405301, 2019.
  * [16] F. S. Roux.  Quantifying entanglement of parametric down-converted states in all degrees of freedom.  Phys. Rev. Research, 2:023137, 2020.
  * [17] F. S. Roux.  Spatiotemporal effects on squeezing measurements.  Phys. Rev. A, 103:013701, 2021.
  * [18] H. J. Groenewold.  On the principles of elementary quantum mechanics.  Physica, 12:405–460, 1946.
  * [19] J. E. Moyal.  Quantum mechanics as a statistical theory.  Math. Proc. Camb. Philos. Soc., 45:99–124, 1949.
  * [20] T. L. Curtright and C. K. Zachos.  Quantum mechanics in phase space.  Asia Pacific Physics Newsletter, 1:37–46, 2012.