Heisenberg-Scaling Measurement Protocol for Analytic Functions with Quantum Sensor Networks

TL;DR

This paper introduces a quantum sensor network protocol that leverages entanglement to significantly improve the precision of estimating analytic functions of multiple parameters, with proven optimality for qubits and potential applications in quantum computing.

Contribution

It generalizes previous work to multiple parameters, demonstrating an entanglement-based protocol that achieves optimal measurement scaling for qubit sensors and proposing extensions for photonic systems.

Findings

Entanglement improves estimation accuracy by a factor proportional to the number of parameters.

The protocol is proven optimal for qubit sensors.

Potential applications include calibration in quantum computing.

Abstract

We generalize past work on quantum sensor networks to show that, for $d$ input parameters, entanglement can yield a factor $\mathcal O(d)$ improvement in mean squared error when estimating an analytic function of these parameters. We show that the protocol is optimal for qubit sensors, and conjecture an optimal protocol for photons passing through interferometers. Our protocol is also applicable to continuous variable measurements, such as one quadrature of a field operator. We outline a few potential applications, including calibration of laser operations in trapped ion quantum computing.

Figures

Equations

$$\hat{H}={\hat{H}}_c(t)+\frac{1}{2}{\theta }_i{\hat{\sigma }}_i^z,$$

$$M=\mathbb{E}[(\tilde{f}-f(\mathbf{\theta }){)}^2]=\mathrm{Var}\tilde{f}+(\mathbb{E}[\tilde{f}]-f(\mathbf{\theta }){)}^2.$$

$$\mathbb{E}[(\tilde{f}-f(\mathbf{\theta }){)}^2]\ge \frac{1}{F}\ge \frac{1}{F_Q},$$

$$\hat{g}=\frac{\partial \hat{H}}{\partial f}=\frac{\partial \hat{H}}{\partial {\theta }_i}\frac{\partial {\theta }_i}{\partial f}=\frac{1}{2}{\hat{\sigma }}_z^i\frac{\partial {\theta }_i}{\partial f}.$$

$$\Vert \hat{g}{\Vert }_s=\sum _{i=1}^N\left|\frac{\partial {\theta }_i}{\partial f}\right|=\sum _{i=1}^N\left|J_{i1}^{-1}\right|.$$

$$1=J_{1i}{J}_{i1}^{-1}\le \left|J_{1i}\right|\left|J_{i1}^{-1}\right|\le \underset{j}{\max }\left|J_{1j}\right|\sum _{i=1}^N\left|J_{i1}^{-1}\right|.$$

$$\Vert \hat{g}{\Vert }_s=\sum _{i=1}^N\left|J_{i1}^{-1}\right|\ge {\left(\underset{j}{\max }\left|\frac{\partial f}{\partial {\theta }_j}\right|\right)}^{-1}.$$

$$M=\mathbb{E}[(\tilde{f}-f(\mathbf{\theta }){)}^2]\ge \frac{1}{F_Q}\ge \underset{j}{\max }\frac{{\left|\frac{\partial f}{\partial {\theta }_j}\right|}^2}{t^2}.$$

$$M_{\mathrm{unentangled}}=f_i(\mathbf{\theta }{)}^2\mathrm{Var}{\tilde{\theta }}_i.$$

$$M_{\mathrm{unentangled}}=\frac{{\Vert \nabla f(\mathbf{\theta })\Vert }^2}{t_{\text{total}}^2},$$

$$M$$

$$=\mathbb{E}[{\mathrm{Var}}_{\tilde{q}}\tilde{q}]+\frac{2f_{ij}(\mathbf{\theta })+f_{ii}(\mathbf{\theta })f_{jj}(\mathbf{\theta })}{4}\mathrm{Var}{\tilde{\theta }}_i\mathrm{Var}{\tilde{\theta }}_j.$$

$$\mathrm{Var}\widetilde{\mathbf{\alpha }\cdot \mathbf{\theta }}\ge \frac{{\max }_i{\alpha }_i^2}{t^2},$$

$$M$$

$$=\frac{\mathbb{E}[{\max }_i{f}_i(\tilde{\mathbf{\theta }}{)}^2]}{t_2^2}+\frac{g_1(\mathbf{\theta })}{t_1^4},$$

$$f_1(\mathbf{\theta }{)}^2+\frac{f_1(\mathbf{\theta })f_{1ii}(\mathbf{\theta })}{t_1^2}+\frac{f_{1i}(\mathbf{\theta }{)}^2}{t_1^2}+\mathcal{O}((\tilde{\mathbf{\theta }}-\mathbf{\theta }{)}^3).$$

$$M=\frac{g_2(\mathbf{\theta })}{t_2^2}+\frac{g_3(\mathbf{\theta })}{t_1^2{t}_2^2}+\frac{g_1(\mathbf{\theta })}{t_1^4}+\mathcal{O}((\tilde{\mathbf{\theta }}-\mathbf{\theta }{)}^3),$$

$$M=\frac{g_2(\mathbf{\theta })}{t_2^2}+\frac{g_3(\mathbf{\theta })}{t_1^2{t}_2^2}+\frac{g_1(\mathbf{\theta })}{t_1^4}.$$

$$\frac{2g_2(\mathbf{\theta })}{t_2^3}+\frac{2g_3(\mathbf{\theta })}{t_2^3{t}_1^2}=\frac{2g_3(\mathbf{\theta })}{t_2^2{t}_1^3}+\frac{4g_1(\mathbf{\theta })}{t_1^5}.$$

$$g_2(\mathbf{\theta })t_1^2=\frac{g_3(\mathbf{\theta })}{r}+\frac{2g_1(\mathbf{\theta })}{r^3}-g_3(\mathbf{\theta }).$$

$$t_1\approx {\left(\frac{2g_1(\mathbf{\theta })}{g_2(\mathbf{\theta })}\right)}^{1/5}{t}_2^{3/5}\approx {\left(\frac{2g_1(\mathbf{\theta })}{g_2(\mathbf{\theta })}\right)}^{1/5}{t}_{\text{total}}^{3/5}.$$

$$t_1=g(\mathbf{\theta })t_{\text{total}}^{3/5},$$

$$\underset{t_{\text{total}}\to \infty }{\lim }M=\underset{t_{\text{total}}\to \infty }{\lim }\frac{g_2(\mathbf{\theta })}{t_{\text{total}}^2}+\frac{g_3(\mathbf{\theta })}{c^2{t}_{\text{total}}^{2+2p}}+\frac{g_1(\mathbf{\theta })}{c^4{t}_{\text{total}}^{4p}}.$$

$$M=\frac{{\max }_i{f}_i(\mathbf{\theta }{)}^2}{t_{\text{total}}^2},$$

$$\mathrm{Var}\widetilde{\mathbf{\alpha }\cdot \mathbf{\theta }}\ge \frac{{\Vert \mathbf{\alpha }\Vert }_1^2}{N^2}.$$

$$M$$

$$M=\frac{{\Vert \nabla f(\tilde{\mathbf{\theta }})\Vert }_1^2}{N_{\text{total}}^2}+\mathcal{O}\left(\frac{1}{N_{\text{total}}^{12/5}}\right).$$

$$M_{\mathrm{unentangled}}=\frac{f_i(\mathbf{\theta }{)}^2}{n_i^2}.$$

$$M_{\mathrm{unentangled}}=\frac{{\Vert \nabla f(\mathbf{\theta })\Vert }_{2/3}^2}{N_{\text{total}}^2}.$$

$$\displaystyle|\psi_{\mathrm{spin-photon}}\rangle=\frac{1}{\sqrt{2}}\big{(}$$

$$+$$

$$f(\tilde{\mathbf{\theta }})=f(\mathbf{\theta })+T_1+\frac{T_2}{2}+\frac{T_3}{6}+\dots .$$

Full text

I Introduction

Entanglement is a valuable resource for quantum technology. In metrology, entangled probes are capable of more accurate measurements than unentangled probes Bollinger et al. (1996); Huelga et al. (1997); Paris (2008); Pezze and Smerzi (2009); Tóth (2012); Zhang and Duan (2014). In addition to using entangled probes to enhance the measurement of a single parameter, using entanglement to estimate many parameters at once, or a function of those parameters, has recently been an area of interest due to potential applications in tasks such as nanoscale nuclear magnetic resonance imaging Genoni et al. (2013); Humphreys et al. (2013); Gao and Lee (2014); Vidrighin et al. (2014); Yue et al. (2014); Zhang and Fan (2014); Kok et al. (2017); Proctor et al. (2018); Eldredge et al. (2018).

In this work, we are interested in generalizing the work of Ref. Eldredge et al. (2018), which demonstrated a lower bound on the variance of an estimator of a linear combination of $d$ parameters coupled to $d$ qubits. We will generalize this approach to measuring an arbitrary real-valued, analytic function of $d$ parameters and show that entanglement can reduce the variance of such an estimate by a factor of $\mathcal{O}(d)$. Finally, we present a protocol which achieves optimal variance asymptotically in the limit of long measurement time. In addition, when the parameters are coupled to $d$ interferometers or to a combination of interferometers and qubits, we propose an analogous Heisenberg-scaling protocol to improve measurement noise. However, in this case, we lack a proof of optimality. We also can use the protocol presented in Ref. Zhuang et al. (2018) to couple the parameters to continuous variables detected by homodyne measurements.

We will also examine the application of such a protocol to field interpolation. Suppose sensors are placed at $d$ spatially separated locations, but we wish to know the field at a point with no sensor. We may pick a reasonable ansatz for the field with no more than $d$ parameters, use our $d$ measurements to fix the degrees of freedom of that ansatz, and compute the field at our desired point. Because the field of interest is a function of the field at $d$ other locations, our protocol offers reduced noise over performing the same procedure without using entanglement.

II Setup

In this work, bold is used to indicate vectors, hats (as in $\hat{H}$) indicate operators, and variables with a tilde (such as $\tilde{f}$) are estimators of the corresponding quantity with no tilde (such as $f$). The notation $\mathbb{E}_{Y}[X]$ means the expected value of $X$ over all possible $Y$. If we merely write $\mathbb{E}[X]$, then we average over all parameters required to define $X$ (e.g. if $Y$ depended on $Z$, then $\mathbb{E}_{Z}[\mathbb{E}_{Y}[X]]$). We define the variance, $\operatorname{Var}_{Y}[X]$, similarly.

We consider a system with $d$ sensor nodes, where node $i$ consists of a single qubit coupled to a real parameter $\theta_{i}$ (see Fig. 1), and suppose that the state evolves under the Hamiltonian

[TABLE]

where $\hat{\sigma}_{i}^{x,y,z}$ are the Pauli operators acting on qubit $i$ and $\hat{H}_{c}(t)$ is a time-dependent control Hamiltonian that we choose, which may include coupling to ancilla qubits. Here, and throughout the paper, repeated indices indicate summation. We want to measure an arbitrary real-valued, analytic function $f(\bm{\theta})$ of $d$ unknown parameters $\bm{\theta}=\langle\theta_{1},\dots,\theta_{d}\rangle$ for time $t_{\text{total}}$. We would like to determine how well the quantity $f(\bm{\theta})$ can be estimated, and find a protocol for doing so. To specify a protocol, we choose an input state, a control Hamiltonian $\hat{H}_{c}(t)$, and a final measurement.

For a general estimator, we use the mean squared error (MSE) $M$ of our estimate $\tilde{f}$ from the true value $f(\bm{\theta})$ as a figure of merit. Explicitly,

[TABLE]

Thus the MSE accounts for both the variance and the bias of the estimator $\tilde{f}$. By proving lower-bounds for $M$ and then showing that these bounds are saturable, we will be demonstrating protocols which are optimal in this combination of bias and variance.

III Lower bound on error

We now identify the minimum possible error of an estimator of $f(\bm{\theta})$ which measures for time $t_{\text{total}}$. For any estimator $\tilde{f}$, biased or otherwise, which uses samples from a probabilistic process (such as physical experiments) to estimate the value $f(\bm{\theta})$, the MSE is bounded by Braunstein and Caves (1994)

[TABLE]

where $F$ is the Fisher information for the parameter $f$ and $F_{Q}$ is the quantum Fisher information evaluated over our input state, with $F_{Q}\geq F$ always Baumgratz and Datta (2015). Bounds on the error of an estimator in terms of the Fisher information are known as Cramér-Rao bounds. The Fisher information measures the sensitivity of the sampled probability distribution to changes in the parameters $\bm{\theta}$. While $F$ tells us something about a particular experimental setup, $F_{Q}$ is maximized over all possible experiments that could be performed on a state.

In order to evaluate the Fisher information for our function of interest $f$, we will use the method presented in Ref. Boixo et al. (2007) and developed for linear functions in Ref. Eldredge et al. (2018). We start by evaluating the generator $\hat{g}=\partial\hat{H}/\partial f$ as defined in Ref. Boixo et al. (2007). By first writing the chain rule, we find that

[TABLE]

Note that $F_{Q}$ can be upper-bounded by the seminorm of this generator, $F_{Q}\leq t^{2}\|\hat{g}\|_{s}^{2}$ Boixo et al. (2007). (The seminorm of an operator is the difference between its maximum and minimum eigenvalues.) However, to evaluate the seminorm, we will need to evaluate the partial derivative in Eq. (4). To do so we must specify a full basis of functions so that the partial derivative can be defined, which requires specifying which variables are held constant during differentiation. We suppose that a set of functions $f_{1},f_{2},f_{3}\dots f_{N}$ are created, with the $f$ of interest equal to $f_{1}$, defining an invertible coordinate transformation on a region $\mathbb{R}^{N}$ around the point $\mathbf{\theta}$. The seminorm is then:

[TABLE]

Here, $J_{ij}^{-1}$ is an element of the Jacobian matrix of the inverse transformation to that defined by the $f$ functions. Depending on which functions are chosen, the value of $\|\hat{g}\|_{s}$ can vary, as can be seen in Ref. Eldredge et al. (2018) for linear functions. We therefore wish to find the smallest possible $\|\hat{g}\|_{s}$, which will provide the tightest possible bound on $F_{Q}$. To do so, we note that $J^{-1}$ and $J$ must obey an inverse relationship, meaning that the following chain of inequalities holds,

[TABLE]

By using the definition of the Jacobian, we can rewrite this as a lower bound on the value of $\|g\|_{s}$ in terms of partial derivatives of $f$:

[TABLE]

All that remains is to note that if we label the $\theta_{i}$ that yields the maximum first derivative as $\theta_{1}$, and then choose $f_{i}=\theta_{i}$ for $i>1$, the lower bound in Eq. (7) is met, since $\partial\theta_{i}/\partial f_{1}$ must be evaluated holding the other $f_{j}$ constant. Invoking the resulting bound on the quantum Fisher information, we find that the quantum Cramér-Rao bound becomes

[TABLE]

Although the quantum Cramér-Rao bound derived in Eq. (8) cannot always be saturated, it can when the generators $\partial\hat{H}/\partial\theta_{i}$ commute, as in Eq. (1) Baumgratz and Datta (2015). We will show later that the inequality in Eq. (8) can be saturated at asymptotic time $t_{\text{total}}$.

From this point forward, to simplify later calculation, we define $f_{i}(\bm{\theta})=\frac{\partial f(\bm{\theta})}{\partial\theta_{i}}$. This definition also generalizes to multiple partial derivatives (i.e. $f_{ij}=\frac{\partial}{\partial\theta_{j}}\frac{\partial f}{\partial\theta_{i}}$).

Before moving on to the optimal protocol, we will consider a protocol which does not use entanglement and does not saturate Eq. (8) as a useful contrast to an entangled strategy. Suppose we estimate each parameter individually, without bias. Then the MSE $\mathbb{E}[(f(\tilde{\bm{\theta}})-f(\bm{\theta}))^{2}]$ can be written as

[TABLE]

Here we assume the measurement of each single parameter can be made in time $t$ with variance $\operatorname{Var}\tilde{\theta}_{i}=\frac{1}{t^{2}}$, the Heisenberg limit for single particles and therefore the best possible measurement for a non-entangled protocol Bollinger et al. (1996). Estimation protocols that allow one to reach a variance proportional to $1/t^{2}$ without entanglement are outlined in detail in Ref. Kimmel et al. (2015); an experimental realization of single phase estimation without entanglement was performed in Ref. Rudinger et al. (2017). While in realistic settings a Heisenberg-limited measurement on one particle may be challenging and include some constant overhead above $1/t^{2}$, this assumption allows us to identify the improvement possible by using entanglement. Our entanglement-free figure of merit is

[TABLE]

where the $\|\cdot\|$ in Eq. (10) denotes the Euclidean norm. More generally, we use $\left\|\bm{v}\right\|_{p}$ to denote the $p$-norm of vector $\bm{v}$. Since Eq. (10) only saturates Eq. (8) in trivial cases where $\nabla f(\bm{\theta})$ is zero in all but one component, the unentangled protocol described is not optimal.

IV Two-step Protocol

We now present a protocol which asymptotically saturates Eq. (8). Our protocol consists of two steps. First, we make an unbiased estimate $\tilde{\bm{\theta}}$ of $\bm{\theta}$ for time $t_{1}$. Second, given our estimates $\tilde{\bm{\theta}}$, we make an unbiased measurement $\tilde{q}$ of the quantity $q=\nabla f(\tilde{\bm{\theta}})\cdot(\bm{\theta}-\tilde{\bm{\theta}})$ using the linear combination protocol in Ref. Eldredge et al. (2018), which takes time $t_{2}$. Our final estimate is $\tilde{f}=f(\tilde{\bm{\theta}})+\tilde{q}$.

It can be shown that our protocol is optimal (in terms of scaling with the total time $t_{1}$ + $t_{2}$) provided that the individual estimations of the parameters satisfy $\mathbb{E}[(\tilde{\theta}_{i}-\theta_{i})^{4}]=\mathcal{O}(t_{1}^{-4})$ and that $t_{1}$ and $t_{2}$ are chosen properly. To simplify our computations, we will make the more concrete assumption that our initial estimates $\tilde{\bm{\theta}}$ are each normally distributed as $\mathcal{N}(\theta_{i},\operatorname{Var}\tilde{\theta}_{i})$. Then as computed in the Appendix, the figure of merit for this protocol is

[TABLE]

In Eq. (12), the first term is the error resulting from the second phase of the protocol, estimating the linear combination. The second term is a residual error remaining from the first phase of the protocol after it is corrected by the linear combination measurement.

For our particular Hamiltonian $\hat{H}=\frac{1}{2}\theta_{i}\hat{\sigma}_{i}^{z}$, as per Ref. Eldredge et al. (2018), we know that the minimum variance of an unbiased estimator of some linear combination $\bm{\alpha}\cdot\bm{\theta}$ given time $t$ is

[TABLE]

which can be achieved with the entangled GHZ state $|\psi_{\mathrm{spin}}\rangle=\frac{1}{\sqrt{2}}(|0\rangle^{\otimes d}+|1\rangle^{\otimes d})$. We can apply this linear combination protocol to the second phase of our protocol by setting $\bm{\alpha}=\nabla f(\tilde{\bm{\theta}})$. For the individual estimators of the first phase, we use the fact that an individual estimation can be made in time $t$ with variance $1/t^{2}$ Bollinger et al. (1996). Using these results, we simplify Eq. (12):

[TABLE]

where we have absorbed the second derivatives into $g_{1}(\bm{\theta})$, which does not depend on time. Without loss of generality, we designate $f_{1}(\tilde{\bm{\theta}})$ as the largest $f_{i}(\tilde{\bm{\theta}})$. We then expand $\mathbb{E}[f_{1}(\tilde{\bm{\theta}})^{2}]$ as

[TABLE]

We may substitute Eq. (16) into Eq. (15) to obtain

[TABLE]

where $g_{2}(\bm{\theta})=f_{1}(\bm{\theta})^{2}$ and $g_{3}(\bm{\theta})$ have been introduced to absorb more time-independent factors.

IV.1 Optimal time allocation

To complete the protocol, we must specify how the total time $t_{\text{total}}$ is to be allocated between $t_{1}$ and $t_{2}$. We want to choose the $t_{1},t_{2}$, under the constraint that $t_{1}+t_{2}=t_{\text{total}}$, which minimize the MSE

[TABLE]

Notice that the $g_{1},g_{2},g_{3}$ functions are only dependent on $\bm{\theta}$ and not $t_{1}$, so we may set the derivative of $M$ with respect to $t_{1}$ equal to [math] and obtain

[TABLE]

Let $r=t_{1}/t_{2}$. Then we may rearrange to obtain

[TABLE]

Since $t_{1}\gg 1$, then $r\ll 1$, so the $r^{-3}$ term dominates the RHS. Thus, $g_{2}(\bm{\theta})t_{1}^{2}\approx\frac{2g_{1}(\bm{\theta})}{r^{3}}$, which implies

[TABLE]

Therefore the best possible allocation satisfies

[TABLE]

where $g$ is a function which depends only on $f$ and $\bm{\theta}$. In particular, $t_{1}=\mathcal{O}(t_{\text{total}}^{3/5})$, so the fraction of time spent on $t_{1}$ vanishes as $t_{\text{total}}\to\infty$. Almost all of the time is spent on $t_{2}$, the linear combination step of the two-step protocol. It can readily be shown that Eq. (17) is asymptotically dominated by the first term when this time allocation is chosen, which (since $t_{2}\to t_{\text{total}}$) is equal to the right-hand-side of the bound in Eq. (8). In other words, this distribution of time asymptotically achieves the optimal MSE.

The two-step protocol exhibits Heisenberg scaling as defined for distributed sensing Eldredge et al. (2018); Ge et al. (2018); Proctor et al. (2018). Comparing Eq. (10) to Eq. (8) shows an improvement of $\mathcal{O}(d)$, maximized when all components of $\nabla f(\bm{\theta})$ are approximately equal. Intuitively, the advantage is maximal when all parameters contribute, but minimal (i.e. no advantage) when only one parameter affects the function value. Similar behavior was noted in the linear combination case Eldredge et al. (2018).

Note that when actually implementing the protocol, the optimal $t_{1}$ is unknown since the function $g$ that determines it depends on the true parameters $\bm{\theta}$. However, we do not need to use the optimal $t_{1}$ to saturate the bound in Eq. (8). If $t_{1}$ is a function $ct_{\text{total}}^{p}$ of the total time where $\frac{1}{2}<p<1$ and some constant $c$, then the protocol will saturate Eq. (12). Suppose that $t_{1}=ct_{\text{total}}^{p}$ for some $\frac{1}{2}<p<1$ and some constant $c$. Since $p<1$, we see that $\lim_{t\to\infty}\frac{t_{2}}{t_{\text{total}}}=1$. Therefore, we may substitute our $t_{1}$ into the MSE formula in Eq. (17) and simplify:

[TABLE]

Since $p>\frac{1}{2}$, the $t_{\text{total}}^{2}$ term is dominant. Thus, as we defined $g_{2}:=f_{1}(\bm{\theta})^{2}=\max_{i}f_{i}(\bm{\theta})^{2}$ under the assumption that $f_{1}(\bm{\theta})^{2}$ was maximal, our asymptotic error is

[TABLE]

which saturates the bound of Eq. (8). Although selecting a non-optimal time allocation does result in a higher MSE, the additional error is $\mathcal{O}\left(t_{\text{total}}^{-4}\right)$, which is insignificant asymptotically. The two-step protocol will therefore be asymptotically optimal for a wide range of time allocations.

V Function Measurement in Other Physical Settings

We now consider a different physical setting for function estimation. Rather than $d$ qubits which accumulate phase for some time $t$, we instead pass $n$ photons through $d$ Mach-Zehnder interferometers and accumulate some fixed phase $\theta_{i}$ encoded into each interferometer (see Fig. 2). For single parameters, the use of entangled states to reduce noise in this setting has been explored in Refs. Holland and Burnett (1993); Kim et al. (1998); Usha Devi and Rajagopal (2009); Demkowicz-Dobrzański et al. (2015); Dinani et al. (2016) with multiparameter cases explored in Refs. Ge et al. (2018); Proctor et al. (2018). In this setting, the relevant limitation is the total number of photons used in the measurement, rather than time. This constraint is particularly relevant when analyzing a biological or chemical sample which is sensitive to light, making it desirable to reduce noise with as few photons as possible. Similar biologically motivated situations are presented in Refs. Kee and Cicerone (2004); Alem et al. (2015); Jensen et al. (2016).

For photons, a two-step protocol with similar structure to the protocol for qubits yields reduced noise compared to any estimate of $f$ derived entirely from local measurements. Suppose we allot $N_{1}$ photons for the first step (individual measurement) and $N_{2}$ photons for the second step (linear combination), for a total of $N_{\text{total}}=N_{1}+N_{2}$ photons. We again begin from the general result of Eq. (12). However, the use of photons which can be apportioned between modes introduces new structure to the problem. We need to partition the $N_{1}$ photons into $N_{1}=n_{1}+\dots+n_{d}$, putting $n_{i}$ photons into the $i$-th interferometer, as some parameters may affect our final result more than others. Thus, in the second term of Eq. (12), we replace $\operatorname{Var}\tilde{\theta}_{i}$ with $\frac{1}{n_{i}^{2}}$ instead of $\frac{1}{t_{1}^{2}}$ Holland and Burnett (1993).

The optimal variance when measuring the linear combination $\bm{\alpha}\cdot\bm{\theta}$ using $N$ total photons is unknown. However, Ref. Proctor et al. (2018) conjectures the optimal variance to be

[TABLE]

Furthermore, Ref. Proctor et al. (2018) provides a protocol achieving the bound in Eq. (25) using a proportionally weighted GHZ state: $|\psi_{\mathrm{photon}}\rangle=\frac{1}{\sqrt{2}}(|n_{1},0,n_{2},0\dots\rangle+|0,n_{1},0,n_{2},\dots\rangle)$, where $n_{i}=N_{\text{total}}\frac{\alpha_{i}}{\sum\alpha_{j}}$ and where, in reference to Fig. 2, the modes are listed from top to bottom. Note that this will only work for $\bm{\alpha}$ proportional to some rational vector as photons are discrete. Since Eq. (25) is saturable, we may simplify the first term of Eq. (12) to obtain

[TABLE]

For fixed $f$ and $\theta$, the $\frac{1}{n_{i}n_{j}}$ terms in Eq. (26) are minimized for the same ratio of $n_{1}:n_{2}:\dots:n_{d}$ regardless of the value of the total number of photons used, $N_{1}$. Each term is proportional to $N_{1}^{-4}$ multiplied by some function of $f,\bm{\theta}$, and $d$. Therefore, the structure of Eq. (26) becomes identical to the structure of Eq. (17), with $N_{1}$ and $N_{2}$ replacing $t_{1}$ and $t_{2}$. As a result, the optimal allocation of photons between $N_{1}$ and $N_{2}$ will yield $N_{1}=\mathcal{O}(N_{\text{total}}^{3/5})$ and $N_{2}=\mathcal{O}(N_{\text{total}})$, meaning that the $N_{2}^{-2}$ term in Eq. (26) is dominant asymptotically. Therefore, for photons, we may asymptotically achieve

[TABLE]

This strategy is optimal if the linear combination estimation strategy presented in Ref. Proctor et al. (2018) is optimal, as conjectured in that work. We stress that our optimality result remains true for spins evolving under Eq. (1) and it is only for photons that our protocol is only conjectured to be optimal.

Eq. (27) also exhibits Heisenberg scaling. Suppose we were to measure each parameter individually and then calculate the function. When measuring the parameters individually, we obtain the same error formula as Eq. (9), except now we set $\operatorname{Var}\tilde{\theta}_{i}=\frac{1}{n_{i}^{2}}$ to get

[TABLE]

The optimal distribution requires an $n_{i}$ proportional to the weight $f_{i}(\bm{\theta})^{2/3}$, yielding an entanglement-free error of

[TABLE]

As with qubits, by comparing Eq. (27) with Eq. (29) in the case where all of the $f_{i}(\bm{\theta})$ are approximately equal, we find that the photonic two-step protocol yields a $\mathcal{O}(d)$ improvement in error over measuring each parameter individually. This improvement when all quantities are equally important can also be seen in Ref. Ge et al. (2018) for the special case of $f$ being a linear combination. As in the qubit case, the improvement in error is lessened when $\nabla f(\bm{\theta})$ is not approximately equal in all components.

In fact, this method can be extended still more generally. Rather than cases where the signal is imprinted on photons by a phase shift, we can consider the protocol developed in Ref. Zhuang et al. (2018), which is capable of entanglement-enhanced distributed sensing of continuous variables by using homodyne measurements. Besides measuring parameters in different physical settings, we may also measure functions of variables coupled to spins, phase-shifts of photons, continuous variables, and any combination of these. In such a hybrid scenario, we can still make use of the two-step protocol. The first step, obtaining initial estimates for the individual parameters, proceeds equivalently, since the measurements of the spins and of the photons can be viewed as occurring in parallel. For the linear combination case, we can assume that the optimal spin and photon input states can be entangled as follows:

[TABLE]

Here, $n_{i}=N_{\text{total}}\frac{\alpha_{i}}{\sum\alpha_{j}}$, where the sum runs over only the $j$ corresponding to photonic modes, denotes the number of photons which pass through the arms of the $i$-th interferometer. The state in Eq. (30) is designed in such a way that the two branches of the overall wavefunction accumulate relative to each other a phase equal to the total linear combination we are interested in. In order to extract this final phase, the state can be unitarily mapped onto a qubit, which contains all of the accumulated phase and is then measured.

One caveat is that the linear combination protocol will accumulate phase proportional to time for the qubits and phase proportional to the number of photons for interferometers. For instance, if $\theta_{1}$ is coupled to a qubit (and therefore has units of frequency) and $\theta_{2}$ is coupled to an interferometer (and is therefore unitless), then the two branches of our state accumulate a relative phase $\theta_{1}t+\theta_{2}n$. Therefore, one may have to adjust $t$ or $n$ in order to get the desired linear combination.

VI Applications

Our protocol is capable of estimating any analytic function of the inputs, allowing for a large variety of potential applications. Essentially, any time multiple sensors are processed into a single signal, our protocol provides enhanced sensitivity using entanglement. In fact, there is no requirement that different $\theta_{i}$ have the same physical origin. For instance, a $\theta_{1}$ representing an electric field and $\theta_{2}$ measuring a magnetic field could be used to measure the Poynting vector.

One potential application of function measurements is the interpolation of non-linear functions. Suppose that an ansatz with $d$ tunable parameters is made for the strength of the field in a region. With readings from $\geq d$ different points, one could determine the parameters of the ansatz and therefore determine the value of the field at other points. Estimations of these ansatz parameters, which are functions of the measured fields, may potentially be improved using entangled states depending on the figure of merit Chiribella et al. (2005); Baumgratz and Datta (2015). Note that this procedure can be carried out even if it is difficult to invert the ansatz in terms of the $d$ measurements. Suppose that $\bm{\theta}=f(\bm{c},\bm{x})$ and that $\bm{c}=f^{-1}(\bm{\theta},\bm{x})$ exists, but has no closed-form solution which can be easily evaluated. First, we make measurements $\hat{\bm{\theta}}$. To create an initial estimate of the values $\bm{c}$, we use a numerical root-finder to find estimates $\tilde{\bm{c}}$. We can now implement the second step of our protocol by finding the first derivatives $\partial c_{i}/\partial\theta_{j}$ using the matrix identity $\frac{\partial\bm{\theta}}{\partial\bm{c}}\cdot\frac{\partial\bm{c}}{\partial\bm{\theta}}=I$. Since $f$ is known, $\partial\bm{\theta}/\partial\bm{c}$ can be inverted to yield the $\partial\bm{c}/\partial\bm{\theta}$ needed to estimate $\hat{\bm{q}}=\partial\bm{c}/\partial\bm{\theta}|_{\bm{\theta}=\hat{\bm{\theta}}}\cdot(\bm{\theta}-\hat{\bm{\theta}})$. Our final estimate is $\hat{\bm{c}}+\hat{\bm{q}}$, which was obtained without having to compute $f^{-1}$ in general.

Interpolation in this manner can proceed by two different schemes. We can either attempt to measure the ansatz parameters themselves, which allows computation of the field at all other points, or we can skip the final computation step by writing the field at a point of interest as a function of all the points that can be measured. This final function can then be directly measured using an entangled protocol, which will be more accurate. However, the first approach has the advantage that knowing the ansatz parameters allows estimation of all points in the space in question.

One particular interpolation of interest arises in ion trap quantum computing. In trapped ion chains, qubits are manipulated using Gaussian laser beams, and two primary sources of error are intensity and beam pointing fluctuations Cirac and Zoller (1995); Häffner et al. (2008); Brown et al. (2011). Our protocol offers better ways to characterize this noise. In order to detect the field error at a qubit’s position without disturbing the qubit, we can perform interpolation by measuring the field’s effect on other ions, possibly of a different atomic species, positioned nearby. Given the ansatz of the Gaussian beam profile, we are able to calculate the field at the qubit of interest and perhaps correct the error. As entanglement of ions is already a key functionality for trapped ion quantum computers, our proposal is immediately applicable in that domain.

VII Outlook

We have presented a Heisenberg-scaling measurement protocol using quantum sensor networks for measuring any multivariate, real-valued, analytic function, and this protocol is consistent with the Heisenberg limit when measuring functions with comparably-sized gradients in each component. Recent advances in the distribution of entanglement, for instance, in satellites distributing entangled photons more than 1000km Yin et al. (2017), strengthen the viability of this scheme over large distances in the near-term. Potential sensing platforms include trapped ions and nitrogen-vacancy defects in diamond, which can also be entangled Leibfried et al. (2004); Bernien et al. (2013); Dolde et al. (2014); Hensen et al. (2015) and are proven platforms for magnetometry and thermometry Taylor et al. (2008); Tzeng et al. (2015). Future work may include proving the optimality of the two-step protocol when constrained by the number of photons, which would require extending the results of Ref. Proctor et al. (2018), as well as further experimental research into quantum networking to explore how entanglement can be reliably distributed for metrological purposes.

We specifically identified field interpolation as a promising application of our work, but we stress that our protocol can assist in the measurement of any analytic function. More work remains to determine when it is optimal to measure the coefficients of interpolation and when it is optimal to directly measure the final function. We are also interested in fleshing out possible intersections between quantum function estimation and machine learning. Supervised machine learning is a type of interpolation: estimating functional outputs for unknown inputs by extracting information from known input-output pairs Russell and Norvig (2010). It is possible our protocol could be used to improve the accuracy of training a machine learning model if the necessary quantity for training was a function of physical measurements. Additionally, the final output of many machine learning algorithms, such as neural networks, is a non-linear but infinitely differentiable function of the inputs Schmidhuber (2015). Our work could aid in computing this complicated function for new input when making predictions.

Acknowledgements.

We would like to thank M. Foss-Feig, S. Rolston, J. Gross, and S. Kimmel for helpful discussions. This work was supported by ARL CDQI, AFOSR, ARO MURI, ARO, NSF PFC at JQI, NSF PFCQC program, DoE ASCR Quantum Testbed Pathfinder program (Award No. DE-SC0019040), and the DoE BES QIS program (Award No. DE-SC0019449). Z.E. was supported in part by the ARCS Foundation.

Appendix A Figure of Merit for Two-Step Protocol

In this section, we derive Eq. (12) in the main text. Specifically, we derive the figure of merit for the two-step protocol in terms of the measurement accuracy of the independent parameters and the measurement accuracy of the linear combination, yielding a general formula which applies to any physical realization.

For the sake of concision, let $\bm{\Delta}=\tilde{\bm{\theta}}-\bm{\theta}$ which satisfies $\mathbb{E}[\bm{\Delta}]=\bm{0}$. Furthermore, let $T_{k}$ be $k!$ times the $k$-th term of the Taylor expansion of $f$ (so $T_{1}=f_{i}(\bm{\theta})\Delta_{i}$, $T_{2}=f_{ij}(\bm{\theta})\Delta_{i}\Delta_{j}$, $T_{3}=f_{ijk}(\bm{\theta})\Delta_{i}\Delta_{j}\Delta_{k}$, etc.). Thus, the Taylor expansion of $f(\tilde{\bm{\theta}})$ would be

[TABLE]

We compute our figure of merit:

[TABLE]

The actual computation of the labeled terms is rather involved and space consuming, so it is presented in Sec. A.1). Notice that we may simplify

[TABLE]

so Eq. (35) evaluates to

[TABLE]

since $\mathbb{E}[T_{2}T_{3}]$ is $\mathcal{O}(\bm{\Delta}^{5})$. Now, this simplifies further as

[TABLE]

since all terms with some $\Delta_{i}$ to a single power will factor out as $\mathbb{E}[\Delta_{i}]=0$. We will assume that $\Delta_{i}\sim\mathcal{N}(0,\frac{1}{t_{1}^{2}})$ is normally distributed. This is not strictly necessary as long as the distribution of errors satisfies $\mathbb{E}[\Delta_{i}^{4}]\leq\mathcal{O}(t_{1}^{-4})$, a condition that is satisfied by phase estimation procedures like those in Ref. Kimmel et al. (2015). However, assuming normality allows the calculation to proceed easily, as we will be able to simplify $\mathbb{E}[\Delta_{i}^{4}]=3\operatorname{Var}\tilde{\theta}_{i}^{2}$. Thus, we arrive at

[TABLE]

A.1 Simplification of labeled terms

In this subsection, we present the simplification of the labeled terms from Eqs. (33-35) in full detail.

Term 2 is simplified by using the definition of $\operatorname{Var}_{\tilde{q}}\tilde{q}$. One needs to be careful as there are two layers of expected values - one for the values of $\tilde{\bm{\theta}}$ and one for the estimator $\tilde{q}$:

[TABLE]

Terms 1 and 3 are simplified by expanding the Taylor series for $f(\tilde{\theta})$ up to $\bm{\Delta}^{4}$ terms; note that $q=\mathcal{O}(\bm{\Delta})$, so we only need to expand the Taylor series up to $\mathcal{O}(\bm{\Delta}^{3})$ terms:

[TABLE]

[TABLE]