Multiparameter transmission estimation at the quantum Cramér-Rao limit on a cloud quantum computer

Consider a single bosonic mode with annihilation operator $\hat{a}$. After impinging on a beam splitter with transmission amplitude $\eta$ along with another “vacuum” mode annihilated by $\hat{v}$, this operator transforms as

If the vacuum mode is initially unpopulated, this leads to the same reduction in intensity for all input states

As such, the input-output relation in Eq. (1) can be used to model a plethora of physical situations that lead to intensity reductions, ranging from loss through a noisy channel to incomplete reflection, depending on what one considers as the output mode. Notably, sequential applications of Eq. (1) lead to another input-output relation of the same form, so we can always model the total loss or transmittance of a physical scenario by considering the parameter $\eta$ alone, encompassing everything from loss in a network to imperfect detectors.

Given this transformation, how can one best estimate $\eta$? Quantum estimation theory dictates that the variance in any estimate of $\eta$ will be lower-bounded by the inverse of the quantum Fisher information (QFI) $\mathsf{F}_{\hat{\rho}}(\eta)$, as encoded in the quantum Cramér-Rao bound (qCRB):

We have here explicitly included the dependence on the probe state $\hat{\rho}$, for the bound can be improved or worsened by considering different probes, and made no reference to the detection procedure because the QFI automatically optimizes over all such protocols. For a pure probe $\hat{\rho}=\ket{\psi}\bra{\psi}$, the QFI takes the form

where $M$ is the number of times the measurement is repeated. This provides a unique figure of merit for choosing an optimal probe state, so we can consider $M=1$ when optimizing such a state.

Classical light with an average number of photons $E$ (i.e., average energy $E$ in units of energy per photon) takes the form of a coherent state with $|\alpha|^{2}=E$:

This state undergoes the simple transformation $|\alpha\rangle\to|\eta\alpha\rangle$ and thus has the QFI

This can be improved upon by using Fock states with definite photon number $E=N$:

Such states lose their purity when photons are lost, transforming as

yet are proven to be optimal for sensing such loss, with QFI

These present a particularly dramatic improvement in sensing with small amounts of loss $\eta\approx 1$.

How can these quantum limits be attained in practice? One method is to herald the presence of Fock states from easier-to-create TMSV states defined in two orthogonal bosonic modes annihilated by $\hat{a}$ and $\hat{b}$:

with average energy $\sinh^{2}r$ in each mode. Then, if one uses a PNRD in mode $b$ and measures some number $N$ photons to be present, this heralds the presence of the Fock state $\ket{N}$ in mode $a$ that can then be used for estimation tasks. Moreover, the optimal measurement scheme to estimate $\eta$ is to use PNRDs in mode $a$, as the probability distribution $\left|\braket{N}{\psi}\right|^{2}$ is the most sensitive possible for estimating $\eta$. The optimal scheme may thus consist of creating a TMSV, sending one mode directly to a PNRD, and having the second mode interact with the absorptive medium before also being sent to a PNRD (Fig. 1a)). Such a setup relies on there being no loss in the mode sent directly to a PNRD, because otherwise the photon numbers in each mode are not always equal, which is the restriction we circumvent in this work.

We note that different parametrizations can be used for estimating loss; some examples include the transmission probability $q=\eta^{2}$ and the parameter $\phi=\arccos\eta$ defined by the master equation $\partial_{\phi}\rho=\tan\phi\left(2\hat{a}^{\dagger}\vphantom{a}\hat{\rho}\hat{a}-\hat{a}^{\dagger}\vphantom{a}\hat{a}\hat{\rho}-\hat{\rho}\hat{a}^{\dagger}\vphantom{a}\hat{a}\right)$, which might arise more naturally in particular physical problems. The QFI for a new parametrization is simply a transformed version of the old, which we write in the suggestive form

so our results can easily be adapted to fit one’s preferred quantification of loss.

Two modes experience loss in parallel versions of Eq. (1):

If these two modes correspond to two different polarization modes, for example, then this process defines a diattenuation; i.e., dichroism [83].

In this multiparameter context, which can readily be expanded to any integer $k$ modes with transmission amplitudes $\eta_{k}$, one is interested in minimizing the covariances between the estimates of all $k$ parameters $\pmb{\eta}$. The covariance matrix is also governed by the qCRB, now bounded by the QFIM:

where the latter has components

and matrix inequalities $\mathbf{A}\succeq\mathbf{B}$ imply that $\mathbf{A}-\mathbf{B}$ is positive semidefinite. One can change parametrizations by generalizing Eq. (11) to include Jacobians with components $J_{ij}=\partial\eta_{i}/\partial\vartheta_{j}$

where ^⊤ denotes the matrix transpose. An important difference from single-parameter estimation is that the qCRB in Eq. (13) is not always guaranteed to be saturable.

Using classical light $\ket{\alpha}_{a}\otimes\ket{\beta}_{b}$ to probe these two parameters leads to a diagonal QFIM (again we take $M=1$)

To choose the optimal probe state, one must make a choice as to how to weigh the relative importance of minimizing $\mathop{\mathrm{Var}}\nolimits(\eta_{1})$ versus $\mathop{\mathrm{Var}}\nolimits(\eta_{2})$. This can often be chosen in a natural way [72] and, indeed, here one simply weighs the fraction of the total energy $|\alpha|^{2}+|\beta|^{2}$ put into each probe-state mode according to the relative weights given to the two aforementioned variances. With no other distinguishing information due to neither mode being prioritized, we will henceforth seek to minimize the total variance $\mathop{\mathrm{Var}}\nolimits(\eta_{1})+\mathop{\mathrm{Var}}\nolimits(\eta_{2})$, which, in this case, is optimized by putting half of the energy into the classical state in each mode.

Since the transformation of Eq. (12) does not mix the two modes, it is clear that the optimal quantum sensing protocol is to repeat the optimal single-mode sensing protocols in parallel: use a probe state $\ket{m}_{a}\otimes\ket{n}_{b}$, impinge it on the sample, then measure the light using a PNRD for each mode. Again, the QFIM is diagonal, taking the form

In this case, the optimal choice for how to apportion the energy between the two modes depends on the actual values of the parameters being estimated. As in the single-parameter case, a larger transmission amplitude leads to a greater quantum enhancement in sensitivity, so this can appropriately be balanced with the amount of energy in each mode. When no a priori information is present, it behooves one to balance the probe states with $m=n$.

It is difficult to produce Fock states with large numbers of photons. Most of the quantum-enhanced transmission estimation experiments have thus relied on TMSV as a source of nonclassical light, especially in the dim regime with $\sinh^{2}r\ll 1$. We next investigate the performance of TMSV for simultaneously estimating two loss parameters. This is especially pertinent when there is a physical relationship between the two modes being probed: for example, if one is interested in estimating the amount of attenuation in two orthogonal polarization modes to measure linear or circular dichroism, TMSV is a natural choice because it directly produces nonclassical (entangled) states within these two modes.

How does TMSV evolve under Eq. (12)? It evolves into the mixed state [84]

where we have defined the (unnormalized) states

and the parameter-dependent function

We accessed the X8 chip on October 14, 2021, which provides four parallel pairs of modes in which we perform our three-parameter sensing protocol, using the RemoteEngine class from Xanadu’s Strawberry Fields library [86]. We created a program to initialize four copies of the TMSV state given in Eq. (10) with $r=1$, which then propagated along the integrated photonic chip and through optical fibers to eight PNRDs that counted the numbers of photons arriving at the detectors at the end of the circuit (see Ref. [22] Fig. 1 for a diagram and picture of the X8 chip and Ref. [22] §Methods for more details about the apparatus, including the microring resonators that generate the squeezed light and the transition-edge sensors that perform the photon-number detection). Since the loss information is encoded in the distribution of photon numbers, we created histograms for the photon-number distributions recorded for each pair of modes after repeating the program $10^{6}$ times, then repeated the entire process $10^{3}$ times. Standard statistical techniques such as maximum likelihood estimation (MLE) could then be performed to estimate the four triads of loss parameters. The source code and recorded data can be found on GitHub at https://github.com/AaronGoldberg9/transmission_TSMV_X8.

The average probability distributions for the first pair of modes for all $10^{9}$ measurements are plotted in Fig. 7 (see also Figs. 16, 17, and 18 in Appendix A for the same plots of the other mode pairs). The average number of photons in each of the modes, labeled from 0 to 7, is found to be 0.442, 0.310, 0.417, 0.292, 0.488, 0.344, 0.429, and 0.323 which should be compared to the expected number $\sinh^{2}r\approx 1.38$ in each mode. These data naively correspond to a transmission of $\approx 20-35\%$ of light in each mode, with loss being introduced anywhere from the generation of the TMSV states to the fiber coupling to the fiber propagation to the detection events. By locating these upper limits on the transmission parameters of $0.45-0.6$ on Fig. 4, not shown, we anticipate that the amount of loss in the present setup will preclude us from reaching a quantum advantage in sensitivity. However, setting $r=1$ using Strawberry Fields only serves to turn on the pump laser that then gets split into one pump for each of the four mode pairs; the actual pump powers lead to much larger, yet unspecified, squeezing parameters $r$, which also imply the presence of higher loss [87].

For each of the $10^{3}$ trials and each pair of modes, we perform MLE to find which values of $\eta_{1}$, $\eta_{2}$, $r$ (where $r$ encompasses both the initial pump power and all losses or inefficiencies prior to TMSV generation, instead of keeping the bulky notation $\eta_{0}r$), $\nu_{1}$, and $\nu_{2}$ in the conditional photon-number-probability distribution [c.f. Eq. (18)]

best reproduce the measured photon-number probability distributions that we will refer to as $q(m,n)$. In the limit of large numbers of experiments, MLE provides an unbiased estimate of the true underlying parameters with the lowest variances possible (saturating the Cramér-Rao lower bound). Our source code for this procedure is also available on GitHub.

What is the likelihood function required for MLE? For a discrete probability distribution, we aim to find the set of estimated parameters $\hat{\pmb{\theta}}\equiv(\hat{\eta}_{1},\hat{\eta}_{2},\hat{r},\hat{\nu}_{1},\hat{\nu}_{2})$⁴⁴4We use the symbol $\pmb{\theta}$ instead of $\pmb{\eta}$ to repeatedly represent the five parameters for a given pair of modes, which constitute a single set of independent parameters being estimated, in order to avoid confusion with the eight loss parameters in the eight different modes across the four parallel experiments. that minimize the Kullback-Leibler divergence $D_{\mathrm{KL}}$ between the measured probability distribution and our parametrized probability distribution $p(m,n|\hat{\pmb{\theta}})$. This can be seen by assuming each of the $\mu=10^{6}$ shots to be statistically independent, with probability distribution depending on the underlying values of the parameters $p(m,n|\pmb{\theta})$, such that the overall likelihood function is given by the product

where $m_{k}$ and $n_{k}$ are the number of photons recorded at the two detectors for the $k$th measurement run. Noticing that we can rearrange this likelihood product by counting the number of times each of the distinct measurement events $(m,n)$ arises, which we can refer to by $\mu_{mn}=q(m,n)\mu$, the likelihood becomes

We can now directly find the parameters that maximize the likelihood function, which also maximize the log-likelihood function. We write

where we repeatedly added constant values to the quantity being optimized and reversed the extrema by inverting the arguments of the logarithms. Numerically minimizing this Kullback-Leibler divergence $D_{\mathrm{KL}}$ over all values of $\eta_{1}$, $\eta_{2}$, $r$, $\nu_{1}$, and $\nu_{2}$ directly gives us MLE’s estimator for the underlying parameters.

The one numerical caveat is that there may be events $(m,n)$ that do not arise, such that the probabilities $p(m,n)$ that contribute to $D_{\mathrm{KL}}$ do not quite sum to unity. To surmount this, we use the standard method from Refs. [88, 89] (there attributed to Fermi under the title “extended maximum likelihood”) of normalizing each $p(m,n)$ by the total of the probabilities for the contributing events with $q(m,n)>0$. In the present case, the difference between the results with this normalization and without it are negligible, because the PNRDs capture the vast majority of the nonnegligible values of $p(m,n)$.

Another method for obtaining precision estimates is to combine all of the $10^{9}$ measurements into a single MLE procedure, ensuring that the likelihood function is gradientless at the estimated parameters, and using the concavity of the likelihood function evaluated at $\pmb{\theta}=\hat{\pmb{\theta}}$ to determine the measured/estimated Fisher information’s components [89] (see Appendix C for a brief derivation):

We perform this gold standard procedure and compare it to the aggregated estimates from performing MLE $10^{3}$ times on the respective trials.

We repeat the estimation procedure for all $M=10^{3}$ trials and analyze the set of estimated transmission parameters. These results for the TMSV in modes 0 and 4 are plotted in Figs. 8 and 9 and the rest of the distributions are displayed in Appendix A in Figs. 19, 20, 21, 22, 23, and 24. The MLE estimates are approximately Gaussian for all of the parameters, qualitatively demonstrating convergence of the estimation protocol. We find that the squeezing parameters are rather large, all approximately $1.3$, leading to overall transmission probabilities of about $8\%$ to $15\%$ in each mode (first data column in Table 1). We also find the dark count rates to range from 0.03 to 0.08 in each mode, which is much higher than one would expect from the detectors alone [90], so we suspect that some of the pump light is contributing to the spurious detection counts. This is confirmed by running the same circuit with the pump turned off, which reduces the dark count rates to $0.01$, 0.01, 0.03, 0.02, 0.02, 0.02, 0.03, and 0.04.

Can we confirm that the estimators have converged to the “true” values of the underlying parameters? From the histograms, one can appreciate that the results are approximately normally distributed. A quick test for the convergence properties is to use the equivariance property of MLE: we perform MLE using a different parametrization, such as using the transmission probabilities $q_{i}=\eta_{i}^{2}$, and inspect whether $\sqrt{\left\langle\hat{q}_{i}\right\rangle}=\left\langle\hat{\eta}_{i}\right\rangle$, where the expectation value is taken over all MLE trials. We compare all of these data in Table 1; near-identical results are visible with all other reparametrizations that we verified, including the aforementioned $\phi=\arccos\eta$ parameter. In Table 1, we also include the results for performing the MLE on the aggregate data from all $10^{9}$ trials, whose covariances are determined by the concavity of the likelihood function at the estimated parameter values through Eq. (31). This significant agreement with the means and standard deviations for our estimators $\hat{\eta}_{i}$ as quoted in the abstract confirms the convergence of our estimators. The most significant outlier from ideal convergence comes from the mode pair 2 and 6 with only a 1% deviation between methods, with all of the rest of the means agreeing to even higher precision.

How do these results compare to the ultimate limit? Taking the estimated parameters to represent the true underlying ones, we can inspect the lowest variances that might be obtained using coherent and Fock states with the same total input energies $E=2\sinh^{2}\hat{r}$. We use Eqs. (16) and (17) with $|\alpha|^{2}=|\beta|^{2}=E/2$ and $m=n=E/2$, respectively to find the $M=1$ results ranging from approximately 0.15 to 0.17. If we now increase to $M=10^{3}$ to represent the different trials performed, we see the variances in our estimates of $\eta_{1}$ and $\eta_{2}$ to differ by 1-2 orders of magnitude with respect to the limits imposed by coherent and Fock states. Remarkably, the variances obtained are smaller than the ultimate limit provided by Fock states with comparable energies (Table 1), implying that more care must be taken in these statistics.

If we instead increase $M$ to $10^{9}$ to account for the total number of experiments performed, the variances of the $10^{3}$ estimates easily obey the quantum Cramér-Rao lower bound. We can then use the techniques of Refs. [88, 89] to analyze all of the data in aggregate, discarding the trial identifications. We then find the estimated parameters quoted in the abstract and variances closer to but always greater than the variances one would obtain from coherent and Fock states with $M=10^{9}$ repetitions (Table 1).

One can finally compare the variances obtained thus with the ultimate limit one could ever expect to obtain using TMSV states that take into account the possible variations in $r$, as in Eq. (23). Then, the measured variances are approximately twice to thrice the ultimate bound: in Table 1, we observe that the estimates approximately saturate the qCRB, with variances coming within a factor of three of the ultimate quantum limit using TMSV as a probe state for simultaneously estimating all three parameters, even in the presence of nonnegligible dark counts that deteriorate the measurement.⁵⁵5The “factor of two” phrase in our abstract for the uncertainties is really $\lesssim\sqrt{3}$. This means that they could certainly beat classical bounds if the transmission parameters were larger and is why we refer to “at the quantum limit” in our title.

Overall, how well do our estimates account for the underlying data? We can measure the root-mean-square (RMS) error between the probability distributions that we measured and the probability distributions that we can construct using our estimators. We find the overall RMS values 0.003 for each of the four joint probability distributions for the four pairs of modes, so the measured histograms are certainly close to what one would expect from our loss model.