Ultimate limits for quantum magnetometry via time-continuous measurements

Recent developments in the field of quantum metrology have shown how quantum probes and quantum measurements allow one to achieve parameter estimation with precision beyond that obtainable by any classical scheme [1, 2]. The estimation of the strength of a magnetic field is a paradigmatic example in this respect, as it can be mapped to the problem of estimating the Larmor frequency for an atomic spin ensemble [3–9].

As a matter of fact, if the system is initially prepared in a spin coherent state, the mean-squared error of the field estimate scales, in terms of the total spin number J, as $1/J$, which is usually referred to as the standard quantum limit (SQL) to precision. If quantum resources, such as spin squeezing or entanglement between the atoms of the spin ensemble, are exploited, one observes a quadratic enhancement and achieves the so-called Heisenberg scaling, i.e. $1/{J}^{2}$ [10, 11]. On the other hand, it has been proved that such ultimate quantum limit may be easily lost in the presence of noise [12] and that typically a SQL-like scaling is observed, with the quantum enhancement reduced to a constant factor. These observations have been rigorously translated into a set of no-go theorems [13, 14], which fostered several attempts to circumvent them. In particular, it has been shown how one can restore a super-classical scaling in the context of frequency estimation, for specific noisy evolution and/or by optimizing the strategy over the interrogation time [15–19], or by exploiting techniques borrowed from the field of quantum error-correction [20–22].

In this manuscript, we put forward an alternative approach: we assume to start the dynamics with a classical state that is monitored continuously in time via the interacting environment [23, 24]. The goal is to recover the information on the parameter leaking into the environment and simultaneously to exploit the back action of the measurement to drive the system into more sensitive conditional states [25–34]. This approach has received much attention recently [35–42] also in the context of quantum magnetometry [43–48].

Here we rigorously address the performance of these protocols, depicted in figure 1, taking into account the information obtained via the time-continuous non-demolition measurements on the environment, as well as the information obtainable via a strong (destructive) measurement on the conditional state of the system. In particular, in the limit of large spin, we derive an analytical formula for the ultimate bound on the mean-squared error of any unbiased estimator, and conclusively show that, for experimentally relevant values of the dynamical parameters, one can observe a Heisenberg-like scaling.

Zoom In Zoom Out Reset image size

Remarkably, at variance with most of the protocols proposed for quantum magnetometry, and in general for frequency estimation, one does not need to prepare an initial spin-squeezed state. The Heisenberg scaling is in fact obtained also for an initial classical spin coherent state, thanks to the spin squeezing generated by time-continuous measurements' back-action. Finally, we analytically prove that the ultimate quantum limit for noisy magnetometry in the presence of collective transverse noise [36] is in fact saturated by our strategy, i.e. one does not need to implement more involved strategies, e.g. jointly measuring the conditional state of the system and the output modes of the environment at different times.

The paper is organized as follows. In section 2, we present the quantum Cramér–Rao bounds that hold for noisy metrology, with emphasis on estimation strategies based on time-continuous, non-demolition measurements and a final strong measurement on the corresponding conditional quantum states. In section 3, we introduce the physical setting for the estimation of a magnetic field via a continuously monitored atomic ensemble. In particular, we focus on the case of large total spin, where a Gaussian picture is able to describe the whole dynamics. In section 4, we present the main results: we first calculate the classical Fisher information (FI) corresponding to the photoccurent obtained via the time-continuous monitoring of the environment, and we discuss how to attain the corresponding bound via Bayesian estimation. We then address the possibility of performing also a strong measurement on the conditional state of the atomic ensemble, and derive the ultimate limit on this kind of estimation strategy, quantified by an effective quantum Fisher information (QFI). Upon studying this quantity, we observe how, in the relevant parameters' regime, the Heisenberg limit can be effectively restored, also discussing the effects of non-unit monitoring efficiency, corresponding to the loss of photons before the detection. Finally, we also prove the optimality of our measurement strategy in the case of ideal detectors. Section 5 closes the paper with some concluding remarks.

A classical estimation problem consists in inferring the value of a parameter θ from a number M of measurement outcomes $\chi =\{{x}_{1},{x}_{2},\,\ldots ,\,{x}_{M}\}$ and their conditional distribution $p(x| \theta )$. We define an estimator $\hat{\theta }(\chi )$ a function from the measurement outcomes to the possible values of θ and we dub it asymptotically unbiased when, in the limit of large number of repetitions of the experiment M, its average is equal to the true value, i.e. $\int {\rm{d}}\chi \,p(\chi | \lambda )\hat{\theta }(\chi )=\theta$, where $p(\chi | \lambda )={{\rm{\Pi }}}_{j=1}^{M}p({x}_{j}| \theta )$. The Cramér–Rao theorem states that the variance of any unbiased estimator is lower bounded as ${\mathrm{Var}}_{\hat{\theta }}(\theta )={(M{ \mathcal F }[p(x| \theta )])}^{-1}$, where

$$\begin{eqnarray}&&{ \mathcal F }[p(x| \theta )]=\int {\rm{d}}{x}\,p(x| \lambda ){({\partial }_{\theta }\mathrm{log}p(x| \lambda ))}^{2}\end{eqnarray} \tag{ 1 }$$

denotes the classical FI.

In the quantum realm, the conditional probability distribution reads $p(x| \theta )=\mathrm{Tr}[{\varrho }_{\theta }{{\rm{\Pi }}}_{x}]$, where ${\varrho }_{\theta }$ is the quantum state of the system labeled by the parameter θ, and ${{\rm{\Pi }}}_{x}$ is a POVM operator describing the quantum measurement. One can prove that the FI corresponding to any POVM is upper bounded ${ \mathcal F }[p(x| \theta )]\leqslant { \mathcal Q }[{\varrho }_{\theta }]$, where ${ \mathcal Q }[{\varrho }_{\theta }]=\mathrm{Tr}[{\varrho }_{\theta }{L}_{\theta }^{2}]$ is the QFI, and L_θ is the so-called symmetric logarithmic derivative, which can be obtained by solving the equation $2{\partial }_{\theta }{\varrho }_{\theta }={L}_{\theta }{\varrho }_{\theta }+{\varrho }_{\theta }{L}_{\theta }$ [49–51]. The QFI depends on the quantum state ${\varrho }_{\theta }$ only, and thus poses the ultimate bound on the precision of the estimation of θ. Moreover, in the single parameter case the bound is always achievable, that is, there exists a (projective) POVM such that the corresponding classical FI equals the QFI.

In this manuscript we consider a quantum system evolving according to a given Hamiltonian ${\hat{H}}_{\theta }$ characterized by the parameter we want to estimate, and coupled to a bosonic environment at zero temperature described by a train of input operators ${\hat{a}}_{{\rm{in}}}(t)$, satisfying the commutation relation $[{\hat{a}}_{{\rm{in}}}(t),{\hat{a}}_{{\rm{in}}}^{\dagger }(t^{\prime} )]=\delta (t-t^{\prime} )$, via an interaction Hamiltonian ${\hat{H}}_{{\rm{int}}}(t)=\hat{c}{\hat{a}}_{{\rm{in}}}^{\dagger }(t)+{\hat{c}}^{\dagger }{\hat{a}}_{{\rm{in}}}(t)$ ($\hat{c}$ being a generic operator in the system Hilbert space) [23]. By tracing out the environment, the unconditional dynamics of the system is described by the Lindblad master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\varrho }{{\rm{d}}{t}}={ \mathcal L }\varrho =-{\rm{i}}[{\hat{H}}_{\theta },\varrho ]+{ \mathcal D }[\hat{c}]\varrho ,\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal D }[c]\varrho =\hat{c}\varrho {\hat{c}}^{\dagger }-({\hat{c}}^{\dagger }\hat{c}\varrho +\varrho {\hat{c}}^{\dagger }\hat{c})/2$.

If one performs a homodyne detection of a quadrature ${\hat{x}}_{{\rm{out}}}(t)={\hat{a}}_{{\rm{out}}}(t)+{\hat{a}}_{{\rm{out}}}^{\dagger }(t)$ on the output operators, i.e. on the environment just after the interaction with the system, one obtains that the dynamics of the system quantum state ${\varrho }^{(c)}$ conditioned on the measurement results (we will omit the dependence of the measured photocurrent ${{\bf{y}}}_{t}$), is described by the stochastic master equation [23]

$$\begin{eqnarray}&&{\rm{d}}{\varrho }^{(c)}=-{\rm{i}}[{\hat{H}}_{\theta },{\varrho }^{(c)}]\,{\rm{d}}{t}+{ \mathcal D }[\hat{c}]{\varrho }^{(c)}\,{\rm{d}}{t}+\sqrt{\eta }{ \mathcal H }[\hat{c}]{\varrho }^{(c)}\,{{\rm{d}}{w}}_{t}.\end{eqnarray} \tag{ 3 }$$

Here η denotes the efficiency of the detection, dw_t is a stochastic Wiener increment (s.t. ${{\rm{d}}{w}}_{t}^{2}={\rm{d}}{t}$), and ${ \mathcal H }[c]{\varrho }^{(c)}=\hat{c}{\varrho }^{(c)}+{\varrho }^{(c)}{\hat{c}}^{\dagger }-\mathrm{Tr}[{\varrho }^{(c)}(\hat{c}+{\hat{c}}^{\dagger })]{\varrho }^{(c)}$ (notice that in principle one could consider other measurement strategies different from homodyne, yielding a different superoperator). The corresponding measurement record during a time step $t\to t+{\rm{d}}{t}$ is given by the infinitesimal current

$$\begin{eqnarray}&&{{\rm{d}}{y}}_{t}=\sqrt{\eta }\,\mathrm{Tr}[{\varrho }^{(c)}(\hat{c}+{\hat{c}}^{\dagger })]\,{\rm{d}}{t}+{{\rm{d}}{w}}_{t}.\end{eqnarray} \tag{ 4 }$$

With the help of such measurement strategies, one can estimate the value of the parameter θ both from the measured photocurrent ${{\bf{y}}}_{T}={\int }_{0}^{T}{{\rm{d}}{y}}_{t}$, and from a final strong (destructive) measurement on the conditional state ${\varrho }^{(c)}$. In this case, as we explicitly show in appendix A (in general both for the classical and quantum case), the proper quantum Cramér–Rao bound reads

$$\begin{eqnarray}&&{\mathrm{Var}}_{\hat{\theta }}(\theta )\geqslant \displaystyle \frac{1}{M({ \mathcal F }[p({{\bf{y}}}_{T})]+{{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal Q }[{\varrho }^{(c)}]])},\end{eqnarray} \tag{ 5 }$$

where the first term at the denominator ${ \mathcal F }[p({{\bf{y}}}_{T})]$ is the FI corresponding to the classical photocurrent ${{\bf{y}}}_{T}$, while the second term is the average of the QFI for the conditional state ${ \mathcal Q }[{\varrho }^{(c)}]$ over all the possible trajectories, i.e. on all the possible measurement outcomes for the photocurrent.

The classical FI ${ \mathcal F }[p({{\bf{y}}}_{T})]$ can be calculated as described in [35] by evaluating

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{T})]={{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[\mathrm{Tr}{[\tau ]}^{2}],\end{eqnarray} \tag{ 6 }$$

where the operator τ evolves according to the stochastic master equation

$$\begin{eqnarray}&&{\rm{d}}\tau =-{\rm{i}}[{\hat{H}}_{\theta },\tau ]\,{\rm{d}}{t}-{\rm{i}}[({\partial }_{\theta }{\hat{H}}_{\theta }),\varrho ]\,{\rm{d}}{t}+{ \mathcal D }[\hat{c}]\tau \,{\rm{d}}{t}+(\hat{c}\tau +\tau {\hat{c}}^{\dagger }){{\rm{d}}{w}}_{t}.\end{eqnarray} \tag{ 7 }$$

The conditional states ${\varrho }^{(c)}$ at time T can be obtained by integrating (3), for a certain stream of outcomes ${{\bf{y}}}_{T};$ then one can first calculate the corresponding QFI ${ \mathcal Q }[{\varrho }^{(c)}]$, and, numerically or when possible analytically, its average over all the possible trajectories explored by the quantum system due to the homodyne monitoring.

A more fundamental quantum Cramér–Rao bound that applies in this physical setting has been derived in [36], by considering the QFI obtained from the unitary dynamics of the global pure state of system and environment. This QFI is obtained by optimizing over all possible POVMs, i.e. one also considers the possibility of performing non-separable (entangled) measurements over the system and all the output modes ${\hat{a}}_{{\rm{out}}}(t)$ at different times. On the other hand, in the previous setting the estimation strategies were restricted to the more experimentally friendly case of sequential/separable measurements on the output modes and on the final conditional state of the system.

The QFI expressing this ultimate QCRB is by definition

$$\begin{eqnarray}&&{\overline{{ \mathcal Q }}}_{{ \mathcal L }}(\theta )=4{\partial }_{{\theta }_{1}}{\partial }_{{\theta }_{2}}\mathrm{log}(| \langle \psi ({\theta }_{1})| \psi ({\theta }_{2})\rangle | ){| }_{{\theta }_{1}={\theta }_{2}=\theta },\end{eqnarray} \tag{ 8 }$$

where $\langle \psi ({\theta }_{1})| \psi ({\theta }_{2})\rangle$ is the fidelity between the global state of system and environment for two different values of the parameter, and where we have highlighted its dependence on the superoperator ${ \mathcal L }$ that defines the unconditional master equation (2). The key insight is that this fidelity can be determined by using operators acting on the system only [36, 52] and it can be expressed as the trace of an operator $\mathrm{Tr}[\bar{\rho }]=\langle \psi ({\theta }_{1})| \psi ({\theta }_{2})\rangle$, which obeys the following generalized master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{\rho }}{{\rm{d}}{t}}=-{\rm{i}}({\hat{H}}_{{\theta }_{1}}\bar{\rho }-\bar{\rho }{\hat{H}}_{{\theta }_{2}})+{ \mathcal D }[\hat{c}]\bar{\rho }.\end{eqnarray} \tag{ 9 }$$

As before, we already assumed that the dependence on the parameter lies only in the system Hamiltonian ${\hat{H}}_{\theta }$ and that we have a single jump operator $\hat{c}$. We remark that the operator $\bar{\rho }$ is not a proper density operator representing a quantum state, except in the limit case ${\theta }_{1}\to {\theta }_{2}$, where we recover the standard master equation (2).

We address the estimation of the intensity of a static and constant magnetic field B acting on a ensemble of N two-level atoms that are continuously monitored [43–46], as depicted in figure 1. The atomic ensemble can be described as a system with total spin $J=N/2$ with collective spin operators defined as ${\hat{J}}_{\alpha }=\tfrac{1}{2}{\sum }_{i=0}^{N}{\sigma }_{i\alpha }$, where $\alpha =x,y,z$ and ${\sigma }_{i\alpha }$ denotes the Pauli matrices acting on the ith spin. The collective operators obey the same angular momentum commutation rules $[{\hat{J}}_{{\rm{i}}},{\hat{J}}_{j}]={\rm{i}}{\varepsilon }_{{ijk}}{\hat{J}}_{k}$, where ${\varepsilon }_{{ijk}}$ is the Levi-Civita symbol. We remark that in the present manuscript we choose units such that ${\hslash }=1$.

We assume that the atomic sample is coupled to a electromagnetic mode a_in(t) corresponding either to a cavity mode in a strongly driven and heavily damped cavity [26], or analogously to a far-detuned traveling mode passing through the ensemble [46]. By considering an interaction Hamiltonian ${\hat{H}}_{{\rm{int}}}=\sqrt{\kappa }{\hat{J}}_{z}({\hat{a}}_{{\rm{in}}}(t)+{\hat{a}}_{{\rm{in}}}^{\dagger }(t))$ and if these environmental light modes are left unmeasured, the evolution of the system is expressed by (2), which in this case corresponds to a collective transverse noise on the atomic sample,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\varrho }{{\rm{d}}{t}}={{ \mathcal L }}_{{tn}}\varrho =-{\rm{i}}\gamma B[{\hat{J}}_{y},\varrho ]+\kappa { \mathcal D }[{\hat{J}}_{z}]\varrho ,\end{eqnarray} \tag{ 10 }$$

where the constants κ and γ represent respectively the strength of the coupling with the noise and with the magnetic field, that is directed on the y-axis and thus perpendicular to the noise generator. At t = 0 we consider the system prepared in a spin coherent state, i.e. a tensor product of single spin states (qubits) directed in the positive x direction,

$$\begin{eqnarray}&&| \psi (0)\rangle =\underset{k=0}{\overset{N}{\displaystyle \bigotimes }}| +{\rangle }_{k}=| J,J{\rangle }_{x},\end{eqnarray} \tag{ 11 }$$

where $| +\rangle$ is the eigenstate of ${\sigma }_{x}$ with eigenvalue +1. We thus have that the spin component on the x direction attains the macroscopic value $\langle {\hat{J}}_{x}(0)\rangle =J$. The unconditional dynamics of $\langle {\hat{J}}_{x}\rangle$ is obtained by applying the operator ${\hat{J}}_{x}$ to both sides of equation (10) and then taking the trace. The result is the following equation describing damped oscillations

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {\hat{J}}_{x}(t)\rangle }{{\rm{d}}{t}}=\gamma B\langle {\hat{J}}_{z}(t)\rangle -\displaystyle \frac{\kappa }{2}\langle \,{\hat{J}}_{x}(t)\rangle ,\end{eqnarray} \tag{ 12 }$$

where we observe how the the dissipative and unitary parts of the dynamics are respectively shrinking the spin vector $\langle \vec{\hat{J}}\rangle$ and causing its Larmor precession around the y-axis. In the following we will assume to measure small magnetic fields, such that $\gamma {Bt}\ll 1$ and we can approximate the solution of the previous equation as

$$\begin{eqnarray}&&\langle {\hat{J}}_{x}(t)\rangle \approx \langle {\hat{J}}_{x}(0)\rangle {{\rm{e}}}^{-\kappa t/2}={{J}{\rm{e}}}^{-\kappa t/2}.\end{eqnarray} \tag{ 13 }$$

If the light modes are continuously monitored via homodyne measurements at the appropriate phase, one allows a continuous 'weak' measurement of ${\hat{J}}_{z};$ the corresponding stochastic master equation (3) for finite monitoring efficiency η reads

$$\begin{eqnarray}&&{\rm{d}}{\varrho }^{(c)}=-{\rm{i}}\gamma B[{\hat{J}}_{y},{\varrho }^{(c)}]{\rm{d}}t+\kappa { \mathcal D }[{\hat{J}}_{z}]{\varrho }^{(c)}{\rm{d}}{t}+\sqrt{\eta \kappa }{ \mathcal H }[{\hat{J}}_{z}]{\varrho }^{(c)}{{\rm{d}}{w}}_{t},\end{eqnarray} \tag{ 14 }$$

while the measurement result at time t corresponds to an infinitesimal photocurrent ${{\rm{d}}{y}}_{t}=2\sqrt{\eta \kappa }\,\mathrm{Tr}[{\varrho }^{(c)}{\hat{J}}_{z}]{\rm{d}}{t}+{{\rm{d}}{w}}_{t}$. It is important to remark how the collective noise characterizing the master equation (10) describes the dynamics also in experimental situations where no additional coupling to the atomic ensemble, with the purpose of performing continuous monitoring, is engineered [53–55]. In this respect, assuming a non-unit efficiency η corresponds to considering both homodyne detectors that are not able to capture all the photons that have interacted with the spin, and environmental degrees of freedom, causing the same kind of noisy dynamics, that cannot be measured during the experiment.

Let us now consider the limit of large spin $J\gg 1$. In this case, the dynamics may be effectively described with the Gaussian formalism as long as $\langle {\hat{J}}_{x}(t)\rangle \approx J$, i.e. for times t small enough to guarantee that $\kappa t\lesssim 1$. We define the effective quadrature operators of the atomic sample, satisfying the canonical commutation relation $[\hat{X},\hat{P}]=i$, as [46, 47]

$$\begin{eqnarray}&&\hat{X}={\hat{J}}_{y}/\sqrt{{\bar{J}}_{t}}\qquad \hat{P}={\hat{J}}_{z}/\sqrt{{\bar{J}}_{t}},\end{eqnarray} \tag{ 15 }$$

where ${\bar{J}}_{t}\equiv | \langle {\hat{J}}_{x}(t)\rangle |$ (notice that in the limit of large spin J we can safely consider the unconditional average value $\langle {\hat{J}}_{x}(t)\rangle$, as the stochastic correction obtained via (14) would be negligible). In the Gaussian description the initial state $| \psi (0)\rangle$ corresponds to the vacuum state $(\hat{X}+{\rm{i}}\hat{P})| 0\rangle =| 0\rangle$, which is Gaussian. As the stochastic master equation (14) becomes quadratic in the canonical operators (and thus preserves the Gaussian character of states)

$$\begin{eqnarray}&&{\rm{d}}{\varrho }^{(c)}=-{\rm{i}}\gamma B\sqrt{{\bar{J}}_{t}}[\hat{X},{\varrho }^{(c)}]\,{\rm{d}}{t}+\kappa {\bar{J}}_{t}{ \mathcal D }[\hat{P}]\,{\varrho }^{(c)}{\rm{d}}{t}+\sqrt{{\bar{J}}_{t}\eta \kappa }\,{ \mathcal H }[\hat{P}]\,{\varrho }^{(c)}{{\rm{d}}{w}}_{t},\end{eqnarray} \tag{ 16 }$$

the whole dynamics can be equivalently rewritten in terms of first and second moments only [56, 57] (see appendix B for the equations describing the whole dynamics in the Gaussian picture). As it will be clear in the following, due to the nature of the coupling, in order to address the estimation of B, we only need the behavior of the mean and the variance of the atomic momentum quadrature $\hat{P}$ calculated on the conditional state ${\varrho }^{(c)}$, which follows the equations

$$\begin{eqnarray}&&{\rm{d}}\langle \hat{P}(t){\rangle }_{c}=-B\gamma \sqrt{{{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}{\rm{d}}{t}+2{\mathrm{Var}}_{c}[\hat{P}(t)]\sqrt{\eta \kappa {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}{{\rm{d}}{w}}_{t},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\mathrm{Var}}_{c}[\hat{P}(t)]}{{\rm{d}}{t}}=-4\eta \kappa {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}{({\mathrm{Var}}_{c}[\hat{P}(t)])}^{2}.\end{eqnarray} \tag{ 18 }$$

The differential equation for the conditional second moment is deterministic and can be solved analytically. For an initial vacuum state, i.e. with $\mathrm{Var}[\hat{P}(0)]=\tfrac{1}{2}$, we obtain the following solution

$$\begin{eqnarray}&&{\mathrm{Var}}_{c}[\hat{P}(t)]=\displaystyle \frac{1}{8\eta J(1-{{\rm{e}}}^{-\tfrac{\kappa t}{2}})+2},\end{eqnarray} \tag{ 19 }$$

that shows how the conditional state of the atomic sample is deterministically driven by the dynamics into a spin-squeezed state.

Here we will present our main results, that is the derivation of ultimate quantum limits on noisy magnetometry via time-continuous measurements of the atomic sample. We will first evaluate the classical FI ${ \mathcal F }[{{\bf{y}}}_{t}]$ corresponding to the information obtainable from the photocurrent, and we will also show how the corresponding bound can be achieved via Bayesian estimation. We will then evaluate the second term appearing in the bound, corresponding to the information obtainable via a strong measurement on the conditional state of the atomic sample. This will allow us to discuss the ultimate limit on the estimation strategy via the effective QFI: we will focus on the scaling with the relevant parameters of the experiment, i.e. with the total spin number J and the monitoring time characterizing each experimental run t, and we will address the role of the detector efficiency η.

4.1. Analytical FI corresponding to the time-continuous photocurrent

As discussed before, the measured photocurrent ${{\bf{y}}}_{t}$ obtained via continuous homodyne detection can be used to extract information about the system and to estimate parameters which appear in the dynamics. The ultimate limit on the precision of this estimate is quantified by the FI ${ \mathcal F }[p({{\bf{y}}}_{t})]$. Given the Gaussian nature and the simple dynamics of the problem we can compute it analytically in closed form, by applying the results of [58]. As we describe in more detail in appendix B, one obtains the formula

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{t})]=2\eta \kappa {{J}{\rm{e}}}^{-\kappa t/2}\,{{\mathbb{E}}}_{p({{\bf{y}}}_{t})}[{({\partial }_{B}\langle \hat{P}(t){\rangle }_{c})}^{2}].\end{eqnarray} \tag{ 20 }$$

By considering (17) and remembering that ${\rm{d}}{w}_{t}={\rm{d}}{y}_{t}-\sqrt{2\eta \kappa {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}\langle \hat{P}(t){\rangle }_{c}{\rm{d}}t$, one obtains that the time evolution of the derivative of the conditional first moment $\langle \hat{P}(t){\rangle }_{c}$ w.r.t. to the parameter B, can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}({\partial }_{B}\langle \hat{P}(t){\rangle }_{c})}{{\rm{d}}{t}}=-\gamma \sqrt{{{J}{\rm{e}}}^{-\kappa t/2}}-4{\mathrm{Var}}_{c}[\hat{P}(t)]\eta \kappa {{J}{\rm{e}}}^{-\kappa t/2}({\partial }_{B}\langle \hat{P}(t){\rangle }_{c}),\end{eqnarray} \tag{ 21 }$$

where ${\mathrm{Var}}_{c}[\hat{P}(t)]$ is obtained from equation (19). We thus observe that the evolution is deterministic and one can easily derive its analytical solution. By applying equation (20), as the average over the trajectories is not needed, we readily obtain the following analytical formula for the FI

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{t})]=\displaystyle \frac{64{\gamma }^{2}\eta {J}^{2}{{\rm{e}}}^{-\kappa t}{({{\rm{e}}}^{\tfrac{\kappa t}{4}}-1)}^{3}}{9{\kappa }^{2}[(4\eta J+1){{\rm{e}}}^{\tfrac{\kappa t}{2}}-4\eta J]}\cdot [-4\eta J-12\eta {{J}{\rm{e}}}^{\tfrac{\kappa t}{4}}+3(4\eta J+3){{\rm{e}}}^{\tfrac{\kappa t}{2}}+(4\eta J+3){{\rm{e}}}^{\tfrac{3\kappa t}{4}}].\end{eqnarray} \tag{ 22 }$$

As intuitively expected, this is a monotonically increasing function of t, since the partial derivative is always positive. To get some insight into this expression we first report the leading term for $t\to 0$

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{t})]\approx \displaystyle \frac{4}{3}{J}^{2}{\gamma }^{2}\kappa {t}^{3},\end{eqnarray} \tag{ 23 }$$

where we explicitly see both Heisenberg scaling J² and a monitoring-enhanced time scaling t³. We can get further intuition about this expression by expanding it around $J=\infty$, the limit in which the Gaussian approximation becomes exact. The leading order in this other expansion is quadratic in J, thus showing again Heisenberg scaling, irregardless of t:

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{t})]\approx \displaystyle \frac{64{\gamma }^{2}\eta {J}^{2}{{\rm{e}}}^{-\kappa t}{({{\rm{e}}}^{\tfrac{\kappa t}{4}}-1)}^{3}(4{{\rm{e}}}^{\tfrac{\kappa t}{4}}+{{\rm{e}}}^{\tfrac{\kappa t}{2}}+1)}{9{\kappa }^{2}({{\rm{e}}}^{\tfrac{\kappa t}{4}}+1)};\end{eqnarray} \tag{ 24 }$$

this last approximations actually reproduces the behavior of the function quite well in the range of parameters we will consider in the following.

We now want to show that one can achieve this classical Cramér–Rao bound from the time-continuous measurement outcomes obtained via an appropriate estimator. In figure 2 we indeed show the posterior distribution as a function of time for a single experimental run, obtained after a Bayesian analysis (see appendix C for details). We observe how the distribution gets narrower in time around the true value and we also explicitly show that its standard deviation ${\sigma }_{\mathrm{est}}$ converges to the one predicted by the Cramér–Rao bound ${\sigma }_{\mathrm{CR}}(t)={ \mathcal F }{[p({{\bf{y}}}_{t})]}^{-1/2}$. In the initial part of the dynamics the values of ${\sigma }_{\mathrm{est}}$ are smaller than the corresponding ${\sigma }_{\mathrm{CR}}$: this is due to the choice of the prior distribution, being narrower than the likelihood and thus implying some initial knowledge on the parameter which is larger than the one obtainable for small monitoring time.

Zoom In Zoom Out Reset image size

4.2. Quantum Cramér–Rao bound for noisy magnetometry via time-continuous measurements

In order to evaluate the quantum Cramér–Rao bound in equation (5) we now need to consider the second term ${{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal Q }[{\varrho }^{(c)}]]$, corresponding to the information obtainable via strong quantum measurement on the conditional state of the system. The conditional state ${\varrho }^{(c)}$ is Gaussian and has a dependence on the parameter B only in the first moments. Therefore the corresponding QFI can be evaluated as prescribed in [59] (see appendix B for more details) obtaining,

$$\begin{eqnarray}&&{ \mathcal Q }[{\varrho }^{(c)}]=\displaystyle \frac{{({\partial }_{B}\langle \hat{P}(t){\rangle }_{c})}^{2}}{{\mathrm{Var}}_{c}[\hat{P}(t)]}.\end{eqnarray} \tag{ 25 }$$

Since, as we proved before, the evolution of both ${\partial }_{B}\langle \hat{P}(t){\rangle }_{c}$ and ${\mathrm{Var}}_{c}[\hat{P}(t)]$ is deterministic, the average over all possible trajectories is also in this case trivial and we have ${{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal Q }[{\varrho }^{(c)}]]={ \mathcal Q }[{\varrho }^{(c)}]$. By exploiting the analytical solution for both quantities, the QFI reads

$$\begin{eqnarray}&&{ \mathcal Q }[{\varrho }^{(c)}]=\displaystyle \frac{32{\gamma }^{2}J{(12\eta J-4\eta {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}-(8\eta J+3){{\rm{e}}}^{\tfrac{\kappa t}{4}}+3)}^{2}}{9{\kappa }^{2}[(4\eta J+1){{\rm{e}}}^{\tfrac{\kappa t}{2}}-4\eta J]}.\end{eqnarray} \tag{ 26 }$$

As expected, for no monitoring of the environment ($\eta =0$), one obtains that ${ \mathcal Q }[{\varrho }^{(c)}]\sim J$, i.e. corresponding to the SQL scaling. This function is also monotonically increasing with t and we can expand it around $J=\infty$ to study the leading term, which shows again a quadratic scaling in J

$$\begin{eqnarray}&&{ \mathcal Q }[{\varrho }^{(c)}]\approx \displaystyle \frac{128{\gamma }^{2}\eta {J}^{2}{{\rm{e}}}^{-\kappa t}{(-3{{\rm{e}}}^{\tfrac{\kappa t}{2}}+2{{\rm{e}}}^{\tfrac{3\kappa t}{4}}+1)}^{2}}{9{\kappa }^{2}({{\rm{e}}}^{\tfrac{\kappa t}{2}}-1)}.\end{eqnarray} \tag{ 27 }$$

We also remark that the QFI is equal to the classical FI for a measurement of the quadrature $\hat{P}$, thus showing that a strong measurement of the operator ${\hat{J}}_{z}$ on the conditional state of the atomic sample is the optimal measurement saturating the corresponding quantum Cramér–Rao bound.

By combining equations (22) and (26), we can now define the effective QFI

$$\begin{eqnarray}&&\widetilde{{ \mathcal Q }}={ \mathcal F }[p({{\bf{y}}}_{t})]+{{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal Q }[{\varrho }^{(c)}]]={ \mathcal F }[p({{\bf{y}}}_{t})]+{ \mathcal Q }[{\varrho }^{(c)}],\end{eqnarray} \tag{ 28 }$$

which represent the inverse of the best achievable variance according to the quantum Cramér–Rao bound (5). The resulting expression can be simplified to get the following simple analytical formula

$$\begin{eqnarray}&&\tilde{{ \mathcal Q }}={K}_{1}J+\eta {K}_{2}{J}^{2},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{K}_{1}=32\displaystyle \frac{{\gamma }^{2}}{{\kappa }^{2}}{(1-{{\rm{e}}}^{-\kappa t/4})}^{2},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{K}_{2}=64\displaystyle \frac{{\gamma }^{2}}{{\kappa }^{2}}\left(1-\displaystyle \frac{8}{3}{{\rm{e}}}^{-\kappa t/4}+2{{\rm{e}}}^{-\kappa t/2}-\displaystyle \frac{1}{3}{{\rm{e}}}^{-\kappa t}\right).\end{eqnarray} \tag{ 31 }$$

We start by studying how this quantity scales with the total spin: in figure 3 we plot $\widetilde{{ \mathcal Q }}$ as a function of J in the appropriate regions of parameters. We remark that the plots will be presented by using $1/\kappa$ as a time unit so that the strength of the interaction becomes $\gamma /\kappa$ and is always fixed to $1\,{{\rm{G}}}^{-1}$ in the following. We observe that, within the validity of our approximation ($\kappa t\lesssim 1$), it is possible to obtain the Heisenberg-like scaling J² for the effective QFI. There is a transition between SQL-like scaling and Heisenberg scaling depending on the relationship between J and $\kappa t$ showing how the quantum enhancement is observed for $J\gg 1/\kappa t$.

Zoom In Zoom Out Reset image size

The same conclusions are drawn if we look at the behavior of $\widetilde{{ \mathcal Q }}$ as a function of the interrogation time t, plotted in figure 4: a transition from a t²-scaling to a monitoring-enhanced t³-scaling is observed for $J\gg 1/\kappa t$. We remark here that the typical scaling obtained in quantum metrology for unitary parameters is of order t². The observed t³-scaling is due to the continuos monitoring of the system. A similar scaling of the FI would be in fact obtained for an equivalent classical estimation problem, where a continuously monitored classical system is estimated via a the Kalman filter [58]. Notice that there are also few recent examples in the literature where a t⁴-scaling can be observed. This is obtained in noiseless quantum metrology problems with time-dependent Hamiltonian and by exploiting open-loop control [60–63]. In particular in [63], it was also shown that a t³-scaling can be achieved without additional control, but by performing repeated (stroboscopic) measurement on the system, analogously to our strategy.

Zoom In Zoom Out Reset image size

The previous results were both shown by considering perfect monitoring of the environment, i.e. for detectors with unit efficiency η. In figure 5 we plot the behaviors of $\widetilde{{ \mathcal Q }}$ as a function of J and t, varying the detector efficiency η; we observe how the quantum enhancements can be obtained for all non-zero values of η. The effect of having a non-unit monitoring efficiency is simply to imply larger values of J to witness the transition between SQL to Heisenberg-scaling, as one can also understand by looking at the role played by η and J in equation (29).

Zoom In Zoom Out Reset image size

We remind that if we consider only the classical FI ${ \mathcal F }[p({{\bf{y}}}_{t})]$, the Heisenberg scaling in terms of J and t³-scaling are always obtained for $\kappa t\lesssim 1$ and for every η, as shown by the expansion (23). However, if the contribution of this term, as well as the contribution of conditioning to the QFI, are too small then the QFI of the unconditional state, i.e. (26) with $\eta =0$, dominates (the term $\eta {K}_{2}{J}^{2}$ in (29) is negligible) and we observe SQL scaling for $\tilde{{ \mathcal Q }}$. We finally mention that the in the regimes where we observe Heisenberg scaling of $\tilde{{ \mathcal Q }}$, the classical FI ${ \mathcal F }[p({{\bf{y}}}_{t})]$ amounts to a relevant part of the total, namely around 25% .

4.3. Optimality of time-continuous measurement strategy for noisy quantum magnetometry

As explained before, the ultimate limit for quantum magnetometry, in the presence of Markovian transversal noise as the one described by the master equation (10), is given by the QFI ${\overline{{ \mathcal Q }}}_{{ \mathcal L }}$ in equation (8). The generalized master equation (9) in this case (considering the large-spin approximation) reads

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{\varrho }}{{\rm{d}}{t}}=-{\rm{i}}\gamma \sqrt{{\bar{J}}_{t}}({B}_{1}\hat{X}\bar{\varrho }-{B}_{2}\bar{\varrho }\hat{X})+\kappa {\bar{J}}_{t}{ \mathcal D }[\hat{P}]\bar{\varrho }.\end{eqnarray} \tag{ 32 }$$

In appendix D we show how this equation can be solved in a phase space picture, since the equation contains at most quadratic terms in $\hat{X}$ and $\hat{P}$ and thus preserves the Gaussian character of the operator $\bar{\varrho }$.

The final result is

$$\begin{eqnarray}&&{\overline{{ \mathcal Q }}}_{{{ \mathcal L }}_{{tn}}}=\tilde{{ \mathcal Q }}(\eta =1)={K}_{1}J+{K}_{2}{J}^{2},\end{eqnarray} \tag{ 33 }$$

i.e., we exactly obtain the effective QFI $\tilde{{ \mathcal Q }}$ defined in equation (28) in the limit of unit efficiency $\eta =1$. This result remarkably proves that our strategy, not only allows to obtain the Heisenberg limit, but also corresponds to the optimal one, given a collective transversal noise master equation (10) and in the presence of perfectly efficient detectors. Indeed, any other more experimentally complicated strategy, based on entangled and non-local in time measurements of the output modes and the system, would not give better results in the estimation of the magnetic field B.

We have addressed in detail estimation strategies for a static and constant magnetic field acting on an atomic ensemble of two-level atoms also subject to transverse noise. In particular, we have evaluated the ultimate quantum limits to precision for strategies based on time-continuous monitoring of the light coupled the atomic ensemble.

After deriving the appropriate quantum Cramér–Rao bound, we have calculated the corresponding effective QFI in the limit of large spin, posing the ultimate limit on the mean-square error of any unbiased estimator. Our results conclusively show that both Heisenberg J²-scaling in terms of spin, and a monitoring-enhanced t³-scaling in terms of the interrogation time, are obtained for $J\gg 1/\kappa t$, confirming what was discussed in [43, 46]. We have remarkably demonstrated that these quantum enhancements are also obtained for not unit monitoring efficiency, i.e. even if one cannot measure all the environmental modes or for not perfectly efficient detectors. Finally we have analytically proven the optimality of our strategy, i.e. that given the master equation describing the unconditional dynamics of the system and ideal detectors, no other measurement strategy would give better results in estimating the magnetic field.

We remark that Heisenberg scaling, or at least a super-classical scaling, can be obtained in the presence of collective or individual (independent) transversal noise, by preparing a highly entangled or spin-squeezed state at the beginning of the dynamics and, for individual noise, by optimizing on the interrogation time t [17–19]. In this respect, the advantage of our protocol lies in the fact that it achieves the Heisenberg scaling even for an initial classical spin-coherent state, exploiting the dynamical spin squeezing that is generated by the weak measurement.

In conclusion, we have shown that time-continuous measurements represent a resource for noisy quantum magnetometry [43, 46, 47]. Indeed, the information leaking into the environment, here represented by light modes coupled to the atomic sample, obtained via homodyne detection, and the corresponding measurement back-action on the atomic sample, may be efficiently (and optimally) exploited in order to obtain the promised quantum enhanced estimation precision.

MGG would like to thank A Doherty and A Serafini for discussions and acknowledges support from Marie Skłodowska-Curie Action H2020-MSCA-IF-2015 (project ConAQuMe, grant nr. 701154). This work has been supported by EU through the collaborative H2020 project QuProCS (Grant Agreement 641277) and by UniMI through the H2020 Transition Grant.

Here we will show how to derive the quantum Cramér–Rao bound for time-continuous homodyne monitoring reported in equation (5).

We start by considering a (classical) estimation problem of a parameter θ described by a conditional probability $p(z,{{\bf{y}}}_{T}| \theta )$, where the vector ${{\bf{y}}}_{T}={({y}_{1},{y}_{2},\ldots ,{y}_{T})}^{{\mathsf{T}}}$ contains the outcomes of sequential measurements performed up to time T, while z corresponds to a final measurement performed on the state of the system that has been conditioned on the previous measurement results ${{\bf{y}}}_{T}$. The corresponding classical FI can be evaluated as

$$\begin{eqnarray}{ \mathcal F }[p(z,{{\bf{y}}}_{T}| \theta )] & = & \displaystyle \int {\rm{d}}{\bf{y}}\,{\rm{d}}{z}\,p(z,{{\bf{y}}}_{T}| \theta ){({\partial }_{\theta }\mathrm{log}p(z,{{\bf{y}}}_{T}| \theta ))}^{2}\\ & = & \displaystyle \int {\rm{d}}{\bf{y}}\,{\rm{d}}{z}\,p(z| {{\bf{y}}}_{T},\theta )p({{\bf{y}}}_{T}| \theta )[{({\partial }_{\theta }\mathrm{log}p(z| {{\bf{y}}}_{T},\theta ))}^{2}\\ & & +2({\partial }_{\theta }\mathrm{log}p(z| {{\bf{y}}}_{T},\theta ))({\partial }_{\theta }\mathrm{log}p({{\bf{y}}}_{T}| \theta ))\\ & & +{({\partial }_{\theta }\mathrm{log}p({{\bf{y}}}_{T}| \theta ))}^{2}]\,,\end{eqnarray} \tag{ A.1 }$$

where the second expression has been obtained by means of the Bayes rule

$$\begin{eqnarray*}&&p(z,{{\bf{y}}}_{T}| \theta )=p(z| {{\bf{y}}}_{T},\theta )p({{\bf{y}}}_{T}| \theta ).\end{eqnarray*}$$

In the following, we will omit the dependence on the parameter θ and we will denote by ${{\mathbb{E}}}_{p(x)}[\cdot ]$ the average over a probability distribution p(x). By considering each term inside the integral separately one obtains

$$\begin{eqnarray}&&{{\mathbb{E}}}_{p(z,{{\bf{y}}}_{T})}[{({\partial }_{\theta }\mathrm{log}p(z| {{\bf{y}}}_{T}))}^{2}]={{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal F }[p(z| {{\bf{y}}}_{T})]],\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&2\,{{\mathbb{E}}}_{p(z,{{\bf{y}}}_{T})}[{\partial }_{\theta }\mathrm{log}p(z| {{\bf{y}}}_{T})\,{\partial }_{\theta }\mathrm{log}p({{\bf{y}}}_{T})]=2\,\int {\rm{d}}{\bf{y}}({\partial }_{\theta }p({{\bf{y}}}_{T}))\int {\rm{d}}{z}({\partial }_{\theta }p(z| {{\bf{y}}}_{T}))=0,\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{{\mathbb{E}}}_{p(z,{{\bf{y}}}_{T})}[{({\partial }_{\theta }\mathrm{log}p({{\bf{y}}}_{T}))}^{2}]={ \mathcal F }[p({{\bf{y}}}_{T})],\end{eqnarray} \tag{ A.4 }$$

where we have used the property $\int {\rm{d}}{z}({\partial }_{\theta }p(z| {{\bf{y}}}_{T}))={\partial }_{\theta }\int {\rm{d}}{z}\,p(z| {{\bf{y}}}_{T})={\partial }_{\theta }(1)=0$. As a consequence, any unbiased estimator $\hat{\theta }$ based on M experiments, i.e. obtained collecting M series of measurement outcomes $({{\bf{y}}}_{T},z)$, satisfies the generalized Cramér–Rao bound

$$\begin{eqnarray}&&{\mathrm{Var}}_{\hat{\theta }}(\theta )\geqslant \displaystyle \frac{1}{M({ \mathcal F }[p({{\bf{y}}}_{T})]+{{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal F }[p(z| {{\bf{y}}}_{T})]])}\,,\end{eqnarray} \tag{ A.5 }$$

where the first term ${ \mathcal F }[p({{\bf{y}}}_{T})]$ is the FI corresponding to the sequential measurements with outcomes ${{\bf{y}}}_{T}$, while the second term is the average of the FI ${ \mathcal F }[p(z| {{\bf{y}}}_{T})]$, corresponding to the final measurement over all the possible trajectories conditioned on the previous measurement results ${{\bf{y}}}_{T}$. The bound in equation (A.5) bears some formal similarity to the Van Tree's inequality [64], which however applies in a quite different situation, i.e. the case where the parameter to be estimated θ is a random variable distributed according to a given probability distribution $p(\theta )$.

The estimation strategy here described is of particular interest when we deal with quantum systems, given the back-action of quantum measurement on the state of the system itself. We can in fact associate each measurement outcome ${y}_{k}$ to a Kraus operator ${M}_{{y}_{k}}$ such that the conditional quantum state, for the system initially prepared in a state ϱ₀ and after obtaining the stream of outcomes ${{\bf{y}}}_{T}$, reads

$$\begin{eqnarray}&&{\varrho }_{{{\bf{y}}}_{T}}^{(c)}=\displaystyle \frac{{\tilde{M}}_{{{\bf{y}}}_{T}}{\varrho }_{0}{\tilde{M}}_{{{\bf{y}}}_{T}}^{\dagger }}{\mathrm{Tr}[{\tilde{M}}_{{{\bf{y}}}_{T}}{\varrho }_{0}{\tilde{M}}_{{{\bf{y}}}_{T}}^{\dagger }]},\end{eqnarray} \tag{ A.6 }$$

where ${\tilde{M}}_{{{\bf{y}}}_{T}}={M}_{{y}_{T}}...{M}_{{y}_{2}}{M}_{{y}_{1}}$ and the probability of obtaining the outcomes ${{\bf{y}}}_{T}$ reads $p({{\bf{y}}}_{T}| \theta )=\mathrm{Tr}[{\tilde{M}}_{{{\bf{y}}}_{T}}{\varrho }_{0}{\tilde{M}}_{{{\bf{y}}}_{T}}^{\dagger }]$ ³ . One can then also perform a strong (destructive) measurement described by POVM operators $\{{{\rm{\Pi }}}_{z}\}$ on the conditional state, and the whole measurement strategy is described by the conditional probabilities

$$\begin{eqnarray}p(z| {{\bf{y}}}_{T},\theta ) & = & \mathrm{Tr}[{\varrho }_{{{\bf{y}}}_{T}}^{(c)}{{\rm{\Pi }}}_{z}],\\ p(z,{{\bf{y}}}_{T}| \theta ) & = & p(z| {{\bf{y}}}_{T},\theta )\,p({{\bf{y}}}_{T}| \theta )\\ & = & \mathrm{Tr}[{\tilde{M}}_{{{\bf{y}}}_{T}}{\varrho }_{0}{\tilde{M}}_{{{\bf{y}}}_{T}}^{\dagger }{{\rm{\Pi }}}_{z}].\end{eqnarray} \tag{ A.7 }$$

Typically the parameter to be estimated θ enters in the the dynamics described by the Kraus operators ${M}_{{y}_{k}}$. For this reason we will start by considering these operators fixed, while we suppose we can optimize over the final measurement $\{{{\rm{\Pi }}}_{z}\}$. We can then apply the quantum Cramér–Rao bound for the conditional states ${\varrho }_{{{\bf{y}}}_{T}}^{(c)}$, stating that ${ \mathcal F }[p(z| {{\bf{y}}}_{T})]\leqslant { \mathcal Q }[{\varrho }_{{{\bf{y}}}_{T}}^{(c)}]$. One then obtains a more fundamental quantum Cramér–Rao bound for our estimation strategy

$$\begin{eqnarray}&&{\mathrm{Var}}_{\hat{\theta }}(\theta )\geqslant \displaystyle \frac{1}{M({ \mathcal F }[p({{\bf{y}}}_{T})]+{{\mathbb{E}}}_{p({{\bf{y}}}_{T})}[{ \mathcal Q }[{\varrho }_{{{\bf{y}}}_{T}}^{(c)}]])}.\end{eqnarray} \tag{ A.8 }$$

Clearly this bound can be readily applied to the time-continuous case discussed in the main text, where the vector of outcomes ${{\bf{y}}}_{T}$ corresponds to a measured homodyne photocurrent, and where the conditional state ${\varrho }_{{{\bf{y}}}_{T}}^{(c)}$ can be obtained via a stochastic master equation as the one in equation (3).

We should also remark that a bound of this kind has already been considered in [39], in a similar physical situation where n probes, that may be prepared in a quantum correlated initial state, are coupled to n independent environments and one performs sequentially n measurement on the respective environments and a final measurement on the conditional state of the probes.

Here we will provide the formulas for the dynamics of the atomic ensemble described by the stochastic master equation (14). As we mentioned in the text, the whole dynamics preserves the Gaussian character of the quantum state and thus can be fully described in terms of the first moments vector $\langle \hat{{\bf{r}}}{\rangle }_{c}$ and of the covariance matrix ${\boldsymbol{\sigma }}$ of the quantum state ${\varrho }^{(c)}$. These are defined in components as $\langle {\hat{r}}_{j}{\rangle }_{c}=\mathrm{Tr}[{\hat{r}}_{j}{\varrho }^{(c)}]$ and ${\sigma }_{{jk}}=\mathrm{Tr}[\{{\hat{r}}_{j}-\langle {\hat{r}}_{j}{\rangle }_{c},{\hat{r}}_{k}-\langle {\hat{r}}_{k}{\rangle }_{c}\}{\varrho }^{(c)}]$ for the operator vector $\hat{{\bf{r}}}={(\hat{X},\hat{P})}^{{\mathsf{T}}}$. In formulae one obtains [56, 57]:

$$\begin{eqnarray}&&{\rm{d}}\langle \hat{{\bf{r}}}{\rangle }_{c}={\bf{u}}\,{\rm{d}}{t}+\displaystyle \frac{{\boldsymbol{\sigma }}M\,{\rm{d}}{\bf{w}}}{\sqrt{2}},\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{\sigma }}}{{\rm{d}}{t}}=D-{\boldsymbol{\sigma }}{{MM}}^{{\mathsf{T}}}{\boldsymbol{\sigma }},\end{eqnarray} \tag{ B.2 }$$

where

$$\begin{eqnarray}D=\left(\begin{array}{cc}2\kappa {{J}{\rm{e}}}^{-\kappa t/2} & 0\\ 0 & 0\end{array}\right),\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}M=\left(\begin{array}{cc}0 & 0\\ \sqrt{2\eta \kappa {{J}{\rm{e}}}^{-\kappa t/2}} & 0\end{array}\right),\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{\bf{u}}={(0,-\gamma B\sqrt{{{J}{\rm{e}}}^{-\kappa t/2}})}^{{\mathsf{T}}},\end{eqnarray} \tag{ B.5 }$$

and ${\rm{d}}{\bf{w}}$ is a vector of Wiener increments such that ${{\rm{d}}{w}}_{j}\,{{\rm{d}}{w}}_{k}={\delta }_{{jk}}{\rm{d}}{t}$, related to the photocurrent via the equation

$$\begin{eqnarray}&&{\rm{d}}{{\bf{y}}}_{t}=-{M}^{{\mathsf{T}}}\langle \hat{{\bf{r}}}{\rangle }_{c}\,{\rm{d}}{t}+\displaystyle \frac{{\rm{d}}{\bf{w}}}{\sqrt{2}}.\end{eqnarray} \tag{ B.6 }$$

The equations (17), (18) and (21) can be obtained from the ones above, remembering that for our definitions ${\sigma }_{22}=2{\mathrm{Var}}_{c}[\hat{P}(t)]$.

The method to calculate the FI corresponding to the time-continuous measurement in the case of linear Gaussian system has been described in [58]. One has to evaluate the formula

$$\begin{eqnarray}&&{ \mathcal F }[p({{\bf{y}}}_{t})]={{\mathbb{E}}}_{p({{\bf{y}}}_{t})}[2({\partial }_{B}\langle \hat{{\bf{r}}}{\rangle }_{c}^{{\mathsf{T}}}){{MM}}^{{\mathsf{T}}}({\partial }_{B}\langle \hat{{\bf{r}}}{\rangle }_{c})],\end{eqnarray} \tag{ B.7 }$$

that, by plugging in the matrices describing our problem, is easily simplified to equation (20).

As the conditional state is Gaussian, also the calculation of the corresponding QFI can be easily obtained, in this case by applying the results presented in [59]. Moreover, as only the first moments of the state depend on the parameter B, the calculation is further simplified and one has

$$\begin{eqnarray}&&{ \mathcal Q }[{\varrho }^{(c)}]=2\,({\partial }_{B}\langle \hat{{\bf{r}}}{\rangle }_{c}^{{\mathsf{T}}})\,{{\boldsymbol{\sigma }}}^{-1}\,({\partial }_{B}\langle \hat{{\bf{r}}}{\rangle }_{c}).\end{eqnarray} \tag{ B.8 }$$

By noticing that the only non-zero entry of the vector ${\partial }_{B}\langle \hat{{\bf{r}}}{\rangle }_{c}$ is the one corresponding to $\langle \hat{P}(t){\rangle }_{c}$, one easily obtain equation (25).

Bayesian analysis has proven to be an efficient tool for estimation in continuously monitored quantum systems [35, 38, 42, 48]. The goal is to reconstruct the posterior distribution of B given the observed current ${{\bf{y}}}_{t}$, by Bayes rule:

$$\begin{eqnarray}&&p(B| {{\bf{y}}}_{t})=\displaystyle \frac{L(B| {{\bf{y}}}_{t})p(B)}{p({{\bf{y}}}_{t})},\end{eqnarray} \tag{ C.1 }$$

where $p(B)$ is the prior distribution, $L(B| {{\bf{y}}}_{t})\equiv p({{\bf{y}}}_{t}| B)$ is the likelihood and $p({{\bf{y}}}_{t})$ serves as a normalization factor. The Bayesian estimator is the mean of the posterior distribution $\hat{B}({{\bf{y}}}_{t})={{\mathbb{E}}}_{p(B| {{\bf{y}}}_{t})}[B]$ and it is proven that the corresponding variance ${\mathrm{Var}}_{\hat{B}}{(B)={{\mathbb{E}}}_{p(B| {{\bf{y}}}_{t})}[{B}^{2}]-({{\mathbb{E}}}_{p(B| {{\bf{y}}}_{t})}[B])}^{2}$ is asymptotically optimal, i.e. tends to saturate the Cramér–Rao bound when the length of the vector ${{\bf{y}}}_{T}$ is large.

The simulated experimental run is obtained by numerically integrating the stochastic differential equation (17) with the Euler–Maruyama method for the 'true' value of the parameter ${B}_{\mathrm{true}}$. Time is discretized with steps of length ${\rm{\Delta }}t$, i.e. to get from time 0 to time T we perform ${n}_{T}=T/{\rm{\Delta }}t$ steps. Experimental data is represented by the observed measurement current ${{\bf{y}}}_{T}={({\rm{\Delta }}{y}_{{t}_{1}},\ldots ,{\rm{\Delta }}{y}_{{t}_{{n}_{T}}})}^{{\mathsf{T}}}$, which corresponds to an n_T-dimensional vector. The outcome at every time step ${\rm{\Delta }}{y}_{{t}_{i}}$ is sampled from a Gaussian distribution with variance ${\rm{\Delta }}t$ and mean $\overline{{\rm{\Delta }}{y}_{{t}_{i}}}(B)=\sqrt{2\eta \kappa {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}\langle \hat{P}({t}_{i}){\rangle }_{c}{\rm{\Delta }}t$. Notice that $\overline{{\rm{\Delta }}{y}_{{t}_{i}}}(B)$ depends explicitly on the parameter B via the quantum expectation value $\langle \hat{P}({t}_{i}){\rangle }_{c}$ on the conditional state.

Since we are estimating only one parameter the posterior can be obtained on a grid on the parameter space, while for more complicated problems Markov chain Monte Carlo methods might be needed to sample from the posterior [35]. In practical terms we need to solve equations (17) and (18) for every value of the parameter B on the grid, assuming to perfectly know all the other parameters; then we need to calculate the likelihood for each value via

$$\begin{eqnarray}&&L(B| {{\bf{y}}}_{T})\propto \displaystyle \prod _{i=0}^{{n}_{T}}\exp \left[-\displaystyle \frac{{({\rm{\Delta }}{y}_{{t}_{i}}-\overline{{\rm{\Delta }}{y}_{{t}_{i}}}(B))}^{2}}{2\,{\rm{\Delta }}t}\right],\end{eqnarray} \tag{ C.2 }$$

by considering the outcomes as independent random variables, i.e. multiplying the corresponding probabilities. We then apply Bayes rule, equation (C.1), assuming a flat prior distribution p(B) on a finite interval. The same analysis is trivially applied to more than one experiment by simply multiplying the likelihood obtained for every different observed measurement current.

Here we explicitly show how to solve equation (32). The characteristic function for a generic operator $\hat{O}$ is defined as

$$\begin{eqnarray}&&\chi [\hat{O}]({\bf{s}})=\mathrm{Tr}[{\hat{D}}_{-{\bf{s}}}\hat{O}],\end{eqnarray} \tag{ D.1 }$$

where the displacement operator is defined as

$$\begin{eqnarray}&&{\hat{D}}_{-{\bf{s}}}=\exp ({\rm{i}}{{\bf{s}}}^{\top }{\rm{\Omega }}\hat{{\bf{r}}}).\end{eqnarray} \tag{ D.2 }$$

In particular we will work in the phase space of a single mode system, so that ${\hat{{\bf{r}}}}^{\top }=(\hat{X},\hat{P})$ is the vector of quadrature operators and ${{\bf{s}}}^{\top }=(x,p)$ is the vector of phase space coordinates.

The action of operators in the Hilbert space corresponds to differential operators acting on the characteristic function via the following mapping [56, 65]

$$\begin{eqnarray}&&\hat{X}\rho \leftrightarrow \left(-{\rm{i}}{\partial }_{p}-\displaystyle \frac{x}{2}\right)\chi ({\bf{s}}),\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&\rho \hat{X}\leftrightarrow \left(-{\rm{i}}{\partial }_{p}+\displaystyle \frac{x}{2}\right)\chi ({\bf{s}}),\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}&&\hat{P}\rho \leftrightarrow \left({\rm{i}}{\partial }_{x}-\displaystyle \frac{p}{2}\right)\chi ({\bf{s}}),\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray}&&\rho \hat{P}\leftrightarrow \left({\rm{i}}{\partial }_{x}+\displaystyle \frac{p}{2}\right)\chi ({\bf{s}}).\end{eqnarray} \tag{ D.6 }$$

If we now define the characteristic function associated to the operator $\bar{\varrho }$ introduced in equation (9)

$$\begin{eqnarray}&&\bar{\chi }({\bf{s}},t)\equiv \chi [\bar{\varrho }]({\bf{s}}),\end{eqnarray} \tag{ D.7 }$$

the quantity of interest in order to compute the QFI is then $\mathrm{Tr}\bar{\varrho }=\bar{\chi }(0,t)$, as evident from the definition (D.1).

By applying the phase space mapping, from the generalized master equation (32) we get to the following partial differential equation for the characteristic function

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{\chi }({\bf{s}},t)}{{\rm{d}}{t}}=\left[{\rm{i}}\gamma \sqrt{{\bar{J}}_{t}}\displaystyle \frac{{B}_{1}+{B}_{2}}{2}x-\displaystyle \frac{\kappa {\bar{J}}_{t}}{2}{p}^{2}-\gamma \sqrt{{\bar{J}}_{t}}({B}_{1}-{B}_{2}){\partial }_{p}\right]\bar{\chi }({\bf{s}},t).\end{eqnarray} \tag{ D.8 }$$

This equation can be solved by performing a Gaussian ansatz, similarly to [66], i.e. assuming that at every time the characteristic function can be written in the following form

$$\begin{eqnarray}&&\bar{\chi }({\bf{s}},t)=C(t)\exp \left[-\displaystyle \frac{1}{4}{{\bf{s}}}^{\top }{{\rm{\Omega }}}^{\top }{\boldsymbol{\sigma }}(t){\rm{\Omega }}\,{\bf{s}}\,+{\rm{i}}{{\bf{s}}}^{\top }{{\rm{\Omega }}}^{\top }{{\bf{s}}}_{m}(t)\right].\end{eqnarray} \tag{ D.9 }$$

The dependence on time and on the parameters ${B}_{1/2}$ is completely contained in the covariance matrix ${\boldsymbol{\sigma }}(t)$, in the first moment vector ${{\bf{s}}}_{m}{(t)}^{\top }=({x}_{m}(t),{p}_{m}(t))$ and in the function $C(t)=\bar{\chi }(0,t)$, which is the final result we are seeking.

By plugging (D.9) into (D.8) and equating the coefficients for different powers of x and p, one obtains a system of differential equations; the relevant ones are the equations coming from the coefficients of order one, and multiplying p and p²:

$$\begin{eqnarray}&&{\dot{{\boldsymbol{\sigma }}}}_{\mathrm{1,1}}(t)=2\kappa {{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}},\end{eqnarray} \tag{ D.10 }$$

$$\begin{eqnarray}&&{\dot{x}}_{m}(t)=-{\rm{i}}\displaystyle \frac{\gamma }{2}\sqrt{{{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}({B}_{1}-{B}_{2}){{\boldsymbol{\sigma }}}_{\mathrm{1,1}}(t),\end{eqnarray} \tag{ D.11 }$$

$$\begin{eqnarray}&&\dot{C}(t)=-{\rm{i}}\gamma \sqrt{{{J}{\rm{e}}}^{-\tfrac{\kappa t}{2}}}({B}_{1}-{B}_{2}){x}_{m}(t)C(t).\end{eqnarray} \tag{ D.12 }$$

These equations are solved analytically with the initial conditions ${{\boldsymbol{\sigma }}}_{\mathrm{1,1}}(0)=1$, ${x}_{m}(0)=0$ and $C(0)=1$ (since for t = 0 the operator $\bar{\varrho }$ corresponds to the initial state of the system $| 0\rangle \langle 0|$ ), yielding

$$\begin{eqnarray}&&C(t)=\exp \left[-\displaystyle \frac{4{\gamma }^{2}}{3{\kappa }^{2}}J{({B}_{1}-{B}_{2})}^{2}{{\rm{e}}}^{-\kappa t}{({{\rm{e}}}^{\tfrac{\kappa t}{4}}-1)}^{2}\cdot (-4{{J}{\rm{e}}}^{\tfrac{\kappa t}{4}}+(6J+3){{\rm{e}}}^{\tfrac{\kappa t}{2}}-2J)\right].\end{eqnarray} \tag{ D.13 }$$

By plugging this term into equation (8), we finally obtain the ultimate QFI ${\overline{{ \mathcal Q }}}_{{{ \mathcal L }}_{{tn}}}$ reported in equation (33).