Ultimate limits for quantum magnetometry via time-continuous measurements

We address the estimation of the magnetic field B acting on an ensemble of atoms with total spin J subjected to collective transverse noise. By preparing an initial spin coherent state, for any measurement performed after the evolution, the mean-square error of the estimate is known to scale as $1/J$, i.e. no quantum enhancement is obtained. Here, we consider the possibility of continuously monitoring the atomic environment, and conclusively show that strategies based on time-continuous non-demolition measurements followed by a final strong measurement may achieve Heisenberg-limited scaling $1/{J^2}$ and also a monitoring-enhanced scaling in terms of the interrogation time. We also find that time-continuous schemes are robust against detection losses, as we prove that the quantum enhancement can be recovered also for finite measurement efficiency. Finally, we analytically prove the optimality of our strategy.

We address the estimation of the magnetic field B acting on an ensemble of atoms with total spin J subjected to collective transverse noise. By preparing an initial spin coherent state, for any measurement performed after the evolution, the mean-square error of the estimate is known to scale as 1/J, i.e. no quantum enhancement is obtained. Here, we consider the possibility of continuously monitoring the atomic environment, and conclusively show that strategies based on time-continuous non-demolition measurements followed by a final strong measurement may achieve Heisenberg-limited scaling 1/J 2 and also a monitoring-enhanced scaling in terms of the interrogation time. We also find that time-continuous schemes are robust against detection losses, as we prove that the quantum enhancement can be recovered also for finite measurement efficiency. Finally, we analytically prove the optimality of our strategy.

I. INTRODUCTION
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][4][5][6][7][8][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 nogo 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][16][17][18][19], or by exploiting techniques borrowed from the field of quantum errorcorrection [20][21][22].
Here we rigorously address the performance of these protocols, depicted in Fig. 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 Heisenberglike scaling. 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 Sec. II, 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 Sec. III, 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 Sec. IV, we present the main results: we first calculate the classical Fisher information 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. 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 V closes the paper with some concluding remarks.

II. QUANTUM CRAMÉR-RAO BOUNDS FOR TIME-CONTINUOUS HOMODYNE MONITORING
A classical estimation problem consists in inferring the value of a parameter θ from a number M of measurement outcomes χ = {x 1 , x 2 , . . . , x M } and their conditional distribution p(x|θ). We define an estimatorθ(χ) 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. dχ p(χ|λ)θ(χ) = θ, where p(χ|λ) = Π M j=1 p(x j |θ). The Cramér-Rao theorem states that the variance of any unbiased estimator is lower bounded as denotes the classical Fisher information (FI).
In the quantum realm, the conditional probability distribution reads p(x|θ) = Tr[ θ Π x ], where θ is the quantum state of the system labeled by the parameter θ, and Π x is a POVM operator describing the quantum measurement. One can prove that the FI corresponding to any POVM is upper bounded is the quantum Fisher information (QFI), and L θ is the so-called symmetric logarithmic derivative, which can be obtained by solving the equation 2∂ θ θ = L θ θ + θ L θ [49][50][51]. The QFI depends on the quantum state θ 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Ĥ θ characterised by the parameter we want to estimate, and coupled to a bosonic environment at zero temperature described by a train of input operatorsâ in (t), satisfying the commutation relation in (t) +ĉ †â in (t) (ĉ 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 where D[c] =ĉ ĉ † − (ĉ †ĉ + ĉ †ĉ )/2. If one performs a homodyne detection of a quadraturê x out (t) =â out (t) +â † out (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 (c) conditioned on the measurement results (we will omit the dependence of the measured photocurrent y t ), is described by the stochastic master equation [23] Here η denotes the efficiency of the detection, dw t is a stochastic Wiener increment (s.t. dw 2 t = dt), and H[c] (c) = c (c) + (c)ĉ † − Tr[ (c) (ĉ +ĉ † )] (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 → t + dt is given by the infinitesimal current With the help of such measurement strategies, one can estimate the value of the parameter θ both from the measured photocurrent y T = T 0 dy t , and from a final strong (destructive) measurement on the conditional state (c) . In this case, as we explicitly show in A (in general both for the classical and quantum case), the proper quantum Cramér-Rao bound reads where the first term at the denominator F[p(y T )] is the FI corresponding to the classical photocurrent y T , while the second term is the average of the QFI for the conditional state Q[ (c) ] over all the possible trajectories, i.e. on all the possible measurement outcomes for the photocurrent. The classical FI F[p(y T )] can be calculated as described in [35] by evaluating where the operator τ evolves according to the stochastic master equation The conditional states (c) at time T can be obtained by integrating (3), for a certain stream of outcomes y T ; then one can first calculate the corresponding quantum Fisher information Q[ (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 modeŝ a 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 where ψ(θ 1 )|ψ(θ 2 ) 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 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 Tr [ρ] = ψ(θ 1 )|ψ(θ 2 ) , which obeys the following generalized master equation As before, we already assumed that the dependence on the parameter lies only in the system HamiltonianĤ θ and that we have a single jump operatorĉ. We remark that the operatorρ is not a proper density operator representing a quantum state, except in the limit case θ 1 → θ 2 , where we recover the standard master equation (2).

III. QUANTUM MAGNETOMETRY: THE PHYSICAL SETTING
We address the estimation of the intensity of a static and constant magnetic field B acting on a ensemble of N twolevel atoms that are continuously monitored [43][44][45][46], as depicted in Fig. 1. The atomic ensemble can be described as a system with total spin J = N/2 with collective spin operators defined asĴ α = 1 2 N i=0 σ iα , where α = x, y, z and σ iα denotes the Pauli matrices acting on the i-th spin. The collective operators obey the same angular momentum commutation rules [Ĵ i ,Ĵ j ] = iε ijkĴk , where ε ijk is the Levi-Civita symbol. We remark that in the present manuscript we choose units such that = 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 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, 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, where |+ is the eigenstate of σ x with eigenvalue +1. We thus have that the spin component on the x direction attains the macroscopic value Ĵ x (0) = J. The unconditional dynamics of Ĵ x is obtained by applying the operatorĴ x to both sides of Eq. (11) and then taking the trace. The result is the following equation describing damped oscillations where we observe how the the dissipative and unitary parts of the dynamics are respectively shrinking the spin vector Ĵ and causing its Larmor precession around the y-axis. In the following we will assume to measure small magnetic fields, such that γBt 1 and we can approximate the solution of the previous equation as If the light modes are continuously monitored via homodyne measurements at the appropriate phase, one allows a continuous "weak" measurement ofĴ z ; the corresponding stochastic master equation (3) for finite monitoring efficiency η reads while the measurement result at time t corresponds to an in- is important to remark how the collective noise characterizing the master equation (11) 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][54][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 1. In this case, the dynamics may be effectively described with the Gaussian formalism as long as Ĵ x (t) ≈ J, i.e. for times t small enough to guarantee that κt 1. We define the effective quadrature operators of the atomic sample, satisfying the canonical commutation relation [X,P ] = i, as [46,47] whereJ t ≡ | Ĵ x (t) | (notice that in the limit of large spin J we can safely consider the unconditional average value Ĵ x (t) , as the stochastic correction obtained via (15) would be negligible). In the Gaussian description the initial state |ψ(0) corresponds to the vacuum state (X + iP )|0 = |0 , which is Gaussian. As the stochastic master equation (15) becomes quadratic in the canonical operators (and thus preserves the Gaussian character of states) the whole dynamics can be equivalently rewritten in terms of first and second moments only [56,57] (see 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 behaviour of the mean and the variance of the atomic momentum quadratureP calculated on the conditional state (c) , which follows the equations The differential equation for the conditional second moment is deterministic and can be solved analytically. For an initial vacuum state, i.e. with Var[P (0)] = 1 2 , we obtain the following solution that shows how the conditional state of the atomic sample is deterministically driven by the dynamics into a spin-squeezed state.

IV. RESULTS
Here we will present our main results, that is the derivation of ultimate quantum limits on noisy magnetometry via timecontinuous measurements of the atomic sample. We will first evaluate the classical Fisher information F[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 quantum Fisher information: 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 η.

A. Analytical FI corresponding to the time-continuous photocurrent
As discussed before, the measured photocurrent 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 F[p(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 B, one obtains the formula By considering (18) and remembering that dw t = dy t − 2ηκJe − κt 2 P (t) c dt, one obtains that the time evolution of the derivative of the conditional first moment P (t) c w.r.t. to the parameter B, can be written as where Var c [P (t)] is obtained from Eq. (20). We thus observe that the evolution is deterministic and one can easily derive its analytical solution. By applying Eq. (21), as the average over the trajectories is not needed, we readily obtain the following analytical formula for the FI 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 → 0 where we explicitly see both Heisenberg scaling J 2 and a monitoring-enhanced time scaling t 3 . We can get further intuition about this expression by expanding it around J = ∞, 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: (25) 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 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 σ est converges to the one predicted by the Cramér-Rao bound σ CR (t) = F[p(y t )] −1/2 . In the initial part of the dynamics the values of σ est are smaller than the corresponding σ 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. prescribed in [59] (see B for more details) obtaining, Since, as we proved before, the evolution of both ∂ B P (t) c and We also remark that the QFI is equal to the classical FI for a measurement of the quadratureP , thus showing that a strong measurement of the operatorĴ z on the conditional state of the atomic sample is the optimal measurement saturating the corresponding quantum Cramér-Rao bound.
By combining Eqs. (23) and (27), we can now define the effective quantum Fisher information 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 formulaQ where We start by studying how this quantity scales with the total spin: in Fig. 3 we plot Q as a function of J in the appropriate regions of parameters. We remark that the plots will be presented by using 1/κ as a time unit so that the strength of the interaction becomes γ/κ and is always fixed to 1 G −1 in the following. We observe that, within the validity of our approximation (κt 1), it is possible to obtain the Heisenberglike scaling J 2 for the effective QFI. There is a transition between SQL-like scaling and Heisenberg scaling depending on the relationship between J and κt showing how the quantum enhancement is observed for J 1/κt. The same conclusions are drawn if we look at the behaviour of Q as a function of the interrogation time t, plotted in Fig. 4: a transition from a t 2 -scaling to a monitoring-enhanced t 3scaling is observed for J 1/κt. We remark here that the typical scaling obtained in quantum metrology for unitary parameters is of order t 2 . The observed t 3 -scaling is due to the continuos monitoring of the system. A similar scaling of the Fisher information 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 4 -scaling can be observed. This is obtained in noiseless quantum metrology problems with time-dependent Hamiltonian and by exploiting open-loop control [60][61][62][63]. In particular in [63], it was also shown that a t 3 -scaling can be achieved without additional control, but by performing repeated (stroboscopic) measurement on the system, analogously to our strategy.
The previous results were both shown by considering perfect monitoring of the environment, i.e. for detectors with unit efficiency η. In Fig. 5 we plot the behaviours of 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 Eq. (30). We remind that if we consider only the classical FI F[p(y t )], the Heisenberg scaling in terms of J and t 3 -scaling are always obtained for κt 1 and for every η, as shown by the expansion (24). 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. (27) with η = 0, dominates (the term ηK 2 J 2 in (30) is negligible) and we observe SQL scaling forQ. We finally mention that the in the regimes where we observe Heisenberg scaling ofQ, the classical FI F[p(y t )] amounts to a relevant part of the total, namely around 25% .

C. 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 (11), is given by the QFI Q L in Eq. (9). The generalized master equation (10) in this case (considering the large-spin approximation) reads In D we show how this equation can be solved in a phase space picture, since the equation contains at most quadratic terms inX andP and thus preserves the Gaussian character of the operator¯ . The final result is i.e., we exactly obtain the effective QFIQ defined in Eq. (29) in the limit of unit efficiency η = 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 (11) 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.

V. CONCLUSION AND DISCUSSION
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 quantum Fisher information 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 2 -scaling in terms of spin, and a monitoring-enhanced t 3scaling in terms of the interrogation time, are obtained for J 1/κ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 superclassical 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][18][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.

ACKNOWLEDGMENTS
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.

Appendix A: Classical and quantum Cramér-Rao bounds for sequential non-demolition measurements
Here we will show how to derive the quantum Cramér-Rao bound for time-continuous homodyne monitoring reported in Eq. (5). We start by considering a (classical) estimation problem of a parameter θ described by a conditional probability p(z, y T |θ), where the vector y T = (y 1 , y 2 , . . . , y T ) 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 y T . The corresponding classical Fisher information can be evaluated as F[p(z, y T |θ)] = dy dz p(z, y T |θ) (∂ θ log p(z, y T |θ)) 2 = dy dz p(z|y T , θ)p(y T |θ) (∂ θ log p(z|y T , θ)) 2 + 2(∂ θ log p(z|y T , θ)) (∂ θ log p(y T |θ)) where the second expression has been obtained by means of the Bayes rule p(z, y T |θ) = p(z|y T , θ)p(y T |θ) .
In the following, we will omit the dependence on the parameter θ and we will denote by E p(x) [·] the average over a probability distribution p(x). By considering each term inside the integral separately one obtains where we have used the property dz (∂ θ p(z|y T )) = ∂ θ dz p(z|y T ) = ∂ θ (1) = 0. As a consequence, any unbiased estimatorθ based on M experiments, i.e. obtained collecting M series of measurement outcomes (y T , z), satisfies the generalized Cramér-Rao bound where the first term F[p(y T )] is the Fisher information corresponding to the sequential measurements with outcomes y T , while the second term is the average of the Fisher information F[p(z|y T )], corresponding to the final measurement over all the possible trajectories conditioned on the previous measurement results y T . The bound in Eq. (A5) 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(θ).
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 0 and after obtaining the stream of outcomes y T , reads whereM y T = M y T . . . M y2 M y1 and the probability of obtaining the outcomes y T reads p(y T |θ) = Tr[M y T 0M † y T ] [65]. One can then also perform a strong (destructive) measurement described by POVM operators {Π z } on the conditional state, and the whole measurement strategy is described by the conditional probabilities 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 {Π z }. We can then apply the quantum Cramér-Rao bound for the conditional states y T ]. One then obtains a more fundamental quantum Cramér-Rao bound for our estimation strategy Clearly this bound can be readily applied to the timecontinuous case discussed in the main text, where the vector of outcomes y T corresponds to a measured homodyne photocurrent, and where the conditional state (c) y T can be obtained via a stochastic master equation as the one in Eq. (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.

Appendix B: Gaussian dynamics and Gaussian Fisher information
Here we will provide the formulas for the dynamics of the atomic ensemble described by the stochastic master equation (15). 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 r c and of the covariance matrix σ of the quantum state (c) . These are defined in components as r j c = Tr r j (c) and σ jk = Tr {r j − r j c ,r k − r k c } (c) for the operator vectorr = (X,P ) T . In formulae one obtains [56,57]: where and dw is a vector of Wiener increments such that dw j dw k = δ jk dt, related to the photocurrent via the equation The Eqs. (18), (19) and (22) can be obtained from the ones above, remembering that for our definitions The method to calculate the Fisher information corresponding to the time-continuous measurement in the case of linear Gaussian system has been described in [58]. One has to evaluate the formula that, by plugging in the matrices describing our problem, is easily simplified to Eq. (21). 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 By noticing that the only non-zero entry of the vector ∂ B r c is the one corresponding to P (t) c , one easily obtain Eq. (26). 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 y t , by Bayes rule: where p(B) is the prior distribution, L(B|y t ) ≡ p(y t |B) is the likelihood and p(y t ) serves as a normalization factor. The Bayesian estimator is the mean of the posterior distribution and it is proven that the correspond- asymptotically optimal, i.e. tends to saturate the Cramér-Rao bound when the length of the vector y T is large. The simulated experimental run is obtained by numerically integrating the stochastic differential equation (18) with the Euler-Maruyama method for the "true" value of the parameter B true . Time is discretized with steps of length ∆t, i.e. to get from time 0 to time T we perform n T = T /∆t steps. Experimental data is represented by the observed measurement current y T = (∆y t1 , . . . , ∆y tn T ) T , which corresponds to an n T -dimensional vector. The outcome at every time step ∆y ti is sampled from a Gaussian distribution with variance ∆t and mean ∆y ti (B) = 2ηκJe − κt 2 P (t i ) c ∆t. Notice that ∆y ti (B) depends explicitly on the parameter B via the quantum expectation value P (t i ) 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 Eqs. (18) and (19) 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 by considering the outcomes as independent random variables, i.e. multiplying the corresponding probabilities. We then apply Bayes rule, Eq. (C1), 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.
where the displacement operator is defined aŝ In particular we will work in the phase space of a single mode system, so thatr = (X,P ) is the vector of quadrature operators and s = (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,66] Xρ If we now define the characteristic function associated to the operator¯ introduced in Eq. (10) the quantity of interest in order to compute the QFI is then Tr¯ =χ(0, t), as evident from the definition (D1). By applying the phase space mapping, from the generalized master equation (33) we get to the following partial differential equation for the characteristic function This equation can be solved by performing a Gaussian ansatz, similarly to [67], i.e. assuming that at every time the characteristic function can be written in the following form χ(s, t) = C(t) exp − 1 4 s Ω σ(t)Ω s + + is Ω s m (t) .
The dependence on time and on the parameters B 1/2 is completely contained in the covariance matrix σ(t), in the first moment vector s m (t) = (x m (t), p m (t)) and in the function C(t) =χ(0, t), which is the final result we are seeking. By plugging (D9) into (D8) 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 2 : (D12) These equations are solved analytically with the initial conditions σ 1,1 (0) = 1, x m (0) = 0 and C(0) = 1 (since for t = 0 the operator¯ corresponds to the initial state of the system |0 0|), yielding By plugging this term into Eq. (9), we finally obtain the ultimate QFI Q Ltn reported in Eq. (34).