# Quantifying nonclassicality of mixed Fock states

Spencer Rogers<EMAIL_ADDRESS>Department of Physics, University of
Rhode Island, Kingston, RI 02881, USA  Tommy Muth  Department of Physics,
University of Rhode Island, Kingston, RI 02881, USA  Wenchao Ge  Department of
Physics, University of Rhode Island, Kingston, RI 02881, USA

###### Abstract

Nonclassical states of bosonic modes are important resources for quantum-
enhanced technologies. Yet, quantifying nonclassicality of these states, in
particular mixed states, can be a challenge. Here we present results of
quantifying the nonclassicality of a bosonic mode in a mixed Fock state via
the operational resource theory (ORT) measure [W. Ge, K. Jacobs, S. Asiri, M.
Foss-Feig, and M. S. Zubairy, Phys. Rev. Res. 2, 023400 (2020)], which relates
nonclassicality to metrological advantage. Generally speaking, evaluating the
ORT measure for mixed states is challenging, since it involves finding a
convex roof. However, we show that our problem can be reduced to a linear
programming problem. By analyzing the results of numerical optimization, we
are able to extract exact, analytical results for the case where three or four
neighboring Fock states have nonzero population. Interestingly, we find that
such a mode can be in distinct phases, depending on the populations. Lastly,
we demonstrate how our method is generalizable to density matrices of higher
ranks. Our findings suggests a viable method for evaluating nonclassicality of
arbitrary mixed bosonic states and potentially for solving other convex roof
optimization problems.

##  I Introduction

For bosonic modes in quantum optics and other analgous fields, the coherent
states, defined as eigenstates of the annihilation operator, are regarded as
the most classical [1]. Not only do they have minimal uncertainty and maintain
their shape under the harmonic oscillator Hamiltonian, but also they remain
separable from each other when subject to linear interactions (as occur in
beam-splitters). By contrast, superpositions of coherent states (e.g., cat and
squeezed states) may have highly nonclassical features, and with these the
potential to revolutionize such technologies as communication [2], computation
[3, 4], and metrology [5]. Indeed, on the metrological side, squeezed states
have already become an important tool in advanced gravitational wave
interferometry [6, 7].

Given the applications of nonclassicality, quantifying it, i.e., ranking the
nonclassicality of different bosonic states, is an important problem. Various
measures of nonclassicality have been proposed, including the nonclassicality
depth [8], nonclassical distance [9, 10], quantifications in terms of the
Schmidt rank [11, 12], and resource-theoretic measures [13, 14, 15]. These
measures satisfy the basic criteria that they are non-negative, and zero only
for classical states (i.e., coherent states or, if mixed, classical mixtures
of coherent states). They also characterize meaningful aspects of
nonclassicality, such as the minimum number of coherent superpositions in a
quantum state (in the case of the Schmidt rank) [11].

Here, we focus on the operational resource theory (ORT) measure of
nonclassicality [16, 17], which has two defining features. First, it is
operational, since it is related to the usefulness of a state for certain
tasks. For pure states, the ORT measure directly quantifies the metrological
power of a state in terms of the quantum Fisher information of quadrature
sensing [18], while for mixed states, the ORT measure is a tight upper bound
on this power; this is the closest possible relationship because some
nonclassical mixed states have zero metrological advantage over coherent
states. This feature endows the ORT measure an operational meaning to evaluate
different nonclassical states on the same footing. Second, the ORT measure is
resource-theoretic. A quantum resource theory [19] defines free states and
resource states via the operations, called free operations, under which it is
impossible to increase that resource. In the ORT measure of nonclassicality,
coherent states, and classical mixtures of them, are considered the free
states, while nonclassical states are considered resource states. The free
operations are classical operations (meaning operations that cannot create
superpositions of coherent states from mixtures of coherent states) — Ref.
[16] rigorously showed that the ORT measure does not increase under such
operations.

Previous works computed the ORT measure, from here on also referred to as the
nonclassicality $\mathcal{N}$, for various pure bosonic states of interest
[16, 17, 20]. It has been shown, for example, that squeezed vacuum has the
highest nonclassicality per unit average energy (which speaks to the
metrological usefulness of squeezed vacuum) [16]. Cat states whose
nonclassicality per unit energy approach that of squeezed vacuum in the
asymptotic limit of large energy have also been discovered [17], thus
extending our knowledge of optimal sensing with Mach-Zehnder interferometry
[21, 16].

While pure states are ideal for applications, mixed states are inevitable in
practice, due to coupling to the environment. For example, mixtures of Fock
states can be generated under a dephasing channel in various bosonic systems
[22, 23, 24]. Computing the ORT measure for a mixed state $\hat{\rho}$,
however, poses a significant challenge. This is because the ORT measure for
mixed states is a convex roof construction, similar to entanglement measures
for mixed states [25, 26] — computing it involves extremizing a function over
all the decompositions of $\hat{\rho}$. Previous calculations only address
very simple classes of states [17]. For more general states, only a lower
bound has been given: $\mathcal{N}(\hat{\rho})\geq\mathcal{W}(\hat{\rho})$,
where $\mathcal{W}(\hat{\rho})$ is the metrological power (to be defined
later) [17]. In this work, we present a numerical method of linear programming
for calculating the ORT measure for mixed Fock states. Based on the numerical
results, we extract analytical ansatz of the ORT measure for rank-3 and rank-4
mixed states and identify distinct phases depending on the populations. Our
method provides a viable way to evaluate arbitray mixed states of bosonic
systems.

Our paper is structured as follows: in section II, we review the ORT measure
of nonclassicality for pure and mixed states. In section III, we explain our
method, based on linear programming, for calculating the ORT measure for mixed
states. We focus on one class of mixed states in particular: those that are
diagonal in the Fock basis. In section IV, we apply our method to rank-3 and
rank-4 states, and show that such states can be in various phases, depending
on their populations. The results here are largely analytical, as we were able
to infer exact values from the numerical results. For higher-rank states,
numerical results are possible, as we demonstrate on truncated thermal states.
Section V contains our concluding remarks.

##  II The ORT Measure

###  II.1 Pure states

A pure bosonic state $\ket{\psi}$ is defined as nonclassical if and only if it
is not a coherent state, i.e., $\hat{a}\ket{\psi}\neq\alpha\ket{\psi}$ for all
complex numbers $\alpha$ [1]. The ORT measure $\mathcal{N}$ assigns a
nonclassicality value to each pure state as [16, 17]:

$\begin{split}\mathcal{N}(\ket{\psi})&=\max_{\mu}\bra{\psi}(\Delta\hat{X}_{\mu})^{2}\ket{\psi}-\frac{1}{2}\\\
&=\langle\hat{a}^{\dagger}\hat{a}\rangle-\left|\langle\hat{a}\rangle\right|^{2}+\left|\langle\hat{a}^{2}\rangle-\langle\hat{a}\rangle^{2}\right|.\end{split}$
(1)

Here, $(\Delta\hat{X}_{\mu})^{2}$ is the variance of a quadrature
$\hat{X}_{\mu}=i(e^{-i\mu}\hat{a}^{\dagger}-e^{i\mu}\hat{a})/\sqrt{2}$,
$\hat{a}$ is the annihilation operator for the mode and $\mu\in[0,2\pi]$.

Whereas the second line in Eq. (1) is generally more convenient for
calculations, the first line makes the connection to metrological power more
apparent. In quantum metrology with unitary encoding, one measures a classical
parameter $\theta$ by applying a transformation
$\hat{U}(\theta)=e^{-i\theta\hat{G}}$ to a quantum sensor, where the generator
$\hat{G}$ is some Hermitian observable of the quantum sensor.
$\hat{U}(\theta)$ may represent the effect of coupling the quantum sensor to a
classical system, and may or may not describe some interaction picture.
Regardless, the effective Hamiltonian for the sensor is $\hat{G}$, while
$\theta$ acts analogously to time. Assuming the sensor is prepared in a pure
state $\ket{\psi}$, the energy-time uncertainty principle tells us that, the
larger the variance $\bra{\psi}(\Delta\hat{G})^{2}\ket{\psi}$, the faster the
sensor will evolve with respect to $\theta$. So, the larger the variance
$\bra{\psi}(\Delta\hat{G})^{2}\ket{\psi}$, the more responsive the sensor will
be to small variations in $\theta$. Indeed, the quantum Fisher information
(QFI) [27] for this protocol, which quantifies the precision with which
$\theta$ can be measured, is simply this variance
$F_{\hat{G}}(\ket{\psi})=\bra{\psi}(\Delta\hat{G})^{2}\ket{\psi}$.111In this
paper, we drop a factor of $4$ from the standard definition of the QFI. From
the first line of Eq. (1), we see that the ORT measure for pure states is
based on the maximum quadrature variance — in other words, the Fisher
information when using the sensor’s optimal quadrature as the generator. Since
the maximum quadrature variance of a coherent state is always $\frac{1}{2}$,
$\mathcal{N}(\ket{\alpha})=0$ for coherent states $\ket{\alpha}$. For
nonclassical states, the maximum quadrature variance is always greater than
$\frac{1}{2}$, and $\mathcal{N}>0$ describes the sensing enhancement over
coherent states. We refer to this sensing enhancement as the metrological
power $\mathcal{W}$ [16, 18]. For pure states,
$\mathcal{N}(\ket{\psi})=\mathcal{W}(\ket{\psi})$.

One may wonder why the ORT measure for pure states is based only on the
maximum variance quadrature, as opposed to the maximum variance Hermitian
generator $\hat{G}$ in general (quadrature or not). This has to do with the
resource-theoretic nature of the ORT measure. A unitary
$e^{-i\theta\hat{X}_{\mu}}$ is a free operation for classical states [13, 16]—
when it acts on a coherent state, it displaces it in phase space, yielding a
new coherent state. Another free operator for evaluating nonclassicality is to
use a phase shifter, i.e., $\hat{G}=\hat{a}^{\dagger}\hat{a}$. It is shown
that the generator turns into an effective quadrature operator $\hat{X}_{\mu}$
in a multi-port phase sensing scheme with a single-mode nonclassical input in
combination with classical light [18].

###  II.2 Mixed states

Before introducing the ORT measure for mixed states, one should first clarify
what it means for a mixed state to be classical or nonclassical. A mixed
bosonic state is defined as classical if and only if its density matrix
$\hat{\rho}$ can be regarded as a classical mixture of coherent states.
Otherwise, the state is considered nonclassical. Equivalently, one may
determine if a state is classical by examining its Glauber-Sudarshan
$P$-function. Any bosonic state $\hat{\rho}$ can be represented using its
Glauber-Sudarshan $P$-function [28, 29] as

$\hat{\rho}=\int P(\alpha,\alpha^{*})\ket{\alpha}\bra{\alpha}d^{2}\alpha.$ (2)

The state $\hat{\rho}$ is defined as classical if $P(\alpha,\alpha^{*})$ is
positive semidefinite, in which case $P(\alpha,\alpha^{*})$ acts as a
classical probability density over the coherent states $\ket{\alpha}$.
Otherwise, the state is nonclassical [30, 31].

The ORT measure of nonclassicality for mixed states is:

$\begin{split}\mathcal{N}(\hat{\rho})&=\min_{\\{q_{j},\ket{\phi_{j}}\\}}\left[\max_{\mu}\sum_{j}q_{j}\bra{\phi_{j}}(\Delta\hat{X}_{\mu})^{2}\ket{\phi_{j}}\right]-\frac{1}{2}\\\
&=\min_{\\{q_{j},\ket{\phi_{j}}\\}}\left[\sum_{j}q_{j}\left(\bar{n}_{j}-|\bar{\alpha}_{j}|^{2}\right)+\left|\sum_{j}q_{j}\left(\bar{\xi}_{j}-\bar{\alpha}_{j}^{2}\right)\right|\right].\end{split}$
(3)

Here, we have defined
$\bar{n}_{j}\equiv\bra{\phi_{j}}\hat{a}^{\dagger}\hat{a}\ket{\phi_{j}}$,
$\bar{\alpha}_{j}\equiv\bra{\phi_{j}}\hat{a}\ket{\phi_{j}}$, and
$\bar{\xi}_{j}=\bra{\phi_{j}}\hat{a}^{2}\ket{\phi_{j}}$. The minimization is
over all possible decompositions (ensembles) $\\{q_{j},\ket{\phi_{j}}\\}$ of
$\hat{\rho}$, satisfying
$\hat{\rho}=\sum_{j}q_{j}\ket{\phi_{j}}\bra{\phi_{j}}$ with $q_{j}>0$. If
$\hat{\rho}$ is mixed (i.e., $\textrm{rank}(\hat{\rho})>1$), then the number
of states in a decomposition of $\hat{\rho}$ may be any number larger than
$\textrm{rank}(\hat{\rho})$, making the calculation of
$\mathcal{N}(\hat{\rho})$ a highly complex optimization problem in general,
hence this paper.

The definition of $\mathcal{N}(\hat{\rho})$ satisfies several important
conditions, as proved in Ref. [16]. (i) Non-negativity:
$\mathcal{N}(\hat{\rho})\geq 0$, where equality holds if and only if
$\hat{\rho}$ is classical. Non-negativity is the minimum requirement for a
nonclassicality measure. (ii) Weak monotonicity: $\mathcal{N}$ cannot increase
under any classical operation $\Lambda$:
$\mathcal{N}(\Lambda[\hat{\rho}])\leq\mathcal{N}(\hat{\rho})$. A classical
operation is an operation that cannot create superpositions of coherent states
from mixtures of coherent states (such as the use of passive linear optical
operations and displacements). For details, see Ref. [16, 17]. Weak-
monotonicity makes the ORT measure resource-theoretic; classical operations
are the free operations of the resource theory. (iii) Convexity:
$\sum_{j}p_{j}\mathcal{N}(\hat{\rho}_{j})\geq\mathcal{N}\left(\sum_{j}p_{j}\hat{\rho}_{j}\right)$
for any quantum states $\hat{\rho}_{j}$ and probabilities $p_{j}$. (iv) Lower-
bounded by metrological power:
$\mathcal{N}(\hat{\rho})\geq\mathcal{W}(\hat{\rho})$, where equality holds for
pure states. The metrological power for mixed states is defined as
$\mathcal{W}(\hat{\rho})\equiv\max[F_{X}(\hat{\rho})-1/2,0]$, where [16, 17]

$F_{X}(\hat{\rho})=\max_{\mu}\left[\min_{\\{q_{j},\ket{\phi_{j}}\\}}\left(\sum_{j}q_{j}\bra{\phi_{j}}(\Delta\hat{X}_{\mu})^{2}\ket{\phi_{j}}\right)\right]$
(4)

is the QFI for the optimal quadrature angle. Notice that for mixed states, the
QFI also involves a minimization over all possible ensembles (and maximization
over quadratures), although the order of the maximization and minimization in
Eq. (3) and (4) makes a subtle yet important difference (in fact, the QFI can
be computed given an eigenstate decomposition of $\hat{\rho}$ [32]). Again,
the metrological power $\mathcal{W}$ describes the sensing enhancement over
classical states: for all classical states, $\mathcal{W}=0$. Notably, there
exist nonclassical mixed states for which $\mathcal{W}=0$ [16, 17]. Thus, the
tight inequality $\mathcal{N}(\hat{\rho})\geq\mathcal{W}(\hat{\rho})$ is the
closest possible relationship between a resource-theoretic nonclassicality
measure and the metrological power.

##  III Method

The ORT measure for mixed states, Eq. (3), is an example of a convex roof
construction. Such constructions commonly arise in other nonclassicality
measures [11, 15, 14], as well as entanglement measures, for mixed states [25,
26].

Evaluating the convex roof is a matter of optimizing an objective function
over the set of all convex decompositions of a state $\hat{\rho}$. Optimizing
the objective function is a challenging problem, a primary difficulty being
that the number of states in a decomposition of a mixed state may be any
number larger than its rank, although some clever techniques that bound or
approximate the answer have been proposed [25, 33].

Rather than adapt previous techniques to our problem, here we introduce a new
method based on linear programming. Our method is similar to variational
methods for ground state energy calculations, in that we essentially make a
variational ansatz over the space of decompositions. It has the pedagogical
advantage that it is highly intuitive and easy to employ.

###  III.1 Fock-diagonal states

In particular, we utilize our method to evaluate the ORT measure for mixed
states that are diagonal in the Fock basis:
$\hat{\rho}=\sum_{n}p_{n}\ket{n}\bra{n}$, where $\ket{n}$ is the Fock state
with photon number $n$ and $\\{p_{n},\ket{n}\\}$ is one possible
decomposition. Such states arise often due to dephasing effects [22, 23, 24],
which eliminate the coherence between different Fock states. We refer to the
probabilities $p_{n}$ as populations. For such states, the ORT measure
simplifies to:

$\begin{split}\mathcal{N}(\hat{\rho})&=\langle{\hat{a}^{\dagger}\hat{a}}\rangle-\max_{\\{q_{j},\ket{\phi_{j}}\\}}\left(\sum_{j}q_{j}\left|\bar{\alpha}_{j}\right|^{2}-\left|\sum_{j}q_{j}\bar{\alpha}_{j}^{2}\right|\right)\\\
&\leq\langle{\hat{a}^{\dagger}\hat{a}}\rangle.\end{split}$ (5)

The inequality is saturated by states in which all nonzero populations have
associated photon numbers that differ by at least $2$ from one another [17].
An example is $p\ket{0}\bra{0}+(1-p)\ket{2}\bra{2}$, whose population
constraints imply that any decomposition state has the form
$\ket{\phi_{j}}=\ket{0}\bra{0}\ket{\phi_{j}}+\ket{2}\bra{2}\ket{\phi_{j}}$ —
it follows that $\alpha_{j}=\bra{\phi_{j}}\hat{a}\ket{\phi_{j}}=0$ for all
$\ket{\phi_{j}}$. It is thus more worthwhile to evaluate the ORT measure for
Fock-diagonal states where neighboring photon numbers have nonzero
populations.

Prior ORT measure calculations [17] only addressed Fock-diagonal states with
two neighboring populations:
$\hat{\rho}_{2F}=p_{n+1}\ket{n+1}\bra{n+1}+(1-p_{n+1})\ket{n}\bra{n}$. It was
found that the optimal decomposition for such states was a “four-prong” set of
superposition states:
$\ket{\phi_{j}}=\sqrt{p_{n+1}}\ket{n+1}+\sqrt{1-p_{n+1}}e^{ij\pi/2}\ket{n}$
with equal weights $q_{j}=1/4$ ($j=0,1,2,3$).222Actually, there are multiple
optimal decompositions (for example a “three-prong set” with relative phases
that are the cubic roots of unity), but they share the property that
$\left|\bra{n+1}\ket{\phi_{j}}\right|=\sqrt{p_{n+1}}$ for all
$\ket{\phi_{j}}$. The phases $e^{ij\pi/2}$, which are the quartic roots of
unity, cause $\sum_{j}q_{j}\bar{\alpha}_{j}^{2}=0$ — thus the absolute value
term in Eq. (5) vanishes with this decomposition. Plugging this decomposition
into Eq. (5) gives:
$\mathcal{N}(\hat{\rho}_{2F})=n+p_{n+1}-(n+1)p_{n+1}(1-p_{n+1})$.333That this
value is optimal may be further justified via its convexity as a function of
$p_{n+1}$. This result serves as a basic check of our method, and provides
some insight into the types of optimal decompositions one may expect in more
complicated cases.

###  III.2 Linear-programming method

Now we will describe our method for calculating the ORT measure. For
simplicity, we will first explain it in the context of calculating
$\mathcal{N}(\hat{\rho}_{2F})$. Then we will explain how to apply the method
to Fock-diagonal states with more nonzero populations.

Any state occurring in a decomposition of $\hat{\rho}$ must be in the support
of $\hat{\rho}$. Thus, for $\hat{\rho}_{2F}$, each decomposition state must be
of the form
$\ket{\phi(x,\theta)}=x\ket{n+1}+\sqrt{1-x^{2}}e^{i\theta}\ket{n}$, where
$x\in[0,1]$ and $\theta\in[0,2\pi]$. Just as there are infinitely many
decompositions of $\hat{\rho}_{2F}$, there are infinitely many probability
distributions $q(x,\theta)$ such that $\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
q(x,\theta)\ket{\phi(x,\theta)}\bra{\phi(x,\theta)}=\hat{\rho}_{2F}$. We
prefer to think of the problem of finding the optimal decomposition as a
problem of finding the optimal probability distribution $q(x,\theta)$. Each
$q(x,\theta)$ satisfies the constraints:

$\displaystyle\int_{0}^{1}dx\int_{0}^{2\pi}d\theta q(x,\theta)$
$\displaystyle=1$ (6) $\displaystyle\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
x^{2}q(x,\theta)$ $\displaystyle=p_{n+1}$ (7)
$\displaystyle\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
x\sqrt{1-x^{2}}e^{i\theta}q(x,\theta)$ $\displaystyle=0$ (8) $\displaystyle
q(x,\theta)$ $\displaystyle\geq 0.$ (9)

Respectively, the above constraints address normalization, population
constraints, lack of off-diagonal elements, and positivity of all
contributions (i.e., that the distribution describes a convex combination of
states). Although we express $q(x,\theta)$ as if it were a continuous
probability distribution here, the optimal distribution (decomposition) may be
discrete. In terms of $q(x,\theta)$, Eq. (3) reads:

$\begin{split}\mathcal{N}(\hat{\rho}_{2F})-\langle{\hat{a}^{\dagger}\hat{a}}\rangle=-\max_{q(x,\theta)}\Bigg{(}&\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
q(x,\theta)\left|\bar{\alpha}_{x,\theta}\right|^{2}\\\
-&\left|\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
q(x,\theta)\bar{\alpha}_{x,\theta}^{2}\right|\Bigg{)}\end{split}$ (10)

where $\bar{\alpha}_{x,\theta}=x\sqrt{1-x^{2}}\sqrt{n+1}e^{-i\theta}$. The
quantity to be maximized on the right hand side of 10 is called the objective
function $\mathcal{L}[q(x,\theta)]$:

$\begin{split}\mathcal{L}[q(x,\theta)]&=\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
q(x,\theta)\left|\bar{\alpha}_{x,\theta}\right|^{2}\\\
&-\left|\int_{0}^{1}dx\int_{0}^{2\pi}d\theta
q(x,\theta)\bar{\alpha}_{x,\theta}^{2}\right|.\end{split}$ (11)

The second term in the objective function may appear problematic, due to the
absolute value. However, without loss of generality, it can be taken to vanish
(for our class of states). This is because, given any $q(x,\theta)$, it is
always possible to find a distribution $\tilde{q}(x,\theta)$ such that
$\mathcal{L}[\tilde{q}(x,\theta)]\geq\mathcal{L}[q(x,\theta)]$ and
$\int_{0}^{1}dx\int_{0}^{2\pi}d\theta\tilde{q}(x,\theta)\bar{\alpha}_{x,\theta}^{2}=0$.
The idea is to reallocate the probability from $q(x,\theta)$ in a four-prong
formation for each $x$:
$\tilde{q}(x,\theta)=\frac{1}{4}\left(\int_{0}^{2\pi}d\theta^{\prime}q(x,\theta^{\prime})\right)\left(\sum_{j=0}^{3}\delta(\theta-\theta_{0}-\frac{j\pi}{2})\right)$,
where $\theta_{0}\in(0,2\pi)$ is arbitrary. It is straightforward to check
that $\tilde{q}(x,\theta)$ satisfies the same constraints (Eq. (6)- (9)) as
the $q(x,\theta)$ from which it was constructed; meanwhile, it causes the
absolute term in the objective function $\mathcal{L}$ to vanish, while
maintaining the value of the remaining term.

Since the $\theta$-dependence of $\tilde{q}(x,\theta)$ is removed, we may
rename it $\mathcal{Q}(x)\equiv\tilde{q}(x,\theta)$. For our purposes, it is
sufficient to optimize the $x$-dependence of $\mathcal{Q}(x)$, subject to the
constraints $\int_{0}^{1}dx\mathcal{Q}(x)=1$,
$\int_{0}^{1}dxx^{2}\mathcal{Q}(x)=p_{n+1}$, and $\mathcal{Q}(x)\geq 0$. In
terms of $\mathcal{Q}(x)$, Eq. (10) becomes:

$\mathcal{N}(\hat{\rho}_{2F})-\langle{\hat{a}^{\dagger}\hat{a}}\rangle=-\max_{\mathcal{Q}(x)}\left(\int_{0}^{1}dx\mathcal{Q}(x)x^{2}(1-x^{2})(n+1)\right).$
(12)

The optimization problem (including constraints) resembles a linear
programming problem [34] for the vector $\mathcal{Q}(x)$, the only caveat
being that $\mathcal{Q}(x)$ does not have a discrete number of components.
Nevertheless, the solution to the problem can be approximated by discretizing
the set of possible $x$ (i.e., binning):

$\begin{split}-\max_{\mathcal{Q}_{x}}\left(\sum_{x\in\mathcal{X}}\mathcal{Q}_{x}x^{2}(1-x^{2})(n+1)\right)&\gtrsim\mathcal{N}(\hat{\rho}_{2F})-\langle{\hat{a}^{\dagger}\hat{a}}\rangle\\\
\end{split}$ (13) $\displaystyle\sum_{x\in\mathcal{X}}\mathcal{Q}_{x}$
$\displaystyle=1$ (14)
$\displaystyle\sum_{x\in\mathcal{X}}x^{2}\mathcal{Q}_{x}$
$\displaystyle=p_{n+1}$ (15) $\displaystyle\mathcal{Q}_{x}$ $\displaystyle\geq
0.$ (16)

Here, $\mathcal{X}$ is a discrete set of points in the range $[0,1]$. In this
work, we take $\mathcal{X}$ to be the set of real numbers $m\Delta\leq 1$,
where $m$ is a non-negative integer and $\Delta$ is the spacing of the binning
with $0<\Delta\ll 1$). The approximation, Eq. (13), should converge as
$\Delta\rightarrow 0$ (here $\gtrsim$ means approximately equal to, but never
less than).

Figure 1: Optimal discrete probability distribution $\mathcal{Q}_{x}$ for
$\hat{\rho}_{2F}=0.16\ket{n+1}\bra{n+1}+0.84\ket{n}\bra{n}$ (the value of $n$
has no effect here). The spacing between $x$ bins is $\Delta=0.01$. One may
notice that $\mathcal{Q}_{\sqrt{p_{n+1}}}=1$ and $\mathcal{Q}_{x}=0$ for all
other $x$, thus reproducing the expected result from Ref. [17].

Here $\mathcal{Q}_{x}$ is a discrete probability distribution (“histogram”)
over the set of points in $\mathcal{X}$. As such, it is essentially a
variational ansatz for the optimal decomposition, which can be optimized by
linear programming. We solve the linear program using the LinearOptimization
function in Mathematica. An example is demonstrated in Fig. 1.
$\mathcal{Q}_{x}$ is spiked at $x=\sqrt{p_{n+1}}$ and zero elsewhere. This
holds for all choices of $p_{n+1}$ (though one should be mindful and bin
appropriately in the extremes $p_{n+1}\ll 1$ and $1-p_{n+1}\ll 1$). Thus, the
optimal decomposition utilizes the four-prong set of states
$\ket{\phi_{j}}=\sqrt{p_{n+1}}\ket{n+1}+\sqrt{1-p_{n+1}}e^{i(j\pi/2+\theta_{0})}$
(where $j=0,1,2,3$), with probability $1/4$ for each. This corroborates our
method, since it reproduces the results of Ref. [17].

Note that, due to a careful choice of parameters (and binning) for Fig. 1, the
discrete distribution matches the true optimal decomposition perfectly. This
will not always be the case (generally, the discrete distribution will be
slightly suboptimal, acting as a bound). However, one can work backwards from
the numerical evidence to conjecture the true optimal decomposition (i.e.
realize that $0.4=\sqrt{0.16}$, in Fig. 1). The conjecture of the ansatz can
then be corroborated by checking it against finer binnings. In this way, we
obtain analytic expressions for $\mathcal{N}(\hat{\rho}_{3F})$ and
$\mathcal{N}(\hat{\rho}_{4F})$, in section IV.

###  III.3 Higher-rank Fock-diagonal states

A more general Fock-diagonal state, with neighboring populations, is:

$\hat{\rho}_{MF}=\sum_{k=0}^{M-1}p_{n+k}\ket{n+k}\bra{n+k}.$ (17)

Any decomposition state for $\hat{\rho}_{MF}$ must be of the form:

$\ket{\phi\left(\vec{x},\vec{\theta}\right)}=\sum_{k=0}^{M-1}x_{k}e^{i\theta_{k}}\ket{n+k},$
(18)

where $x_{k}\in[0,1]$, $\theta_{k}\in[0,2\pi]$, and $\sum_{k}x_{k}^{2}=1$.
Without loss of generality, one can take $\theta_{M-1}=0$, and
$x_{0}=\sqrt{1-\sum_{k=1}^{M-1}x_{k}^{2}}$. $\hat{\rho}_{MF}$ may then be
expressed as a summation over these states with probabilities
$q(\vec{x},\vec{\theta})$ satisfying certain constraints. The nonclassicality
calculation then becomes a matter of optimizing $q(\vec{x},\vec{\theta})$.
Similar to the case of $\hat{\rho}_{2F}$, it is sufficient to optimize a
probability distribution $\mathcal{Q}(\vec{x})$ over the amplitudes $x_{k}$
with symmetric $\vec{\theta}$-dependence. The optimization problem can be
solved approximately by discretizing the set of possible $\vec{x}$:

$\begin{split}-\max_{\mathcal{Q}_{\vec{x}}}\left(\sum_{\vec{x}\in\mathcal{X}}\mathcal{Q}_{\vec{x}}\left(\sum_{k=0}^{M-2}x_{k+1}x_{k}\sqrt{n+k+1}\right)^{2}\right)\\\
\gtrsim\mathcal{N}(\hat{\rho}_{MF})-\langle\hat{a}^{\dagger}\hat{a}\rangle\end{split}$
(19)

$\displaystyle\sum_{\vec{x}\in\mathcal{X}}\mathcal{Q}_{\vec{x}}$
$\displaystyle=1$ (20)
$\displaystyle\sum_{\vec{x}\in\mathcal{X}}x_{k}^{2}\mathcal{Q}_{\vec{x}}$
$\displaystyle=p_{n+k}$ (21) $\displaystyle\mathcal{Q}_{\vec{x}}$
$\displaystyle\geq 0.$ (22)

In this work, we take
$\left(x_{1},x_{2},...,x_{M-1}\right)=\left(l_{1}\Delta,l_{2}\Delta,...,l_{M-1}\Delta\right)$,
where $\sum_{k=1}^{M-1}l_{k}^{2}\Delta^{2}\leq 1$, the $l_{k}$ are non-
negative integers, and $0<\Delta\ll 1$. The approximation, Eq. (19) should
converge as $\Delta\rightarrow 0$. Note that $\mathcal{Q}_{\vec{x}}$ is
essentially a vector in $\mathbb{R}^{d}$, where $d$ is the number of elements
in $\mathcal{X}$. For fixed $\Delta$, $d$ scales as $\Delta^{-M+1}$. The
optimization may be carried out using linear programming, as before, although
the computational complexity will increase with $M$.

Note that an upper bound on the ORT measure can easily be found by plugging in
the simple decomposition
$\mathcal{Q}_{\vec{x}}=\prod_{k=1}^{M-1}\delta_{x_{k},\sqrt{p_{k}}}$:

$\mathcal{N}(\hat{\rho}_{MF})\leq\langle\hat{a}^{\dagger}\hat{a}\rangle-\left(\sum_{k=0}^{M-2}\sqrt{p_{k+1}p_{k}(n+k+1)}\right)^{2}.$
(23)

Equality holds in the $M=2$ case (where this simple decomposition is optimal),
as demonstrated previously. For $M>2$, we find (see section IV) that equality
holds in some regimes (i.e., some choices of the populations $p_{k}$), but not
all.

##  IV Numerical Results and Analytical ansatz

###  IV.1 Rank-3 case

For fixed $n$ in
$\hat{\rho}_{3F}=p_{n+2}\ket{n+2}\bra{n+2}+p_{n+1}\ket{n+1}\bra{n+1}+(1-p_{n+2}-p_{n+1})\ket{n}\bra{n}$,
there are two free parameters: $p_{n+2}$ and $p_{n+1}$. The possible choices
$(p_{n+2},p_{n+1})$ of these parameters form the 45-45-90 triangle defined by
$p_{n+2}+p_{n+1}\leq 1$, $p_{n+2}\geq 0$, and $p_{n+1}\geq 0$. From our
numerical investigations, we deduced that there are three distinct regimes, or
phases, within this triangular phase space (see Fig. 2).

Figure 2: Phase space for $\hat{\rho}_{3F}$ assuming the lowest populated
photon number is $n=0$. The expression for the ORT measure
$\mathcal{N}(\hat{\rho}_{3F})$ depends on which of three distinct phases the
state is in (although $\mathcal{N}(\hat{\rho}_{3F})$ is still piecewise
continuous as shown in Fig. 3). These phases are named according to what their
optimal decompositions look like. For $n>0$, the general shape of these phases
is maintained, although the point where the curves intersect the bottom axis
shifts left as $p_{n+2}=\frac{2+n}{3+2n}$.

####  IV.1.1 Triplet phase

In Fig. 2, the triplet phase is the upper, gray regime (the lower boundary of
which is given in ensuing subsections). Each state in the optimal
decomposition of $\hat{\rho}_{3F}$ is of the form:

$\begin{split}\ket{\phi_{T}(\theta_{j},\gamma_{j})}&=\sqrt{p_{n+2}}\ket{n+2}+\sqrt{p_{n+1}}e^{i\theta_{j}}\ket{n+1}\\\
&+\sqrt{1-p_{n+2}-p_{n+1}}e^{i\gamma_{j}}\ket{n}.\end{split}$ (24)

As stated before, the decomposition states come in a four-prong set
($j=0,1,2,3$) with equal weights and the relative phases $\theta_{j}$ and
$\gamma_{j}$ is determined by the choice of the phase of the coherence
$\alpha_{j}$. For context, this decomposition is a straightforward
generalization of the optimal decomposition for density matrices that are
diagonal in two neighboring Fock states,
$\hat{\rho}_{2F}=p_{n+1}\ket{n+1}\bra{n+1}+p_{n}\ket{n}\bra{n}$ (see section
III and [17]). The ORT measure for states in this regime is:

$\begin{split}&\mathcal{N}(\hat{\rho}_{3F})=2p_{n+2}+p_{n+1}+n-\\\
&\left(\sqrt{p_{n+2}}\sqrt{n+2}+\sqrt{1-p_{n+2}-p_{n+1}}\sqrt{n+1}\right)^{2}p_{n+1}.\end{split}$
(25)

In this case, the bound in Eq. (23) is saturated.

What stands out about this decomposition is the fact that
$\left|\bra{n+2}\ket{\phi_{T}(\theta_{j},\gamma_{j})}\right|^{2}=p_{n+2}$,
$\left|\bra{n+1}\ket{\phi_{T}(\theta_{j},\gamma_{j})}\right|^{2}=p_{n+1}$, and
$\left|\bra{n}\ket{\phi_{T}(\theta_{j},\gamma_{j})}\right|^{2}=1-p_{n+2}-p_{n+1}$.
We call this decomposition the “triplet” decomposition (and the phase where it
is optimal the triplet phase) because of this symmetric treatment of the three
photon numbers.

Of course, it is always possible to decompose
$\hat{\rho}_{3F}=p_{n+2}\ket{n+2}\bra{n+2}+p_{n+1}\ket{n+1}\bra{n+1}+p_{n}\ket{n}\bra{n}$
into the four states $\ket{\phi_{T}(\theta_{j},\gamma_{j})}$. However, it
turns out that the triplet decomposition is not optimal for all choices of
$(p_{n+2},p_{n+1})$. For some values of $(p_{n+2},p_{n+1})$, a more optimal
decomposition is possible.

####  IV.1.2 Upper-pair phase

In Fig. 2, the upper-pair phase is the lower-left, blue regime. Five states
are a part of the optimal decomposition of $\hat{\rho}_{3F}$. One is the Fock
state $\ket{n}$. The other four states are of the form:

$\begin{split}&\ket{\phi_{U}(\theta_{j},\gamma_{j})}=\sqrt{1-f}e^{i\gamma_{j}}\ket{n}\\\
&+\sqrt{f}\left(\sqrt{\frac{p_{n+2}}{p_{n+2}+p_{n+1}}}\ket{n+2}+\sqrt{\frac{p_{n+1}}{p_{n+2}+p_{n+1}}}e^{i\theta_{j}}\ket{n+1}\right),\end{split}$
(26)

where:

$f=\frac{(2+n)p_{n+2}}{(1+n)p_{n+1}+(3+2n)p_{n+2}}.$ (27)

The states $\ket{\phi_{U}(\theta_{j},\gamma_{j})}$ contribute a total
probability $\left(p_{n+2}+p_{n+1}\right)f^{-1}$ (each contributes $1/4$ of
this). The ORT measure is:

$\begin{split}\mathcal{N}(\hat{\rho}_{3F})&=2p_{n+2}+p_{n+1}+n\\\
&-(p_{n+2}+p_{n+1})f^{-1}\left|\bra{\phi_{U}(0,0)}\hat{a}\ket{\phi_{U}(0,0)}\right|^{2}.\end{split}$
(28)

What stands out about this decomposition is the fact that
$\left|\bra{n+2}\ket{\phi_{U}(\theta_{j},\gamma_{j})}\right|^{2}/\left|\bra{n+1}\ket{\phi_{U}(\theta_{j},\gamma_{j})}\right|^{2}=p_{n+2}/p_{n+1}$,
the population ratio between $\ket{n+2}$ and $\ket{n+1}$ in the original
density matrix $\hat{\rho}_{3F}$. By contrast,
$\left|\bra{n}\ket{\phi_{U}(\theta_{j},\gamma_{j})}\right|^{2}$ does not have
a neat connection to the population ratios in the original density matrix.
Moreover, $\ket{n}$ is also a state in the optimal decomposition. Due to this
asymmetric treatment of the photon numbers, we call this the upper-pair phase.

One may wonder how the value of $f$ was obtained, given its unintuitive
appearance. First, we deduced from the linear optimization that, for some
values of $(p_{n+2},p_{n+1})$, the optimal decomposition involved the Fock
state $\ket{n}$ and states of the form $\ket{\phi_{U}(\theta_{j},\gamma_{j})}$
(see Eq. (26)). We then optimized the value of $f$ in Eq. (26) analytically,
yielding Eq. (27).

Such a decomposition is not possible everywhere, though. As previously
mentioned , the states $\ket{\phi_{U}(\theta_{j},\gamma_{j})}$ contribute a
total probability $(p_{n+2}+p_{n+1})f^{-1}$. This cannot exceed $1$. Thus:

$p_{n+2}+p_{n+1}\leq\frac{(2+n)p_{n+2}}{(1+n)p_{n+1}+(3+2n)p_{n+2}}.$ (29)

This sets the upper bound on this phase (blue dashed curve in Fig. 2): note
that, at the upper bound,
$\ket{\phi_{U}(\theta_{j},\gamma_{j}}=\ket{\phi_{T}(\theta_{j},\gamma_{j})}$.
To find the highest possible value of $p_{n+2}$, we can set $p_{n+1}=0$ and
saturate the inequality:

$p_{n+2}^{\text{max}}=\frac{2+n}{3+2n}.$ (30)

####  IV.1.3 Lower-pair phase

In Fig. 2, the lower-pair phase is the lower-right, orange regime. It is
analogous to the upper-pair phase, except that in the optimal decomposition it
is the lower-two Fock states $\ket{n+1}$ and $\ket{n}$ that are “paired.” Five
states are a part of the optimal decomposition of $\hat{\rho}_{3F}$. One is
the Fock state $\ket{n+2}$. The other four are states of the form:

$\begin{split}&\ket{\phi_{L}(\theta_{j},\gamma_{j})}=\sqrt{1-g}\ket{n+2}\\\
&+\sqrt{g}\left(\sqrt{\frac{p_{n+1}}{1-p_{n+2}}}e^{i\theta_{j}}\ket{n+1}+\sqrt{\frac{1-p_{n+2}-p_{n+1}}{1-p_{n+2}}}e^{i\gamma_{j}}\ket{n}\right),\end{split}$
(31)

where:

$g=\frac{(1+n)(-1+p_{n+1}+p_{n+2})}{-3-2n+(1+n)p_{n+1}+(3+2n)p_{n+2}}.$ (32)

The states $\ket{\phi_{L}(\theta_{j},\gamma_{j})}$ contribute a total
probability $(1-p_{n+2})g^{-1}$ (each contributes $1/4$ of this). The ORT
measure is:

$\begin{split}\mathcal{N}(\hat{\rho}_{3F})&=2p_{n+2}+p_{n+1}+n\\\
&-(1-p_{n+2})g^{-1}\left|\bra{\phi_{L}(0,0)}\hat{a}\ket{\phi_{L}(0,0)}\right|^{2}.\end{split}$
(33)

We call this the “lower-pair phase” due to the asymmetric treatment of the
photon numbers.
$\left|\bra{n+1}\ket{\phi_{L}(\theta_{j},\gamma_{j})}\right|^{2}/\left|\bra{n}\ket{\phi_{L}(\theta_{j},\gamma_{j})}\right|^{2}=p_{n+1}/(1-p_{n+2}-p_{n+1})$,
the population ratio between $\ket{n+1}$ and $\ket{n}$ in the original state.
By contrast, $\left|\bra{n+2}\ket{\phi_{L}(\theta_{j},\gamma_{j})}\right|^{2}$
does not have a simple connection to the population ratios in the original
density matrix; also, $\ket{n+2}$ is itself a state in the optimal
decomposition.

The value of $g$ was obtained in the same way that $f$ in Eq. (27) was. We
deduced from the linear optimization that, for some values of
$(p_{n+2},p_{n+1})$, the optimal decomposition involved the Fock state
$\ket{n+2}$ and states of the form $\ket{\phi_{L}(\theta_{j},\gamma_{j})}$
(see Eq. (31)). We then optimized $g$ in Eq. (31) analytically.

Similarly, the total probability $(1-p_{n+2})g^{-1}$ conttibuted by the states
$\ket{\phi_{L}(\theta_{j},\gamma_{j})}$ cannot exceed $1$ such that

$1-p_{n+2}\leq\frac{(1+n)(-1+p_{n+1}+p_{n+2})}{-3-2n+(1+n)p_{n+1}+(3+2n)p_{n+2}}.$
(34)

This sets the upper bound on this phase (orange dotted curve in Fig. 2): note
that, at the upper bound,
$\ket{\phi_{L}(\theta_{j},\gamma_{j}}=\ket{\phi_{T}(\theta_{j},\gamma_{j})}$.
To find the lowest possible value of $p_{n+2}$, we can set $p_{n+1}=0$ and
saturate the inequality:

$p_{n+2}^{\text{min}}=\frac{2+n}{3+2n}.$ (35)

Having scanned over various values of $(p_{n+2},p_{n+1})$, and numerically
approximated $\mathcal{N}(\hat{\rho}_{3F})$ for each, we are reasonably
confident (although we lack analytical proof) that these three phases
(triplet, upper-pair, and lower-pair) are the only three, and that we have
characterized the nonclassicality exactly over the entire phase space of
$\hat{\rho}_{3F}$. Assuming the ansatz, the nonclassicality is plotted in Fig.
3

Figure 3: Nonclassicality $\mathcal{N}(\hat{\rho}_{3F})$ over the entire phase
space of $\hat{\rho}_{3F}$. Here, the lowest populated photon number is
assumed to be $n=0$. At the corners, which represent Fock states, the
nonclassicality takes the values $0$, $1$, and $2$. The nonclassicality is
continuous despite there being three distinct phases (see Fig. 2). The points
in each phase are colored differently here for reference.

###  IV.2 Rank-4 case

For fixed $n$ in
$\hat{\rho}_{4F}=p_{n+3}\ket{n+3}\bra{n+3}+p_{n+2}\ket{n+2}\bra{n+2}+p_{n+1}\ket{n+1}\bra{n+1}+(1-p_{n+3}-p_{n+2}-p_{n+1})\ket{n}\bra{n}$,
there are three free parameters: $p_{n+3}$, $p_{n+2}$, and $p_{n+1}$. The
possible choices $(p_{n+3},p_{n+2},p_{n+1})$ of these parameters form the
triangular pyramid defined by $p_{n+3}+p_{n+2}+p_{n+1}\leq 1$, $p_{n+3}\geq
0$, $p_{n+2}\geq 0$, and $p_{n+1}\geq 0$. From our numerical investigations,
we found six distinct regimes, or phases, within this pyramidal phase space.
These fall into three categories: which we call quartet, triplet, and pair.

Figure 4: Nonclassicality $\mathcal{N}(\hat{\rho}_{4F})$ over the phase space
of $\hat{\rho}_{4F}$, assuming the lowest populated photon number of $n=0$. To
our knowledge, there are six distinct regimes (phases), in which there are
distinct ways of optimally decomposing $\hat{\rho}_{4F}$. The complete phase
space composed of these regimes is the triangular pyramid defined by
$p_{n+3}+p_{n+2}+p_{n+1}\leq 1$. This plot was created by calculating the ORT
measure at discrete points within each phase, and some phases may be larger
than they appear. At the vertices of the pyramid (which represent Fock
states), the ORT measure takes the values $0$, $1$, $2$, and $3$.

####  IV.2.1 Quartet Phase

The quartet phase is shown in Fig. 4 (upper-left). This phase comprises the
majority of the phase space, in the same way that the triplet phase comprises
the majority of the phase space in the rank-3 case. This phase has the most
simple optimal decomposition, utilizing states of the form:

$\begin{split}\ket{\phi_{Q}(\theta_{j},\gamma_{j},\zeta_{j})}&=\sqrt{p_{n+3}}\ket{n+3}+\sqrt{p_{n+2}}e^{i\theta_{j}}\ket{n+2}\\\
&+\sqrt{p_{n+1}}e^{i\gamma_{j}}\ket{n+1}+\sqrt{p_{n}}e^{i\zeta_{j}}\ket{n},\end{split}$
(36)

where $p_{n}=1-\sum_{k=1}^{3}p_{n+k}$. These states are a straightforward
generalization of the states in Eq. (24). As before, the decomposition states
come in a symmetric set with equal weights. We call this the quartet phase
because all four photon numbers are treated symmetrically, with amplitudes
equal to the square root of the populations. The ORT measure for states in
this regime is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-\left(\sum_{k=0}^{2}\sqrt{p_{n+k+1}p_{n+k}(n+k+1)}\right)^{2}.\end{split}$
(37)

In this case, the bound in Eq. (23) is saturated.

It is always possible to decompose $\hat{\rho}_{4F}$ into the states in Eq.
(36). However, this decomposition is not optimal for all choices of
$(p_{n+3},p_{n+2},p_{n+1})$. The regions where the quartet decomposition is
not optimal (and, consequently, the boundaries of the quartet phase) are
described in the following.

####  IV.2.2 Triplet Phases

The triplet phases are natural generalizations of the pair phases in the
rank-3 case. There are four triplet phases: the $n$-triplet, $(n+1)$-triplet,
$(n+2)$-triplet, and $(n+3)$-triplet.

$n$-triplet—The $n$-triplet phase is shown in the upper-middle panel of Fig.
4. The optimal decomposition consists of one Fock state, $\ket{n}$, as well as
states of the form:

$\begin{split}&\ket{\phi_{0T}(\theta_{j},\gamma_{j},\zeta_{j})}=\sqrt{1-f_{0}}e^{i\zeta_{j}}\ket{n}+\sqrt{f_{0}}\Bigg{(}\sqrt{\frac{p_{n+3}}{1-p_{n}}}\ket{n+3}\\\
&+\sqrt{\frac{p_{n+2}}{1-p_{n}}}e^{i\theta_{j}}\ket{n+2}+\sqrt{\frac{p_{n+1}}{1-p_{n}}}e^{i\gamma_{j}}\ket{n+1}\Bigg{)}.\end{split}$
(38)

Such states contribute a total probability $(1-p_{n})f_{0}^{-1}$, which cannot
exceed $1$, and so $1\geq f_{0}\geq 1-p_{n}$. Here, $f_{0}$ is chosen to
maximize the objective, i.e., make
$(1-p_{n})f_{0}^{-1}\left|\bra{\phi_{0T}(0,0,0)}\hat{a}\ket{\phi_{0T}(0,0,0)}\right|^{2}$
as large as possible. $f_{0}$ is a quite complicated function of $n$ and the
populations, so we do not show its algebraic form here. The ORT measure for
states in this phase is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-(1-p_{n})f_{0}^{-1}\left|\bra{\phi_{0T}(0,0,0)}\hat{a}\ket{\phi_{0T}(0,0,0)}\right|^{2}.\end{split}$
(39)

$f_{0}=1-p_{n}$ determines a boundary on this phase: note that
$\ket{\phi_{0T}(\theta_{j},\gamma_{j},\zeta_{j})}\rightarrow\ket{\phi_{Q}(\theta_{j},\gamma_{j},\zeta_{j})}$
as $1-p_{n}\rightarrow f_{0}$. Interestingly, this condition is not sufficient
to determine the entire boundary of this phase, because there are points
$(p_{n+3},p_{n+2},p_{n+1})$ satisfying $f_{0}\geq 1-p_{n}$ for which an even
better decomposition exists, hence the $(n+2,n+1)$-pair phase (see lower right
panel of Fig. 4). Ultimately, the $n$-triplet phase consists of the points
satisfying $f_{0}\geq 1-p_{n}$ that are not part of the $(n+2,n+1)$-pair
phase.

Lastly, it is worth noting that the $p_{n+3}=0$ face of the pyramidal phase
space of $\hat{\rho}_{4F}$ corresponds to the rank-$3$ case with populated
photon numbers $n+2$, $n+1$, and $n$; the $n$-triplet phase is a smooth
extension of the upper-pair phase on this rank-$3$ face.

$(n+3)$-triplet—The $(n+3)$-triplet phase is shown in the upper-right panel of
Fig. 4. The optimal decomposition consists of one Fock state, $\ket{n+3}$, as
well as states of the form:

$\begin{split}&\ket{\phi_{3T}(\theta_{j},\gamma_{j},\zeta_{j})}=\sqrt{1-f_{3}}\ket{n+3}\\\
&+\sqrt{f_{3}}\Bigg{(}\sqrt{\frac{p_{n+2}}{1-p_{n+3}}}e^{i\theta_{j}}\ket{n+2}+\sqrt{\frac{p_{n+1}}{1-p_{n+3}}}e^{i\gamma_{j}}\ket{n+1}\\\
&+\sqrt{\frac{p_{n}}{1-p_{n+3}}}e^{i\zeta_{j}}\ket{n}\Bigg{)}.\end{split}$
(40)

Such states contribute a total probability $(1-p_{n+3})f_{3}^{-1}$, which
cannot exceed $1$, and so $1\geq f_{3}\geq 1-p_{n+3}$. Here, $f_{3}$ is chosen
to maximize the objective, i.e., make
$(1-p_{n+3})f_{3}^{-1}\left|\bra{\phi_{3T}(0,0,0)}\hat{a}\ket{\phi_{3T}(0,0,0)}\right|^{2}$
as large as possible. Since $f_{3}$ is a quite complicated function of $n$ and
the populations, we do not show its algebraic form here. The ORT measure for
states in this phase is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-(1-p_{n+3})f_{3}^{-1}\left|\bra{\phi_{3T}(0,0,0)}\hat{a}\ket{\phi_{3T}(0,0,0)}\right|^{2}.\end{split}$
(41)

$f_{3}=1-p_{n+3}$ determines a boundary on this phase: note that
$\ket{\phi_{3T}(\theta_{j},\gamma_{j},\zeta_{j})}\rightarrow\ket{\phi_{Q}(\theta_{j},\gamma_{j},\zeta_{j})}$
as $1-p_{n+3}\rightarrow f_{3}$. This condition, however, is not sufficient to
determine the entire boundary of this phase, because there are points
$(p_{n+3},p_{n+2},p_{n+1})$ satisfying $f_{3}\geq 1-p_{n+3}$ for which an even
better decomposition exists (such points belong to the $(n+2,n+1)$-pair
phase). In fact, there are even points satisfying both $f_{3}\geq 1-p_{n+3}$
and $f_{0}\geq 1-p_{n}$, all of which belong to the $(n+2,n+1)$-pair phase.

The $p_{n+3}+p_{n+2}+p_{n+1}=1$ face of the pyramidal phase space of
$\hat{\rho}_{4F}$ corresponds to the rank-$3$ case with populated photon
numbers $n+3$, $n+2$, and $n+1$; the $(n+3)$-triplet phase is a smooth
extension of the lower-pair phase on this rank-$3$ face.

$(n+1)$-triplet— The $(n+1)$-triplet phase is shown in the lower-left panel of
Fig. 4. The optimal decomposition consists of one Fock state, $\ket{n+1}$, as
well as states of the form:

$\begin{split}&\ket{\phi_{1T}(\theta_{j},\gamma_{j},\zeta_{j})}=\sqrt{1-f_{1}}e^{i\gamma_{j}}\ket{n+1}\\\
&+\sqrt{f_{1}}\Bigg{(}\sqrt{\frac{p_{n+3}}{1-p_{n+1}}}\ket{n+3}+\sqrt{\frac{p_{n+2}}{1-p_{n+1}}}e^{i\theta_{j}}\ket{n+2}\\\
&+\sqrt{\frac{p_{n}}{1-p_{n+1}}}e^{i\zeta_{j}}\ket{n}\Bigg{)}.\end{split}$
(42)

Such states contribute a total probability $(1-p_{n+1})f_{1}^{-1}$, which
cannot exceed $1$, and so $1\geq f_{1}\geq 1-p_{n+1}$. $f_{1}$ maximizes the
objective, by making
$(1-p_{n+1})f_{1}^{-1}\left|\bra{\phi_{1T}(0,0,0)}\hat{a}\ket{\phi_{1T}(0,0,0)}\right|^{2}$
as large as possible. The ORT measure for states in this phase is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-(1-p_{n+1})f_{1}^{-1}\left|\bra{\phi_{1T}(0,0,0)}\hat{a}\ket{\phi_{1T}(0,0,0)}\right|^{2}.\end{split}$
(43)

$f_{1}=1-p_{n+1}$ determines a boundary on this phase: note that
$\ket{\phi_{1T}(\theta_{j},\gamma_{j},\zeta_{j})}\rightarrow\ket{\phi_{Q}(\theta_{j},\gamma_{j},\zeta_{j})}$
as $1-p_{n+1}\rightarrow f_{1}$. As far as we could determine from our
numerical investigations, the $(n+1)$-triplet decomposition is optimal for all
points $(p_{n+3},p_{n+2},p_{n+1})$ in $1\geq f_{1}\geq 1-p_{n+1}$.

The $p_{n+3}+p_{n+2}+p_{n+1}=1$ face of the pyramidal phase space of
$\hat{\rho}_{4F}$ corresponds to the rank-$3$ case with populated photon
numbers $n+3$, $n+2$, and $n+1$; the $(n+1)$-triplet phase is a smooth
extension of the upper-pair phase on this rank-$3$ face.

$(n+2)$-triplet— The $(n+2)$-triplet phase is shown in the lower-middle panel
of Fig. 4. The optimal decomposition consists of one Fock state, $\ket{n+2}$,
as well as states of the form:

$\begin{split}&\ket{\phi_{2T}(\theta_{j},\gamma_{j},\zeta_{j})}=\sqrt{1-f_{2}}e^{i\theta_{j}}\ket{n+2}\\\
&+\sqrt{f_{2}}\Bigg{(}\sqrt{\frac{p_{n+3}}{1-p_{n+2}}}\ket{n+3}+\sqrt{\frac{p_{n+1}}{1-p_{n+2}}}e^{i\gamma_{j}}\ket{n+1}\\\
&+\sqrt{\frac{p_{n}}{1-p_{n+2}}}e^{i\zeta_{j}}\ket{n}\Bigg{)}.\end{split}$
(44)

Such states contribute a total probability $(1-p_{n+2})f_{2}^{-1}$, which
cannot exceed $1$, and so $1\geq f_{2}\geq 1-p_{n+2}$. $f_{2}$ maximizes the
objective, by making
$(1-p_{n+2})f_{2}^{-1}\left|\bra{\phi_{2T}(0,0,0)}\hat{a}\ket{\phi_{2T}(0,0,0)}\right|^{2}$
as large as possible. The ORT measure for states in this phase is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-(1-p_{n+2})f_{2}^{-1}\left|\bra{\phi_{2T}(0,0,0)}\hat{a}\ket{\phi_{2T}(0,0,0)}\right|^{2}.\end{split}$
(45)

$f_{2}=1-p_{n+2}$ determines a boundary on this phase: note that
$\ket{\phi_{2T}(\theta_{j},\gamma_{j},\zeta_{j})}\rightarrow\ket{\phi_{Q}(\theta_{j},\gamma_{j},\zeta_{j})}$
as $1-p_{n+2}\rightarrow f_{2}$. As far as we could determine from our
numerical investigations, the $(n+2)$-triplet decomposition is optimal for all
points $(p_{n+3},p_{n+2},p_{n+1})$ in $1\geq f_{2}\geq 1-p_{n+2}$.

The $p_{n+3}=0$ face of the pyramidal phase space of $\hat{\rho}_{4F}$
corresponds to the rank-$3$ case with populated photon numbers $n+2$, $n+1$,
and $n$; the $(n+2)$-triplet phase is a smooth extension of the lower-pair
phase on this rank-$3$ face.

####  IV.2.3 Pair Phase

The $(n+2,n+1)$-pair phase is shown in the lower right panel of Fig. 4, and to
our knowledge, is the only phase of its kind. This phase is in surface contact
with the $n$-triplet and $(n+3)$-triplet phases. The optimal decomposition
breaks $\hat{\rho}_{4F}$ into two parts: one part has support only for
$\ket{n+3}$ and $\ket{n}$ (any decomposition is optimal for this part), while
the other part is optimally decomposed into states of the form:

$\begin{split}&\ket{\phi_{P}(\theta_{j},\gamma_{j},\zeta_{j})}=\sqrt{1-f}\left(\sqrt{g}\ket{n+3}+\sqrt{1-g}e^{i\zeta_{j}}\ket{n}\right)+\\\
&\sqrt{f}\left(\sqrt{\frac{p_{n+2}}{p_{n+2}+p_{n+1}}}e^{i\theta_{j}}\ket{n+2}+\sqrt{\frac{p_{n+1}}{p_{n+2}+p_{n+1}}}e^{i\gamma_{j}}\ket{n+1}\right),\end{split}$
(46)

where:

$g=\frac{(3+n)p_{n+2}}{(1+n)p_{n+1}+(3+n)p_{n+2}}$ (47)

and $f=f(n,p_{n+2},p_{n+1})$ maximizes
$(p_{n+2}+p_{n+1})f^{-1}|\bra{\phi_{P}(0,0,0)}\hat{a}\ket{\phi_{P}(0,0,0)}|^{2}$.
These states contribute a total probability $(p_{n+2}+p_{n+1})f^{-1}$.
Interestingly, $f$ and $g$ depend only on $n$, $p_{n+2}$, and $p_{n+1}$ (not
on $p_{n+3}$ or $p_{n}$). Consequently,
$(p_{n+2}+p_{n+1})f^{-1}|\bra{\phi_{P}(0,0,0)}\hat{a}\ket{\phi_{P}(0,0,0)}|^{2}$
depends only on $n$, $p_{n+2}$, and $p_{n+1}$. The ORT measure in this phase
is:

$\begin{split}\mathcal{N}(\hat{\rho}_{4F})&=3p_{n+3}+2p_{n+2}+p_{n+1}+n\\\
&-(p_{n+2}+p_{n+1})f^{-1}\left|\bra{\phi_{P}(0,0,0)}\hat{a}\ket{\phi_{P}(0,0,0)}\right|^{2}.\end{split}$
(48)

This regime is bounded by the inequalities:

$\displaystyle p_{n+2}+p_{n+1}$ $\displaystyle\leq f$ (49)
$\displaystyle(1-f)f^{-1}g(p_{n+2}+p_{n+1})$ $\displaystyle\leq p_{n+3}$ (50)
$\displaystyle(1-f)f^{-1}(1-g)(p_{n+2}+p_{n+1})$ $\displaystyle\leq p_{n}.$
(51)

Respectively, these inequalities ensure the
$\ket{\phi_{P}(\theta_{j},\gamma_{j},\zeta_{j})}$ states do not over-fill the
total population, $\ket{n+3}$ population, and $\ket{n}$ population.

Lastly, we note that the $(n+2,n+1)$-pair phase may be regarded as a further
refinement of the $n$-triplet and $(n+3)$-triplet phases. It lies within the
union of $f_{1}\geq 1-p_{n+1}$ and $f_{3}\geq 1-p_{n+3}$ and covers the
intersection of $f_{1}\geq 1-p_{n+1}$ and $f_{3}\geq 1-p_{n+3}$. Within its
domain of validity, the $(n+2,n+1)$-pair decomposition is superior to that of
the corresponding triplet phases.

As far as we were able to determine from our numerical investigations, these
six phases are the only six in the rank-4 case. We do not have an analytical
proof, and of course could only check a finite number of parameter choices
$(p_{n+3},p_{n+2},p_{n+1})$, so these ansatz may be considered conjectures. In
support of our results, however, we note that there is a high degree of
consistency to them, particularly in regards to how they map down to the
rank-3 case.

###  IV.3 Higher-rank cases

To demonstrate the scope of our linear programming numerical method, we apply
it to higher-rank mixed Fock states of the form given by Eq. (17). As an
example, we consider truncated thermal states, where the populations of
$\ket{n}$ $(n\geq M)$ are reduced to zero:

$\hat{\rho}_{MF}=\sum_{k=0}^{M-1}p_{k}\ket{k}\bra{k},$ (52)

where the populations are:

$p_{k}=\mathcal{A}\frac{n_{\text{th}}^{k}}{(1+n_{\text{th}})^{k+1}}$ (53)

with $n_{\text{th}}$ the mean photon number of the field before the truncation
and:

$\mathcal{A}=\frac{1}{1-\left(\frac{n_{\text{th}}}{1+n_{\text{th}}}\right)^{M}}$
(54)

is the normalization factor after truncating the thermal states. Even though
true (i.e., untruncated) thermal states are diagonal in the Fock basis, they
are classical, meaning $\mathcal{N}=0$ and their $P$-function is positive.
This classicality is only possible because their expansion in terms of Fock
states has infinitely many terms, which can be equivalently represented by a
mixture of coherent states [30]. For truncated thermal states, they should
evaluate nonclassicality.

Figure 5: Nonclassicality per unit average energy, $\mathcal{N}/n_{M}$, along
with average energy $n_{M}$, of truncated thermal states. The original thermal
state has average energy $n_{th}=0.5$ and the truncated thermal states are
obtained by projecting out all photon numbers greater than or equal to $M$.
Such truncated thermal states should resemble the original thermal state
(which has nonclassicality $0$) more and more as $M$ increases, and indeed we
see $\mathcal{N}/n_{M}\rightarrow 0$ as $M$ increases.

Increasing the number of terms in a truncated thermal state should cause the
nonclassicality to approach zero (i.e., recovers the thermal state’s
classicality). We demonstrate this in Fig. 5, where we evaluate the ORT
measure, per unit average energy, using our numerical method in the rank-5 and
rank-6 cases. We choose a low $n_{\text{th}}=0.5$ such that the expected
photon number after truncation at $M=6$ is close to $n_{\text{th}}$:
$n_{6}=0.491758$. Nonclassicality per unit energy $\mathcal{N}/n_{M}$ indeed
decreases from its $M=2$ value of $0.25$ to a $M=6$ value of $\lesssim.011$.

###  IV.4 Optimal Decompositions and Convexity

Before concluding, we would like to remark on a mathematical observation about
our results. The insights here not only support our claims (about the optimal
decompositions we found), but may be helpful for solving or checking the
solution of other convex roof optimization problems.

The optimal decompositions we discovered may be divided into two categories:
“simple” and “composite.” Examples of “simple” decompositions are the triplet
decomposition in the rank-3 case and the quartet decomposition in the rank-4
case. Recall that these decompositions are always valid regardless of the
specific populations, although they may not be optimal (in the ORT measure
sense). A state is “simply-decomposed” if its optimal decomposition is a
simple decomposition: for example,
$0.4\ket{2}\bra{2}+0.4\ket{1}\bra{1}+0.2\ket{0}\bra{0}$ is simply-decomposed
since it lies in the triplet phase (see Fig. 2). Equivalently, a state is
simply-decomposed if Eq. (23) is saturated. By contrast,
$\hat{\rho}_{\text{comp}}=0.2\ket{2}\bra{2}+0.2\ket{1}\bra{1}+0.6\ket{0}\bra{0}$
is not simply-decomposed, since it lies in the upper-pair phase (Fig. 2).
Nevertheless, $\hat{\rho}_{\text{comp}}$ may be rewritten as a convex
combination of density matrices that are simply-decomposed (it is
“compositely-decomposed”):
$\hat{\rho}_{\text{comp.}}=p_{1}\hat{\rho}_{\text{simp},1}+p_{2}\hat{\rho}_{\text{simp},2}$.
In this case $\hat{\rho}_{\text{simp},1}=\ket{0}\bra{0}$,
$\hat{\rho}_{\text{simp},2}=0.25\ket{2}\bra{2}+0.25\ket{1}\bra{1}+0.5\ket{0}\bra{0}$,
$p_{1}=0.2$, and $p_{2}=0.8$ — interestingly, $\hat{\rho}_{\text{simp},2}$ is
precisely at the boundary of the triplet and upper-pair phases (Fig. 2). The
upper-pair decomposition of $\hat{\rho}_{\text{comp.}}$ makes precise use of
this convex combination, see Eq. (26) and (27). Recall that the ORT measure
satisfies convexity, meaning
$\sum_{j}p_{j}\mathcal{N}(\hat{\rho}_{j})\geq\mathcal{N}\left(\sum_{j}p_{j}\hat{\rho}_{j}\right)$
for any quantum states $\hat{\rho}_{j}$ and probabilities $p_{j}$. The
equality can hold because the second term in Eq. (3) vanishes for the mixed
Fock states. In this case, the equality holds if and only if the decomposition
into $\hat{\rho}_{j}$ and probabilities $p_{j}$ is optimal. A simple and
important test for our proposed decompositions is that they saturate this
inequality. In section III, we implied (however indirectly) that
$\mathcal{\hat{\rho}_{\text{comp.}}}=p_{1}\mathcal{N}(\hat{\rho}_{\text{simp},1})+p_{2}\mathcal{N}(\hat{\rho}_{\text{simp},2})$,
and so our analytical results pass this test.

These observations hold quite generally across all of our analytical results
in the rank-3 and rank-4 cases. Rank-3 states in the triplet phase are simply-
decomposed, and rank-4 states in the quartet phase are simply-decomposed. All
other rank-3 and rank-4 states are compositely-decomposed, and may be thought
of as convex combinations of simply-decomposed states (particularly states
which lie at the boundaries of the simple phases). We fully expect these
patterns to hold in higher-rank cases as well, although we did not derive any
analytical results for them.

##  V Conclusion

In this work, we introduced a method, based on linear programming, for
evaluating the ORT measure of Fock-diagonal mixed states. In doing so, we have
further advanced the ORT measure as a useful tool for characterizing bosonic
nonclassicality. Our work also suggests that linear programming may be useful
for solving convex roof optimization problems in other contexts, such as
entanglement quantification.

We have provided in-depth results for the case where three or four neighboring
Fock states are populated. Full solutions were obtained by carefully analyzing
the optimal histograms given by the linear programming method. We found that,
depending on the populations, Fock-diagonal states can be in various different
phases, wherein different types of decompositions are optimal (for ORT measure
purposes). We also noted that the optimal decompositions we found fall into
two categories: “simple” and “composite,” and that this distinction may be
useful for solving other convex roof optimization problems.

We gave some numerical results for rank-$5$ and rank-$6$ states, showing that
the method remains applicable. As the rank increases, the computational
complexity does increase, as more bins are needed for the optimal histogram to
be accurate. However, we suspect that such a challenge can be met by employing
an iterative, adaptive approach.

Since the ORT measure is resource-theoretic, and related to metrological
power, our results help facilitate deeper understanding of the resources
behind quantum-enhanced sensing. Moving forward, we would like to use our
results to study the tradeoff between entanglement generation and individual
mode nonclassicality in a beam-splitter [12, 35, 36, 37]. We would also like
to calculate the nonclassicality of mixed states that have coherences in the
Fock basis.

## Acknowledgements

We thank Luke Ellert-Beck and Garrett Jepson for helpful discussions.

This research is supported by NSF Award 2243591.

## References

  * Hillery [1985] M. Hillery, Classical pure states are coherent states, Physics Letters A 111, 409 (1985).
  * Braunstein and Kimble [1998] S. L. Braunstein and H. J. Kimble, Teleportation of continuous quantum variables, Phys. Rev. Lett. 80, 869 (1998).
  * Braunstein and Van Loock [2005] S. L. Braunstein and P. Van Loock, Quantum information with continuous variables, Rev. Mod. Phys. 77, 513 (2005).
  * Mirrahimi  _et al._ [2014] M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, Dynamically protected cat-qubits: a new paradigm for universal quantum computation, New J. Phys. 16, 045014 (2014).
  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum metrology, Phys. Rev. Lett. 96, 010401 (2006).
  * Aasi  _et al._ [2013] J. Aasi, J. Abadie, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Enhanced sensitivity of the ligo gravitational wave detector by using squeezed states of light, Nat. Photonics 7, 613 (2013).
  * Tse  _et al._ [2019] M. Tse, H. Yu, N. Kijbunchoo, A. Fernandez-Galiana, P. Dupej, L. Barsotti, C. Blair, D. Brown, S. Dwyer, A. Effler, _et al._ , Quantum-enhanced advanced ligo detectors in the era of gravitational-wave astronomy, Phys. Rev. Lett. 123, 231107 (2019).
  * Lee [1991a] C. T. Lee, Measure of the nonclassicality of nonclassical states, Phys. Rev. A 44, R2775 (1991a).
  * Marian  _et al._ [2002] P. Marian, T. A. Marian, and H. Scutaru, Quantifying nonclassicality of one-mode gaussian states of the radiation field, Physical Rev. Lett. 88, 153601 (2002).
  * Hillery [1987] M. Hillery, Nonclassical distance in quantum optics, Phys. Rev. A 35, 725 (1987).
  * Gehrke  _et al._ [2012] C. Gehrke, J. Sperling, and W. Vogel, Quantification of nonclassicality, Phys. Rev. A 86, 052118 (2012).
  * Vogel and Sperling [2014] W. Vogel and J. Sperling, Unified quantification of nonclassicality and entanglement, Phys. Rev. A 89, 052302 (2014).
  * Tan  _et al._ [2017] K. C. Tan, T. Volkoff, H. Kwon, and H. Jeong, Quantifying the coherence between coherent states, Physical Rev. Lett. 119, 190405 (2017).
  * Yadin  _et al._ [2018] B. Yadin, F. C. Binder, J. Thompson, V. Narasimhachar, M. Gu, and M. S. Kim, Operational resource theory of continuous-variable nonclassicality, Phys. Rev. X 8, 041038 (2018).
  * Kwon  _et al._ [2019] H. Kwon, K. C. Tan, T. Volkoff, and H. Jeong, Nonclassicality as a quantifiable resource for quantum metrology, Phys. Rev. Lett. 122, 040503 (2019).
  * Ge  _et al._ [2020] W. Ge, K. Jacobs, S. Asiri, M. Foss-Feig, and M. S. Zubairy, Operational resource theory of nonclassicality via quantum metrology, Phys. Rev. Research 2, 023400 (2020).
  * Ge and Zubairy [2020] W. Ge and M. S. Zubairy, Evaluating single-mode nonclassicality, Phys. Rev. A 102, 043703 (2020).
  * Ge  _et al._ [2023] W. Ge, K. Jacobs, and M. S. Zubairy, The power of microscopic nonclassical states to amplify the precision of macroscopic optical metrology, npj Quantum Information 9, 5 (2023).
  * Chitambar and Gour [2019] E. Chitambar and G. Gour, Quantum resource theories, Rev. Mod. Phys. 91, 025001 (2019).
  * Merlin  _et al._ [2022] J. Merlin, E. Devibala, and A. B. M. Ahmed, Operational resource measure of nonclassicality for number states filtered coherent states, Eur. Phys. J. D 76, 47 (2022).
  * Lang and Caves [2013] M. D. Lang and C. M. Caves, Optimal quantum-enhanced interferometry using a laser power source, Phys. Rev. Lett. 111, 173601 (2013).
  * Ferrini  _et al._ [2010] G. Ferrini, D. Spehner, A. Minguzzi, and F. W. J. Hekking, Noise in bose josephson junctions: Decoherence and phase relaxation, Phys. Rev. A 82, 033621 (2010).
  * Arqand  _et al._ [2020] A. Arqand, L. Memarzadeh, and S. Mancini, Quantum capacity of a bosonic dephasing channel, Phys. Rev. A 102, 042413 (2020).
  * Strzałka  _et al._ [2024] M. Strzałka, R. Filip, and K. Roszak, Qubit-environment entanglement in time-dependent pure dephasing, Phys. Rev. A 109, 032412 (2024).
  * Tóth  _et al._ [2015] G. Tóth, T. Moroder, and O. Gühne, Evaluating convex roof entanglement measures, Phys. Rev. Lett. 114, 160501 (2015).
  * Chen  _et al._ [2005] K. Chen, S. Albeverio, and S.-M. Fei, Entanglement of formation of bipartite quantum states, Phys. Rev. Lett. 95, 210501 (2005).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * Glauber [1963] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev. 130, 2529 (1963).
  * Sudarshan [1963] E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, New York, 1997).
  * Lee [1991b] C. T. Lee, Measure of the nonclassicality of nonclassical states, Phys. Rev. A 44, R2775 (1991b).
  * Tóth and Apellaniz [2014] G. Tóth and I. Apellaniz, Quantum metrology from a quantum information science perspective, J. Phys. A: Math. Theor. 47, 424006 (2014).
  * Audenaert  _et al._ [2001] K. Audenaert, F. Verstraete, and B. De Moor, Variational characterizations of separability and entanglement of formation, Phys. Rev. A 64, 052304 (2001).
  * Boyd and Vandenberghe [2004] S. P. Boyd and L. Vandenberghe,  _Convex Optimization_ (Cambridge University Press, 2004).
  * Ge  _et al._ [2015] W. Ge, M. E. Tasgin, and M. S. Zubairy, Conservation relation of nonclassicality and entanglement for gaussian states in a beam splitter, Phys. Rev. A 92, 052328 (2015).
  * Arkhipov  _et al._ [2016] I. I. Arkhipov, J. Peřina Jr., J. Svozilík, and A. Miranowicz, Nonclassicality invariant of general two-mode gaussian states, Scientific Reports 6, 26523 EP (2016).
  * Liu  _et al._ [2023] J. Liu, W. Ge, and M. S. Zubairy, Classical-nonclassical polarity of gaussian states, arXiv preprint arXiv:2310.12104 (2023).