Estimation of classical parameters via continuous probing of complementary quantum observables

In our formulation of the estimation problem we will treat the classical parameter $\vec{\gamma}$ as unknown with an assigned probability distribution $P(\vec {\gamma })$. The measurements then cause an update of the probability distribution, which is governed by Bayes rule for conditional probabilities. Until the actual value of $\vec {\gamma }$ is known, we thus have to treat each candidate value with a probability factor, and for each possible value of $\vec {\gamma }$, the corresponding state of the quantum system evolves under the quantum filtering equation with the corresponding dependence of $\vec {\gamma }$.

To this end, it is convenient [7, 9–11, 18] to consider the augmented Hilbert space $\mathfrak {H}=\mathfrak {H}_{\mathrm {s}}\otimes \mathfrak {H}_{\gamma }$, where $\mathfrak {H}_{\mathrm {s}}$ and $\mathfrak {H}_{\gamma }$ refer to the quantum system Hilbert space and the space of classical states for the variable $\vec{\gamma}$, respectively. The latter space describes states with definite values of the parameters, and superposition states are not populated. The quantum mechanical notation, however, still applies and, e.g. describes a probability distribution for a set of values $\vec {\gamma }_k$ as a diagonal density matrix $\hat \varrho _{\gamma }=\sum _k P_k |{\vec{\gamma} _k}\rangle \langle {\vec{\gamma} _k}|$. The classical variables $\vec {\gamma }$ are equivalent to QND variables of an ancillary quantum system that interacts with our probe system, and for which the Bayesian probability update is fully equivalent to the quantum measurement back-action. When we incorporate the parameters $\vec {\gamma }$ in this way we can directly apply the filtering equation on the augmented space. We note that the equivalence of classical unknown parameters and QND variables of ancillary quantum systems goes both ways. This permits certain types of quantum evolution to be described exactly by classical probability distributions, and has recently been applied, e.g. in the context of decoherence of two-level quantum systems driven by coherent laser light [19–21].

The observer who does not know the value of $\vec {\gamma }$ describes the combined quantum and classical system by the augmented density matrix

$$\begin{eqnarray} &&{\hat{\varrho}} = \sum_{k=1}^N P_k |{\vec{\gamma}_k}\rangle\langle{\vec{\gamma}_k}|\otimes{\hat{\varrho}}_k^{\mathrm{s}}, \end{eqnarray} \tag{ 4 }$$

where $\hat \varrho ^{\mathrm {s}}_k$ is the normalized system state density matrix associated with the specific value $\vec {\gamma }=\vec {\gamma }_k$.

The combined system evolves according to the quantum filter equation

$$\begin{eqnarray}&&\fl {\rm d} {\hat{\varrho}} = \left(\frac{\mathrm{i}}{\hbar} [{\hat{\varrho}},\skew3\hat H] + \sum_j\Gamma_j \mathcal{D}[\hat{\cal O}_j]{\hat{\varrho}} + \sum_{n=1}^{M}M_n \mathcal{D}[\hat{\cal M}_n]{\hat{\varrho}}\right){\rm d} t \nonumber\\ &&+ \sum_{n=1}^{M}\!\sqrt{M_n}\mathcal{H}[\hat{\cal M}_n]{\hat{\varrho}} \left(\mathrm{d}Y_n^{\mathrm{D}}(t)-\eta_n\sqrt{M_n}\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}}\,\mathrm{d}t \right), \end{eqnarray} \tag{ 5 }$$

where the operator $\mathcal {H}$ is defined as

$$\begin{eqnarray} \mathcal{H}[\skew4\hat{f}]{\hat{\varrho}} &= \skew4\hat{f} {\hat{\varrho}} + {\hat{\varrho}} \skew4\hat{f}^{\dagger} - \langle{\skew4\hat{f} + \skew4\hat{f}^{\dagger}}\rangle_{\mathrm{E}} {\hat{\varrho}}. \end{eqnarray} \tag{ 6 }$$

Here we use the notation $\langle {\hat {f}}\rangle _{\mathrm {E}}=\mathrm {Tr}(\hat {f}\hat {\varrho })=\sum _{k=1}^NP_k\mathrm {Tr}_{\mathrm {S}}(\skew4\hat{f}\hat \varrho _k^{\mathrm {s}})$ to explicitly recall that the expectation value of the signal should be determined by the full augmented quantum state, equivalent to a weighted average over the ensemble of states of the quantum system governed by the different values of $\vec {\gamma }$.

If the dependence on $\vec {\gamma }$ only enters via the Hamiltonian $\hat {H}(\vec {\gamma })$, the different terms in equation (5) are implemented as the following product operators on $\mathfrak {H}=\mathfrak {H}_{\mathrm {s}} \otimes \mathfrak {H}_{\gamma }$:

The quantum filter equation (5) is established as an augmented version of equation (1) to efficiently include the probabilistic representation of the unknown variables and their correlations with the quantum system. As the main objective of our analysis is the estimation of the parameter $\vec {\gamma }$, we shall in the following section reformulate equation (5) to explicitly extract the conditioned dynamics of the probability distribution P_k in equation (4).

The stochastic process appearing in equation (5),

$$\begin{eqnarray} &&\sqrt{\eta_n}\mathrm{d}V_n(t):= \mathrm{d}Y_n^{\mathrm{D}}(t)-\eta_n\sqrt{M_n}\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}}\,\mathrm{d}t. \end{eqnarray} \tag{ 8 }$$

is not a standard Wiener process. This is because we subtract from the measured signal a weighted average based on our probabilistic description of $\vec {\gamma }$, while the measured photocurrent dY^D_n in the experiment is governed by the actual value of the unknown parameter $\vec {\gamma }=\vec {\gamma }_{k_0}$.

Equation (2) characterizes the properties of such a realistic detection record, and when inserted in equation (5), we find that the stochastic process in equation (8) can be rewritten as

$$\begin{eqnarray} &&\mathrm{d}V_n \!= \!\mathrm{d}W_n - \sqrt{\eta_n M_n} [ \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}}-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{0} ] \mathrm{d}t. \end{eqnarray} \tag{ 9 }$$

Hence, dV_n(t) is a stochastic Gaussian process with variance dt, and with (statistical) mean value $\langle \langle \mathrm {d}V_n(t)\rangle \rangle =\sqrt {\eta _n M_n}(\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _{0}-\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _{\mathrm {E}})\mathrm {d}t$, reflecting precisely the difference in the expectation value assumed by the weighted average over different $\vec {\gamma }_k$ and by the correct value $\vec {\gamma }_{k_0}$.

As shown in detail in appendix A.3, equation (5) leads to two separate equations for $\hat {\varrho }_k^{\mathrm {s}}$ and P_k that permit a full simulation of the detection process

$$\begin{eqnarray}&&\fl {\rm d} {\hat{\varrho}}_{k}^{\mathrm{s}} = \frac{\mathrm{i}}{\hbar} [{\hat{\varrho}}_{k}^{\mathrm{s}},\skew3\hat H(\vec{\gamma}_{k})]{\rm d} t + \sum_j\Gamma_j \mathcal{D}[\hat{\cal O}_j]{\hat{\varrho}}_{k}^{\mathrm{s}}\,{\rm d} t\nonumber\\ &&+ \sum_{n=1}^{M}M_n \mathcal{D}[\hat{\cal M}_n]{\hat{\varrho}}_{k}^{\mathrm{s}}\,{\rm d} t +\sqrt{M_n} \mathcal{H}[\hat{\mathcal{M}}_n]{\hat{\varrho}}_{k}^{\mathrm{s}} \,\mathrm{d}\tilde{V}^{(k)}_n(t) \end{eqnarray} \tag{ 10 }$$

with ${\rm d} \skew3\tilde {V}^{(k)}_n(t)=\mathrm {d}Y_n^{\mathrm {D}}(t)-\eta _n\sqrt {M_n}\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _k\,\mathrm {d}t$ and

$$\begin{equation} \fl \mathrm{d}P_k=\mathrm{Tr}_{\mathrm{S}}({\rm d} \hat\rho_k^{\mathrm{s}}) =P_k\sum_{n=1}^{M}\sqrt{\eta_nM_n}[ \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} ]\mathrm{d}V_n(t). \end{equation} \tag{ 11 }$$

In both equations the photocurrent dY^D_n(t), observed or simulated according to equation (1), appears, and equation (11) agrees with the expression given in [10]. We note, however, that the stochastic term in our equation (10) contains the expectation value $\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _k$ corresponding to the parameter value $\vec {\gamma }_k$, while in [10] the ensemble average $\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _{\mathrm {E}}$ has been used.

By inserting the expression (9) for dV_n in equation (11), we see that the change of dP_k due to the measurements is given by

$$\begin{equation}\fl ( \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} ) [ \mathrm{d}W_n(t)+\sqrt{\eta_nM_n}( \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{0}-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} )\mathrm{d}t ]. \end{equation} \tag{ 12 }$$

This equation has a natural interpretation: for parameter values $\vec{\gamma} _k$ where the expected mean current $\propto \langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _k$ differs in the same (opposite) direction from the ensemble mean as the one expected for the actual value $\langle {\skew3\hat {\cal {M}}_n+\skew3\hat {\cal {M}}_n^{\dagger }}\rangle _{0}$, the two parentheses will typically have the same (opposite) sign, and P_k will increase (decrease). Due to the random contribution dW_n(t), however, the probabilities will show fluctuations, and their increase (decrease) with time will appear only as an average trend, leading, in particular, to a typically increasing value for the probability of the correct value P_k0.

Following [10], we have first investigated the case in which the initial state of the atom is represented by the Bloch vector $\vec {r}_0=(0,1,0)$ with M_1,2 = 0 but M₃ ≠ 0 (i.e. we probe only the component of the spin along z and M ≡ M₃), and where the magnetic field has a known strength while it has equal prior probabilities to point in the positive and in the negative x-axis directions. In figure 1 (left) we show the behavior of 〈〈cos θ〉〉(t) as function of time in the two cases of a weak $\vert \vec {b}_{\mathrm {u}}\vert = 0.1$ (lower, black curve) and a strong $\vert \vec {b}_{\mathrm {u}}\vert = 1.5$ field (upper, red curve). The curves are obtained by averaging over 104 simulated detection records. The stronger the field amplitude, the better is the directional estimate, but, as also observed in [10], the quantum filter does not unambiguously identify the direction $\vec {b}_{\mathrm {u}}$ of the unknown field.

Zoom In Zoom Out Reset image size

We emphasize that the inability to determine the direction of the field unambiguously is not due to the method of analysis of the simulated measurement data. The filter equation, indeed, provides a full Bayesian analysis, but the method of detection restricts the information content of the measurement data and hence the sensitivity of the quantum system to the Hamiltonian unknown parameter(s). For example, the measured z-component of the spin is sensitive to a precession around the y-axis, but only the first transient evolution of the Bloch vector away from the equatorial plane yields information about the sense of rotation. This is reflected by the curves in figure 1 (left), which reach constant levels after finite time, and hereafter the detection scheme cannot teach us more about the direction of the magnetic field. In the case of the weak field, we learn only very little, because the measurement processes forces the z-component of the spin to acquire definite values on a time scale too short for the result to be sensitive to the small rotation angle toward the positive or negative z-axis caused by the magnetic field. Hereafter, further evolution of the z-component may be suppressed by a quantum Zeno-like effect, and the sense of rotation away from the z-axis is anyway not discernible by measurements of σ_z, only.

This situation changes if the experiment is changed and all the three spin components are detected. As illustrated in figure 1 (right), the accumulation of results from all three detectors causes an effective update of the probability distribution that proceeds much longer in time than when only one detector is applied. For high strength of the field we obtain a probability distribution which is well converged to the correct $\vec {b}_{\mathrm {u}}$ within the measurement time in the figure, and for the weaker field, we show a more slowly, but steadily, growing probability of the correct value. In these numerical experiments 104 trials have been performed in order to accumulate sufficient data for the statistical averages within a reasonable computational time. A small fraction (∼ 1%) of the simulated trajectories are unstable in the case of the three detectors and they have been rejected from the statistical averages. The probing of several components on the one hand allows discrimination between left and right circular precession, and, on the other hand prevents (Zeno-) locking of the system to eigenstates of the probed quantities as they do not commute and do not have common eigenstates.

Finally, we have analyzed the scenario in which the magnetic field has an isotropic prior probability distribution, represented by an ensemble $\mathcal {V}_{B}$ with N = 98 directions on the unit sphere. Here, we assume the stronger field $\vert \vec {b}_{\mathrm {u}}\vert = 1.5$, and we assume equal strength probing of all three spin components.

The initial condition for the probabilities with all P_k = 1/N is illustrated by the (green) sphere displayed in figure 2 at time Mt = 0. We study the convergence of the probability distribution as a function of the probing time, and for long times (MT = 25), the filter converges steadily toward a single direction. We note that the spheres are cut on the top because in the simulation we have considered an interval θ = (0,π) equally spaced with N_θ = 7 grid points and an interval ϕ = [0,2π) with N_ϕ = 14 grid points for the azimuthal angle.

Zoom In Zoom Out Reset image size

The equal strength probing of the three cartesian spin components σ_x, σ_y and σ_z is equivalent [10, 22] to probing of any other cartesian set, including for example one parallel component and two perpendicular components to the applied magnetic field, and we find that the monitoring of all three spin components lead to unambiguous identification of the direction of the applied field.

The spin vector may, inadvertently, align, parallel or anti-parallel with the applied magnetic field axis, but the random back-action of the probing of the non-commuting orthogonal spin observables will kick the system away from these directions and give rise to renewed observable precession.

We consider an atomic quantum system with degenerate ground states |g_m〉 and excited states |e_m〉, where m = ± 1/2 denotes the azimuthal quantum number with respect to the quantization z-axis. The atom interacts through its magnetic moment $\hat {\vec {\mu }}$ with a classical magnetic field $\skew3\vec {B}=(B_x,B_y,B_z)$, which drives a Larmor precession of the ground state atomic spin. The atom is probed by an off-resonant laser field, which is linearly polarized along the x-axis. Due to the dipole selection rules, the two circularly polarized field components with annihilation operators $\hat a_{\pm } = (\mp \hat a_x + \mathrm {i} \hat a_y)/\sqrt {2}$ couple individually to the two different ground state populations.

During off-resonant probing (i.e. Δ ≫ g), the atomic excited state can be adiabatically eliminated, and the quantized field–atom Hamiltonian is given by [27]

$$\begin{eqnarray} &&\skew3\hat H = \frac{\hbar g^2}{\Delta}\sum_{\ell = \pm 1} \hat a_{\ell}^{\dagger} \hat a_{\ell}|{g_{\ell/2}}\rangle\langle{g_{\ell/2}}| + \hat{\vec{\mu}}\cdot\vec{B}, \end{eqnarray} \tag{ A.1 }$$

where Δ = ω_L − ω_A is the laser atom detuning, $g = \vec {d}\cdot \vec {E}_0/\hbar$, where $\vec {d}$ is the atomic electric dipole moment and $|\vec {E}_0| = \sqrt {\hbar \omega _{\mathrm {L}}/(V\epsilon _0)}$ is the electric field per photon in the quantization volume V and ₀ is the (electric) vacuum permeability. As result of equation (A.1), the two circularly polarized field components experience phase shifts that depend on the atomic occupation of the two ground states. This implies a (Faraday) rotation of the field polarization, which is proportional to the population difference between the ground states |g_±1/2〉.

Because of the strong linearly polarized probe field with photon number N_ph ≫ 1, the Stokes operators of the field can be written as

$$\begin{eqnarray} \skew3\hat J_x &=& \frac{\hat a_x^{\dagger} \hat a_x - \hat a_y^{\dagger} \hat a_y}{2} \approx \frac{N_{\mathrm{ph}}}{2},\nonumber\\ \skew3\hat J_y &=& \frac{\hat a_x^{\dagger} \hat a_y + \hat a_y^{\dagger} \hat a_x}{2} \simeq \frac{\sqrt{N_{\mathrm{ph}}}}{2} (\hat a_y + \hat a_y^{\dagger}) = \sqrt{N_{\mathrm{ph}}}\,\hat y,\\ \skew3\hat J_z &=& \frac{\hat a_x^{\dagger} \hat a_y - \hat a_y^{\dagger} \hat a_x}{2\mathrm{i}} \simeq \frac{\sqrt{N_{\mathrm{ph}}}}{2\mathrm{i}} (\hat a_y - \hat a_y^{\dagger}) = \sqrt{N_{\mathrm{ph}}}\,\hat p_y,\nonumber \end{eqnarray} \tag{ A.2 }$$

defining the canonical conjugate operators $\hat {y}$ and $\hat {p}_y$. The polarization rotation of the field is conveniently measured by subtracting the intensities of polarization components linearly polarized at ±45° with respect to the incident field, and when this quantity is expressed in terms of the Stokes observables we recover an expression proportional to the operator $\hat {y}$.

By writing $\hat {\vec {\mu }}=\mu \hat {\vec {\sigma }}$, where $\hat {\vec {\sigma }}=(\hat \sigma _x,\hat \sigma _y,\hat \sigma _z)$ are the Pauli matrices, and using the definitions in equation (A.2), the total Hamiltonian can be written

$$\begin{eqnarray} &&\skew3\hat H = \frac{\hbar g^2}{\Delta}\sqrt{N_{\mathrm{ph}}}\hat p_y\hat{\sigma}_z + \mu \sum_{\alpha = x,y,z}\hat\sigma_{\alpha} B_{\alpha}. \end{eqnarray} \tag{ A.3 }$$

A term $\frac {\hbar g^2}{2\Delta }(\hat a_x^{\dagger } \hat a_x + \hat a_y^{\dagger } \hat a_y)$ has been omitted from equation (A.3) as it gives rise to a common phase shift, but no polarization rotation of the probe laser field.

The continuous interaction between the atom and the light beam can then be represented as a sequence of interactions of the atom with one beam segment of length L = cτ after the other. Each segment, in turn, is described as a single harmonic oscillator mode described by the operators in equation (A.2). In the time interval τ the atom interacts with N_ph = Φτ photons, where Φ is the photon flux, and the dynamics is well described by the coarse-grained propagator

$$\begin{eqnarray} &&\hat U = \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\skew3\hat H_{\tau}} \simeq \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\kappa_{\tau}\hat{p}_y\,\hat{\sigma}_{\vec{m}}} \,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mu_{\tau}\vert\vec{B}\vert\hat\sigma_{\vec{n}_B}}, \end{eqnarray} \tag{ A.4 }$$

where

$$\begin{eqnarray} &&\kappa_{\tau} = \frac{\hbar g^2\tau}{\Delta}\sqrt{N_{\mathrm{ph}}}, \quad \mu_{\tau} = \mu \tau, \end{eqnarray} \tag{ A.5 }$$

and $\vec {n}_B\equiv (n_B^x,n_B^y,n_B^z)$ is the unit vector pointing in the magnetic field direction. We assume that with probe beams of appropriate polarization and propagation direction we can probe the atomic ground state spin projection $\hat {\sigma }_{\vec {m}} = \hat {\vec {\sigma }}\cdot \vec {m}$ along any direction $\vec {m}$, and we have introduced $\hat {\sigma }_{\vec {n}_B} = \hat {\vec {\sigma }}\cdot \vec {n}_B$. Since N_ph∝τ and g∝τ^−1/2 (through the volume V =Aτc), the dimensionless coupling constant κ_τ is proportional to τ^1/2, while μ_τ is linear in τ, and we have thus neglected terms of order τ^3/2 and higher in equation (A.4).

The incident laser beam is in a coherent state of linearly polarized light described by a Gaussian wave function π^−1/4 e^−p2_y/2 , associated with the observable $\hat {p}_y$ introduced in equation (A.2). The joint state of an atomic ground superposition state $|{\psi _{\mathrm {s}}(t)}\rangle =\sum _{\ell = \pm 1}c_{\ell /2}|{g_{\ell /2}}\rangle$ and the incident y-polarization component of the quantized probe field can thus be written in the product basis |p_y,g_ℓ/2〉 ≡ |p_y〉⊗|g_ℓ/2〉,

$$\begin{eqnarray} &&|{\Psi_{\mathrm{ps}}(t)}\rangle = \frac{1}{\pi^{1/4}}\sum_{\ell = \pm 1} c_{\ell/2} \int\mathrm{d}p_y \,\mathrm{e}^{-p^2_y/2} |{p_y,g_{\ell/2}}\rangle. \end{eqnarray} \tag{ A.6 }$$

Under the action of (A.4), this state evolves into the entangled state

$$\begin{eqnarray} |{\Psi_{\mathrm{ps}}(t+\tau)}\rangle &=& \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mu_{\tau}\vert\vec{B}\vert\hat\sigma_{\vec{n}_B}} \,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\kappa_{\tau}\hat{p}_y\,\hat{\sigma}_{\vec{m}}}|{\Psi_{\mathrm{ps}}(t)}\rangle\nonumber\\ &=&\frac{1}{\pi^{1/4}} \sum_{\ell = \pm 1} \int\mathrm{d}p_y \,c^{\prime}_{\ell/2}(p_y)\, \mathrm{e}^{-p^2_y/2 } |{p_y,g_{\ell/2}}\rangle, \end{eqnarray} \tag{ A.7 }$$

where, to first order in τ (second order in $\sqrt {\tau }$), the new expansion coefficients c'_ℓ/2(p_y) are

$$\begin{eqnarray}&&\fl c^{\prime}_{\ell/2} = c_{\ell/2}\!\!\left[ 1 - \frac{1}{2}\left(\frac{\kappa_{\tau}}{\hbar}\right)^2 \!\!p_y^2 - \mathrm{i}\ell\left( m_z\frac{\kappa_{\tau}}{\hbar} p_y + n_B^z \frac{\mu_{\tau}\vert\skew3\vec{B}\vert}{\hbar} \right) \right]\nonumber\\ &&-\mathrm{i}c_{-\ell/2}\left[ \frac{\kappa_{\tau}}{\hbar} p_y (m_x-\mathrm{i}\ell m_y) +\frac{\mu_{\tau}\vert\skew3\vec{B}\vert}{\hbar} (n_B^x-\mathrm{i}\ell n_B^y) \right]. \end{eqnarray} \tag{ A.8 }$$

To describe the back-action due to the field measurement, it is convenient to transform the entangled state (A.7) to the y rather than p_y representation of the field, using the relation

$$\begin{eqnarray} &&|{p_y}\rangle=\frac{1}{\sqrt{2\pi}}\int\mathrm{d}y \,\mathrm{e}^{-\mathrm{i}yp_y}|{y}\rangle. \end{eqnarray} \tag{ A.9 }$$

Hence, the state (A.7) can be rewritten as

$$\begin{eqnarray} &&|{\Psi_{\mathrm{ps}}(t+\tau)}\rangle = \frac{1}{\pi^{1/4}}\sum_{\ell = \pm 1} \int\mathrm{d}y \,\tilde{c}_{\ell/2}(y)\, \mathrm{e}^{-y^2/2 } |{y,g_{\ell/2}}\rangle \end{eqnarray} \tag{ A.10 }$$

with the new coefficient $\tilde{c}_{\ell /2}(y)$ given by

$$\begin{eqnarray}&&\fl \tilde{c}_{\ell/2} = \mathrm{i}c_{-\ell/2}\left[ \frac{\mu_{\tau}\vert\skew3\vec{B}\vert}{\hbar} (\mathrm{i}\ell n_B^y-n_B^x) +\frac{\kappa_{\tau}}{\hbar} y (\ell m_y+\mathrm{i}m_x) \right] \nonumber\\ &&+c_{\ell/2}\!\!\left[ 1 - \frac{1}{2}\left(\frac{\kappa_{\tau}}{\hbar}\right)^2 (1 - y^2) - \mathrm{i}\ell\!\left( \!n_B^z \frac{\mu_{\tau}\vert\skew3\vec{B}\vert}{\hbar} -\mathrm{i}m_z\frac{\kappa_{\tau}}{\hbar} y \right) \right]. \end{eqnarray} \tag{ A.11 }$$

A measurement of the light probe observable $\hat y$ with outcome y^D projects the state of the system (A.10) onto the state component with that definite value, i.e. the atomic part of the system becomes

$$\begin{eqnarray} &&|{\psi_{\mathrm s}(t+\tau)}\rangle = \frac{1}{\pi^{1/4}}\sum_{\ell = \pm 1}\,\skew3\tilde{c}_{\ell/2}(y^{\mathrm{D}}) \,\mathrm{e}^{-\frac{(y^{\mathrm{D}})^2}{2}} |{g_{\ell/2}}\rangle, \end{eqnarray} \tag{ A.12 }$$

where the explicit dependence of $\tilde{c}_{\ell /2}$ on the measurement outcome y^D causes an (infinitesimal) transfer of amplitude among the two atomic states.

Given the quantum state (A.10), the probability to measure a given value y^D is

$$\begin{eqnarray} &&\mathcal{P}(y^{\mathrm{D}}) = \frac{\mathrm{e}^{-(y^{\mathrm{D}})^2}}{\sqrt{\pi}}\sum_{\ell = \pm 1} \vert\tilde{c}_{\ell/2}(y^{\mathrm{D}})\vert^2 \simeq \pi^{-1/2}\, \mathrm{e}^{-(y^{\mathrm{D}}-y_0)^2}, \end{eqnarray} \tag{ A.13 }$$

where $y_0 = \frac {\kappa _{\tau }}{\hbar }\langle {\hat \sigma _{\vec {m}}}\rangle$. This explicitly shows how the optical probing yields a signal proportional to the desired mean value $\langle \hat {\sigma }_{\vec {m}}\rangle$, and it enables us to model the measurement outcome y^D as a stochastic variables

$$\begin{eqnarray} &&y^{\mathrm{D}} = \frac{\kappa_{\tau}}{\hbar}\langle{\hat\sigma_{\vec{m}}}\rangle + \frac{\Delta W}{\sqrt{2\tau}}, \end{eqnarray} \tag{ A.14 }$$

where ΔW is a (finite) Gaussian Wiener increment with mean zero and variance τ.

By replacing y with the expression (A.14) for y^D in equation (A.12), and by expanding the expressions for the state amplitudes to lowest order in τ, we obtain, in the continuous limit, the quantum filtering equation for the state of the atomic system

$$\begin{eqnarray} &&{\rm d} {\hat{\varrho}}_{\mathrm s} = -\mathrm{i}\frac{\mu_B}{\hbar} [\hat{\sigma}_{\vec{B}},{\hat{\varrho}}_{\mathrm s}] {\rm d} t + M \mathcal{D}[\hat\sigma_{\vec{m}}]{\hat{\varrho}}_{\mathrm s}\,{\rm d} t + \sqrt{M}\mathcal{H}[\hat\sigma_{\vec{m}}]{\hat{\varrho}}_{\mathrm s}\,{\rm d} W(t). \end{eqnarray} \tag{ A.15 }$$

Here the interaction parameter M is given explicitly by

$$\begin{eqnarray} &&M = \frac{g^4\tau^2}{4(\omega_{\mathrm{L}} - \omega_{\mathrm{A}})^2}\Phi, \end{eqnarray} \tag{ A.16 }$$

and an infinitesimal Wiener process models the noise in the detected signal, ${\rm d} Y^{\mathrm {D}}(t)=2\sqrt {M}\langle \hat {\sigma }_m\rangle \mathrm {d}t + \mathrm {d}W$.

We note that several probe fields may be used to simultaneously probe different spin components, and, to lowest order in τ, their effects on the quantum state commute. They may hence be included as separate terms with independent Wiener processes dW_n.

In order to derive equations (10) and (11) we first need the stochastic master equation for the (unnormalized) density operator $\hat \rho _k^{\mathrm {s}}$, which is given by

$$\begin{eqnarray}&&\fl \frac{{\rm d} \hat\rho_k^{\mathrm{s}}}{P_k} = \!\left(\frac{\mathrm{i}}{\hbar} [{\hat{\varrho}}_k^{\mathrm{s}},\skew3\hat H(\vec{\gamma}_k)] + \sum_j\Gamma_{\!\!j} \mathcal{D}[\hat{\cal O}_j]{\hat{\varrho}}_k^{\mathrm{s}} + \sum_{n=1}^{M}M_n \mathcal{D}[\hat{\cal M}_n]{\hat{\varrho}}_k^{\mathrm{s}}\right){\rm d} t \nonumber\\ &&+ \sum_{n=1}^{M}\!\sqrt{\eta_nM_n}\!\left(\hat{\mathcal{M}}_n{\hat{\varrho}}_k^{\mathrm{s}} +{\hat{\varrho}}_k^{\mathrm{s}}\hat{\mathcal{M}}_n^{\dagger} -\langle{\hat{\mathcal{M}}_n+\hat{\mathcal{M}}_n^{\dagger}}\rangle_E{\hat{\varrho}}_k^{\mathrm{s}} \right)\!\mathrm{d}V_n(t), \end{eqnarray} \tag{ A.17 }$$

where we note that ${\mathrm { d}} \hat \rho _k^{\mathrm {s}}\equiv {\mathrm { d}} (\langle {\vec {\gamma }_k}|\hat \varrho |{\vec {\gamma }_k}\rangle )$. Given this, the equation of motion for the probability P_k is easily obtained by computing the trace of equation (A.17), which provides precisely equation (11). Now, since $\hat \varrho _k^{\mathrm {s}}=\hat \rho _k^{\mathrm {s}}/P_k$ we have

$$\begin{eqnarray} &&\mathrm{d}{\hat{\varrho}}_k^{\mathrm{s}}= \frac{1}{P_k}\cdot\mathrm{d}\hat\rho_k^{\mathrm{s}}+\hat\rho_k^{\mathrm{s}}\cdot\mathrm{d}\left(\frac{1}{P_k}\right) +\mathrm{d}\hat\rho_k^{\mathrm{s}}\cdot\mathrm{d}\left(\frac{1}{P_k}\right), \end{eqnarray} \tag{ A.18 }$$

where (classical Itô formula [28])

$$\begin{eqnarray} &&\mathrm{d}\left(\frac{1}{P_k}\right)=-\frac{1}{P_k^2}\mathrm{d}P_k +\frac{1}{P_k^3}(\mathrm{d}P_k)^2. \end{eqnarray} \tag{ A.19 }$$

Thus, we have

$$\begin{eqnarray} &&\frac{(\mathrm{d}P_k)^2}{P_k^2} = \sum_{n=1}^{M}\eta_nM_n\left( \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} \right)^2\!\mathrm{d}t, \end{eqnarray} \tag{ A.20 }$$

where we used the fact that dW_n(t)dt = 0 and dW_n(t) dW_n'(t) = δ_n,n' dt. Consequently, we obtain

$$\begin{eqnarray}&&\fl \mathrm{d}\!\!\left(\!\frac{1}{P_k}\!\right)\!=\!\frac{1}{P_k} \!\sum_{n\!=\!1}^{M}\left\{\eta_nM_n( \langle{\skew4\hat{\cal{M}}_n\!+\!\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k\!-\!\langle{\skew4\hat{\cal{M}}_n\!+\!\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} ) ( \langle{\skew4\hat{\cal{M}}_n\!+\!\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k\!-\!\langle{\skew4\hat{\cal{M}}_n\!+\!\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{k_0} )\mathrm{d}t\!-\!\sqrt{\eta_nM_n}\right. \nonumber\\ &&\times\left. ( \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k-\langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} )\mathrm{d}W_n(t)\right\}, \end{eqnarray} \tag{ A.21 }$$

and therefore

$$\begin{eqnarray} \mathrm{d}\hat\rho_k^{\mathrm{s}}\cdot\mathrm{d}\!\!\left(\!\frac{1}{P_k}\!\right)&=& \sum_{n=1}^{M}\eta_nM_n\!( \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_{\mathrm{E}} - \langle{\skew4\hat{\cal{M}}_n+\skew4\hat{\cal{M}}_n^{\dagger}}\rangle_k )\nonumber\\ &&\times \left( \hat{\mathcal{M}}_n{\hat{\varrho}}_k^{\mathrm{s}} +{\hat{\varrho}}_k^{\mathrm{s}}\hat{\mathcal{M}}_n^{\dagger} -\langle{\hat{\mathcal{M}}_n+\hat{\mathcal{M}}_n^{\dagger}}\rangle_{\mathrm{E}}{\hat{\varrho}}_k^{\mathrm{s}} \right)\mathrm{d}t. \end{eqnarray} \tag{ A.22 }$$

Putting all together into equation (A.18) and by using the definition (2) we derive (10).

We can represent the (normalized) density matrices associated with each value $\vec {\gamma }_k$ as a Bloch vector, and propagate the collection of Bloch vectors as subject to the noisy detection signals—governed by equation (2). These equations have the form

$$\begin{eqnarray}&&\fl \mathrm{d}\vec{r}_{k} = 2\{ [(\vec{b}_{k}\times\vec{r}_{k})-(\alpha_1+\alpha_2+\alpha_3)\vec{r}_{k}+\vec{\alpha}*\vec{r}_k+ (1-\delta_{k_0,k})( \vec{{I}}_k(1-r_{k}^2)-\vec{r}_{k}\times(\vec{r}_{k}\times\vec{{I}}_k) ) ]\mathrm{d}t \nonumber\\ &&+ \mathrm{d}\vec{\Omega}(1-r_{k}^2)-\vec{r}_{k}\times(\vec{r}_{k}\times\mathrm{d}\vec{\Omega}) \}, \end{eqnarray} \tag{ A.23 }$$

where we have passed to dimensionless units by performing the replacement t → Mt with $M=\max \{M_1,M_2,M_3\}$. Besides this, we have defined the vectors $\vec \alpha$, with components α_n = M_n/M, $\vec {\eta }=(\eta _1,\eta _2,\eta _3)$, $\vec {b}_k=\mu _{\!B}\vec {B}_k/(\hbar M)$ and $\mathrm {d}\vec {\Omega }$, with components ${\rm d} \Omega _n = \sqrt {\eta _n\alpha _n}\,\mathrm {d}W_n$. The symbol * in equation (A.23) indicates the pointwise product, $(\vec {p}*\vec {q})_n \equiv p_nq_n$.

Instead, the probabilities represented as the column vector $\vec {P}=(P_1,P_2,\ldots ,P_N)^{\mathsf{T}}$ obey the following matrix equation:

$$\begin{eqnarray}&&\fl \mathrm{d}\vec{P} = 4\{\vec{P}*(\vec{\upsilon}^{{\sf T}}\cdot \mathbf{C})^{{\sf T}}- \vec{P}[\vec{P}^{{\sf T}}\cdot(\vec{\upsilon}^{{\sf T}}\cdot {\bf C})^{{\sf T}})]\nonumber\\ &&-\vec{P}*[(\mathbf{C}\,\vec{P})^{{\sf T}}\tilde{\mathbf{C}}]^{{\sf T}}+\vec{P}[({\bf C}\,\vec{P})^{{\sf T}}\cdot(\tilde{\mathbf{C}}\,\vec{P})] \}\mathrm{d}t +G\,\mathrm{d}\vec{\Omega}. \end{eqnarray} \tag{ A.24 }$$

Here, C and $\skew3\tilde {C}$ are 3 × N matrices and G is an N × 3 matrix containing the Bloch vector solutions of equation (A.23), C_n,j = x⁽ⁿ⁾_j, $\skew3\tilde {C}_{n,j}=\eta _n\alpha _nC_{n,j}$ and $G_{j,n}=2P_j(x_j^{(n)}-\vec {P}^{\mathsf {T}}\cdot \vec {X}_n)$ where $\vec {X}_n=(x_1^{(n)},x_2^{(n)},\ldots ,x_N^{(n)})^{\mathsf {T}}$, and $\vec {\upsilon }=\vec {\eta }*\vec {\alpha }*\vec {r}_{k_0}$.

For the numerical simulation of both (A.23) and (A.24) we employed an Itô–Euler integrator [28] with a time step Δt ranging from 2 × 10⁻⁷ to 10⁻⁵ M⁻¹ depending on the size of the set $\mathcal {V}_B$ and on the number of switched off detectors. Such a choice enabled us to have an efficient integrator, even though some of the quantum trajectories might have been unstable. To solve the instability problem, we have first tried to apply an implicit Milstein method [28, 29], but since both (A.23) and (A.24) are nonlinear, one has to solve numerically at each time step (e.g. by means of the Nelder–Mead method [25]) implicit equations like $\vec {r}_k(t+\Delta t)=\vec {r}_k(t)+f(\vec {r}_k(t+\Delta t))$, where f is the rhs of equation (A.23) plus some additional term due to the Milstein routine. While such a strategy might solve the problem, we have numerically observed that such an approach is significantly more time consuming than the Itô–Euler integrator. Thus, we also applied a derivative free order 2.0 weak predictor corrector method [29], which turns out to be quite efficient in the case of a single Wiener noise process, but in the case of three detectors we noticed, as for other predictor–corrector methods, that the instability could not be fixed. Hence, we employed the simple Itô–Euler integrator with (rather) small time steps (up to Δt = 2 × 10⁻⁷ M⁻¹). We noticed that with such a simple strategy the unstable quantum trajectories could have been reduced or even suppressed, but at the expenses of a very long numerical computation.