Noisy atomic magnetometry in real time

The optical magnetometry scheme we consider, see figure 1(a), consists of an ensemble of N spin-1/2 atoms prepared in a coherent spin state (CSS) polarised along the x-direction [28], so that the initial mean and variance of the ensemble angular momentum operators, denoted by the vector $\hat{\boldsymbol{J}}(t)={({\hat{J}}_{x}(t),{\hat{J}}_{y}(t),{\hat{J}}_{z}(t))}^{\mathrm{T}}$ in the Heisenberg picture, read ${\langle \hat{\boldsymbol{J}}(0)\rangle }_{\text{CSS}}={(J,0,0)}^{\mathrm{T}}$ and ${\langle {{\Delta}}^{2}\hat{\boldsymbol{J}}(0)\rangle }_{\text{CSS}}={(0,J/2,J/2)}^{\mathrm{T}}$, respectively, where J = N/2. As shown in figure 1(b), one may then naturally visualise the distribution of the angular momentum with help of the Bloch sphere representation. The aim of the scheme is to measure and estimate in real-time the magnetic field B_t being directed along y, which fluctuates according to the Ornstein–Uhlenbeck (OU) stochastic process [55]:

$$\begin{equation}\mathrm{d}{B}_{t}=-\chi {B}_{t}\,\mathrm{d}t+\mathrm{d}{W}_{B},\end{equation} \tag{ 3 }$$

where by dW_B we denote the Wiener noise of zero mean, $\mathbb{E}\!\left[\mathrm{d}{W}_{B}\right]=0$, variance $\mathbb{E}\!\left[\mathrm{d}{W}_{B}^{2}\right]={q}_{B}\,\mathrm{d}t$. The magnetic field B_t induces a Larmor precession around the y-axis at the frequency ω_L(t) = γB_t , with γ being the effective (constant) gyromagnetic ratio. However, as discussed in more detail below, in our analysis we restrict to dynamics at small times in the presence of small magnetic fields, under which the angular momentum operator $\hat{\boldsymbol{J}}(t)$ deviates only slightly from pointing along the x-direction—the direction along which the atoms are initially pumped. In principle, this regime could also be achieved over larger timescales in experiments involving stroboscopic probing [56–59]. However, such implementations would then require the atoms to be initially pumped in the direction perpendicular to the magnetic field, as depicted in figure 1.

Zoom In Zoom Out Reset image size

The atomic ensemble is continuously monitored by exploiting the paramagnetic Faraday rotation effect [16], which twists the linear polarisation of the (off-resonant) light propagating through the ensemble in the z-direction—the light probe, see figure 1(a). As a result, a quantum non-demolition measurement is realised [33, 34] that is undertaken continuously in time [35–37], which can be effectively described by means of the homodyne-like continuous measurement [60] with the shot-noise in the measured signal taking the form of a Wiener process [43, 61]. Note that, even if higher-spin atoms are considered, the following measurement model still applies as long as the contribution from the atomic polarisability tensor component can be suppressed [43, 62, 63]. In particular, the photocurrent being measured at time t is proportional to the mean value of the collective atomic spin-component along the probe, ${\langle {\hat{J}}_{z}(t)\rangle }_{(\text{c})}$, i.e.:

$$\begin{equation}y(t)\mathrm{d}t=2\eta \sqrt{M}{\langle {\hat{J}}_{z}(t)\rangle }_{(\text{c})}\,\mathrm{d}t+\sqrt{\eta }\,\mathrm{d}W,\end{equation} \tag{ 4 }$$

where η is the detection efficiency, M is the measurement strength, and dW is the Wiener differential fulfilling $\mathbb{E}\!\left[\mathrm{d}{W}^{2}\right]=\mathrm{d}t$ according to the Itô lemma [55].

Importantly, due to the active influence the continuous measurement (4) exerts on the atoms, the ensemble dynamics becomes conditional—denoted by the subscript (c)—and the atoms evolve differently depending on a particular trajectory of the measurement outcomes registered up to (and including) time t, y _⩽t = {y(τ)}_0⩽τ⩽t [61]. In particular, the ensemble is characterized at time t by the quantum state conditioned on the past measurement record, ρ_(c)(t) ≡ ρ(t| y _⩽t ), which undergoes further conditional dynamics described by a stochastic master equation:

$$\begin{equation}\mathrm{d}{\rho }_{(\text{c})}\!(t)=-\text{i}{\gamma}{B}_{t}[{\hat{J}}_{y},{\rho }_{(\text{c})}\!(t)]\,\mathrm{d}t+\sum\limits _{\alpha =x,y,z}{\gamma }_{\alpha }\mathcal{D}[{\hat{J}}_{\alpha }]{\rho }_{(\text{c})}\!(t)\,\mathrm{d}t\end{equation} \tag{ 5 }$$

$$\begin{equation}\quad +M\mathcal{D}[{\hat{J}}_{z}]{\rho }_{(\text{c})}\!(t)\,\mathrm{d}t+\sqrt{M\eta }\mathcal{H}[{\hat{J}}_{z}]{\rho }_{(\text{c})}\!(t)\,\mathrm{d}W,\end{equation} \tag{ 6 }$$

where the superoperators $\mathcal{D}$ and $\mathcal{H}$ are defined as $\mathcal{D}[\mathcal{O}]\rho =\mathcal{O}\rho {\mathcal{O}}^{{\dagger}}-\frac{1}{2}({\mathcal{O}}^{{\dagger}}\mathcal{O}\rho +\rho {\mathcal{O}}^{{\dagger}}\mathcal{O})$ and $\mathcal{H}[\mathcal{O}]\rho =\mathcal{O}\rho +\rho {\mathcal{O}}^{{\dagger}}-\mathrm{T}\mathrm{r}[(\mathcal{O}+{\mathcal{O}}^{{\dagger}})\rho ]\rho$ for any (also non-Hermitian) operator $\mathcal{O}$ [30].

The first term in (5) arises simply from the free Hamiltonian, $\hat{H}={\gamma}{B}_{t}{\hat{J}}_{y}$, responsible for the Larmor precession of the ensemble spin. In contrast, both terms in (6) describe the continuous measurement of ${\hat{J}}_{z}$ introduced above [31, 32], of which the former is responsible for measurement-induced decoherence, or the backaction, while the latter constitutes the non-linear information gain provided when observing a particular measurement signal in (4) [30]. Furthermore, in our analysis we model all possible, e.g. environment-induced, decoherence mechanisms affecting the ensemble (as a whole) by introducing the second term in (5) with three dissipative components being parametrised by distinct effective rates, γ_α , in the three directions α ∈ {x, y, z}. Let us also emphasise that the Wiener differential, dW, appearing in (6) is the same as in the measurement dynamics (4)—a particular fluctuation of the photocurrent in (4) after being registered drives the conditional state of the ensemble according to (6).

Finally, let us highlight the difference between the fluctuations of the estimated magnetic field B_t , as dictated by the OU process (3), and the collective decoherence introduced above in (5), of which only the latter should be interpreted as a form of 'noise'—in the statistical sense, forcing the quantum state of ensemble to lose its purity over time. In particular, due to its presence, the expectation value of ${\hat{J}}_{\alpha }\ (\alpha =x,y,z)$ exponentially decreases with a rate of γ_α , whose inverse can be identified as the phenomenological (ensemble) spin-decoherence time ${T}_{2}^{\ast }$ [12]. Within the theory of open quantum systems [64], this may be associated with tracing out some environmental degrees of freedom or the environment itself monitoring continuously each component ${\hat{J}}_{\alpha }$, whose stochastic trajectory of outcomes is inaccessible and, hence, must be averaged out. On the contrary, the above dynamical model (3)–(6) describes the evolution along a single trajectory of the fluctuating magnetic field. As a result, the impact of field fluctuations on the performance in magnetometry should not be treated on equal grounds with the decoherence of the atomic ensemble, but rather associated with the inability to perfectly estimate the field in real time due to its stochastic nature. In fact, if one on purpose decided to ignore the field fluctuations and rather consider the effective dynamics averaged over all possible field trajectories (see appendix F), one would recover the above model of collective decoherence in the direction of the field (here, ${\hat{J}}_{y}$), but with a decoherence rate that accumulates as γ_y ∼ t² over time.

In order to define the regime of interest in which the atomic sensor operates, let us first consider the unconditional dynamics of the atomic ensemble,

$$\begin{equation}\mathrm{d}\rho (t)=-\text{i}{\gamma}{B}_{t}[{\hat{J}}_{y},\rho (t)]\,\mathrm{d}t+\sum\limits _{\alpha =x,y,z}{\gamma }_{\alpha }\mathcal{D}[{\hat{J}}_{\alpha }]\rho (t)\,\mathrm{d}t+M\mathcal{D}[{\hat{J}}_{z}]\rho (t)\,\mathrm{d}t,\end{equation} \tag{ 7 }$$

which can be simply obtained by dropping the stochastic term in (5) and (6), so that it now describes the evolution of the atomic state at time t, averaged over all the stochastic trajectories of the measured outcomes, i.e. [55]:

$$\begin{equation}\rho (t)={\mathbb{E}}_{p({\boldsymbol{y}}_{\leqslant t})}\!\left[{\rho }_{(\text{c})}(t)\right]=\int \!\mathrm{d}{\boldsymbol{y}}_{\leqslant t}\,p({\boldsymbol{y}}_{\leqslant t})\,\rho (t\vert {\boldsymbol{y}}_{\leqslant t}),\end{equation} \tag{ 8 }$$

where by ∫d y _⩽t we denote the integral over all possible measurement records y _⩽t = {y(τ)}_0⩽τ⩽t . Using equation (7) to define the dynamics of any observable $\hat{O}(t)$ in the Heisenberg picture, whose mean must then obey $\langle \hat{O}(t)\rangle =\mathrm{T}\mathrm{r}\!\left\{\hat{O}(t)\rho \right\}=\mathrm{T}\mathrm{r}\!\left\{\hat{O}\rho (t)\right\}$, we observe that the (unconditional) evolution of the mean spin-component along the direction of pumping, $\langle {\hat{J}}_{x}(t)\rangle$, is governed by the solution to the following set of coupled differential equations:

$$\begin{equation}\mathrm{d}\langle {\hat{J}}_{x}(t)\rangle ={\gamma}{B}_{t}\langle {\hat{J}}_{z}(t)\rangle \mathrm{d}t-\frac{1}{2}(M+{\gamma }_{y}+{\gamma }_{z})\langle {\hat{J}}_{x}(t)\rangle \mathrm{d}t,\end{equation} \tag{ 9 }$$

$$\begin{equation}\mathrm{d}\langle {\hat{J}}_{z}(t)\rangle =-{\gamma}{B}_{t}\langle {\hat{J}}_{x}(t)\rangle \mathrm{d}t-\frac{1}{2}({\gamma }_{x}+{\gamma }_{y})\langle {\hat{J}}_{z}(t)\rangle \mathrm{d}t,\end{equation} \tag{ 10 }$$

$$\begin{equation}\mathrm{d}{B}_{t}=-\chi {B}_{t}\,\mathrm{d}t+\mathrm{d}{W}_{B},\end{equation} \tag{ 11 }$$

of which the last (stochastic) one just corresponds to the OU process with $\mathbb{E}\!\left[\mathrm{d}{W}_{B}^{2}\right]={q}_{B}$, describing the magnetic-field fluctuations in (3).

Although this implies that the dynamics of $\langle {\hat{J}}_{x}(t)\rangle$ is stochastic, we show in appendix A that, if one focuses on short timescales such that $\bar{{\omega }_{\text{L}}(t)}\,t\ll 1$, where $\bar{{\omega }_{\text{L}}(t)}={\gamma}\vert \,{\bar{B}}_{t}\vert$ is the time-averaged Larmor frequency, then the mean spin-component along the direction of the pump must decay exponentially as follows:

$$\begin{equation}\langle {\hat{J}}_{x}(t)\rangle \approx J{\text{e}}^{-(M+{\gamma }_{y}+{\gamma }_{z})\,t/2}=J{\text{e}}^{-r\,t/2}\end{equation} \tag{ 12 }$$

with an effective decay rate: r = M + γ_y + γ_z . For this to be true, not only the time-average value of the magnetic field, ${\bar{B}}_{t}=\frac{1}{t}{\int }_{0}^{t}\mathrm{d}\tau {B}_{\tau }$, must be small $(\vert \,{\bar{B}}_{t}\vert \ll \frac{1}{{\gamma}t})$, but also the distribution of ${\bar{B}}_{t}$ should be narrow enough, which we ensure by verifying that $\left\vert \mathbb{E}\!\left[{\bar{B}}_{t}\right]\pm 2\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\!\left[\,{\bar{B}}_{t}\right]}\right\vert \lesssim \frac{1}{{\gamma}t}$ according to the 68–95–99.7 rule for normal distributions (see also figure 3(a)). We explicitly prove in appendix A that for approximation (12) to hold it is sufficient to consider short timescales such that rt ≲ 1, while also requiring the parameters of the OU process in (11) (or (3)) to obey:

$$\begin{equation}\chi \ll \frac{4\,r}{3}\quad \text{and}\quad {q}_{B}\lesssim \frac{3\,{r}^{3}}{4{{\gamma}}^{2}}.\end{equation} \tag{ 13 }$$

Furthermore, by restricting to short timescales at which the approximation (12) is valid, we effectively deal with spin dynamics in the regime in which $\hat{\boldsymbol{J}}(t)$ only slightly deviates from its initial x-polarisation. This allows us then to perform the Gaussian approximation [45, 65] and introduce the effective (time-dependent) quadrature operators:

$$\begin{equation}\hat{X}(t)=\frac{{\hat{J}}_{y}}{\sqrt{\vert \langle {\hat{J}}_{x}(t)\rangle \vert }}\approx \frac{{\hat{J}}_{y}}{\sqrt{J{\text{e}}^{-rt/2}}}\quad \text{and}\quad \hat{P}(t)=\frac{{\hat{J}}_{z}}{\sqrt{\vert \langle {\hat{J}}_{x}(t)\rangle \vert }}\approx \frac{{\hat{J}}_{z}}{\sqrt{J{\text{e}}^{-rt/2}}},\end{equation} \tag{ 14 }$$

which satisfy the canonical commutation relationship $[\hat{X}(t),\hat{P}(t)]=\text{i}\frac{{\hat{J}}_{x}}{\left\vert \langle {\hat{J}}_{x}(t)\rangle \right\vert }\approx \text{i}$, as long as $\vert \langle {\hat{J}}_{x}(t)\rangle \vert \gg 1$ [45, 65], which is assured within the approximation (12) whenever J ≫ 1—recall J = N/2  ≳  10⁵ in optical-magnetometry experiments [25–27]. As a consequence, it is enough to consider the evolution of the first and second moments of the quadratures, so that when focussing on the probed spin-component in the z-direction that fulfils ${\hat{J}}_{z}\approx \sqrt{J{\text{e}}^{-rt/2}}\hat{P}(t)$, its conditional dynamics (5) and (6) is completely described by the following set of equations:

$$\begin{equation}y(t)\mathrm{d}t=2\eta \sqrt{M}{\langle {\hat{J}}_{z}(t)\rangle }_{(\text{c})}\mathrm{d}t+\sqrt{\eta }\,\mathrm{d}W,\end{equation} \tag{ 15 }$$

$$\begin{equation}\mathrm{d}{B}_{t}=-\chi {B}_{t}\,\mathrm{d}t+\mathrm{d}{W}_{B}\quad (\text{with}\enspace \mathbb{E}\!\left[\mathrm{d}{W}_{B}^{2}\right]={q}_{B})\end{equation} \tag{ 16 }$$

$$\begin{equation}\mathrm{d}{\langle {\hat{J}}_{z}(t)\rangle }_{(\text{c})}=-{\gamma}{B}_{t}J{\text{e}}^{-rt/2}\,\mathrm{d}t+2\sqrt{\eta M}{\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}\mathrm{d}W,\end{equation} \tag{ 17 }$$

$$\begin{equation}\mathrm{d}{\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}=-4M\eta {\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}^{2}\mathrm{d}t+{\gamma }_{y}{J}^{2}{\text{e}}^{-rt}\,\mathrm{d}t,\end{equation} \tag{ 18 }$$

which assume $[\hat{X}(t),\hat{P}(t)]\approx \text{i}$ to hold—satisfied whenever J ≫ 1, rt ≲ 1 and the conditions (13) are valid. Note that equations (15)–(17) are linear with respect to one another, while equation (18) can be solved independently. Moreover, they involve only Gaussian (in particular, Wiener) stochastic-noise terms, and are derived based on the Gaussian approximation (14). Hence, we refer in short to the evolution described by (15)–(18) as the (conditional) dynamics of ${\hat{J}}_{z}$ in the linear-Gaussian regime.

Firstly, let us note that within the linear-Gaussian regime, the equations of motion (15)–(17) are independent of the decoherence rate γ_x , which is thus redundant. This is a consequence of the approximation (12) (see also appendix A) and can be intuitively explained—deviations of $\hat{\boldsymbol{J}}(t)$ from the x-direction are then too small for the collective noise generated by the $\mathcal{D}[{\hat{J}}_{x}]$-term in equation (5) to have any effect on the quadratures (14). Furthermore, the decoherence rate γ_z is also unnecessary, since by transforming the parameters of the continuous measurement (4) as follows: M → M − γ_z , η → ηM/(M − γ_z ) and $y\to y\sqrt{M/(M-{\gamma }_{z})}$; we retrieve the conditional dynamics (15)–(18) with γ_z = 0. Hence, the impact of the collective noise generated by the $\mathcal{D}[{\hat{J}}_{z}]$-term in equation (5) instead, can always be interpreted and incorporated into a modified form of the continuous measurement (4).

As a result, without loss of generality, we restrict the effective decay rate of $\langle {\hat{J}}_{x}(t)\rangle$ introduced in (12) to take the form r = M + γ_y , when considering the (conditional) dynamics of ${\hat{J}}_{z}$ in (15)–(17). We then solve for the dynamics of ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}$ in (18), which constitutes a fully decoupled differential equation. For the complete analytical solution we refer the reader to appendix B, where we also show that the time-dependence of the (conditional) ${\hat{J}}_{z}$-variance (satisfying ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(0)\rangle }_{(\text{c})}={\langle {{\Delta}}^{2}{\hat{J}}_{z}(0)\rangle }_{\text{CSS}}=J/2$ at t = 0) can be described by two distinct short-time and long-time behaviours:

where ${t}^{\ast }={(2J\sqrt{M{\gamma }_{y}\eta })}^{-1}$ is the transition time between the two regimes.

Importantly, note that t ≪ t* implies $2Jt{\gamma }_{y}\ll \sqrt{{\gamma }_{y}/(\eta M)}$. Then, if also γ_y < ηM, we may infer 2Jtγ_y ≪ 1 and approximate 1 + 2Jtγ_y ≈ 1 in (19a). As a result, we then recover the noiseless (γ_y = 0) solution for the variance within the short-time regime, previously found by Geremia et al [44], i.e.:

$$\begin{equation}{\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}=\frac{J}{2}\,\frac{1+2Jt{\gamma }_{y}}{1+2JtM\eta }{\text{e}}^{-(M+{\gamma }_{y})t/2}\approx \frac{J}{2+4JtM\eta },\end{equation} \tag{ 20 }$$

despite the actual presence of the noise (γ_y > 0). In fact, we prove also in appendix B that ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}$ is a non-decreasing function at t ≈ 0 if γ_y ⩾ ηM. Hence, the noise may be considered negligible at small times only if γ_y < ηM.

The dynamics of ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}$ enables us to study the phenomenon of spin-squeezing of the atomic ensemble [28] by directly computing the (conditional) squeezing parameter introduced by Wineland et al [66] along the pumping x-direction:

$$\begin{equation}{\xi }^{2}(t)=\frac{{\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}}{{\langle {\hat{J}}_{x}(t)\rangle }^{2}}{\left[\frac{{\langle {{\Delta}}^{2}{\hat{J}}_{z}(0)\rangle }_{\text{CSS}}}{{\langle {\hat{J}}_{x}(0)\rangle }_{\text{CSS}}^{2}}\right]}^{-1}=\frac{2J{\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}}{{\langle {\hat{J}}_{x}(t)\rangle }^{2}}\end{equation} \tag{ 21 }$$

which when satisfying ξ²(t) < 1 indicates a gain in interferometric precision over the CSS [66]. In particular, just like we decomposed ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}$ according to (19), we can similarly split the behaviour of the squeezing parameter in the linear-Gaussian regime, in which $\langle {\hat{J}}_{x}(t)\rangle \approx J{\text{e}}^{-r\,t/2}$, as follows

From now on in all plots, in order to confirm that apart from (13) the condition rt ≲ 1 (with now r = M + γ_y ) is satisfied and, hence, the linear-Gaussian approximation (12) is valid, we introduce the rescaled time t_S = (M + γ_y )t and limit its range to t_S ≲ 1, what assures the timescale of the linear-Gaussian regime to always be consistently considered.

In figure 2 we present explicitly the exact dynamics of the squeezing parameter (21) for two important cases: (a) when γ_y < ηM and the spin-squeezing (ξ²(t) < 1) occurs within a finite-time window (see figure 2(a)); and (b) when γ_y ⩾ ηM for which spin-squeezing is forbidden (see figure 2(b)). In both cases, it is evident that the exact solution for ξ²(t) very closely follows the two-regime behaviour in (22), with a clear transition at t ≈ t*. Moreover, as seen explicitly from the two-regime solution (22)—and the exact solution that can be straightforwardly derived from the exact evolution of the ${\hat{J}}_{z}$-variance (19) described in appendix B—the dynamics of the squeezing parameter is specified solely by the properties of the continuous measurement (η and M), collective decoherence (γ_y ), and the number of atoms (J = N/2).

Zoom In Zoom Out Reset image size

In what follows, we turn to the problem of sensing the magnetic field, B_t , in real time by constructing the optimal estimator, ${\tilde{B}}_{t}$, of its fluctuating value in the form of a KF [51, 52]. However, we have already included in figure 3 the evolution of the corresponding minimal estimation error, ${{\Delta}}^{2}{\tilde{B}}_{t}$, in order to emphasise the fact that the stochasticity of the B_t -signal affects the estimation procedure, but not the spin-squeezing phenomenon. In particular, for a given magnetic-field trajectory generated according to the OU process (3)—of a magnitude that does not invalidate the linear-Gaussian regime, see figure 3(a)—we obtain a measurement signal y(t) from which we can infer via (4) the conditional dynamics of the mean spin-component ${\langle {\hat{J}}_{z}(t)\rangle }_{(\text{c})}$, see figure 3(b). Now, the evolution of its variance, ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}$, allows us to compute the squeezing parameter ξ²(t), which, as described above, evolves in the linear-Gaussian regime in the same manner for any particular stochastic trajectory of the magnetic field—see figure 3(c). In contrast, the stochastic fluctuations of B_t affect the estimation procedure, so that the KF reaches a constant error (the SS) after a certain time despite the atomic ensemble still being spin-squeezed (ξ²(t) < 1). This can be directly seen by comparing the evolution of ${{\Delta}}^{2}{\tilde{B}}_{t}$ in figure 3(d) against the one of ξ²(t) in figure 3(c) for t_S   ≳  10⁻³. The transition time when this occurs depends on the parameters of the OU process (3)—as later shown in section 4.2.

Zoom In Zoom Out Reset image size

In the absence of collective noise and field fluctuations (γ_y = 0 and q_B = χ = 0), the Riccati equation (33) can be explicitly solved [44]. Since γ_y = 0, the conditional dynamics of the ${\hat{J}}_{z}$-variance is then given by equation (20). Moreover, since no fluctuations of the field are being considered, the noise correlation matrix just reads Q_t = diag(1, 0). As a result, the decoupled system of differential equations introduced in (34) and (35) can be analytically solved, providing the solution for the minimal aMSE in estimating the B-field (23) as

$$\begin{equation}{{\Delta}}^{2}{\tilde{B}}_{t}^{\text{HL}}({\sigma }_{0})={[{\mathbf{\Sigma }}_{t}]}_{22}=\frac{{M}^{2}}{16\eta {{\gamma}}^{2}{J}^{2}}\,\frac{(1+2JMt\eta )}{a(t){\text{e}}^{-Mt}+4\,(1+4J\eta ){\text{e}}^{-Mt/2}+b(t)},\end{equation} \tag{ 36 }$$

where we have emphasised its dependence on the width of the prior, σ₀, and defined:

$$\begin{equation}a(t)=-(1+2\eta J(4+Mt)),\end{equation} \tag{ 37 }$$

$$\begin{equation}b(t)=\frac{{M}^{2}}{16\eta {{\gamma}}^{2}{J}^{2}{\sigma }_{0}^{2}}+\frac{{M}^{3}t}{8{{\gamma}}^{2}J{\sigma }_{0}^{2}}+c(t).\end{equation} \tag{ 38 }$$

with c(t) = (Mt − 3) + 2ηJ(Mt − 4). Note that for an infinitely wide prior, σ₀ → ∞, b(t) = c(t) and the aMSE (36) matches consistently the solution obtained in [44] with ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{HL}}(\infty )$ exhibiting the non-classical scaling in both t and J—following the HL ∼1/N² with N = J/2—as mentioned already in equation (2). This becomes clear after making the approximation summarised in figure 4, i.e.:

where the terms (39a) and (39b) are obtained by Taylor-expanding ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{HL}}(\infty )$ in equation (36) to leading order in time t and (JMt)⁻¹, respectively.

Zoom In Zoom Out Reset image size

Within the Bayesian approach to estimation theory there exist various general lower bounds on the minimal aMSE, which are commonly referred to as Bayesian (or global) Cramér–Rao Bounds (BCRBs) [72]. Still, in the case of linear and Gaussian stochastic processes many of these simplify to a single form [40, 41]. When interested in estimating the process at its last step—here, the magnetic field B_t at the final time t—the applicable bound is the so-called marginal unconditional BCRB [41]:

$$\begin{equation}{{\Delta}}^{2}{\tilde{B}}_{t}\geqslant {({J}_{B})}^{-1},\end{equation} \tag{ 40 }$$

where J_B is the Bayesian information (BI) [40] evaluated at time t:

$$\begin{equation}{J}_{B}={\mathbb{E}}_{p({B}_{t},{\boldsymbol{y}}_{\leqslant t})}\!\left[{\left({\partial }_{{B}_{t}}\,\mathrm{log}\,p({B}_{t},{\boldsymbol{y}}_{\leqslant t})\right)}^{2}\right].\end{equation} \tag{ 41 }$$

The BI can be conveniently split into two terms:

$$\begin{equation}{J}_{B}={J}_{P}+{J}_{M},\end{equation} \tag{ 42 }$$

where J_P represents our knowledge about B_t prior to any estimation, and J_M accounts for the contribution of the measurement record y _⩽t . Formally, the BI associated with our prior knowledge, J_P , corresponds to the Fisher information (FI) [13] evaluated with respect to the estimated field value B_t for the distribution p(B_t ) defined in (24), i.e.:

$$\begin{equation}{J}_{P}=F[p({B}_{t})]={\mathbb{E}}_{p({B}_{t})}\left[{\left({\partial }_{{B}_{t}}\,\mathrm{log}\,p({B}_{t})\right)}^{2}\right].\end{equation} \tag{ 43 }$$

In contrast, J_M reads

$$\begin{equation}{J}_{M}={\mathbb{E}}_{p({B}_{t},{\boldsymbol{y}}_{\leqslant t})}\left[{\left({\partial }_{{B}_{t}}\,\mathrm{log}\,p({\boldsymbol{y}}_{\leqslant t},{B}_{t})\right)}^{2}\right]=\int \mathrm{d}{B}_{t}\,p({B}_{t})F[p({\boldsymbol{y}}_{\leqslant t}\vert {B}_{t})],\end{equation} \tag{ 44 }$$

where

$$\begin{equation}F[p({\boldsymbol{y}}_{\leqslant t}\vert {B}_{t})]={\mathbb{E}}_{p({\boldsymbol{y}}_{\leqslant t}\vert {B}_{t})}\!\left[{\left({\partial }_{{B}_{t}}\,\mathrm{log}\,p({\boldsymbol{y}}_{\leqslant t}\vert {B}_{t})\right)}^{2}\right]\end{equation} \tag{ 45 }$$

is now the FI of the distribution p( y _⩽t |B_t ) evaluated again with respect to B_t .

The prior contribution to the BI (43) can always be ignored by letting the prior distribution describing our knowledge about the field at t = 0, p(B₀), be of infinite width [40]—see appendix E.1 where we demonstrate that J_P = 0 for ${\sigma }_{0}^{2}\to \infty$. In contrast, the contribution of the measurement record to the BI (44) in general requires one to explicitly evaluate the FI of p( y _⩽t |B_t ) for every value B_t may take. In this work, we avoid doing so by constructing an upper-bound on F[p( y _⩽t |B_t )] that is determined by the noise in the atomic magnetometry scenario, which turns out to be independent of the form of the continuous measurement and the initial state of the atoms.

In order to do so, we discretize the time evolution such that t = kδt—this does not prevent us from obtaining the results in the continuous-time limit, which can be recovered by finally letting δt → 0 [30]. As a result, the continuous-measurement record y _⩽t can be described by a finite set of (time-ordered) outcomes, y _k = {y₀, y₁, ..., y_k }, collected at the end of each time interval indexed by j = 0, 1, ..., k. Similarly, the evolution of the magnetic field corresponds now to a discrete-time stochastic process, B _k = {B₀, B₁, ..., B_k }, with B_j representing the value taken by the B-field throughout the jth time interval. Moreover, the conditional quantum state of the ensemble at time t = kδt, previously described by ρ_(c)(t) ≡ ρ(t| y _⩽t ) as the solution to the conditional dynamics (6), can now be explicitly written as

$$\begin{equation}{\rho }_{k}\equiv \rho (t\vert {\boldsymbol{y}}_{k})=\frac{{E}_{{y}_{k}}\,{{\Lambda}}_{{B}_{k}}\!\left[\dots {E}_{{y}_{1}}\,{{\Lambda}}_{{B}_{1}}\!\left[{E}_{{y}_{0}}\,{{\Lambda}}_{{B}_{0}}\!\left[{\rho }_{0}\right]{E}_{{y}_{0}}^{{\dagger}}\right]{E}_{{y}_{1}}^{{\dagger}}\dots \,\right]{E}_{{y}_{k}}^{{\dagger}}}{\mathrm{T}\mathrm{r}\!\left\{{E}_{{y}_{k}}\,{{\Lambda}}_{{B}_{k}}\!\left[\dots {E}_{{y}_{1}}\,{{\Lambda}}_{{B}_{1}}\!\left[{E}_{{y}_{0}}\,{{\Lambda}}_{{B}_{0}}\!\left[{\rho }_{0}\right]{E}_{{y}_{0}}^{{\dagger}}\right]{E}_{{y}_{1}}^{{\dagger}}\dots \,\right]{E}_{{y}_{k}}^{{\dagger}}\right\}},\end{equation} \tag{ 46 }$$

where the continuous measurement is formally described by a set of measurement operators ⁵ ${E}_{{y}_{j}}$ that form a positive-operator valued measure (POVM), i.e. ${\left\{{E}_{{y}_{j}}^{{\dagger}}{E}_{{y}_{j}}\right\}}_{{y}_{j}}$ with ${\sum }_{{y}_{j}}\;{E}_{{y}_{j}}^{{\dagger}}{E}_{{y}_{j}}=1$ for every j. In the expression (46), every ${{\Lambda}}_{{B}_{j}}$ denotes the quantum channel (a completely positive trace-preserving map) describing the evolution of the atomic ensemble during the jth interval in between measurements. Such evolution will then solely incorporate the Larmor precession under the magnetic field B_j (constant in that time interval), and the collective decoherence. A scheme describing this sequence of channels interspersed with measurements is depicted in figure 5(a) (top).

Zoom In Zoom Out Reset image size

The denominator of the expression (46) should be interpreted as the probability of registering the measurement record y _k given a particular (discrete) trajectory of the B-field B _k , i.e.:

$$\begin{equation}p({\boldsymbol{y}}_{k}\vert {\boldsymbol{B}}_{k})=\mathrm{T}\mathrm{r}\!\left\{{E}_{{y}_{k}}\,{{\Lambda}}_{{B}_{k}}\!\left[\dots {E}_{{y}_{1}}\,{{\Lambda}}_{{B}_{1}}\!\left[{E}_{{y}_{0}}\,{{\Lambda}}_{{B}_{0}}\!\left[{\rho }_{0}\right]{E}_{{y}_{0}}^{{\dagger}}\right]{E}_{{y}_{1}}^{{\dagger}}\dots \,\right]{E}_{{y}_{k}}^{{\dagger}}\right\}.\end{equation} \tag{ 47 }$$

Moreover, its marginal distribution can then be determined using the Bayes' rule as

$$\begin{equation}p({\boldsymbol{y}}_{k}\vert {B}_{k})=\frac{p({\boldsymbol{y}}_{k},{B}_{k})}{p({B}_{k})}=\frac{1}{p({B}_{k})}\int \mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,p({\boldsymbol{y}}_{k}\vert {\boldsymbol{B}}_{k}),\end{equation} \tag{ 48 }$$

constituting the discrete-time equivalent of p( y _⩽t |B_t ) appearing in equation (44).

The crucial step in our construction is the observation that the dynamics of the atomic ensemble in between measurements can be CS [22, 53, 54]. In particular, each quantum channel ${{\Lambda}}_{{B}_{j}}$ in equation (46) can be equivalently replaced by a probabilistic mixture of unitary evolutions such that only the classical mixing probability depends on the instantaneous B-field value B_j —see figure 5(b). In other words, we 'shift' the encoding of the parameter B_t from a quantum channel to a classical probability, which we show to be Gaussian. Namely, for any j = 0, 1, ..., k:

$$\begin{equation}{{\Lambda}}_{{B}_{j}}[\,\cdot \,]=\int \!\mathrm{d}{\omega }_{j}\,q({\omega }_{j}\vert {B}_{j})\,{\mathcal{U}}_{{\omega }_{j},\delta t}[\,\cdot \,],\end{equation} \tag{ 49 }$$

where ω_j is an auxiliary frequency-like random variable distributed according to the Gaussian distribution:

$$\begin{equation}{\omega }_{j}\,\sim \,q({\omega }_{j}\vert {B}_{j})=\mathcal{N}\!\left({\gamma}{B}_{j},\frac{{\gamma }_{y}}{\delta t}\right)\!,\end{equation} \tag{ 50 }$$

which parametrises now a noiseless Larmor precession of the ensemble within the jth time interval described by the unitary transformation:

$$\begin{equation}{\mathcal{U}}_{{\omega }_{j},\delta t}[\,\cdot \,]={\text{e}}^{-\text{i}{\omega }_{j}\hat{H}\delta t}\,\cdot \,{\text{e}}^{\text{i}{\omega }_{j}\hat{H}\delta t}\end{equation} \tag{ 51 }$$

with the Hamiltonian $\hat{H}={\hat{J}}_{y}$ consistent with the dynamics (5). A step by step demonstration of this statement is presented in appendix E.2. Moreover, within the setting here considered, the random variable B_j follows a (discrete-time) OU process (3), while the initial B₀ is drawn from a Gaussian prior distribution of variance ${\sigma }_{0}^{2}$.

As a consequence, by now defining the set of particular frequencies ω_j chosen in each time interval, ω _k = {ω₀, ω₁, ..., ω_k }, we can interpret the conditional distribution of the measurement record (47) as a probabilistic average over the frequency-like parameters:

$$\begin{equation}p({\boldsymbol{y}}_{k}\vert {\boldsymbol{B}}_{k})={\mathbb{E}}_{q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})}\!\left[p({\boldsymbol{y}}_{k}\vert {\boldsymbol{\omega }}_{k})\right]=\int \mathrm{d}{\boldsymbol{\omega }}_{k}\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})\,p({\boldsymbol{y}}_{k}\vert {\boldsymbol{\omega }}_{k}),\end{equation} \tag{ 52 }$$

where $q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})={\prod }_{j=0}^{k}\;q({\omega }_{j}\vert {B}_{j})$ is just a product distribution and

$$\begin{equation}p({\boldsymbol{y}}_{k}\vert {\boldsymbol{\omega }}_{k})=\mathrm{T}\mathrm{r}\!\left\{{E}_{{y}_{k}}{\mathcal{U}}_{{\omega }_{k}\delta t}\left[\dots {E}_{{y}_{1}}{\mathcal{U}}_{{\omega }_{1}\delta t}\left[{E}_{{y}_{0}}{\mathcal{U}}_{{\omega }_{0}\delta t}\left[{\rho }_{0}\right]{E}_{{y}_{0}}^{{\dagger}}\right]{E}_{{y}_{1}}^{{\dagger}}\dots \,\right]{E}_{{y}_{k}}^{{\dagger}}\right\}\end{equation} \tag{ 53 }$$

is now the conditional probability of detecting a set of measurement records y _k , given a particular set of frequencies ω _k dictating subsequent unitary evolutions of the atomic ensemble within each time interval.

The decomposition (52) proves the equivalence between a circuit of non-unitary quantum channels describing the noisy dynamics of the atomic ensemble, and the average (${\mathbb{E}}_{q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})}\!\left[\dots \right]$ denoted in figure 5 with ⟨...⟩) of circuits involving noiseless unitary evolutions, whose frequencies ω _k are drawn from a distribution q( ω _k | B _k ) containing all the information about the field trajectory B _k . We represent this equivalence schematically in figure 5(a).

Focussing on the marginal distribution of interest (48), we substitute the decomposition (52) into it, in order to obtain

$$\begin{equation}p({\boldsymbol{y}}_{k}\vert {B}_{k})=\int \!\mathrm{d}{\boldsymbol{\omega }}_{k}\,p({\boldsymbol{y}}_{k}\vert \,{\boldsymbol{\omega }}_{k})\left[\frac{1}{p({B}_{k})}\int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})\right]={\mathcal{S}}_{{\boldsymbol{\omega }}_{k}\to {\boldsymbol{y}}_{k}}\!\left[{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right],\end{equation} \tag{ 54 }$$

where we can now identify ${\mathcal{S}}_{{\boldsymbol{\omega }}_{k}\to {\boldsymbol{y}}_{k}}[\,\cdot \,]=\int \!\mathrm{d}{\boldsymbol{\omega }}_{k}\,p({\boldsymbol{y}}_{k}\vert \,{\boldsymbol{\omega }}_{k})\,\cdot$ as a stochastic map independent of the estimated B_k , and define the probability distribution function:

$$\begin{equation}{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})=\frac{1}{p({B}_{k})}\int \mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k}),\end{equation} \tag{ 55 }$$

which effectively describes the information about B_k contained within q( ω _k | B _k ).

Crucially, as the FI is generally contractive (monotonic) under the action of stochastic maps [13, 73], we can now upper-bound F[p( y _⩽t |B_t )] in (45) as follows

$$\begin{equation}F[p({\boldsymbol{y}}_{k}\vert {B}_{k})]=F\!\left[{\mathcal{S}}_{{\boldsymbol{\omega }}_{k}\to {\boldsymbol{y}}_{k}}\!\left[{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right]\right]\leqslant F\left[{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right],\end{equation} \tag{ 56 }$$

and evaluate the FI of ${\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})$, which we denote as ${\mathbb{P}}_{{B}_{t}}({\boldsymbol{\omega }}_{< t})$ upon letting δt → 0 (with t = kδt) to recover the continuous-time limit. In appendix E.3, we explicitly show that for the case of infinitely wide (σ₀ → ∞) Gaussian prior distribution p(B₀), the FI reads

$$\begin{equation}F\left[{\mathbb{P}}_{{B}_{t}}({\boldsymbol{\omega }}_{< t})\right]=\sqrt{\frac{{{\gamma}}^{2}}{{\gamma }_{y}{q}_{B}}}\,\mathrm{tanh}\!\left(t\sqrt{\frac{{q}_{B}{{\gamma}}^{2}}{{\gamma }_{y}}}\right),\end{equation} \tag{ 57 }$$

and being independent of the estimated B_t directly sets an upper bound in J_M , i.e.:

$$\begin{equation}{J}_{M}=\int \mathrm{d}{B}_{t}\,p({B}_{t})F[p({\boldsymbol{y}}_{\leqslant t}\vert {B}_{t})]\end{equation} \tag{ 58 }$$

$$\begin{equation}\leqslant \int \mathrm{d}{B}_{t}\,p({B}_{t})F\left[{\mathbb{P}}_{{B}_{t}}({\boldsymbol{\omega }}_{< t})\right]=\sqrt{\frac{{{\gamma}}^{2}}{{\gamma }_{y}{q}_{B}}}\,\mathrm{tanh}\!\left(t\sqrt{\frac{{q}_{B}{{\gamma}}^{2}}{{\gamma }_{y}}}\right).\end{equation} \tag{ 59 }$$

Hence, by using also the fact that J_P = 0 when σ₀ → ∞, we conclude that

$$\begin{equation}{{\Delta}}^{2}{\tilde{B}}_{t}\geqslant {\left({J}_{P}+{J}_{M}\right)}^{-1}={J}_{M}^{-1}\end{equation} \tag{ 60 }$$

$$\begin{equation}\geqslant \sqrt{\frac{{\gamma }_{y}\,{q}_{B}}{{{\gamma}}^{2}}}\,\mathrm{coth}\left(t\sqrt{\frac{{q}_{B}{{\gamma}}^{2}}{{\gamma }_{y}}}\right)\equiv {{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}({q}_{B})\end{equation} \tag{ 61 }$$

$$\begin{equation}\geqslant \frac{1}{{{\gamma}}^{2}}\frac{{\gamma }_{y}}{t}={{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}(0),\end{equation} \tag{ 62 }$$

where the last inequality follows from x coth(x) ⩾ 1 for any x > 0 with the CS limit, ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}({q}_{B})$, monotonically increasing with q_B —as expected, the error is predicted to deteriorate with an increase in the strength of the field fluctuations.

Although we differ the full proof of the upper bound (59) to appendix E.3, let us note that in the simplest scenario with the magnetic field being time-invariant, the CS limit can be straightforwardly derived and takes the form ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}(0)$ as in equation (62). For a constant B-field, the probability distribution (55) simplifies to the product ${\mathbb{P}}_{B}({\boldsymbol{\omega }}_{k})=q({\boldsymbol{\omega }}_{k}\vert B)={\prod }_{j=0}^{k}q({\omega }_{j}\vert B)$, whose FI is now evaluated with respect to B and reads:

$$\begin{equation}F\left[{\mathbb{P}}_{B}({\boldsymbol{\omega }}_{k})\right]=F\left[\prod\limits _{j=0}^{k}q({\omega }_{j}\vert B)\right]=\sum\limits _{j=0}^{k}F\left[q({\omega }_{j}\vert B)\right]=k\,F\left[q({\omega }_{j}\vert B)\right]={{\gamma}}^{2}\frac{k\delta t}{{\gamma }_{y}}={{\gamma}}^{2}\frac{t}{{\gamma }_{y}},\end{equation} \tag{ 63 }$$

so that the CS limit (62) then directly follows.

In general, the CS limit ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}({q}_{B})$ defined in (61) can always be approximated independently of the value of q_B by splitting the timescales into two regimes, as follows:

where the expression (64a) implies that the CS limit (61) always simplifies at short times to its form applicable in absence of field fluctuations, i.e. ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}(0)$ in (62).

Finally, let us highlight that the derivation of the CS limit is independent of the measurement dynamics and the initial state, and depends only on the form of the noisy quantum channel describing the evolution of the system in each discretised step between subsequent measurements.

Returning to the particular measurement model and the estimation strategy considered, we determine the SS solution of the KF derived in section 3, as it may then be directly compared with the fundamental CS limit determined above. Still, for simplicity, we focus here only on the SS solution when the magnetic field fluctuates according to the Wiener process [55], i.e. the special case of the OU process (3) with χ = 0. The complete solution for χ > 0 is presented in appendix D.

In general, the SS is attained by the KF when its covariance matrix Σ_t no longer changes in time, so that $\frac{\mathrm{d}{\mathbf{\Sigma }}_{t}}{\mathrm{d}t}=0$ in equation (33). This corresponds to the solution of equating the rhs of (33) to zero—solving the corresponding (non-linear) algebraic Riccati equation [70]. In our case, finding the SS solution becomes much easier after noting that for t ≫ t* (guaranteed as t → ∞ for SS) ${\langle {{\Delta}}^{2}{\hat{J}}_{z}(t)\rangle }_{(\text{c})}\approx {V}_{ > {t}^{\ast }}(t)$, as in (19b). Then, the SS solution for ${{\Delta}}^{2}{\tilde{B}}_{t}={[{\mathbf{\Sigma }}_{t}]}_{22}$ can be shown (see appendix D) to read

$$\begin{equation}{{\Delta}}^{2}{\tilde{B}}_{t}^{\text{SS}}={\left(\frac{{q}_{B}{\gamma }_{y}}{{{\gamma}}^{2}}+\frac{1}{{\gamma}J}\sqrt{\frac{{q}_{B}^{3}}{M\eta }}{\text{e}}^{(M+{\gamma }_{y})t/2}\right)}^{1/2},\end{equation} \tag{ 65 }$$

whose second term in the parenthesis—the one surviving in the absence of noise (γ_y = 0)—has been derived previously in [46]. In contrast, it is the first term arising due to the decoherence considered here (γ_y > 0) that always dominates for large J.

In particular, at the relevant timescales $t\lesssim {(M+{\gamma }_{y})}^{-1}$ whenever $J\gtrsim \frac{{\gamma}}{{\gamma }_{y}}\sqrt{{q}_{B}/(\eta M)}$, the second term within the brackets in (65) becomes negligible and the SS solution (65) coincides with the CS limit (64b) valid at long timescales, i.e.:

$$\begin{equation}{{\Delta}}^{2}{\tilde{B}}_{t}^{\text{SS}}\,\underset{J\gg 1}{=}\,{{\Delta}}^{2}{\tilde{B}}_{\infty }^{\text{CS}}({q}_{B})=\frac{\sqrt{{q}_{B}{\gamma }_{y}}}{{\gamma}},\end{equation} \tag{ 66 }$$

which moreover implies that—see also the main plot of figure 6—${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{SS}}\approx {{\Delta}}^{2}{\tilde{B}}_{t}^{\text{CS}}({q}_{B})$, as long as further $t\;\gtrsim \;\frac{1}{{\gamma}}\sqrt{{\gamma }_{y}/{q}_{B}}$ for which the CS limit takes the form (64b).

Zoom In Zoom Out Reset image size

Note that this proves that in presence of decoherence (γ_y > 0) for large enough ensemble size (J = N/2 ≫ 1) the chosen type of continuous measurement (homodyne-like model) and the initial state of the atomic ensemble (CSS) always give the best possible precision in estimating the fluctuating magnetic field (q_B > 0) in the SS regime—bearing in mind that the linear-Gaussian approximation made throughout our work must hold. Crucially, this conclusion holds for large atomic ensembles also at short timescales at which the SS solution does not apply, as we will now explicitly show.

Once the CS limit and the SS solution of the KF have been explicitly derived, we finally have all the information we need to fully understand the evolution of the estimation error in figure 3(d), which we supplement now with additional plots of the aMSE with respect to time and ensemble size (J = N/2) in figures 6 and 7, respectively.

Zoom In Zoom Out Reset image size

4.3.1. Estimation error as a function of time

We depict the time evolution of the estimation error, i.e. the aMSE of the KF, in figure 6. Importantly, as shown in the main plot, for 'large' ensembles (we quantify 'large' below) the aMSE attains the fundamental CS limit (61) that holds for any potential measurement and estimation strategy, hence, proving that the magnetometry scheme achieves then the ultimate precision dictated by the decoherence.

In such a setting, the behaviour of error as a function of time can be intuitively explained. Firstly, at very short timescales for which the atomic decoherence can be ignored, the aMSE follows the noiseless-like (or supra-classical) scaling of 1/t³. Once the impact of decoherence 'kicks in', its behaviour transitions to a 'classical'-like regime exhibiting 1/t scaling, as determined by the CS limit in (64a). At long times, at which the field fluctuations are optimally compensated by the KF, it reaches its SS and the error saturates at a constant value $\sqrt{{q}_{B}{\gamma }_{y}}/{\gamma}$. However, let us emphasise again that this value arises due to the presence of decoherence (γ_y > 0), which in this case dominates the SS solution of the KF in equation (66) and, in fact, assures the error to attain the ultimate CS limit, in particular, its long-time behaviour (64b).

In contrast, for 'small' ensembles, as presented in the inset, the aMSE does not reach the CS limit. It is so, as the impact of the atomic decoherence can then be actually neglected in comparison to the field fluctuations. As a consequence, the KF reaches quickly with time its SS solution, which is now effectively independent of the decoherence, i.e. ${{\Delta}}^{2}{\tilde{B}}_{t}^{\text{SS}}\approx \sqrt[4]{{q}_{B}^{3}/(M\eta {{\gamma}}^{2}{J}^{2})}$ in equation (65). One can then interpret the magnetometer to operate in a 'perfect' manner, with its precision being dictated fully by just the characteristics of the field.

Note that the above two settings are formally distinguished by whether the SS solution in (65) attains or not the long-time CS limit (64b)—or, in other words, whether the green line can cross the black line in figure 6. Hence, this condition defines naturally what constitutes a 'large' or 'small' ensemble above, i.e. $J\;\gtrsim \;{J}_{\text{CS}}^{\prime }$ or $J\lesssim {J}_{\text{CS}}^{\prime }$, respectively, where

$$\begin{equation}{J}_{\text{CS}}^{\prime }=\frac{{\gamma}}{{\gamma }_{y}}\sqrt{\frac{{q}_{B}}{\eta M}}\end{equation} \tag{ 67 }$$

follows from the derivation of equation (66).

Moreover, by evaluating the times at which the dominant behaviours of the error appearing in figure 6 cross one another, we explicitly determine the transition times between all the different regimes:

$$\begin{equation}{t}_{\text{CS}}=\frac{1}{J}\sqrt{\frac{3}{\eta M{\gamma }_{y}}},\hspace{25.0pt}{t}_{\text{CS}}^{\prime }=\frac{1}{{\gamma}}\sqrt{\frac{{\gamma }_{y}}{{q}_{B}}},\hspace{25.0pt}{t}_{\text{SS}}=\frac{{3}^{1/3}}{{{\gamma}}^{2}{J}^{2}\eta M{q}_{B}},\end{equation} \tag{ 68 }$$

which are valid as long as $t\lesssim {(M+{\gamma }_{y})}^{-1}$ (and the linear-Gaussian approximation holds). In particular, t_CS marks the transition time of the aMSE from the noiseless-like to the 'classical'-like regime, whose later transition to the steady-state regime occurs then at ${t}_{\text{CS}}^{\prime }$. Similarly, for 'small' ensembles t_SS indicates when the aMSE reaches its SS.

4.3.2. Estimation error as a function of the number of atoms

The aMSE as a function of the ensemble size (J = N/2) is presented in figure 7, where we show two distinct scenarios differentiated by the timescales considered. Still, we observe that the ultimate CS limit dictated by the decoherence is attained in both cases as long as a sufficiently large ensemble is considered—proving then the measurement scheme and the estimation procedure to be optimal.

In particular, for 'short' timescales (main plot), the error as a function of J initially follows the Heisenberg-like scaling 1/J² before attaining the short-time CS limit (64a) once the effect of the collective noise becomes significant. For 'long' timescales (inset), the KF optimally compensates for the field fluctuations and its SS solution applies as long as J is sufficiently large. Still, the long-time CS limit (64b) is achieved as J → ∞, as it dictates in this case the SS solution (66) affected by the collective noise (γ_y > 0). The time threshold differentiating between these two cases, ${t}_{\text{CS}}^{\prime }$, is then given by (64)—coinciding consistently with the definition in (68).

Importantly, both when the aMSE (solid blue line) attains the CS limit at 'short' timescales $t\lesssim {t}_{\text{CS}}^{\prime }$ (the main plot of figure 7), or at 'long' timescales $t\;\gtrsim \;{t}_{\text{CS}}^{\prime }$ (the inset of figure 7), it saturates at a constant value as J is increased. In each case, that occurs at different ensemble sizes: J_CS and ${J}_{\text{CS}}^{\prime }$, respectively. In the latter case, one can also identify a threshold value J_SS, above which the aMSE exhibits a $1/\sqrt{J}$-scaling before saturating at ${J}_{\text{CS}}^{\prime }$. All these can be determined by comparing the dominant behaviours of the error within each regime, and read:

$$\begin{equation}{J}_{\text{CS}}=\frac{1}{t}\sqrt{\frac{3}{\eta M{\gamma }_{y}}},\hspace{25.0pt}{J}_{\text{SS}}=\frac{{3}^{2/3}}{{\gamma}{t}^{2}\sqrt{\eta M{q}_{B}}},\end{equation} \tag{ 69 }$$

while ${J}_{\text{CS}}^{\prime }$ coincides consistently with the definition (67).

The marginal probability density function of the magnetic field at time t—or equivalently after the time-step k = t/δt—can be written as:

$$\begin{equation}p({B}_{k})=\int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})=\int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\,\prod\limits _{j=1}^{k}p({B}_{j}\vert {B}_{j-1})\,p({B}_{0}),\end{equation} \tag{ E.3 }$$

where the initial field B₀ is drawn from a Gaussian distribution which effectively specifies our a priori knowledge about the field at t = 0, i.e.:

$$\begin{equation}{B}_{0}\sim p({B}_{0})=\frac{1}{\sqrt{2\pi {\sigma }_{0}^{2}}}\,\mathrm{exp}\!\left(-\frac{{B}_{0}^{2}}{2{\sigma }_{0}^{2}}\right),\end{equation} \tag{ E.4 }$$

while each transition probability p(B_j |B_j−1) for all j = 1, 2, ..., k is given by the OU process (3) as a Gaussian distribution [55]:

$$\begin{equation}p({B}_{j}\vert {B}_{j-1})=\sqrt{\frac{1}{2\pi {V}_{P}}}\,\mathrm{exp}\!\left(-\frac{{({B}_{j}-{B}_{j-1}{\text{e}}^{-\chi \delta t})}^{2}}{2{V}_{P}}\right)\end{equation} \tag{ E.5 }$$

with variance

$$\begin{equation}{V}_{P}=\frac{{q}_{B}}{2\chi }(1-{\text{e}}^{-2\chi \delta t}).\end{equation} \tag{ E.6 }$$

Hence, after explicitly evaluating the multiple integrals in (E.3), we arrive at the marginal Gaussian distribution for B_k :

$$\begin{equation}p({B}_{k})=\sqrt{\frac{1}{2\pi {V}_{P}^{(k)}}}\,\mathrm{exp}\!\left(-\frac{{B}_{k}^{2}}{2{V}_{P}^{(k)}}\right),\end{equation} \tag{ E.7 }$$

whose mean is zero and its variance

$$\begin{equation}{V}_{P}^{(k)}={\sigma }_{0}^{2}{\text{e}}^{-2k\chi \delta t}+\frac{{q}_{B}}{2\chi }(1-{\text{e}}^{-2k\chi \delta t}),\end{equation} \tag{ E.8 }$$

which in the case of a magnetic field following a Wiener process—a special case of the OU process (3) with χ = 0—simplifies to ${\mathrm{lim}}_{\chi \to 0}\,{V}_{P}^{(k)}=k\delta t{q}_{B}+{\sigma }_{0}^{2}$.

In general, the FI evaluated with respect to a Gaussian distribution, e.g. the marginal distribution (E.7), corresponds to the inverse of its variance. Hence, the prior contribution to the BI defined in equation (43) of the main text simply reads

$$\begin{equation}{J}_{P}=F[p({B}_{k})]=\frac{1}{{V}_{P}^{(k)}}=\frac{1}{{\sigma }_{0}^{2}{\text{e}}^{-2k\chi \delta t}+\frac{{q}_{B}}{2\chi }(1-{\text{e}}^{-2k\chi \delta t})},\end{equation} \tag{ E.9 }$$

and in the continuous-time limit of δt → 0 with k = t/δt, it becomes

$$\begin{equation}{J}_{P}=\frac{1}{{\sigma }_{0}^{2}{\text{e}}^{-2\chi t}+\frac{{q}_{B}}{2\chi }(1-{\text{e}}^{-2\chi t})}.\end{equation} \tag{ E.10 }$$

In the special case of a Wiener process (χ → 0), equation (E.10) reduces to ${\mathrm{lim}}_{\chi \to 0}\,{J}_{P}={({q}_{B}t+{\sigma }_{0}^{2})}^{-1}$. On the other hand, for any χ, q_B , t ⩾ 0, expression (E.10) vanishes if we do not possess any prior knowledge about the initial field B₀, i.e. when we consider σ₀ → ∞ in equation (E.4), for which ${\mathrm{lim}}_{{\sigma }_{0}\to \infty }\,{J}_{P}=0$.

Consider a unitary evolution, which is generated by $\xi \hat{H}$ with $\hat{H}$ being a time-independent Hamiltonian and $\xi \in \mathbb{R}$, being applied to a state ρ₀ for a time interval τ, i.e.:

$$\begin{equation}{\mathcal{U}}_{\xi ,\delta t}[{\rho }_{0}]={\text{e}}^{-\text{i}\xi \hat{H}\tau }{\rho }_{0}{\text{e}}^{\text{i}\xi \hat{H}\tau }\end{equation} \tag{ E.11 }$$

where the frequency-like parameter ξ is randomly distributed according to a Gaussian probability density:

$$\begin{equation}\xi \sim {p}_{{\mu }_{\tau },{\sigma }_{\tau }}(\xi )=\frac{1}{\sqrt{2\pi {\sigma }_{\tau }^{2}}}{\text{e}}^{-\frac{{(\xi -{\mu }_{\tau })}^{2}}{2{\sigma }_{\tau }^{2}}},\end{equation} \tag{ E.12 }$$

whose mean μ_τ and standard deviation σ_τ are some smooth time-dependent functions.

Theorem 1. The quantum map Λ_τ obtained by averaging over ξ after a time τ, i.e.:

$$\begin{equation}{\rho }_{\tau }={{\Lambda}}_{\tau }[{\rho }_{0}]={\mathbb{E}}_{p(\xi )}\!\left[{\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]=\int \mathrm{d}\xi \,{p}_{{\mu }_{\tau },{\sigma }_{\tau }}(\xi )\,{\text{e}}^{-\text{i}\xi \hat{H}\tau }{\rho }_{0}{\text{e}}^{\text{i}\xi \hat{H}\tau },\end{equation} \tag{ E.13 }$$

corresponds to the solution of the master equation:

$$\begin{equation}\frac{\mathrm{d}{\rho }_{\tau }}{\mathrm{d}\tau }=-\text{i}\omega (\tau )[\hat{H},{\rho }_{\tau }]+{\Gamma}(\tau )\left(\hat{H}{\rho }_{\tau }\hat{H}-\frac{1}{2}\left\{{\hat{H}}^{2},{\rho }_{\tau }\right\}\right)\end{equation} \tag{ E.14 }$$

$$\begin{equation}=-\text{i}\omega (\tau )[\hat{H},{\rho }_{\tau }]-\frac{1}{2}{\Gamma}(\tau )\left[\hat{H},[\hat{H},{\rho }_{\tau }]\right],\end{equation} \tag{ E.15 }$$

where the effective time-dependent frequency and decay parameters read, respectively:

$$\begin{equation}\omega (\tau )={\mu }_{\tau }+\tau {\dot{\mu }}_{\tau }\quad \text{and}\quad {\Gamma}(\tau )=2{\sigma }_{\tau }^{2}\,\tau \left(1+\frac{\dot{{\sigma }_{\tau }}}{{\sigma }_{\tau }}\tau \right).\end{equation} \tag{ E.16 }$$

Proof. By explicitly differentiating ρ_τ defined in (E.13) with respect to τ, we obtain:

$$\begin{align}\hfill \frac{\mathrm{d}{\rho }_{\tau }}{\mathrm{d}\tau }& =\frac{\mathrm{d}}{\mathrm{d}\tau }{{\Lambda}}_{\tau }[{\rho }_{0}]=\frac{\mathrm{d}}{\mathrm{d}\tau }{\mathbb{E}}_{p(\xi )}\!\left[{\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]=\frac{\mathrm{d}}{\mathrm{d}\tau }\left(\int \mathrm{d}\xi \,{p}_{{\mu }_{\tau },{\sigma }_{\tau }}(\xi )\ {\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right)\hfill \\ \hfill & =\int \mathrm{d}\xi \left(\frac{\mathrm{d}}{\mathrm{d}\tau }{p}_{{\mu }_{\tau },{\sigma }_{\tau }}(\xi )\right){\mathcal{U}}_{\xi \tau }[{\rho }_{0}]+\int \mathrm{d}\xi \,{p}_{{\mu }_{\tau },{\sigma }_{\tau }}(\xi )\frac{\mathrm{d}}{\mathrm{d}\tau }{\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\hfill \\ \hfill & =\frac{1}{{\sigma }_{\tau }^{3}}\,{\mathbb{E}}_{p(\xi )}\!\left[\left((\xi -{\mu }_{\tau }){\sigma }_{\tau }{\dot{\mu }}_{\tau }+({(\xi -{\mu }_{\tau })}^{2}-{\sigma }_{\tau }^{2}){\dot{\sigma }}_{\tau }\right){\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]-\text{i}[\hat{H},{\mathbb{E}}_{p(\xi )}\!\left[\xi {\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]]\hfill \\ \hfill & =\frac{{\dot{\mu }}_{\tau }}{{\sigma }_{\tau }^{2}}\,{\mathbb{E}}_{p(\xi )}\!\left[\xi {\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]+\frac{{\dot{\sigma }}_{\tau }}{{\sigma }_{\tau }^{3}}\,{\mathbb{E}}_{p(\xi )}\!\left[{(\xi -{\mu }_{\tau })}^{2}{\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]-\left(\frac{{\mu }_{\tau }{\dot{\mu }}_{\tau }}{{\sigma }_{\tau }^{2}}+\frac{\dot{{\sigma }_{\tau }}}{{\sigma }_{\tau }}\right){\rho }_{\tau }\hfill \\ \hfill & \quad -\text{i}[\hat{H},{\mathbb{E}}_{p(\xi )}\!\left[\xi {\mathcal{U}}_{\xi \tau }[{\rho }_{0}]\right]],\hfill \end{align} \tag{ E.17 }$$

where by using relations for the moments of a Gaussian distribution we can further simplify the following expressions:

$$\begin{equation}{\mathbb{E}}_{p(\xi )}\!\left[\xi \,{\mathcal{U}}_{\xi \tau }\!\left[{\rho }_{0}\right]\right]=-\mathrm{i}{\sigma }_{\tau }^{2}\tau \,[\hat{H},{\rho }_{\tau }]+{\mu }_{\tau }{\rho }_{\tau },\end{equation} \tag{ E.18 }$$

$$\begin{equation}{\mathbb{E}}_{p(\xi )}\!\left[{(\xi -{\mu }_{\tau })}^{2}{\mathcal{U}}_{\xi \tau }\!\left[{\rho }_{0}\right]\right]=-{\sigma }_{\tau }^{4}{\tau }^{2}\left[\hat{H},[\hat{H},{\rho }_{\tau }]\right]+{\sigma }_{\tau }^{2}{\rho }_{\tau }.\end{equation} \tag{ E.19 }$$

As a result, we can write the dynamics (E.17) as

$$\begin{equation}\frac{\mathrm{d}{\rho }_{\tau }}{\mathrm{d}\tau }=-\text{i}({\mu }_{\tau }+{\dot{\mu }}_{\tau }\tau )\,[\hat{H},{\rho }_{\tau }]-{\sigma }_{\tau }^{2}\tau \left(1+\frac{{\dot{\sigma }}_{\tau }}{{\sigma }_{\tau }}\tau \right)\left[\hat{H},[\hat{H},{\rho }_{\tau }]\right],\end{equation} \tag{ E.20 }$$

which is the desired form stated above in equation (E.15). □

Corollary. For the case of N spin-1/2 particles evolving for a time τ according to the dynamics (E.15) with $\hat{H}={\hat{J}}_{y}$, ω(τ) = γB and Γ(τ) = γ_y , i.e.:

$$\begin{equation}\frac{\mathrm{d}{\rho }_{\tau }}{\mathrm{d}\tau }=-\text{i}{\gamma}B[{\hat{J}}_{y},{\rho }_{\tau }]-\frac{1}{2}{\gamma }_{y}[{\hat{J}}_{y},[{\hat{J}}_{y},{\rho }_{\tau }]]\end{equation} \tag{ E.21 }$$

with B being constant over the time τ, the effective quantum map describing the evolution Λ_τ in (E.13) can be interpreted as mixture of unitary channels with the Gaussian mixing probability (E.12) of mean μ_τ = γB and standard deviation ${\sigma }_{\tau }=\sqrt{\frac{{\gamma }_{y}}{\tau }}$.

The above statement can be straightforwardly verified by explicitly computing the effective time-dependent frequency and decay parameters in (E.16) for the choice of μ_τ = γB and ${\sigma }_{\tau }=\sqrt{\frac{{\gamma }_{y}}{\tau }}$, which consistently simplify to:

$$\begin{equation}\omega (\tau )={\gamma}B+\tau \,0={\gamma}B,\end{equation} \tag{ E.22 }$$

$$\begin{equation}{\Gamma}(\tau )=2\frac{{\gamma }_{y}}{\tau }\tau \left(1+\frac{1}{2}{\left(\frac{{\gamma }_{y}}{\tau }\right)}^{-1/2}\left(-\frac{{\gamma }_{y}}{{\tau }^{2}}\right){\left(\frac{{\gamma }_{y}}{\tau }\right)}^{-1/2}\tau \right)=2{\gamma }_{y}\left(1-\frac{1}{2}\right)={\gamma }_{y}.\end{equation} \tag{ E.23 }$$

As discussed in the main text, upper-bounding the measurement-record contribution to the BI, i.e. J_M defined in equation (44), corresponds effectively—see inequality (59)—to computing the FI:

$$\begin{equation}F\!\left[{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right]=\int \!\mathrm{d}{\boldsymbol{\omega }}_{k}\,{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\left[-{\partial }_{{B}_{k}}^{2}\,\mathrm{log}\left({\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right)\right]\end{equation} \tag{ E.24 }$$

for the probability distribution defined in equation (55):

$$\begin{equation}{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})=\frac{1}{p({B}_{k})}\int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k}).\end{equation} \tag{ E.25 }$$

Note that ${\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})$ contains a Gaussian marginal distribution p(B_k ) specified in (E.7), as well as products of Gaussian distributions:

$$\begin{equation}p({\boldsymbol{B}}_{k})=\prod\limits _{j=1}^{k}p({B}_{j}\vert {B}_{j-1})p({B}_{0})\quad \text{and}\quad q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})=\prod\limits _{j=0}^{k}q({\omega }_{j}\vert {B}_{j}),\end{equation} \tag{ E.26 }$$

where p(B_j |B_j−1) is the Gaussian transition probability of the OU process defined in (E.5) with variance V_P given by (E.6), p(B₀) is the prior Gaussian distribution with zero mean and variance ${\sigma }_{0}^{2}$, while each q(ω_j |B_j ) is the mixing probability introduced within the CS method in equation (50), which is also a Gaussian with mean γB_j and variance

$$\begin{equation}{V}_{Q}=\frac{{\gamma }_{y}}{\delta t}.\end{equation} \tag{ E.27 }$$

It follows from the definitions in (E.26) that the integral in (E.25) can be rewritten as a set of nested integrals, i.e.:

$$\begin{align}\hfill \int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})& =\int \!\mathrm{d}{\boldsymbol{B}}_{k-1}\prod\limits _{j=1}^{k}p({B}_{j}\vert {B}_{j-1})q({\omega }_{j}\vert {B}_{j})p({B}_{0})q({\omega }_{0}\vert {B}_{0})\hfill \\ \hfill & =q({\omega }_{k}\vert {B}_{k})\int \mathrm{d}{B}_{k-1}\,p({B}_{k}\vert {B}_{k-1})q({\omega }_{k-1}\vert {B}_{k-1})\dots \int \mathrm{d}{B}_{1}\,p({B}_{2}\vert {B}_{1})q({\omega }_{1}\vert {B}_{1})\hfill \\ \hfill & \quad \times \int \mathrm{d}{B}_{0}\,p({B}_{1}\vert {B}_{0})q({\omega }_{0}\vert {B}_{0})p({B}_{0}),\hfill \end{align} \tag{ E.28 }$$

which we would like to simplify. To do so, we first need to prove the following lemma:

Lemma 2. Let us consider a recurrence relation between ${\mathcal{P}}_{j}({B}_{j})$ and ${\mathcal{P}}_{j-1}({B}_{j-1})$ for j = 0, 1, 2, ... as a generalised convolution of Gaussian distributions:

$$\begin{equation}{\mathcal{P}}_{j}({B}_{j})=\int \!\mathrm{d}{B}_{j-1}\,\frac{1}{\sqrt{2\pi {V}_{P}}}{\text{e}}^{-\frac{{({B}_{j}-{B}_{j-1})}^{2}}{2{V}_{P}}}\frac{1}{\sqrt{2\pi {V}_{Q}}}{\text{e}}^{-\frac{{({\omega }_{j-1}-{\gamma}{B}_{j-1})}^{2}}{2{V}_{Q}}}{\mathcal{P}}_{j-1}({B}_{j-1}),\end{equation} \tag{ E.29 }$$

where V_P , V_Q ⩾ 0 and ${\mathcal{P}}_{0}({B}_{0})={C}_{0}{\text{e}}^{-\frac{{({B}_{0}-{\mu }_{0})}^{2}}{2{V}_{0}}}$ with some fixed C₀, V₀ ⩾ 0 and ${\mu }_{0}\in \mathbb{R}$.

Then, for all j ⩾ 1:

$$\begin{equation}{\mathcal{P}}_{j}({B}_{j})={C}_{j}{\text{e}}^{-\frac{{({B}_{j}-{\mu }_{j})}^{2}}{2{V}_{j}}}\end{equation} \tag{ E.30 }$$

where the parameters C_j , μ_j and V_j are given as the solution to the following (coupled) recurrence relations:

$$\begin{equation}{C}_{j}={C}_{j-1}{\left(2\pi \left({{\gamma}}^{2}{V}_{P}+{V}_{Q}+\frac{{V}_{P}{V}_{Q}}{{V}_{j-1}}\right)\right)}^{-1/2}{\text{e}}^{-\frac{{({\omega }_{j-1}-{\gamma}{\mu }_{j-1})}^{2}}{2({V}_{Q}+{{\gamma}}^{2}{V}_{j-1})}}\end{equation} \tag{ E.31 }$$

$$\begin{equation}{\mu }_{j}=\frac{{V}_{Q}{\mu }_{j-1}+{V}_{j-1}{\gamma}\,{\omega }_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}\end{equation} \tag{ E.32 }$$

$$\begin{equation}{V}_{j}={V}_{P}+\frac{{V}_{Q}{V}_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}.\end{equation} \tag{ E.33 }$$

Proof. As for any recurrence problem, it is sufficient to prove that the solution (E.30) holds for j = 0, and that the recurrence relation (E.29) is fulfilled for any j ⩾ 1. The first part is trivially satisfied by definition, while to prove the latter we substitute P_j−1(B_j−1), defined according to (E.30), into (E.29) and explicitly perform the integration, i.e.:

$$\begin{align}\hfill {\mathcal{P}}_{j}({B}_{j})& =\int \!\mathrm{d}{B}_{j-1}\,\frac{1}{\sqrt{2\pi {V}_{P}}}{\text{e}}^{-\frac{{({B}_{j}-{B}_{j-1})}^{2}}{2{V}_{P}}}\frac{1}{\sqrt{2\pi {V}_{Q}}}{\text{e}}^{-\frac{{({\omega }_{j-1}-{\gamma}{B}_{j-1})}^{2}}{2{V}_{Q}}}{C}_{j-1}{\text{e}}^{-\frac{{({B}_{j-1}-{\mu }_{j-1})}^{2}}{2{V}_{j-1}}}\hfill \\ \hfill & =\frac{{C}_{j-1}}{\sqrt{2\pi \left({V}_{Q}+{{\gamma}}^{2}{V}_{P}+\frac{{V}_{P}{V}_{Q}}{{V}_{j-1}}\right)}}\mathrm{exp}\left\{-\frac{{B}_{j}^{2}-2\alpha {B}_{j}+\beta }{2\left({V}_{P}+\frac{{V}_{Q}{V}_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}\right)}\right\},\hfill \end{align} \tag{ E.34 }$$

where α and β are constants independent of B_j and B_j−1, and equal to

$$\begin{equation}\alpha =\frac{{V}_{Q}{\mu }_{j-1}+{V}_{j-1}\,{\gamma}\,{\omega }_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}},\end{equation} \tag{ E.35 }$$

$$\begin{equation}\beta =\frac{{V}_{Q}{\mu }_{j-1}^{2}+{V}_{j-1}{\omega }_{j-1}^{2}+{V}_{P}{({\omega }_{j-1}-{\gamma}{\mu }_{j-1})}^{2}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}.\end{equation} \tag{ E.36 }$$

Hence, by 'completing the square' we rewrite (E.34) as

$$\begin{equation}{\mathcal{P}}_{j}({B}_{j})=\frac{{C}_{j-1}{\text{e}}^{-\frac{-{\alpha }^{2}+\beta }{2\left({V}_{P}+\frac{{V}_{Q}{V}_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}\right)}}}{\sqrt{2\pi \left({V}_{Q}+{{\gamma}}^{2}{V}_{P}+\frac{{V}_{P}{V}_{Q}}{{V}_{j-1}}\right)}}\mathrm{exp}\left\{-\frac{{({B}_{j}-\alpha )}^{2}}{2\left({V}_{P}+\frac{{V}_{Q}{V}_{j-1}}{{V}_{Q}+{{\gamma}}^{2}{V}_{j-1}}\right)}\right\},\end{equation} \tag{ E.37 }$$

and, after substituting for α and β from (E.35) and (E.36), we arrive at the expression (E.30) for P_j (B_j ) with C_j , μ_j and V_j specified by the recurrence relations (E.31)–(E.33). □

Now, using the above lemma we may rewrite equation (E.28) as

$$\begin{equation}\int \mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})\,q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})=q({\omega }_{k}\vert {B}_{k}){\mathcal{P}}_{k}({B}_{k}),\end{equation} \tag{ E.38 }$$

with ${\mathcal{P}}_{k}({B}_{k})$ being now defined according to equation (E.30) with variances V_P and V_Q in (E.31)–(E.33) specified in our case by equations (E.6) and (E.27), respectively. Consequently, the FI of ${\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})$ in equation (E.24) reads

$$\begin{align}\hfill F\!\left[{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\right]& =\int \!\mathrm{d}{\boldsymbol{\omega }}_{k}\,{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\left[-{\partial }_{{B}_{k}}^{2}\,\mathrm{log}\left(\frac{1}{p({B}_{k})}\int \mathrm{d}{\boldsymbol{B}}_{k-1}\,p({\boldsymbol{B}}_{k})q({\boldsymbol{\omega }}_{k}\vert {\boldsymbol{B}}_{k})\right)\right]\hfill \\ \hfill & =\int \!\mathrm{d}{\boldsymbol{\omega }}_{k}{\mathbb{P}}_{{B}_{k}}({\boldsymbol{\omega }}_{k})\left[-{\partial }_{{B}_{k}}^{2}\,\mathrm{log}\left(\frac{q({\omega }_{k}\vert {B}_{k})}{p({B}_{k})}{\mathcal{P}}_{k}({B}_{k})\right)\right]\hfill \end{align} \tag{ E.39 }$$

$$\begin{equation}=\frac{{{\gamma}}^{2}}{{V}_{Q}}-\frac{1}{{V}_{P}^{(k)}}+\frac{1}{{V}_{k}},\end{equation} \tag{ E.40 }$$

where the expression (E.40) follows from the fact that we are dealing with a product (and quotient) of Gaussian distributions within log(...) in (E.39), so that the FI becomes just the sum (and difference) of the inverses of their respective variances. In particular, V_Q /γ² is the variance of q(ω_k |B_k ) when treating B_k as the random variable, ${V}_{P}^{(k)}$ is the variance of p(B_k ) specified in (E.8); while V_k is the variance of ${\mathcal{P}}_{k}({B}_{k})$ given by the recurrence relation (E.33) that must still be solved.

Although the recursive relation (E.33) (with ${V}_{0}={\sigma }_{0}^{2}$) admits a general solution, for simplicity we present only its form for the relevant setting, in which we do not possess any prior knowledge about the initial field B₀, i.e. when σ₀ → ∞ in (E.4). Then,

$$\begin{equation}{\left.{V}_{k}\right\vert }_{{\sigma }_{0}\to \infty }=\frac{1}{2{V}_{P}}+\frac{1}{2{\gamma}}\sqrt{{V}_{P}(4{V}_{Q}+{V}_{P}{{\gamma}}^{2})}\left(1+\frac{2}{-1+{\left(\frac{2{V}_{Q}+{\gamma}\left({V}_{P}{\gamma}+\sqrt{{V}_{P}(4{V}_{Q}+{V}_{P}{{\gamma}}^{2})}\right)}{2{V}_{Q}+{\gamma}\left({V}_{P}{\gamma}-\sqrt{{V}_{P}(4{V}_{Q}+{V}_{P}{{\gamma}}^{2})}\right)}\right)}^{k}}\right),\end{equation} \tag{ E.41 }$$

which—after substituting for ${V}_{P}=\frac{{q}_{B}}{2\chi }(1-{\text{e}}^{-2\chi \delta t})$ and V_Q = γ_y /δt according to expressions (E.6) and (E.27), respectively, and taking the continuous-time limit δt → 0 with k = t/δt—takes the form:

$$\begin{equation}\underset{\delta t\to 0}{\mathrm{lim}}\left\{{\left.{V}_{k=t/\delta t}\right\vert }_{{\sigma }_{0}\to \infty }\right\}=\sqrt{\frac{{\gamma }_{y}\,{q}_{B}}{{{\gamma}}^{2}}}\,\mathrm{coth}\!\left(t\sqrt{\frac{{q}_{B}{{\gamma}}^{2}}{{\gamma }_{y}}}\right).\end{equation} \tag{ E.42 }$$

Finally, by noting that the term γ²/V_Q = γ² δt/γ_y in (E.40) vanishes when letting δt → 0, and so does the term $1/{V}_{P}^{(k)}$ when σ₀ → ∞ (see equation (E.8)), we arrive at the FI of ${\mathbb{P}}_{{B}_{t}}({\boldsymbol{\omega }}_{< t})$ defined within the continuous-time δt → 0 limit as

$$\begin{equation}F\!\left[{\mathbb{P}}_{{B}_{t}}({\boldsymbol{\omega }}_{< t})\right]\,\underset{{\sigma }_{0}\to \infty }{=}\,\frac{1}{{\mathrm{lim}}_{\delta t\to 0}\left\{{\left.{V}_{k=t/\delta t}\right\vert }_{{\sigma }_{0}\to \infty }\right\}}=\sqrt{\frac{{{\gamma}}^{2}}{{\gamma }_{y}{q}_{B}}}\,\mathrm{tanh}\!\left(t\sqrt{\frac{{q}_{B}{{\gamma}}^{2}}{{\gamma }_{y}}}\right),\end{equation} \tag{ E.43 }$$

which is the expression stated in equation (57) of the main text.