Quantifying nonclassicality of mixed Fock 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]:

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) [28] 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}$.¹¹1In 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].

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 [29, 30] as

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 [31, 8].

The ORT measure of nonclassicality for mixed states is [16]:

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]

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.

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 [23, 24, 25], 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:

It is obvious that $\mathcal{N}(\hat{\rho})\leq\langle\hat{a}^{\dagger}\hat{a}\rangle$, which 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$).²²2Actually, 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})$.³³3That 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.

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 an arbitrary number of 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:

Respectively, the above constraints address normalization, population constraints, lack of nonzero 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.  (5) reads:

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 Eq. (7) is called the objective function $\mathcal{L}[q(x,\theta)]$:

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_{k=0}^{1}\sum_{j=0}^{3}\delta(\theta-\theta_{0}-\frac{j\pi}{2}-2\pi k)\right)$, where $\theta_{0}\in[0,\frac{\pi}{2})$ is arbitrary. It is straightforward to check that $\tilde{q}(x,\theta)$ satisfies the same constraints (Eq.  (6a)- (6d)) 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 solved, we may consider instead $\mathcal{Q}(x)\equiv\int_{0}^{2\pi}d\theta\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.  (7) becomes:

The optimization problem (including constraints) resembles a linear programming problem [35] 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):

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.  (10), should converge as $\Delta\rightarrow 0$ (here $\gtrsim$ means approximately equal to, but never less than).

Here $\mathcal{Q}_{x}$ is a discrete probability distribution (“histogram”) over the set of points in $\mathcal{X}$, satisfying constraints which reproduce the original density matrix. As such, $\mathcal{Q}_{x}$ 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 from linear optimization matches the true optimal decomposition perfectly. This will not always be the case, since the chosen binning may skip over the optimal point(s) (generally, the discrete distribution will be slightly suboptimal, acting as a bound). However, one can work backwards from the numerical evidence to make an ansatz for a better decomposition (which may, in fact, be the true optimal one). For example, by realizing that $0.4=\sqrt{0.16}$, in Fig. 1, one may conjecture the general trend that $\mathcal{Q}_{x}$ is spiked at $x=\sqrt{p_{n+1}}$. Such an ansatz can then be corroborated by checking it against finer binnings and other values for the populations. In this way, we obtain analytical ansatzes for $\mathcal{N}(\hat{\rho}_{3F})$ and $\mathcal{N}(\hat{\rho}_{4F})$, in section IV.

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

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

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 mixture of 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})\equiv\sum_{\vec{\theta}}\tilde{q}(\vec{x},\vec{\theta})$ over the vectors $\vec{x}$. This is because, given any $q(\vec{x},\vec{\theta})$, one may always construct (see Appendix B) a $\tilde{q}(\vec{x},\vec{\theta})$ with $\vec{\theta}$-dependence such that:

where $\hat{\rho}_{\vec{x}}=\sum_{k=0}^{M-1}x_{k}^{2}\ket{n+k}\bra{n+k}$ is a normalized Fock-diagonal state indexed by $\vec{x}$, $\sum_{\vec{\theta}}\tilde{q}(\vec{x},\vec{\theta})|\bar{\alpha}_{\vec{x},\vec{\theta}}|^{2}\geq\sum_{\vec{\theta}}q(\vec{x},\vec{\theta})|\bar{\alpha}_{\vec{x},\vec{\theta}}|^{2}$, and $\sum_{\vec{\theta}}\tilde{q}(\vec{x},\vec{\theta})\bar{\alpha}_{\vec{x},\vec{\theta}}^{2}=0$. Constructed in this way, $\tilde{q}(\vec{x},\vec{\theta})$ will be at least as optimal as the originating $q(\vec{x},\vec{\theta})$. Note that the existence of such a $\tilde{q}(\vec{x},\vec{\theta})$ demonstrates that Eq.  (5) may be rewritten:

in general, for Fock-diagonal states.

The optimization problem for $\mathcal{Q}(\vec{x})$ can be solved approximately by discretizing the set of possible $\vec{x}$:

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$, $l_{k}$ are non-negative integers, and $0<\Delta\ll 1$. The approximation, Eq.  (16), 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 for Fock-diagonal states 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}}}$:

Equality holds in the $M=2$ case (where this simple decomposition is optimal), as demonstrated previously. For $M>2$, equality holds in some regimes (i.e., some choices of the populations $p_{n+k}$), although not all.

States that saturate the equality of Eq.  (18) are of some fundamental interest, acting as building blocks to decompose other states. As implied by Eq.  (14) and  (17a)- (17b), for any probability distribution $\mathcal{Q}(\vec{x})$:

We claim that, for the optimal $\tilde{\mathcal{Q}}(\vec{x})$, any $\hat{\rho}_{\vec{x}}$ (which is itself a Fock-diagonal state) for which $\tilde{\mathcal{Q}}(\vec{x})>0$ must saturate the equality of Eq.  (18). To see this, suppose that for some contributing $\vec{x}^{*}$, $\hat{\rho}_{\vec{x}^{*}}$ does not saturate Eq.  (18). This would imply that there exists an optimal decomposition $\tilde{\mathcal{R}}(\vec{y})$ of $\hat{\rho}_{\vec{x}^{*}}$ such that: $\sum_{\vec{y}}\tilde{\mathcal{R}}({\vec{y}})\left(\sum_{k=0}^{M-2}y_{k+1}y_{k}\sqrt{n+k+1}\right)^{2}$ is greater than $\left(\sum_{k=0}^{M-2}x_{k+1}^{*}x_{k}^{*}\sqrt{n+k+1}\right)^{2}.$ One could then construct a better decomposition (than $\tilde{\mathcal{Q}}(\vec{x})$) of $\hat{\rho}_{MF}$ by reallocating the probability assigned to $\hat{\rho}_{\vec{x}^{*}}$: $\tilde{\mathcal{S}}(\vec{x})=\tilde{\mathcal{Q}}(\vec{x})(1-\delta_{\vec{x},\vec{x}^{*}})+\tilde{\mathcal{Q}}(\vec{x}^{*})\tilde{\mathcal{R}}(\vec{x})$. But, by assumption, $\mathcal{Q}(\vec{x})$ was optimal. Thus, there is a contradiction with having $\hat{\rho}_{\vec{x}^{*}}$ not saturate Eq.  (18). Our results reflect this, as discussed in subsection IV.1.4.

We call Fock-diagonal states that saturate the equality of Eq.  (18) simply-decomposed states, and Fock-diagonal states that do not saturate Eq.  (18) compositely-decomposed states. The optimal decomposition of a compositely decomposed state is a convex combination $\hat{\rho}_{MF}=\sum_{\vec{x}}\tilde{\mathcal{Q}}(\vec{x})\hat{\rho}_{\vec{x}}$ of simply-decomposed states. However, each contributing $\hat{\rho}_{\vec{x}}$ could have a rank less than or equal to $M$.

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 identified three distinct regimes, or phases, within this triangular phase space (see Fig. 2). For each phase, we provide an analytical ansatz, based on numerical evidence, for the ORT measure and associated optimal decomposition.