Estimation of classical parameters via continuous probing of complementary quantum observables

We discuss how continuous probing of a quantum system allows estimation of unknown classical parameters embodied in the Hamiltonian of the system. We generalize the stochastic master equation associated with continuous observation processes to a Bayesian filter equation for the probability distribution of the desired parameters, and we illustrate its application by estimating the direction of a magnetic field. In our example, the field causes a ground state spin precession in a two-level atom which is detected by the polarization rotation of off-resonant optical probes, interacting with the atomic spin components.


I. INTRODUCTION
High precision metrology with quantum systems is a research field of high current activity. Through the quantization of their energy levels, elementary quantum systems provide fundamental time and frequency standards and, due to the highly developed means for preparation, control, and detection of these systems, they serve as excellent probes of various perturbations such as applied electric and magnetic fields, and inertial effects associated with rotation, acceleration, and classical or relativistic gravitational effects.
The theoretical research proceeds along different directions according to the different measurement schemes. Thus, for experiments where a quantum system is subject to a perturbation for a given short duration of time, the search for the initial quantum states on which different values of the perturbation leads to the most distinguishable outcomes has promoted the use of concepts such as the Fisher information and the Cramer-Rao bound [1], and has identified squeezed states, Schrödinger cat-like states and generalizations hereof as useful resource states in metrology. Along a different path, measurements that occur continuously in time, such as continuous wave laser spectroscopy, are made subject to analyses, that serve to exhaust the information about the desired parameters from the entire sequence of measurement data. While a simple relationship between the average signal, e.g., an absorption profile, and an unknown physical parameter provides a relatively straightforward method for estimation, the systematic extraction of a reliable error-bar on the estimate is more challenging [2][3][4].
In this work we present such an analysis for the nontrivial case of the continuous probing of a single quan-tum system. We detail a Bayesian analysis, which treats our description of unknown parameter values and the unknown state of the quantum system on an equal footing and where any data received serves to update our prior estimate of the parameters of interest. This idea has been previously applied to the case of quantum-non-demolition (QND) probing of quantum systems [2,5], where it is largely equivalent to a Kalman filter [3]. The formalism presented here, however, is more general, and new physical features appear, when not only a single QND variable is being detected.
The paper is organized as follows: In Sec. II we present our general quantum filter equation and we give an explicit derivation for our specific physical model system. In Sec. III we introduce the unknown classical parameters and we show how their representation as effective incoherent quantum degrees of freedom augments the quantum filtering equation to automatically provide a Bayesian update formula for the unknown classical parameters. In Sec. IV we present numerical simulations and we show how the fluctuating measurement signals of optical probes interacting with an atomic spin gradually filters the probability distribution of a magnetic field with known strength, but unknown direction, and how it ultimately reveals this direction with good precision. In Sec. V we conclude and discuss the outlook and perspectives of our work. Details about the derivation of the augmented stochastic master equation (SME) and on its numerical solution are provided in the appendix.

II. OPTICAL PROBING OF A SINGLE ATOMIC SPIN
A. General quantum filtering equation A quantum system with HamiltonianĤ subject to Markovian damping, described by Lindblad operatorsÔ j and rates Γ j , is described by a reduced system density matrix̺ s , which obeys the master equation where the operator D is defined as [6] Interaction with continuous quantized probe fields cause entanglement of the system with the fields which, if the field degrees of freedom are subsequently traced out, is described by inclusion of further master equation terms D[M n ]̺ with Lindblad operatorsM n and effective interaction parameters M n . The entanglement of the system and the fields, however, implies that measurements of the probe field variables after the interaction lead to a back-action on the state of the quantum system.
Subject to the quantum back-action due to continuous amplitude measurements on the probe field, the reduced density matrix of the quantum system obeys the following stochastic master equation where the operator H is defined as [6] H and where dW n are infinitesimal Wiener processes accounting for the noisy contribution to the field amplitude measurement dY D n (t) = η n M n M n +M † n dt + dW n (t). (5) In the equations (3) and (5) the parameters η n account for the detector efficiencies and transmission losses of the probe beams between the system and the detector. If η n = 0, the corresponding probe field merely contributes extra damping and decoherence to the system, while if η n = 1, and in the absence of other damping or ineffective probing terms, the system may be described by a stochastic wave function rather than a density matrix [7][8][9][10]. The random Wiener noise increments dW n constitute the so-called innovation processes, i.e., they are the difference between the experimentally observed signals dY D n (t) and their quantum-mechanical expectation values η n √ M n M n +M † n dt determined from the current quantum state of the system. This difference amounts to the shot-noise in field amplitude measurements by homodyne detection, and the formalism can also be adapted to treat, e.g., photon counting measurements, where the innovation process is described by an (infinitesimal) Poisson process. In either case, the measurement back-action has only infinitesimal influence in every small time step while for certain measurements accumulated over time it ultimately causes the collapse on a random eigenstate normally attributed to the von Neuman projection postulate.
We have presented the quantum filtering equation in a general form following a Markovian and perturbative treatment of the interaction of the system with its environment. For every concrete physical example one has to validate such a treatment and to explicitly analyze the appropriate interactions and field measurement schemes to obtain the operators and parameters in Eq. (3). In the next subsection we will recall such a derivation for the explicit case of an atom with a degenerate ground state that interacts with off-resonant probe laser fields.

B. Dynamics of a two-state atom and an off-resonant laser field
The detection of optical phase shifts and polarization rotation is at the heart of spin squeezing and quantum entanglement schemes involving the collective spin degrees of freedom associated with large atomic ensembles [11][12][13][14], and it also constitutes the basis for atomic magnetometers which today explore the quantum limits of resolution [15][16][17]. The spins in large atomic ensembles are well approximated as harmonic oscillator degrees of freedom and Gaussian approximations and classical filtering theory apply [12,18,19], while single atoms must be described by the full density matrix formalism, that we will address in the following.
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 axis z. The atom interacts through its magnetic momentˆ µ with a classical magnetic field B = (B x , B y , B z ), which drives a Larmor precession of the ground state atomic spin. The atom is coupled to an off-resonant linearly polarized laser field. This field, which is linearly polarized along the x-axis, can be decomposed into two circular components with annihilation operatorsâ ± = (∓â x + iâ y )/ √ 2. Due to the dipole selection rules, the two circularly polarized field components 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 [20]: In Eq. (6), ∆ = ω L − ω A is the laser atom detuning, g = d· E 0 / , where d is the atomic electric dipole moment and | E 0 | = ω L /(V ǫ 0 ) is the electric field per photon in the quantization volume V , and ǫ 0 is the (electric) vacuum permeability.
As as result of Eq. (6), the two circularly polarized field components experience phase shifts that depend on the atomic occupation of the two ground states. The resulting phase difference between the two field components 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 aŝ defining the canonical conjugate operatorsŷ andp y . The polarization rotation of the field is conveniently measured by subtracting the intensities of polarization components linearly polarized at ±45 degrees 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ŷ. By writingˆ µ = µˆ σ, whereˆ σ = (σ x ,σ y ,σ z ) are the Pauli matrices, and using the definitions in Eq. (7), the total Hamiltonian can be written A term g 2 2∆ (â † xâx +â † yây ) has been omitted from Eq. (8) as it gives rise to a common phase shift, but no polarization rotation of the probe laser field.
To properly account for the entanglement created between the atom and the continuous probe field, we treat the laser beam as a sequence of segments of length L = cτ , area A = V /L and volume V , each initially prepared in a coherent state before the interaction. The continuous interaction between the atom and the light beam can then be represented as a sequence of interactions of the atom with one segment after the other. Each segment, in turn, is described as a single harmonic oscillator mode described by the operators in Eq. (7). In the time interval τ the atom interacts with N ph = Φτ photons, where Φ is the photon flux, and we assume that the interaction is sufficiently weak, that the dynamics can be well described by the coarse-grained propagator where and n B ≡ (n x B , n y B , n z B ) is the unit vector pointing in the magnetic field direction. We have introducedσ nB = σ · n B and, in the part describing the atom-light interaction, we have introducedσ m =ˆ σ · m, for the Pauli operator along an arbitrary unit vector direction m. To this aim we assume that we can probe the atomic ground state spin along any such direction with probe beams of appropriate polarization and propagation direction, or by application of a unitary spin rotation prior to probing of the z-component with a fixed beam set-up. 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 in Eq. (9), we have thus neglected terms of order τ 3/2 and higher.
The incident laser beam is in a coherent state of linearly polarized light described by a Gaussian wave function π −1/4 e −p 2 y /2 in the argument p y , associated with the observablep y introduced in Eq. (7). The joint state of an atomic ground superposition state |ψ s (t) = ℓ=±1 c ℓ/2 |g ℓ/2 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 , |Ψ ps (t) = 1 π 1/4 ℓ=±1 c ℓ/2 dp y e −p 2 y /2 |p y , g ℓ/2 . (11) Under the action of (9), this state evolves into the entangled state where, to first order in τ (second order in √ τ ), the new expansion coefficients c ′ ℓ/2 (p y ) are:

C. A quantum filtering equation for the two-state atom
To describe the back-action due to the field measurement, it is convenient to transform the entangled state (12) to the y rather than p y representation of the field, using the relation and the state (12) can be rewritten as: with the new coefficientc ℓ/2 (y) given bỹ A measurement of the light probe observableŷ with outcome y D projects the state of the system (15) onto the state component with that definite value, i.e., the atomic part of the system becomes where the explicit dependence ofc ℓ/2 on the measurement outcome y D causes an (infinitesimal) transfer of amplitude among the two atomic states.
Given the quantum state (15), the probability to measure a given value y D is where y 0 = κτ σ m . This explicitly shows how the optical probing yields a signal proportional to the desired mean value σ m , and it allows us to to model the measurement outcome y D as a stochastic variable: where ∆W is a (finite) Gaussian Wiener increment with mean zero and variance τ .
By replacing y with the expression (19) for y D in Eq. (17), 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 Here the interaction parameter M is given explicitly by and an infinitesimal Wiener process models the noise in the detected signal, dY D (t) = 2 √ M σ m dt + dW . Thus, the system evolves according to a stochastic equation of the same form as the standard quantum filter equation (3).
We note that several probe fields may be used to 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 .
The modifications in Eq. (3) due to finite detector efficiency can also be understood from first principles in the model system, and lead to similar correction factors in Eq. (20).

III. CONDITIONAL DYNAMICS AND ESTIMATION OF A CLASSICAL PARAMETER
We are interested in the use of quantum systems to estimate a classical physical parameter or a set of parameters γ. Here, in order to keep the discussion as general as possible, we treat the unknown parameter γ as a vector quantity to indicate that it may be a set of parameters such as a damping rates, energy shifts, and coupling strengths, or, as in our example below, the directional components of a vector magnetic field. The experiment is sensitive to the value of these parameters, e.g., if they are coefficients in the HamiltonianĤ =Ĥ( γ) acting on the system, and if this dependence results in a change of the observables probed in the experiment.
We describe the quantum dynamics of the combined system with γ belonging to a finite set of values V γ = { γ k : k = 1, . . . , N }.
For an observer who knows the true value of γ = γ k0 , the system is described by our original reduced system stochastic master equation (3) with H = H( γ k0 ), conditioned by the measurement signals dY D n (t) = √ η n dW n (t) + η n M n M n +M † n 0 dt. (23) where M n +M † n 0 = Tr S {(M n +M † n )̺ s 0 }, and where dW n (t) are standard Wiener processes. Here Tr S (·) is the trace on the system Hilbert space.
In the following we will assume that we probe all three spin components (σ x ,σ y ,σ z ) with measurement strengths (M 1 , M 2 , M 3 ) and efficiencies (η 1 , η 2 , η 3 ). By introducing the three-dimensional real Bloch vector, r = Tr(̺ s 0ˆ σ), the stochastic master equation (22) can be rewritten in the Bloch vector representation, where we have passed to dimensionless units by performing the replacement t → M t with M = max{M 1 , M 2 , M 3 }. Besides this, we have defined the vectors α, with components α n = M n /M , η = (η 1 , η 2 , η 3 ), b k = µ B B k /( M ), and d Ω, with components dΩ n = √ η n α n dW n . The symbol * in Eq. (24) indicates the pointwise product, ( p * q) n ≡ p n q n .
Equation (24) was derived and analyzed in detail in Refs. [21,22], for the special case of identical probing strengths and efficiencies for all three orthogonal spin directions. As shown in that work, in the absence of a magnetic field, the Bloch vector is driven towards a steady state radial (purity) distribution and an isotropic angular distribution within the Bloch sphere. The magnetic field, in turn, provides a torque for the atomic spin vector, and the resulting spin precession, in the plane perpendicular to the field, reveals itself in a modulation of the polarization rotation measurements.

A. Augmented quantum filter equation
In our formulation of the estimation problem we will treat the classical parameter γ as unknown with an assigned probability distribution P ( γ). The measurements then cause an update of the probability distribution, which is governed by Bayes rule for conditional probabilities. Until the actual value of γ is known, we thus have to treat each candidate value with a probability factor, and for each possible value of γ, the corresponding state of the quantum system evolves under the quantum filtering equation with the corresponding dependence of γ.
To this end it is convenient [2,4,[23][24][25] to consider the augmented Hilbert space H = H s ⊗ H γ , where H s and H γ refer to the quantum system Hilbert space and the space of classical states for the variable γ, 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 γ k as a diagonal density matrix ̺ γ = k P k | γ k γ k |. The classical variables γ are equivalent to quantum non-demolition (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 backaction. When we incorporate the parameters γ in this way we can directly apply the filtering equation on the augmented space.
The observer who does not know the value of γ describes the combined quantum and classical system by the augmented density matrix where̺ s k is the normalized system state density matrix associated with the specific value γ = γ k .
The combined system evolves according to the quantum filter equation where the operator H is defined as: Here we use the notation f E = Tr(f̺) = N k=1 P k Tr S (f̺ s k ) 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 γ.
If the dependence on γ only enters via the Hamiltonian H( γ), the different terms in (26) are implemented as the following product operators on

B. Detection signal properties
The stochastic process appearing in Eq. (26), is not a standard Wiener process. This is because we subtract from the measured signal a weighted average based on our probabilistic description of γ, while the measured photocurrent dY D n in the experiment is governed by the actual value of the unknown parameter γ = γ k0 .
Equation (23) characterizes the properties of such a realistic detection record, and, when inserted in Eq. (26), we find that the stochastic process in Eq. (30) can be rewritten as Hence, dV n (t) is a stochastic Gaussian process with variance dt, and with (statistical) mean value dV n (t) = √ η n M n ( M n +M † n 0 − M n +M † n E )dt, reflecting precisely the difference in the expectation value assumed by the weighted average over different γ k and by the correct value γ k0 .
Equation (26) permits a full simulation of the detection process and provides a time dependent solution of the form̺ where each value of γ k is associated with an unnormalized ρ s k . Equation (32) is of course equivalent to Eq. (25) with the normalized system density matrix̺ s k and the probability distribution P k .
As shown in detail in the appendix V A, Eq. (26) leads to two separate equations for̺ s k and P k : with dṼ (k) In both equations the photocurrent dY D n (t), observed or simulated according to Eq. (22), appears, and Eq. (34) agrees with the expression given in Ref. [24]. We note, however, that the stochastic term in our Eq. (33) contains the expectation value M n +M † n k corresponding to the parameter value γ k , while in Ref. [24] the ensemble average M n +M † n E has been used. By inserting the expression (31) for dV n in Eq. (34), we see that the change of dP k due to the measurements is given by This equation has a natural interpretation: For parameter values γ k where the expected mean current ∝ M n +M † n k differs in the same (opposite) direction from the ensemble mean as the one expected for the actual value M n +M † n 0 , the two parentheses will typically have the same 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 .

IV. VECTOR MAGNETOMETRY WITH A TWO-LEVEL ATOM
We now apply the formalism to the two-level atom coupled to a magnetic field B with known magnitude | B|, but unknown direction in space. Since such a field will cause a spin precession around the magnetic field axis, we expect that optical probing of a single spin component will not be sensitive to the magnetic field projection along the spin direction probed, while the simultaneous probing of all three spin components is bound to reveal any motion of the mean spin vector due to the magnetic precession.
The acquisition of data from a real or a simulated experiment cause a continuous update of the probability distribution P k for the magnetic field, represented in the following by the dimensionless vector, b k = µ B B k /( M ). The angular measure, cos , quantifies the scattering of the magnetic field directions b k inferred from a single experimental run around the actual value b u . By carrying out a large number of simulations, we thus quantify the average performance of the method by the (average) scalar product as a function of the measurement time t.
A. Direction of a magnetic field along a given axis Following Ref. [24], we have first investigated the case in which the initial state of the atom is represented by the Bloch vector r 0 = (0, 1, 0) with M 1,2 = 0 but M 3 = 0 (i.e., we probe only the component of the spin along z and M ≡ M 3 ), 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 Fig. 1 we show the behavior of cos θ (t) as function of time in the two cases of a weak | b u | = 0.1 (lower, black curve) and a strong | b u | = 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 Ref. [24], the quantum filter does not unambiguously identify the direction b u of the unknown field.
This situation changes when all the three spin components are detected. As illustrated in Fig. 2, the accumulation of results from all three detectors lead, especially for high strength of the field, to a final probability distribution which is well converged to the correct b u . 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. There is a number of competing effects that may explain the dependence on the quality of our estimate on the field strength and the number of spin components probed: Measuring the z-component of a spin may lead to a (Zeno-effect) suppression of, and, hence, insensitivity to, slow precession of the spin around the x-axis, while even for strong fields, the measurement of a single spin component does not allow discrimination between left and right circular precession around the x-axis. The probing of several components on the one hand allow 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.

B. Estimation of a spherically random direction of a magnetic field
Finally, we have analyzed the scenario in which the magnetic field has an isotropic prior probability distribution, represented by an ensemble V B with N = 98 directions on the unit sphere, as illustrated in Fig. 3. Here, we assume the stronger field | b u | = 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 Fig. 3 at time M t = 0. We study the convergence of the probability distribution as a function of the probing time, and for long times (M T = 15), the filter converges well to a single direction. We note that the spheres are cut on the top because in the simulation we have consid- ered an interval θ = (0, π) equally spaced with N θ = 7 grid points and an interval φ = [0, 2π) with N φ = 14 grid points for the azimuthal angle.
The equal strength probing of the three cartesian spin components σ x , σ y , σ z is equivalent [21,24] 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.

V. CONCLUSION
In this work, we have demonstrated a Bayesian filter for classical parameters, which affect the dynamics of a quantum system. Previous studies along the same lines have focused on quantum non-demolition measurements of typically a single variable, but as shown by our analysis, a non-QND setting may be analyzed by the very same assumptions and methods. Non-QND probing may have specific advantages and provide more decisive results, when the parameters affect different observables of the quantum probe, as illustrated explicitly by our numerical simulations.
The Bayesian filter is derived from a standard quantum filter formulation of conditioned quantum dynam-ics. In this mapping, we model the classical parameters as QND observables of auxiliary quantum systems, and their classical probability distribution thus coincides with the conventional reduced density matrix elements for a quantum system. Since the quantum state description cannot be completed by further knowledge in the form of hidden variables, our formulation of the parameter estimation problem, indeed, provides the tightest and most precise probability distribution for the variable of interest conditioned on the measurement outcome and on the prior probability distribution.
The discretisation of the parameter space and solution of a quantum system master equation associated with each potential parameter value naturally puts limit on the precision of the method and the number of variables that can be realistically determined. A natural next step would be to apply methods that gradually refine the parameter space around the most likely values and suppress the most unlikely ones from the calculation. Such weighted stochastic differential equations are known in statistics [26], and we imagine that they may be used to decide objective means to suppress and to breed new parameter values without enlarging the memory and computational demands of the method. Another interesting approach, put forward in Ref. [27], involves projection of the complete system on a non-linear lower-dimensional manifold on which the integration of the stochastic differential equations of motion is faster. Alternatively, maximum likelihood methods and random searches through the parameter space, e.g., by Markov Chain Monte Carlo methods [28], may be effectively applied to even very large search spaces. We imagine that our simulations may serve as useful reference data for testing such alternative estimation techniques.
where we note that dρ s k ≡ d( γ k |̺| γ k ). Given this, the equation of motion for the probability P k is easily obtained by computing the trace of Eq. (37), which provides precisely Eq. (34). Now, since̺ s k =ρ s k /P k we have: where (classical Itô formula [29]) Thus, we have where we used the fact that dW n (t)dt = 0 and dW n (t)dW n ′ (t) = δ n,n ′ dt. Consequently, we obtain Putting all together into Eq. (38) and by using the definition (23) we derive (33).

B. Numerical simulations
In analogy to Eq. (24), we can represent the (normalized) density matrices associated with each value γ k as a Bloch vector, and propagate the collection of Bloch vectors as subject to the noisy detection signals -governed by Eq. (23). Using the same notation as in Eq. (24), these equations have the form d r k = 2 ( b k × r k ) − (α 1 + α 2 + α 3 ) r k + α * r k Instead, the probabilities represented as the column vector P = (P 1 , P 2 , . . . , P N ) T obey the following matrix equation d P = 4 P * ( υ T · C) T − P P T · ( υ T · C) T ) − P * (C P ) TC T + P (C P ) T · (C P ) dt + G d Ω. (44) Here, C andC are 3 × N matrices and G is an N × 3 matrix containing the Bloch vector solutions of Eq. (43), C n,j = x (n) j ,C n,j = η n α n C n,j , and G j,n = 2P j (x (n) j − P T · X n ) where X n = (x (n) 1 , x (n) 2 , · · · , x (n) N ) T , and υ = η * α * r k0 .
For the numerical simulation of both (43) and (44) we employed an Itô-Euler integrator [29] with a time step ∆t ranging from 2 × 10 −7 · M −1 to 10 −5 · M −1 depending on the size of the set 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 Miltstien method [29,30], but since both (43) and (44) are nonlinear, one has to solve numerically at each time step (e.g., by means of the Nelder-Mead method [28]) implicit equations like r k (t + ∆t) = r k (t) + f ( r k (t + ∆t)), where f is the r.h.s. of Eq. (43) plus some additional term due to the Miltstien 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 [30], which turns out to be quite efficient in the case of a single Wiener noise process, but in the case of three detectors, whose generalization is not straightforward, 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 −7 · M −1 ). We noticed that with such a sim-ple strategy the unstable quantum trajectories could have been reduced or even suppressed, but at the expenses of a very long numerical computation.