This corresponds to the single measurement operator \( L_1 = \mathbf{N} = \sum_{n=0}^{n_{\max}} \, n\, |n\rangle\langle n| \, ,\) as can now be performed routinely in various experimental setups [44, 45, 46]. The most interesting point concerns the correlated evolution of eigenstate populations, where we can obtain various invariant combinations, like:

The simplest setting corresponds to \( L_1 = \mathbf{X} = \int_{\mathbb{R}} x\, |x\rangle\langle x| \, dx\) , representing for instance the position operator, in the complete orthonormal basis \( \{|x\rangle\}\) . We also write the density operator \( \rho = \int_{\mathbb{R}^2} \, \rho(x_1,x_2)\, |x_1\rangle\langle x_2| \, dx_1\, dx_2\) in the position representation (accepting a lack of rigor for the uncountably infinite basis). The off-diagonal amplitudes then decay as

naturally expressing how coherences between positions further apart decay faster. Regarding the diagonal components, i.e. the population distribution \( p(x) = \rho(x,x)\) along \( x \in \mathbb{R}\) , we obtain in particular:

Such expressions indicate how most elements of the “shape” of \( p(x)\) are preserved independently of measurement realizations. In the form (7), the continuous limit of (9) becomes:

Thus, we multiply the distribution \( p_0(x)\) by a Gaussian, whose stochastic mean \( \gamma_{t}\) reflects the integrated measurement output, while the standard deviation \( \sigma_t\) decreases deterministically. Finally, we show in appendix how a heterodyne measurement, with \( L_1 = \mathbf{X}\) and \( L_2 = i \, \mathbf{X}\) , induces in addition a particularly correlated evolution of the off-diagonal phases.

Note that here we have assumed no other dynamics than the position measurement \( X\) , which becomes infinitely precise as \( t\) goes to infinity (at the price of infinite uncertainty on momentum, the conjugate variable). In Section 3.2 we add realistic elements to the dynamics, under which the (un)certainties remain finite, and give the corresponding full state expressions in reduced form.

We finish with an example where it is natural to measure several commuting hermitian \( L_k\) , namely in quantum error correction (QEC), see e.g. [25, chapter 10],[47]. In the most common paradigm of stabilizer codes, each \( L_k\) compares adjacent qubits to monitor an error syndrome, with eigenvalues \( \pm 1\) . The eigenspace associated to +1 on all \( L_k\) corresponds to the valid codespace, other measurement results identify corresponding errors.

The simplest error-correcting code is the 3-qubit repetition code. The Hilbert space is \( \mathcal{H}=\text{span}(\{|0\rangle,|1\rangle\} \otimes \{|0\rangle,|1\rangle\} \otimes \{|0\rangle,|1\rangle\}) \simeq \mathbb{C}^8\) and the measurements each compare whether two qubits have the same or a different value in the canonical basis. Denoting \( P_{\{|v_1\rangle,|v_2\rangle,...\}}\) the orthonormal projector onto the subspace spanned by \( v_1,v_2,...\) , we thus have:

Together, these \( L_k\) distinguish four eigenspaces, each one twofold degenerate:

The twofold degeneracy gives the possibility to encode a qubit in the codespace \( |\psi\rangle = c_0 |000\rangle + c_1 |111\rangle\) ; in the other eigenspaces, it allows to identify a single spurious bit-flip \( 0 \leftrightarrow 1\) , via changes in the measured eigenvalues, without any loss of quantum information \( c_0,c_1\) .

• In the form of Proposition 2, we have \( c^{(a,b)}_t = c^{(a,b)}_0\) when \( (|a\rangle,|b\rangle) \subset V_k\) for some \( k\) (no loss of encoded quantum information) and \( c^{(a,b)}_t = c^{(a,b)}_0 \, e^{-(1-\eta) 8 t}\) otherwise (two of the three syndromes \( L_k\) differentiate \( V_j\) from \( V_{j'\neq j}\) , inducing decay of coherences among them). Regarding diagonal elements, \( z^{(\alpha)}_t = \frac{\rho_t(a,a)}{\rho_t(b,b)} = z^{(\alpha)}_0\) is conserved for \( (|a\rangle,|b\rangle) \subset V_k\) for each \( k=0,1,2,3\) . This just says that within each eigenspace, relative populations do not change. Since we expect stochastic motion along \( 3\) dimensions, there are no further deterministic correlations among diagonal elements: the relative populations between the different eigenspaces \( V_0,\, V_1,\, V_2,\, V_3\) can evolve in all directions as a function of the three measurement signals. In the form of Proposition 3, we can thus write in particular \( \; \rho_t(a,b) = \rho_0(a,b) \gamma_{t,k} \;\) for each \( 2\times2\) matrix block corresponding to a subspace \( (|a\rangle,|b\rangle) \subset V_k\) , with the real numbers \( \gamma_{t,k}\) evolving independently (up to overall normalization) as a function of measurement results. The evolution of the other \( 2\times2\) matrix blocks, corresponding to a coherence between two different subspaces \( V_k,V_j\) , is then deduced from the evolution of these diagonals and the deterministic dynamics of \( c^{(a,b)}_t\) .

• In the above setting, one measurement channel is redundant: applying only e.g. \( L_1\) and \( L_3\) measurements is sufficient towards identifying to which subspace \( V_k\) the state belongs. Then \( \mathcal{M}_t\) has one dimension less and there is another deterministic variable. Solving for \( \alpha\) according to Proposition 2 yields:

  \[ z^{(\alpha*)}_t = \frac{\rho_t(000,000) \rho_t(010,010)}{\rho_t(001,001)\rho_t(100,100)} = z^{(\alpha*)}_0 \; . \]

  The conservation of \( z^{(\alpha*)}\) indicates how, when e.g. the state converges towards \( V_0\) under stochastic measurement results, the population on \( V_2\) must converge to zero quicker than the populations on \( V_1\) and \( V_3\) . This reflects that confusing \( V_0\) with \( V_2\) would correspond to confusing the eigenvalues of both \( L_1\) and \( L_3\) , whereas confusing e.g. \( V_0\) and \( V_1\) implies confusion on a single measurement operator. The \( c^{(a,b)}_t\) similarly follow different decay rates in absence of \( L_2\) .
• The repetition code is meant to counter spurious bit flips. These are added to the model through operators \( L_{3+k}=\sqrt{\Gamma_k} \sigma_{x,k}\) , where \( \sigma_{x,k} = |0\rangle\langle 1|+|1\rangle\langle 0|\) on qubit \( k\) , and with \( \eta_4=\eta_5=\eta_6=0\) . Thus the SDE is driven by no new stochastic processes. However, as in the classical example with the rolling wheel, the commutation of stochastic vector fields with new deterministic dynamics can induce diffusion in many more directions. This turns out to be the case here indeed, as we have quickly generated up to 15 different vector fields in \( \mathfrak{G}_F\) by hand even for the special situation \( \eta_1=\eta_2=\eta_3\) and \( \Gamma_1=\Gamma_2=\Gamma_3\) (see appendix). There is thus no drastic and exact model reduction in this realistic situation, compared to modeling the full state \( \rho_t\) with 64-1 real components. In other words: the 3 values of integrated measurement results \( y_{t,1}, y_{t,2}, y_{t,3}\) would not be enough anymore to identify the state, and in fact the associated memory effect would have to be modeled with at least 15 variables. An idea towards efficient approximate model reduction, requiring little online adaptation [27], could be to typically drop terms corresponding to higher-order commutators. Anyways, any small approximate continuous-time model may already help refine error-correction decisions compared to a discrete-time projective measurement model.

According to Proposition 1, just measuring \( L_1\) and \( L_2\) in presence of \( H_0 = \Delta \bf{ a}^\dagger \bf{ a}\) yields deterministic manifolds of dimension 2. Adding a Hamiltonian proportional to any combination of \( \bf{ X}\) and \( \bf{ P}\) , which is the form of typical control signals, does not change the manifold dimension. When the system decays into an environment of nonzero temperature \( n_{th} > 0\) , an operator \( L_3=\bf{ a}^†\) is added with \( \eta_3=0\) . This situation features deterministic manifolds of dimension 4. When \( \Delta=0\) and not measuring \( L_2\) (homodyne detection), all these dimensions are divided by two.
Thus, for \( n_{th} >0\) , attempting to describe the state with a smart integral of each output signal is bound to fail. Nevertheless, it is sufficient to track two integrals per output channel, in order to capture the corresponding memory effect, and obtain again a compact and fully general description of \( \rho_t\) .

The aim is thus to write the solution in compact form for the SDE:

Here \( u_t\) and \( v_t\) are real control signals, possibly time-dependent, and \( \Delta=0\) without loss of generality by redefining variables in a rotating frame. Inspired by the Gaussian form observed in the QND measurement of \( \bf{ X}\) , the following result is obtained.

The solutions for the evolution of those variables, as well as the full proof, can be found in appendix.

This result concerns the most general initial condition. It may further simplify for particular \( \rho_0\) . For instance, for \( n_{th}=0\) and starting with a so-called coherent state \( \rho_0 = |\alpha\rangle\langle \alpha|\) , \( \alpha \in \mathbb{C}\) , the solution \( \rho_t\) is independent of measurement results. Indeed, the Lindblad equation which reflects the mixture over all possible measurement realizations then remains in a pure state \( \rho_t = |\alpha(t)\rangle\langle \alpha(t)|\) , and this is only possible if all measurement-conditioned trajectories coincide [52, Ch.4,Lem.2].

Despite \( \bf{ X}\) and \( \bf{ P}\) not commuting, \( [G_\bf{ X},\, G_\bf{ P}] = 0\) . In the criterion of Proposition 1, commutation with the deterministic dynamics gives two new vector fields \( G_{i\bf X},\, G_{i\bf P} \in \mathfrak{G}_F\) . Further commutations give no new element. Adding Hamiltonians in \( \bf{ X}\) , \( \bf{ P}\) , \( \bf{ a}^† \bf{ a}\) and dissipation in a thermal environment (\( L_3=\sqrt{(1+n_{th})}\bf{ a}\) , \( L_4=\sqrt{n_{th}}\bf{ a}^†\) with \( \eta_3=\eta_4=0\) ) maintains this property. In all cases, there are thus deterministic manifolds of dimension 4.

A most general setting introduces different rates on the two channels. Then, defining \( \tilde\bf{ X}_t = \cos(\Delta t)\, \bf{ X} + \sin(\Delta t)\, \bf{ P}\) and \( \tilde\bf{ P}_t = \cos(\Delta t)\, \bf{ P} - \sin(\Delta t)\, \bf{ X}\) in a rotating frame with respect to \( H_0 = \Delta \bf{ a}^\dagger \bf{ a}\) , we consider the model:

The block-triangular set of equations governing these variables is provided in appendix, with simplifications in special cases. The proof is also given in appendix.

Putting \( \Gamma_p=0\) yields back a QND measurement of \( L_1=\bf{ X}\) (see Section 3.1), now associated to further dynamics. For a general \( \Delta\) and thermal environment, \( \rho_t\) then remains confined onto \( \mathcal{M}_t\) still of dimension 4, and best described by Proposition 5 ; when \( \Delta=0\) , it reduces to dimension 2; when in addition \( \Gamma_{\ell}=0\) , it reduces to dimension 1 in agreement with Proposition 2.

Consider subsystem B a harmonic oscillator (photon annihilation operator \( \bf{ b}\) ) and subsystem A a finite-dimensional system on \( \mathcal{H}_A \simeq \mathbb{C}^d\) , with the dynamics:

Operator \( Q_A\) is Hermitian and no other dynamics are allowed on subsystem A, while \( u_t+i v_t\) is a control signal on subsystem B whose annihilation channel is monitored. This setting is typical in “dispersive readout” of qubits with \( Q_A=\sigma_z\) , see e.g. [58, 59]. It can be viewed as an indirect QND measurement of \( Q_A = \sum_{s=1}^d \, \lambda_s \, |s\rangle\langle s|\) , considered non-degenerate here for simplicity of notation. Previous work has established compact expressions for solutions of (13) for particular initial conditions, e.g. a separable initial state with subsystem B in a coherent state [59]. The following provides a low-dimensional characterization of the solution for any initial state — see proof in Appendix.

• For \( u_t=v_t=0\) , the manifolds \( \mathcal{M}_t\) are of dimension \( 2d\) and characterized by the Abelian algebra: \( \; \mathfrak{G}_F = \text{span}\{G_{|s\rangle\langle s| \otimes \delta \bf{ b}} : s=1,2,...,d;\; \delta =1,i \} \; .\) For general control signals, \( \mathcal{M}_t\) has dimension \( 4d-2\) , characterized by the Abelian algebra: \( \; \mathfrak{G}_F = \text{span}\{G_{|s\rangle\langle s| \otimes \delta \bf{ b}} : s=1,2,...,d;\; \delta =1,i \} \cup \{ G_{|s\rangle\langle s| \otimes \delta \bf{ I}} : s=1,2,...,d-1;\; \delta =1,i \} \; .\) Furthermore, the set of states reachable by system (13), under all measurement realizations and all control signals \( u_t,v_t\) , is confined to a time-dependent manifold \( \tilde{\mathcal{M}}_t\) of dimension \( 3 d^2 +2d - 1\) . This provides the minimal representation linking input and output signals of (13).

  The presence of drives on the harmonic oscillator here modifies the manifold dimension, unlike when the harmonic oscillator was considered alone. More detrimental modifications — in the sense of increasing the manifold dimension — seem to show up when trying to add thermal photons, or when adding dynamics on subsystem A which does not commute with \( Q_A\) . From a complementary viewpoint, this bounds the controllability of the qudit-harmonic oscillator joint system [60], when seeing measurement-induced evolution as a resource; and thus in principle, including any devisable feedback control protocols.
• The full solution can be expressed in the form \( \rho_t = \sum_{j,k=1}^{d} |j\rangle\langle k| \otimes \rho_{t,(j,k)} \; ,\) where each \( \rho_{t,(j,k)}\) follows a Gaussian kernel, see Appendix. For \( u_t=v_t=0\) , all the variables are deterministically tied to the Gaussian means of diagonal elements \( \rho_{t,(k,k)}\) . Nonzero control signals further induce stochastic components in the relative normalizations of these Gaussians. The algebraic criterion has been an invaluable tool in initiating (and stopping) the search for dependencies among the variables involved.
• Asymptotically, for large times, the full solution characterizes the following behavior:
  The harmonic oscillator state \( \rho_{t,(k,k)}\) converges towards a so-called coherent state \( =|\alpha_k\rangle\langle \alpha_k|\) whose amplitude \( \alpha_k\) is conditioned on qudit level \( k\) , possibly moving under control drives, but independent of the measurement results.
  For \( \eta \ll 1\) (or if the system is initialized with the \( |\alpha_k\rangle\langle \alpha_k|\) ), the qudit level populations still undergo significant dynamics after this initial transient; these dynamics correspond to the standard QND measurement process, with \( k\) -conditioned coherent state amplitudes \( \alpha_k\) replacing the measurement eigenvalues of Section 3.1. See Appendix for details.

Consider subsystem B a harmonic oscillator (photon annihilation operator \( \bf{ b}\) ) and subsystem A another harmonic oscillator (photon annihilation operator \( \bf{ a}\) ), with the dynamics:

Such coupling model is obtained after Rotating Wave Approximation in a dipolar coupling between quasi-resonant oscillators A and B, with \( u_t+i v_t\) a complex drive amplitude on A and \( \Delta\) expressing a small residual resonance mismatch [4, Chap.3.2]. When the B subsystem is subject to heterodyne measurement and a thermal environment, this is again a linear quantum system [42].

• The system (14) with \( n_{th}>0\) features deterministic manifolds \( \mathcal{M}_t\) of dimension \( 8\) , characterized by the algebra:

  \[ \mathfrak{G}_F = \text{span}\{G_{\delta {\bf b}},\, G_{\delta {\bf b}^†},\, G_{\delta {\bf a}},\, G_{\delta {\bf a}^†};\; \delta =1,i \} \; . \]

  For \( n_{th}=0\) , the manifold further reduces to dimension 4 as all the \( ^†\) terms in \( \mathfrak{G}_F\) disappear.

  Proof: Operators acting on B alone generate an algebra: \( \mathfrak{G}_B = \text{span}\{\, G_\bf{ b}\,,\; G_{i\bf{ b}} \, \} \; \) when \( n_{th}=0\) , adding \( G_\bf{{ b}^†},\; G_{i\bf{ b}^†}\) for \( n_{th}>0\) . The commutator of \( G_\bf{ b}\) (respectively \( G_\bf{{ b}^†}\) ) with the coupling Hamiltonian just yields \( G_\bf{ a}\) (respectively \( G_\bf{{ a}^†}\) ). Commutation of the latter two yields \( G_{[\bf{ a},\bf{ a}^†]} = G_\bf{ I} = 0\) . Further commutation of these A-dynamics with any B-dynamics vanishes. Further commutation of these A-dynamics with the coupling Hamiltonian yields back \( G_\bf{ b}\) and \( G_\bf{{ b}^†}\) . Hence the algebra is closed at this point. \( \square\)

• Adding simple dynamics on the B subsystem would not change this manifold structure. First, detuning on B, i.e. a Hamiltonian term in \( \bf{ I} \otimes \tilde\Delta \bf{ b}^† \bf{ b}\) , would be symmetric to detuning on A, by going to another rotating frame; this does not change the dissipation on B, it just implies re-combining the output ports \( dy_{t,1}\) and \( dy_{t,2}\) accordingly with the quadratures’ new rotation speed, and similarly for the control signals on A. Second, adding drives on B, i.e. a Hamiltonian term in \( \bf{ I} \otimes (\tilde{u}_t (\bf{ b}+\bf{ b}^\dagger) - i\tilde{v}_t (\bf{ b}-\bf{ b}^\dagger) )\) , can also be converted into drives on the A subsystem. Indeed, notice that

  \[ g ({\bf a} \otimes {\bf b}^† + {\bf a}^† \otimes {\bf b}) + {\bf I} \otimes (\tilde{u}_t ({\bf b}+{\bf b}^\dagger) - i\tilde{v}_t ({\bf b}-{\bf b}^\dagger) ) = g (({\bf a}+\beta_t) \otimes {\bf b}^† + ({\bf a}+\beta_t)^† \otimes {\bf b})\; , \]

  with \( \beta_t = \frac{\tilde{u}_t + i \tilde{v}_t}{g}\) . Hence, consider the time-dependent unitary change of frame \( \tilde{\rho}_t = \bf{ U}_t \rho_t \bf{ U}_t^\dagger\) with \( \bf{ U}_t = \bf{ D}(\beta_t) \otimes \bf{ I}\) . Here \( \bf{ D}(\alpha) = \exp(\alpha \bf{ a}^† - \alpha^* \bf{ a})\) is the displacement operator [4] on harmonic oscillator mode A, satisfying \( \bf{ D}(\alpha)\, \bf{ a}\, \bf{ D}(\alpha)^\dagger = \bf{ a} - \alpha\) and \( \tfrac{d\bf{ D}(\alpha_t) }{dt}\, \bf{ D}(\alpha_t)^\dagger = \tfrac{d\alpha}{dt} \bf{ a}^\dagger - \tfrac{d\alpha^*}{dt} \bf{ a} + (\tfrac{d\alpha}{dt}\alpha^* - \tfrac{d\alpha^*}{dt} \alpha)\, \bf{ I}\) . With these elements and since \( \bf{ U}_t\) commutes with subsystem B, one checks that \( d\tilde{\rho}\) follows the same equation (14), except \( u_t\) and \( v_t\) are replaced by \( u_t - \tfrac{\Delta \tilde{u}_t}{g} + \tfrac{d\alpha^*/dt-d\alpha/dt}{2i}\) and \( v_t - \tfrac{\Delta \tilde{v}_t}{g} + \tfrac{d\alpha^*/dt+d\alpha/dt}{2}\) respectively.
• This structure is not surprising, since (14) again models a so-called linear quantum system [42], now on two harmonic oscillator modes. The state can again be expressed with a Gaussian kernel, on the Wigner function in joint state space. We leave the derivation of the corresponding expression as an exercise, possibly exploiting linear quantum systems theory which has notably contributed to the study of such networks of modes.

According to the algebraic criterion, in this case the number of stochastic variables is supposed to decrease. This can be confirmed on our Wigner representation, by checking the algebraic criterion now for the classical SDE governing the evolution of the 3-dimensional system \( \xi_t,\theta_t\) and \( t\) . The variable \( t\) must be explicitly included since the expressions are all explicitly time-dependent through \( a_t,s_t,d_t\) . The Stratonovich form is identical to the Itō form, because the Wiener vector field only depends on \( t\) , which itself has no stochastic evolution (\( dt = 1 \cdot dt + 0 \cdot dw^1_t\) ). Then the commutator of the deterministic vector field with the Wiener vector field, turns out indeed to be proportional to the Wiener vector field itself (no new direction), if and only if \( n_{th}=0\) . In other words, the variables \( \xi_t,\theta_t\) remain confined to a one-dimensional manifold. For \( n_{th}>0\) a new component appears in the commutator, confirming that in this case the pair of variables \( (\xi_t,\theta_t)\) driven by a single Wiener process \( dw^1_t\) , can diffuse to span the entire plane at any given time \( t\) , and can thus not be further reduced.

We can parameterize the time-dependent curve to which \( \xi_t,\theta_t\) are confined for \( n_{th}=0\) , by saying that a combination \( z_t(\xi_t,\theta_t,t)\) of the variables (not reduced to \( z_t=t\) ) must evolve deterministically. Canceling the contribution of \( dw^1_t\) in \( dz\) leads to the constraint \( \tfrac{\partial z_t}{\partial \xi} = 4 e^t\, \tfrac{\partial z_t}{\partial \theta}\) . This suggests the solution \( z_t = 4 \xi_t + e^{-t} \theta_t\) . We further check that the evolution of such \( z_t\) is deterministic, i.e. the right-hand side of

depends only on \( z_t\) and on \( t\) . (Of course this solution is far from unique as any function of \( z_t\) and \( t\) would also be a solution.) \( \square\)

We first give the detailed computations concerning the algebraic criterion of Prop.1. After noting that \( [L_1,L_2] = 0\) and hence \( [G_{L_1},\; G_{L_2}]=0\) , we compute further commutation with the deterministic evolution. We compute \( [L_j^\dagger L_j , \; L_1 ] = \frac{1}{2}( \sigma_{-A}\pm\sigma_{zA}\sigma_{-B} + \sigma_{-B}\pm \sigma_{zB}\sigma_{-A})\) with signs depending on \( j \in \{1,2\}\) , and similarly for \( [L_2,\; L_j^\dagger L_j ] \) . Taken separately, this corresponds to two new vector fields \( G_{L'_k}\) with \( L'_1 = \sigma_{zA}\sigma_{-B} + \sigma_{zB}\sigma_{-A}\) and \( L'_2 = i (\sigma_{zA}\sigma_{-B} - \sigma_{zB}\sigma_{-A})\) , whose further commutations happen to generate yet other directions. When the individual channels \( L_1\) and \( L_2\) have arbitrary rates, \( \mathfrak{G}_F\) contains all these vector fields and the model reduction would not be efficient. However, when the rates of \( F_{L_1}\) and \( F_{L_2}\) are equal (and thus, without requiring \( \eta_1=\eta_2\) since they do not appear there), we see instead that the commutator with their sum yields no new diffusion directions; you can see this thanks to the opposite signs in the above individual commutators.

One can further check that when adding nonzero drives in \( \sigma_x\) and/or \( \sigma_y\) on A and/or B, under any conditions, the diffusion spans manifolds of dimension 10 — hence, expect no easy explicit formula for the solution. Adding detuning, i.e. a Hamiltonian in \( \sigma_{z,A}\) or \( \sigma_{z,B}\) , just doubles the manifold to dimension 4.

As stated in the main text, towards writing the deterministic manifold equations, we translate the dynamics (15) to the Pauli basis.

• In fact, as the ground state on any qubit should be invariant, we readily introduce displaced Z variables, replacing \( r(Zq),\; r(qZ),\; r(ZZ)\) respectively by \( r(Zq^*) = r(Zq)-r(Iq),\; r(qZ^*) = r(qZ)-r(qI),\; r(ZZ^*)=r(ZZ)-r(ZI)-r(IZ)\) for \( q\in \{I,X,Y\}\) . The corresponding equations follow systematically, e.g. 

  \[ \begin{eqnarray*} dr(XI)_t &=& -r(XI)_t\, dt + \big( -r(ZI^*)_t + r(XX)_t - (r(XI)+r(IX))_t r(XI)_t \big) \, \sqrt{\eta_1}\, dw_{t,1}\\ && + \big( -r(XY)_t - (r(YI)-r(IY))_t\, r(XI)_t \big) \, \sqrt{\eta_2}\, dw_{t,2} \; . \end{eqnarray*} \]

  The others are similar and we do not write them all down here. Note that the output signals in this basis write:

  \[ dy_{t,1}= \sqrt{\eta_1} (r(XI)+r(IX))_t \; dt + dw_{t,1}\;\; , ~ dy_{t,2}= \sqrt{\eta_2} (r(YI)-r(IY))_t \; dt + dw_{t,2} \; . \]
• We next group the variables in sums and differences, e.g. \( X_s = r(XI)+r(IX)\) , \( XY_d = r(XY)-r(YX)\) , \( YZ_s = r(YZ)+r(ZY)\) and so on. The dynamics are obtained as sums and differences, with no Itō specificities. This decouples the system into a block of 9 variables, not influenced by, yet influencing (block-triangular form), the 6 remaining ones \( X_d,\, Y_s,\, Z_d,\, XY_s, \, XZ_d,\, YZ_s\) .

• We finally write the ratio between each of these variables and \( r(ZZ)_t\) , computing the resulting dynamics with due Itō corrections. Using \( dy_{t,1}\) and \( dy_{t,2}\) to include the stochastic contributions, this yields the following rather compact dynamics:

  \[ \begin{eqnarray*} dB_{1,t} &=& B_{1,t}\, dt - 2 \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_1 = YZ_d / r(ZZ)\\ dB_{2,t} &=& B_{2,t}\, dt - 2 \sqrt{\eta_1}\, dy_{t,1} ~ \text{for } B_2 = XZ_s / r(ZZ)\\ dB_{3,t} &=& 2 B_{3,t}\, dt + B_{1,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_3 = r(YY) / r(ZZ)\\ dB_{4,t} &=& 2 B_{4,t}\, dt - B_{2,t} \sqrt{\eta_1}\, dy_{t,1} ~ \text{for } B_4 = r(XX) / r(ZZ)\\ dB_{5,t} &=& 2 B_{5,t}\, dt + B_{1,t} \sqrt{\eta_1}\, dy_{t,1} + B_{2,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_5 = XY_d / r(ZZ)\\ dB_{6,t} &=& 2 B_{6,t}\, dt + B_{2,t} \sqrt{\eta_1}\, dy_{t,1} + B_{1,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_6 = Z_s / r(ZZ)\\ dB_{7,t} &=& 3 B_{7,t}\, dt + (2B_4-B_6)_t \sqrt{\eta_1}\, dy_{t,1} - B_{5,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_7 = X_s / r(ZZ)\\ dB_{8,t} &=& 3 B_{8,t}\, dt - B_{5,t} \sqrt{\eta_1}\, dy_{t,1} - (2B_3+B_6)_t \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_8 = Y_d / r(ZZ)\\ dB_{0,t} &=& 4 B_{0,t}\, dt + B_{7,t} \sqrt{\eta_1}\, dy_{t,1} + B_{8,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } B_0 = 1/ r(ZZ) \end{eqnarray*} \]

  concerning the first block of 9 variables, and

  \[ \begin{eqnarray*} dR_{1,t} &=& R_{1,t}\, dt ~ \text{for } R_1 = XZ_d / r(ZZ)\\ dR_{2,t} &=& R_{2,t}\, dt ~ \text{for } R_2 = YZ_s / r(ZZ) \\ dR_{3,t} &=& 2 R_{3,t}\, dt - R_{1,t} \sqrt{\eta_1}\, dy_{t,1} - R_{2,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } R_3 = Z_d / r(ZZ)\\ dR_{4,t} &=& 2 R_{4,t}\, dt - R_{2,t} \sqrt{\eta_1}\, dy_{t,1} + R_{1,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } R_4 = XY_s / r(ZZ)\\ dR_{5,t} &=& 3 R_{5,t}\, dt - R_{3,t} \sqrt{\eta_1}\, dy_{t,1} - R_{4,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } R_5 = X_d / r(ZZ)\\ dR_{6,t} &=& 3 R_{6,t}\, dt + R_{4,t} \sqrt{\eta_1}\, dy_{t,1} - R_{3,t} \sqrt{\eta_2}\, dy_{t,2} ~ \text{for } R_6 = Y_s / r(ZZ) \end{eqnarray*} \]

  for the remaining block of 6 variables. This is our basis for identifying deterministically evolving variables (like \( R_1\) and \( R_2\) already are).

We must identify 13 deterministic combinations out of the above 15 equations.

• Since \( B_1\) and \( B_2\) each move by a different the stochastic measurement result \( dy_2\) or \( dy_1\) , we would keep those as representing the 2 stochastic variables.

• We next try to eliminate the stochastic term involving just \( B_1\) in the expression of \( dB_{3,t}\) . By Itō calculus, we note that \( d(\frac{B_1^2}{2}) = \big(2 \frac{B_1^2}{2} + 2 \eta_2\big) \, dt - 2 B_1 \sqrt{\eta_2}\, dy_2\) , such that we can deduce

  \[ d\tilde{B}_{3,t} = (2 \tilde{B}_{3,t} + \eta_2)\, dt ~ \text{for } \tilde{B}_3 = (\tfrac{B_1^2}{4}+B_3) \; . \]

  A similar treatment with \( d(\frac{B_2^2}{2})\) and \( d(\frac{B_1 B_2}{2})\) yields

  \[ \begin{eqnarray*} d\tilde{B}_{4,t} &=& (2 \tilde{B}_{4,t} + \eta_1)\, dt ~ \text{for } \tilde{B}_4 = (\tfrac{B_2^2}{4}-B_4) \; ,\\ d\tilde{B}_{5,t} &=& 2 \tilde{B}_{5,t}\, dt ~ \text{for } \tilde{B}_5 = (\tfrac{B_1 B_2}{2}+B_5) \; . \end{eqnarray*} \]
• The stochastic terms in \( dB_{6,t}\) must be eliminated with a combination of \( d(B_1^2)\) and \( d(B_2^2)\) . For the stochastic terms in \( dB_{7,t}\) and \( dB_{8,t}\) , first we use products like e.g. \( d(B_4 B_2)\) to cancel \( B_{4,t} \sqrt{\eta_1}\, dy_{t,1}\) , then we add cubic terms in \( B_1,\, B_2\) to eliminate the remaining terms. For \( dB_0\) , after \( d(B_7 B_2 + B_8 B_1)\) , adding \( d(B_3^2+B_4^2+B_5^2+B_6^2)\) cancels all stochastic terms.
• A similar procedure is applied to the \( R_k\) , combining them in products with \( B_1,B_2\) . Overall, this yields the following final result.

This is but one way to describe the dynamics, as combinations of these variables yield other deterministic expressions — e.g. \( R_1/R_2\) is constant — which the reader may set up according to preferences.

The measurement vector fields \( G_\bf{{ b}_l}\) all commute. Commutation with the deterministic terms, for \( u_t=v_t=0\) , yield linear combinations of the \( \bf{ a}_k\) , possibly premultiplied by \( |j\rangle\langle j|\) . Explicitly, writing \( \bf{ q} = \sum_j q_{j,k}\; |j\rangle\langle j| \otimes \bf{ a}_k\) , one commutation of \( G_\bf{{ q}}\) with the deterministic vector field amounts to the linear map:

where \( q_{j,k}\) is considered a column vector, in tensor product structure of the \( j\) and \( k\) coordinates; column \( l\) of the matrix \( \beta\) contains the components \( \beta_{l,k}\) ; matrix \( \Delta\) is diagonal with components \( \Delta_k\) ; and matrix \( X\) is diagonal with components \( \chi_{j,k}\) (according to the same tensor product structure).

For generic values of \( \Delta, X, \beta\) , independently of the number of measurement channels \( m\) , the map (28) generates all \( 2 d n\) independent real linear combinations of \( \{\, |j\rangle\langle j| \otimes \bf{ a}_k , \; i|j\rangle\langle j| \otimes \bf{ a}_k : j=1,...,d ; k=1,...,n \, \}\) . By commuting the latter with the drive terms, in \( u_t\) and \( v_t\) , we obtain contributions in \( G_{|j\rangle\langle j| \otimes \bf{ I}}\) and \( G_{i|j\rangle\langle j| \otimes \bf{ I}}\) . Thus, taking into account that \( G_{i\mathbf{I}} = G_{\mathbf{I}}=0\) , in total the dimension of the manifold on which measurement outputs can spread the state is \( 2 d (n+1)-2\) , in presence of control signals and for general initial conditions; for \( u_t=v_t=0\) it reduces to dimension \( 2d n\) .

Like for our other examples with harmonic oscillators, we consider their Wigner function representation, involving a Gaussian kernel to model the low stochastic dimension.

• Thanks to the dispersive coupling with the qudit, we can write

  \[ \begin{equation} \rho_{t} = \sum_{j,j'=1}^d\; |j\rangle\langle j'| \otimes \rho_{t,(j,j')} \; , \end{equation} \]

  (29)

  where each \( \rho_{t,(j,j')}\) is an operator on the Hilbert space of the \( n\) cavities. For \( j=j'\) , these are unnormalized density operators, while for \( j \neq j'\) they need not be hermitian nor positive. Conditioned on the measurement results, the \( \rho_{t,(j,j')}\) evolve independently of each other. Indeed, the only place where one qudit level can influence another, is when taking the trace to define the \( dy_{t,l}\) .

• We represent each \( \rho_{t,(j,j')}\) via its multi-variate Wigner function \( \mathcal{W}^{\{\rho_{t,(j,j')}\}}(q)\) . This is just the straightforward generalization of the single oscillator case, i.e. a pseudo-probability distribution over \( 2n\) variables, corresponding to the column vector \( q=(x_1,p_1,x_2,p_2,...,x_n,p_n)\) . Taking its marginal with respect to e.g. \( (p_1,x_2,p_3,x_4,...)\) gives the (truly real) probability distribution for a measurement of \( x_1 \otimes p_2 \otimes x_3 \otimes ...\) which is indeed a valid observable. By using the translations to Wigner space of annihilation and creation operators acting on \( \rho\) , as recalled right before Section 6.2.2, the stochastic master equation (16),(17) translates into:

  \[ \begin{eqnarray} \nonumber d\mathcal{W}^{(j,j')}_t &=& \sum_{k=1}^n \; \left(u_{t,k} \frac{\partial}{\partial p_k} - v_{t,k} \frac{\partial}{\partial x_k} \right) \, \mathcal{W}^{(j,j')}_t\, dt + \delta_k^{(j,j')}\; \left(x_k \frac{\partial}{\partial p_k} - p_k \frac{\partial}{\partial x_k} \right) \, \mathcal{W}^{(j,j')}_t\,dt \\ && + \sum_{k=1}^n \; - i\, \omega_k^{(j,j')} \; \left( 2\, x_k^2 + 2\, p_k^2 -1 -\frac{1}{8}\frac{\partial^2}{\partial x_k^2} - \frac{1}{8}\frac{\partial^2}{\partial p_k^2} \right) \, \mathcal{W}^{(j,j')}_t\,dt \end{eqnarray} \]

  (30)

  \[ \begin{eqnarray} && + \sum_{l,k,k'} \beta_{l,k} \bar{\beta}_{l,k'} \; \left( \mathbf{1}_{k,k'}+ \tfrac{1}{8}(\tfrac{\partial}{\partial x_k} + i \tfrac{\partial}{\partial p_k})(\tfrac{\partial}{\partial x_{k'}} - i \tfrac{\partial}{\partial p_{k'}})\right) \mathcal{W}^{(j,j')}_t\,dt \\ \nonumber && + \sum_{l,k,k'} \frac{\beta_{l,k} \bar{\beta}_{l,k'}}{4} \; \left(\, (x_k+i p_k )(\tfrac{\partial}{\partial x_{k'}} \text{-} i \tfrac{\partial}{\partial p_{k'}}) +(x_{k'}\text{-}i p_{k'}) (\tfrac{\partial}{\partial x_{k}} + i \tfrac{\partial}{\partial p_{k}}) \, \right) \mathcal{W}^{(j,j')}_t\,dt \\ \nonumber && + \sum_{l=1}^m \, 2 \sqrt{\eta_l} \left(\sum_{k=1}^n \mathfrak{R}(\beta_{l,k})\left(x_k + \tfrac{1}{4}\tfrac{\partial}{\partial x_k}\right) - \mathfrak{I}(\beta_{l,k})\left(p_k + \tfrac{1}{4}\tfrac{\partial}{\partial p_k}\right) - s_l \right)\, \mathcal{W}^{(j,j')}_t \; dw_{t,l} \; . \\\end{eqnarray} \]

  Here \( \delta_k^{(j,j')} = \Delta_k + \frac{\chi_{j,k}+\chi_{j',k}}{2}\) ; \( \omega_k^{(j,j')} = \frac{\chi_{j,k}-\chi_{j',k}}{2}\) ; the indicator function \( \mathbf{1}_{k,k'}\) equals 1 if \( k=k'\) and 0 otherwise; \( \bar{a}\) denotes the complex conjugate of \( a\) (without transpose, for vectors and matrices, i.e. \( \bar{a}^\dagger\) is the transpose of \( a\) ); \( \mathfrak{R}\) and \( \mathfrak{I}\) respectively denote real and imaginary part; and the measurement results are described by

  \[ \begin{equation} dy_{t,l} = 2 \sqrt{\eta_l}\, s_l\, dt + dw_{t,l} \; = \; 2 \sqrt{\eta_l}\, \sum_k (\; \mathfrak{R}(\beta_{l,k}) \text{Tr}(\rho_t \, \mathbf{X}_k) - \mathfrak{I}(\beta_{l,k}) \text{Tr}(\rho_t \, \mathbf{P}_k) \; ) \, dt + dw_{t,l} \; . \end{equation} \]

  (31)

  In a more compact notation, we have:

  \[ \begin{eqnarray} \nonumber d\mathcal{W}^{(j,j')}_t &=& \left(\, u^\dagger (\mathbf{I}_n \otimes J) + q^\dagger (\delta^{(j,j')} \otimes J) \,\right)\, \text{grad}_q \, \mathcal{W}^{(j,j')}_t\,dt \\ \nonumber && \; -i\,\left(\, q^\dagger (2\omega^{(j,j')} \otimes \mathbf{I}_2) q - \text{Tr}(\omega^{(j,j')}) - \tfrac{1}{8}\text{Tr}( (\omega^{(j,j')} \otimes \mathbf{I}_2)\, \text{hess}_{q} )\, \right) \, \mathcal{W}^{(j,j')}_t\,dt \\ \nonumber && + \left( \frac{1}{2}\text{Tr}(\tilde{R}) + \frac{1}{2} q^\dagger \tilde{R}\, \text{grad}_q + \frac{1}{8}\text{Tr}( \tilde{R}\, \text{hess}_q) \right) \, \mathcal{W}^{(j,j')}_t\,dt \\ && + \sum_{l=1}^m \, 2 \sqrt{\eta_l} \left( \beta_{l,A}^\dagger q + \tfrac{1}{4}\beta_{l,A}^\dagger\, \text{grad}_q - s_l \right) \mathcal{W}^{(j,j')}_t \; dw_{t,l} \; ; \end{eqnarray} \]

  (32)

  \[ \begin{eqnarray} dy_{t,l} &=& 2 \sqrt{\eta_l}\, \beta_{l,A}^\dagger \, Q_{t} \; dt + dw_{t,l} \; . \\\end{eqnarray} \]

  (33)

  Here \( J=[0 \; ,\; 1 \; ;\; -1 \; ,\; 0 ] \in \mathbb{R}^{2 \times 2}\) ; the diagonal matrices \( \delta^{(j,j')}\) and \( \omega^{(j,j')}\) \( \in \mathbb{R}^{n \times n}\) contain e.g. \( \delta_k^{(j,j')}\) as component \( k\) ; the vector \( u \in \mathbb{R}^2\) contains the control inputs in natural order; \( \text{grad}_q\) and \( \text{hess}_q\) respectively denote the gradient and the hessian with respect to \( q\) ; and column vectors \( \beta_{l,A} = [\mathfrak{R}(\beta_{l,1}),\; -\mathfrak{I}(\beta_{l,1}),\; \mathfrak{R}(\beta_{l,2}),\; -\mathfrak{I}(\beta_{l,2}),\; ... ]\) , \( Q_{t} = [\,\text{Tr}(\rho_t\, \mathbf{X}_1),\;\text{Tr}(\rho_t\, \mathbf{P}_1),\;\text{Tr}(\rho_t\, \mathbf{X}_2),\;\text{Tr}(\rho_t\, \mathbf{P}_2),\;... \,] \; \in \mathbb{R}^{2n}\) . We have further defined

  \[ \begin{eqnarray*} \tilde{R} &=& \sum_{l=1}^m \; R^{(re)}_l \otimes I_2 + R^{(im)}_l \otimes J ~ \in \mathbb{R}^{2n \times 2n}\\ && \text{with } R^{(re)}_{l,(k,k')} = \mathfrak{R}(\beta_{l,k}\bar{\beta}_{l,k'}) \;\; \text{ and } \;\; R^{(im)}_{l,(k,k')} = \mathfrak{I}(\beta_{l,k}\bar{\beta}_{l,k'}) \; . \end{eqnarray*} \]

  Note that \( \tilde{R}\) is real symmetric and \( \text{Tr}(\tilde{R}) = 2 \sum_l \beta_{l,A}^\dagger \beta_{l,A}\) . For the single cavity setting of Section 6.3.1, we just have \( \tilde{R}=2\, \mathbf{I}_{2n}\) . For further compactness, we will use below the notation \( \tilde{J} = (\mathbf{I}_n \otimes J)\) , \( \tilde\delta^{(j,j')} = (\delta^{(j,j')} \otimes J)\) and \( \tilde\omega^{(j,j')} = -i(\omega^{(j,j')} \otimes \mathbf{I}_2)\) . All these matrices belong to \( \mathbb{C}^{2n \times 2n}\) , the first two are real skew-symmetric and the last one is imaginary diagonal; they are thus all skew-hermitian, i.e. generators of rotations, as one would guess from the parameters they involve.

• Describing a function of \( 2n\) variables involves, in full generality, a number of variables which is exponential in \( n\) . In the present case, thanks to the low-dimensional confinement, we can reduce this complexity essentially to the initial state only: from there, the evolution in time involves a Gaussian and thus a number of variables at most quadratic in \( n\) . Moreover, in agreement with the algebraic criterion, the number of variables influence by the stochastic measurement results is only linear in \( n\) . Explicitly, we write

  \[ \begin{eqnarray} \mathcal{W}^{\{\rho_{t,(j,j')}\}}(q) &=& \frac{\mu_t^{(j,j')}}{\pi^{n}\text{det}(S_t^{(j,j')})^{1/2}}\int \mathcal{W}^{\{\rho_0^{(j,j')}/\mu_0^{(j,j')}\}}(q_0) \, K^{(j,j')}_{t}(q,q_0)\, dq_0 \; , \end{eqnarray} \]

  (34)

  \[ \begin{eqnarray} K^{(j,j')}_{t}(q,q_0) &=& \exp\Bigg[-\left(q\text{-}\bar{\gamma}^{(j,j')}_t\text{-}\bar{M}^{(j,j')}_t q_0\right)^\dagger\; (S_t^{(j,j')})^{-1}\; \left(q\text{-}\gamma^{(j,j')}_t\text{-}M^{(j,j')}_t q_0\right) \\ \nonumber && \phantom{KKKKKKKKK} \; + q_0^\dagger D_t^{(j,j')} q_0 + (\bar{\lambda}^{(j,j')})^\dagger q_0 \; \Bigg] \; . \\\end{eqnarray} \]

  The time-dependent parameters whose equations we seek are scalars \( \mu_t^{(j,j')} \in \mathbb{C}\) ; vectors \( \gamma_t^{(j,j')}\) and \( \lambda_t^{(j,j')}\) \( \in \mathbb{C}^{2n}\) ; and matrices \( S_t^{(j,j')}\) , \( M_t^{(j,j')}\) and \( D_t^{(j,j')}\) \( \in \mathbb{C}^{2n \times 2n}\) . For \( j=j'\) , these parameters are all real. For \( j \neq j'\) , they can be complex and we have introduced some (not transposed) conjugates e.g. \( \bar{\gamma}_t^{(j,j')}\) to simplify expressions below. By analogy with Section 6.3.1, we have anticipated that the variables in latin letter will evolve deterministically, while those in greek letters will involve a stochastic component.

  In order to reach the correct dimension \( 2d(n+1)-2\) for the deterministically evolving manifold containing all trajectories, we can view the \( \gamma^{(j,j)}\) as independent stochastic variables (thus \( (2n)d\) real values), as well as part of the \( \mu^{(j,j')}\) (covering the \( 2d-2\) remaining variables); all the other variables, in particular the \( \lambda^{(j,j')}\) and the \( \gamma^{(j\neq j')}\) , will be deterministically related to them.

• We next express \( d\mathcal{W}^{(j,j')}_t\) in two ways. On the left-hand side, viewing \( \mathcal{W}^{(j,j')}_t\) as a function of the Gaussian parameters, we express it with \( d\mu_{t,(j,j')}, d\gamma_t^{(j,j')}, ...\) , going up to second order Taylor expansion of the Gaussian in order to cover the Itō correction. On the right-hand side, viewing \( \mathcal{W}^{(j,j')}_t\) as a function of \( q\) , we compute the right-hand side of (32) for the particular Gaussian form (34). By separating, under the integral, the terms involving various powers of \( q, q_0\) and the deterministic and stochastic parts, we obtain the following set of equations. To reduce clutter, we drop the indices \( ^{(j,j')}\) and introduce \( Z_t=-S_t^{-1}\) :

  \[ \begin{align} q^\dagger\, & dZ_t\, q + \tfrac{1}{2} q^\dagger (Z_t+\bar{Z}_t^\dagger) d\gamma_t\, d\bar{\gamma}_t^\dagger (Z_t+\bar{Z}_t^\dagger)^\dagger \, q \end{align} \]

  (35)

  \[ \begin{align} \nonumber & \;\; = q^\dagger\, [\tilde{\delta},\, Z_t]\, q \, dt + \tfrac{1}{2} q^\dagger\, \{ \tilde{R},\, Z_t \} \, q \, dt + 2 q^\dagger \, \tilde{\omega} \, q \, dt + \tfrac{1}{8} q^\dagger\, (Z_t+\bar{Z}_t^\dagger)\, (\tilde{R} - \tilde{\omega})\, (Z_t+\bar{Z}_t^\dagger) \, q \, dt \\ \nonumber q_0^\dagger\, & dD_t \, q_0 + q_0^\dagger\, (d\bar{M}^\dagger_t Z_t M_t + \bar{M}^\dagger_t Z_t dM_t) \, q_0 + q_0^\dagger \, (\bar{M}^\dagger_t dZ_t M_t ) \, q_0 \\ + \tfrac{1}{2} & q_0^\dagger\, \left(d\lambda_t + \bar{M}_t^\dagger\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right) \left(d\bar{\lambda}_t^\dagger + d\bar{\gamma}_t^\dagger (Z_t+\bar{Z}_t^\dagger) M_t \right) \, q_0 \end{align} \]

  (36)

  \[ \begin{align} \nonumber & \;\; = \tfrac{1}{8} \; q_0^\dagger\, \bar{M}^\dagger \, (Z_t+\bar{Z}_t^\dagger) \, (\tilde{R} - \tilde{\omega}) \, (Z_t+\bar{Z}_t^\dagger) \, M \, q_0 \, dt\\ \nonumber q_0^\dagger\, & \bar{M}_t^\dagger (dZ_t+d\bar{Z}_t^\dagger) \, q + q_0^\dagger\, d\bar{M}_t^\dagger (Z_t+\bar{Z}_t^\dagger) \, q + q_0^\dagger\, \left(d\lambda_t + \bar{M}_t^\dagger\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right)d\bar{\gamma}_t^\dagger (Z_t+\bar{Z}_t^\dagger) \, q \\ & \;\; = -q_0^\dagger\, \bar{M}_t^\dagger (Z_t+\bar{Z}_t^\dagger)\, (\tilde{\delta}-\tfrac{1}{2}\tilde{R} \, q \, dt + \tfrac{1}{4} q_0^\dagger\, \bar{M}_t^\dagger (Z_t+\bar{Z}_t^\dagger)\,(\tilde{R}-\tilde{\omega})\,(Z_t+\bar{Z}_t^\dagger)\, q \, dt \end{align} \]

  (37)

  \[ \begin{align} \nonumber q^\dagger\, & (Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,\text{stochastic}} \\ & \;\; = -2 q^\dagger \, (\mathbf{I}_{2n}+\tfrac{Z_t+\bar{Z}_t^\dagger}{4})\; {\textstyle \sum_{l=1}^m} \; \sqrt{\eta_l} \, \beta_{l,A} \; dw_{t,l} \end{align} \]

  (38)

  \[ \begin{align} \nonumber q_0^\dagger\, & \left(\, d\lambda_{t,\text{stochastic}}+\bar{M}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,\text{stochastic}}\, \right) \\ & \;\; = -2 q_0^\dagger \, \bar{M}^\dagger \, \tfrac{Z_t+\bar{Z}_t^\dagger}{4} \; {\textstyle \sum_{l=1}^m} \; \sqrt{\eta_l} \, \beta_{l,A} \; dw_{t,l} \end{align} \]

  (39)

  \[ \begin{align} \nonumber & \frac{d\mu_{t,\text{stochastic}}}{\mu_t}\, + \bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,\text{stochastic}} \\ & \;\; = -2 {\textstyle \sum_{l=1}^m} \; \sqrt{\eta_l} \, \left( \bar{\gamma}_t^\dagger \, \tfrac{Z_t+\bar{Z}_t^\dagger}{4} \, \beta_{l,A} + s_l \right) \; dw_{t,l} \end{align} \]

  (40)

  \[ \begin{align} \nonumber q^\dagger\, & (dZ_t+d\bar{Z}_t^\dagger)\, \gamma_t + q^\dagger\,(Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,\text{det}} \\ \nonumber +& q^\dagger \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \frac{d\mu_t}{\mu_t} + q^\dagger\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t\, \bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \\ & \;\; = q^\dagger\,(Z_t+\bar{Z}_t^\dagger)\, \tilde{J} \, u_t \, dt + q^\dagger\, (\tilde{\delta}+ \tfrac{1}{2}\tilde{R})\, (Z_t+\bar{Z}_t^\dagger)\, \gamma_t \, dt \end{align} \]

  (41)

  \[ \begin{align} \nonumber & \;\; \phantom{=} \; + \tfrac{1}{4} q^\dagger\,(Z_t+\bar{Z}_t^\dagger)\, (\tilde{R}-\tilde{\omega}) \, (Z_t+\bar{Z}_t^\dagger)\, \gamma_t \, dt\\ \nonumber q_0^\dagger\, & d\lambda_{t,\text{det}} + q_0^\dagger\, \bar{M}^\dagger_t \, (dZ_t+d\bar{Z}_t^\dagger)\, \gamma_t + q_0^\dagger\, d\bar{M}_t^\dagger \, (Z_t+\bar{Z}_t^\dagger)\, \gamma_t + q_0^\dagger\, \bar{M}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,\text{det}}\\ +& q_0^\dagger \left(\bar{M}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t + d\lambda_t \right)\; \left(\tfrac{d\mu_t}{\mu_t} + \bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right) \end{align} \]

  (42)

  \[ \begin{align} \nonumber & \;\; = q_0^\dagger\, \bar{M}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, \tilde{J}\, u_t \, dt + \tfrac{1}{4}\, q_0^\dagger\,\bar{M}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\,(\tilde{R}-\tilde{\omega}) \, (Z_t+\bar{Z}_t^\dagger)\, \gamma_t\\ \nonumber & \frac{d\mu_{t,\text{det}}}{\mu_t} + \frac{d(\text{det}(Z_t)^{1/2})}{\text{det}(Z_t)^{1/2}} + \bar{\gamma}^\dagger_t\, dZ_t \, \gamma_t + \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_{t,det} \\ \nonumber + & d \bar{\gamma}^\dagger_t\, Z_t \, d\gamma_t + \tfrac{d\mu_t}{\mu_t} \, \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t + \tfrac{1}{2} \left( \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right)^2 \\ & \;\; = \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, \tilde{J}\, u_t\, dt + \tfrac{1}{8} \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, (\tilde{R}-\tilde{\omega})\,(Z_t+\bar{Z}_t^\dagger)\, \gamma_t \, dt \end{align} \]

  (43)

  \[ \begin{align} & \;\; \phantom{=} \; + \tfrac{1}{2}\, \text{Tr}\left(\, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\,(\tilde{R}-\tilde{\omega}) \right)\, dt \; . \\\end{align} \]

• Using (38)-(40), we express all the quadratic terms appearing on the left-hand side of the other equations, i.e. the Itō terms:

  \[ \begin{align} & (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \, d\bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger) = 4 (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, B \, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, dt \end{align} \]

  (44)

  \[ \begin{align} & \left(\, d\lambda_t +\bar{M}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t\, \right) \left( d\bar{\lambda}^\dagger_t + d\bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, M_t \right) = 4 \bar{M}^\dagger_t\, \tfrac{Z_t+\bar{Z}_t^\dagger}{4} \, B \, \tfrac{Z_t+\bar{Z}_t^\dagger}{4} \, M_t \, dt \end{align} \]

  (45)

  \[ \begin{align} & \left(\, d\lambda_t +\bar{M}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t\, \right) \, d\bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger) = \bar{M}^\dagger_t\,(Z_t+\bar{Z}_t^\dagger)\, B \, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4}) \; dt \end{align} \]

  (46)

  \[ \begin{align} & (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \, \left( \tfrac{d\mu_t}{\mu_t} + \bar{\gamma}^\dagger_t (Z_t+\bar{Z}_t^\dagger) d\gamma_t \right) = 4\, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, Z_t\, dt \end{align} \]

  (47)

  \[ \begin{align} \nonumber & \phantom{KKKK} \text{ with } Z_t = B \, \tfrac{Z_t+\bar{Z}_t^\dagger}{4}\, \gamma_t + {\textstyle \sum_l }\; \eta_l s_l \, \beta_{l,A} \\ & \left(\bar{M}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t + d\lambda_t \right)\; \left(\tfrac{d\mu_t}{\mu_t} + \bar{\gamma}^\dagger_t \, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right) = 4 \bar{M}^\dagger_t \, \tfrac{Z_t+\bar{Z}_t^\dagger}{4} \, Z_t \, dt \end{align} \]

  (48)

  \[ \begin{align} & d\bar{\gamma}^\dagger_t \, Z_t \, d\gamma_t = 2\, \text{Tr}\left(\, (Z_t+\bar{Z}_t^\dagger)^{-1} \, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, B\, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, \right) \, dt \end{align} \]

  (49)

  \[ \begin{align} & \tfrac{d\mu_t}{\mu_t} \, \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t + \tfrac{1}{2} \left( \bar{\gamma}^\dagger_t\, (Z_t+\bar{Z}_t^\dagger)\, d\gamma_t \right)^2 \end{align} \]

  (50)

  \[ \begin{align} & \phantom{KKK} = 4 \bar{\gamma}^\dagger_t\, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, Z_t \, dt - 2 \bar{\gamma}^\dagger_t (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, B \, (\mathbf{I}_{2n} + \tfrac{Z_t+\bar{Z}_t^\dagger}{4})\, \gamma_t \, dt \\\end{align} \]

  Here we have introduced \( B = \sum_{l=1}^m \eta_l \beta_{l,A}\beta_{l,A}^\dagger \; \in \mathbb{R}^{2n \times 2n}\) . For the single cavity setting of Section 6.3.1, we just have \( B = \eta\, \mathbf{I}_{2n}\) .
  

We solve these equations sequentially.

• First, from (35) and (44), we obtain an autonomous, deterministic equation for \( Z_t\) , and thus for \( S_t\) in the original formulation. Decomposing it in a symmetric part \( Z_S=(Z+\bar{Z}^\dagger)/2\) and a skew-symmetric one \( Z_A=(Z-\bar{Z}^\dagger)/2\) , we get:

  \[ \begin{eqnarray} \tfrac{d}{dt} Z_S &=& [\tilde{\delta},\, Z_S] + \tfrac{1}{2} \{ \tilde{R},\, Z_S\} + 2 \tilde\omega + \tfrac{1}{2} Z_S (\tilde{R}-\tilde{\omega}) Z_S - 2 (\mathbf{I}_{2n}+\tfrac{Z_S}{2}) B (\mathbf{I}_{2n}+\tfrac{Z_S}{2}) \; , \end{eqnarray} \]

  (51)

  \[ \begin{eqnarray} \tfrac{d}{dt} Z_A &=& [\tilde{\delta},\, Z_A] + \tfrac{1}{2} \{ \tilde{R},\, Z_A\} \; . \end{eqnarray} \]

  (52)

  Starting with \( Z_{A,0}=0\) at \( t=0\) , we will have \( Z_{A,t}=0\) for all times. We can thus identify \( Z_t = Z_{S,t}\) . We are not going to explicitly integrate (51) here, but one can note that it asymptotically converges (modulo generic appropriate observation conditions) towards the stationary point \( Z_t = -2 \, \mathbf{I}_{2n}\) . We report below its translation back to \( S_t\) .

• Next, replacing the previous point in (37) and using (46), we obtain a deterministic equation for \( M_t\) , depending on \( S_t\) :

  \[ \tfrac{d}{dt} M_t = \tilde{\delta}\, M_t - 2 Z_t^{-1} \tilde{\omega}\, M_t - \tfrac{1}{2} \tilde{R} \, M_t + 2 Z_t^{-1} (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, B \, M_t \; . \]

  This is a linear equation, but nontrivial, since \( M_{0} \neq 0\) and \( Z_t\) is time-varying. Once \( Z_t\) has converged to \( Z_t = -2 \, \mathbf{I}_{2n}\) (if we can say so, since there is no reason for such timescale separation), it would become a time-independent linear equation, involving rotation of \( M_t\) with \( \tilde{\delta},\tilde{\omega}\) and exponential decay with rate \( \tilde{R}/2\) .

• Next, we can replace the preceding results in (36), and using (45), we obtain \( D_t\) as an explicit integral on \( M_t\) :

  \[ \tfrac{d}{dt} D_t = 2 \bar{M}^\dagger_t\, \tilde{\omega}\, M_t - 2 \bar{M}^\dagger_t\, B \, M_t \; . \]

  This closes a first set of deterministic variables.

• After inserting the corresponding Itō terms into (38),(41), we get the stochastic equation:

  \[ \begin{eqnarray} d\gamma_t &=& \tilde{J}\, u_t \, dt + \tilde{\delta} \, \gamma_t \, dt - \tfrac{1}{2} \tilde{R} \, \gamma_t \, dt - 2 Z_t^{-1} \tilde{\omega}\, \gamma_t\, dt \end{eqnarray} \]

  (53)

  \[ \begin{eqnarray} && + Z_t^{-1}\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, \left( 2\, B \gamma_t\, dt - {\textstyle \sum_l }\; \sqrt{\eta_l} \, \beta_{l,A} \, dy_{t,l} \right) \; , \\\end{eqnarray} \]

  where thus \( dy_{t,l}\) follows (31). These variables are too numerous to be all independent, according to the manifold dimension expected from the algebraic criterion. In fact, each \( \gamma_t^{(j \neq j')}\) can be expressed as a deterministic function of \( \gamma_t^{(j,j)}\) and \( \gamma_t^{(j',j')}\) \( \in \mathbb{R}^{2n}\) , such that those variables together contribute \( 2dn\) dimensions to the confining manifold. The deterministic link goes as follows.

  First, define \( L_t^{(j,j')} = Z_t^{(j,j')}\, \left(\mathbf{I}_{2n}+\tfrac{Z_t^{(j,j')}}{2}\right)^{-1} \, \gamma_t^{(j,j')}\) . This results in

  \[ \begin{equation} dL_t^{(j,j')} = Z_t^{(j,j')}\, \left(\mathbf{I}_{2n}+\tfrac{Z_t^{(j,j')}}{2} \right)^{-1} \tilde{J}\, u_t \, dt -{\textstyle \sum_l}\; \sqrt{\eta_l} \beta_{l,A}\, dy_{t,l} \; + \left( \tilde{\delta}^{(j,j')} - \tilde{\omega}^{(j,j')} + \tfrac{\tilde{R}}{2} \right) \, L_t^{(j,j')} \; . \end{equation} \]

  (54)

  Any linear combination \( \sum_{j,j'} \; C^{(j,j')}\; L_t^{(j,j')}\) with coefficients satisfying \( \sum_{j,j'} \; C^{(j,j')} = 0\) would cancel the explicit dependence on measurements \( dy_{t,l}\) .

  The remaining issue is to find a linear combination which also provides an autonomous deterministic term, thus avoiding indirect dependence on the measurements.
• The term in \( dy_{t,l}\) will cancel as soon as \( \sum_{j,j'} \; C^{(j,j')} = 0\) . Note that the property also holds if the \( C^{(j,j')} \in \mathbb{C}^{2n \times 2n}\) are time-dependent.
• The term in \( \tilde{J}\, u_t \, dt\) will be an autonomous deterministic contribution, as long as the \( C^{(j,j')}_t\) only depend on deterministic variables.
• The term in \( \tilde{R}\) is independent of \( (j,j')\) . Hence, if we take \( X^{(j,j')}\) a scalar linear combination of the \( L^{(j,j')}\) , then it would lead to a term proportional to \( \tilde{R}\; X^{(j,j')}\) , which is compatible with an autonomous equation for \( X^{(j,j')}\) . If the \( C^{(j,j')}\) are matrix coefficients, then the result of their commutation with \( \tilde{R}\) will have to be investigated.
• Unfortunately, the term in \( \tilde{\delta}^{(j,j')}-\tilde{\omega}^{(j,j')}\) is not of this type, calling for a more complicated solution. However, we observe that (i) these are all rotations and (ii) if we take the triplet \( (j\neq j'),\; (j,j),\; (j',j')\) , then their corresponding coefficients are linearly related.

This last observation is the key towards proposing two things. First, we will search for deterministic variables as a linear combination of said triplet, thus

Second, we will decompose \( \tilde{\delta}^{(j,j')}-\tilde{\omega}^{(j,j')}\) into a mean contribution and a deviation:

The first term on the right hand side is the same for all three cases, it thus plays a role similar to \( \tilde{R}\) . The second term depends on the row of (56), but it generates a rotation; moreover, all these terms commute since \( \delta^{(j,j')}\) and \( \omega^{(j,j')}\) are diagonal, and we recall \( \tilde{\delta}^{(j,j')} = \delta^{(j,j')} \otimes J\) . This suggests to cancel the row-dependent parts of (56) by going to a corresponding rotating frame, which would be transparent for the first, common term. The algebraic criterion, fixing the manifold dimension, encourages us the other terms should behave favorably.

We thus write:

The constraint \( \sum_{j,j'} \; C_t^{(j,j')} = 0\) for all \( t\) then imposes, up to pre-multiplying by the same constant matrix:

One checks that these \( C^{(j,j',k)}_0\) , although not full-rank, commute with \( A \otimes J\) for any \( A\) , and of course with \( A \otimes \mathbf{I}_2\) . They thus maintain the property

Finally, similar commutation properties, as well as the property

lead to the result:

In hindsight, (59) shows that \( \tilde{X}_t = e^{-\tilde{\omega}^{(j,j')} t} \, X_t = \; C^{(j,j',1)}_0\, L^{(j,j)}_t + C^{(j,j',2)}_0 \, L^{(j',j')}_t + C^{(j,j',3)}_0 \, L^{(j,j')}_t \; \) in fact undergoes deterministic dynamics. Taking all things together, we have

as a deterministic variable, which together with the knowledge of \( \gamma^{(j,j)}_t\) and \( \gamma^{(j',j')}_t\) define:

• Next, now again separately for each \( (j,j')\) , consider (39),(42), in which \( d\gamma_t\) can be inserted from the last point. This readily yields

  \[ \begin{equation} d\lambda_t = 4 \bar{M}_t^\dagger\, \tilde{\omega} \, \gamma_t\, dt - 2 \bar{M}_t^\dagger\, \left( 2\, B \gamma_t\, dt - {\textstyle \sum_l }\; \sqrt{\eta_l} \, \beta_{l,A} \, dy_{l,t} \right) \; . \end{equation} \]

  (61)

  We recognize that the measurement-dependent term can be eliminated in a linear combination with \( d\gamma\) : defining

  \[ h_t = \gamma_t + Z_t^{-1}\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, \tfrac{(\bar{M}_t^\dagger)^{-1}}{2} \, \lambda_t \; , \]

  we obtain the autonomous deterministic evolution:

  \[ \tfrac{d}{dt} h_t = \tilde{J} \, u_t + \tilde{\delta}\, h_t - \tfrac{1}{2} \tilde{R} \, h_t + \tilde{\omega} \, h_t \; . \]

  For \( u_t=0\) , this would involve a matrix exponential inducing constant rotation and decay, in fact the same as for \( M_t\) when \( Z_t = -2 \, \mathbf{I}_{2n}\) . In any case, it establishes that \( \lambda_t^{(j,j')}\) is a purely deterministic function of \( \gamma_t^{(j,j')}\) , i.e. once \( \gamma_t^{(j,j')}\) is given, the value of \( \lambda_t^{(j,j')}\) is independent of the measurement results.

• Last but not least, there remains to treat \( \mu_t\) . Inserting the Itō terms and the previous results into (40),(43) yields

  \[ \begin{eqnarray*} \frac{d\mu_{t}}{\mu_t} + \frac{d(\text{det}(Z_t)^{1/2})}{\text{det}(Z_t)^{1/2}} &=& 2 \sum_l \, \sqrt{\eta_l}\, (\beta_{l,A}^\dagger \gamma_t - s_l) \, dw_{t,l} \\ && + \tfrac{1}{4} \text{Tr}((\tilde{R}-\tilde{\omega})Z_t)\, dt + 2 \bar{\gamma}_t^\dagger \tilde{\omega}\, \gamma_t\, dt + \tfrac{1}{2} \text{Tr}(\tilde{R}-\tilde{\omega}) \, dt \\ && - \text{Tr}\left(\, B \, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, Z_t^{-1}\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, \right) \, dt \; . \end{eqnarray*} \]

  Inserting \( \frac{d(\text{det}(Z_t)^{1/2})}{\text{det}(Z_t)^{1/2}} = \tfrac{1}{2}\, \text{Tr}(Z_t^{-1}\, dZ_t)\) with \( dZ_t\) derived above, reduces this equation to:

  \[ \frac{d\mu_{t}}{\mu_t} = 2 \sum_l \, \sqrt{\eta_l}\, (\beta_{l,A}^\dagger \gamma_t - s_l) \, dw_{t,l} - \text{Tr}(\tilde{\omega}\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, Z_t^{-1})\, dt + 2 \bar{\gamma}_t^\dagger \, \tilde{\omega} \, \gamma_t\, dt \; . \]

  Taking into account that \( \mu^{(j,j')} = \bar{\mu}^{(j',j)}\) , these are \( d^2\) stochastic equations, whereas there should remain \( 2d-2\) independent stochastic variables. One way of parameterizing this is to consider the \( \mu^{(j=j')}\) as independent stochastic, thus following (recall that \( \tilde{\omega}^{(j,j')}=0\) for \( j=j'\) ):

  \[ \frac{d\mu^{(j,j)}_{t}}{\mu^{(j,j)}_t} \;=\; 2 \sum_l \, \sqrt{\eta_l}\, (\beta_{l,A}^\dagger \gamma^{(j,j)}_t - s_l) \, dw_{t,l} \; = \; 2 \sum_l \, \sqrt{\eta_l}\, \beta_{l,A}^\dagger (\gamma^{(j,j)}_t - Q_{t}) \, dw_{t,l}\; . \]

  These contribute \( d-1\) independent variables, as they are subject to one constraint fixing the total normalization. Like for the single-oscillator case, we see that \( \mu_t^{(j,j)}\) varies more if its Gaussian position \( \gamma^{(j,j)}\) deviates more from the total state mean positions \( Q_t\) .

  Next, we observe that the stochastic term comports one part that is independent of \( (j,j')\) — and could thus cancel when taking linear combinations of \( \mu^{(j,j')}\) — while the other part involves \( \gamma^{(j,j)}_t\) . To make such a term appear from another stochastic variable, we consider \( \Gamma_t = \bar{\gamma}_t^\dagger \, Z_t\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})^{-1}\, \gamma_t\) . A few computations then indeed yield:

  \[ \begin{eqnarray} d(\log(\mu)+\Gamma) &=& - 2 \sum_l\, \sqrt{\eta_l}\, s_l\, dw_{t,l} - 2 \sum_l \, \eta_l \, s_l^2 \, dt \end{eqnarray} \]

  (62)

  \[ \begin{eqnarray} && + \text{Tr}\left( B\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, Z_t^{-1}\, \right) \, dt - \text{Tr}(\tilde{\omega}\, (\mathbf{I}_{2n}+\tfrac{Z_t}{2})\, Z_t^{-1})\, dt \\ \nonumber && + 2\, L_t^\dagger \, \tilde{J} \, u_t \, dt \; , \\\end{eqnarray} \]

  with \( L_t\) essentially related to \( \gamma_t\) and following the stochastic equation (54). The first row of (62) is now totally independent of \( (j,j')\) . The second row only contains deterministic variables. In the last row, the presence of \( L_t\) still has to be handled. The idea, in order to obtain purely deterministic variables, is to take linear combinations of (62) for various \( (j,j')\) : its first row would disappear, while the last row would involve the deterministic variables \( \tilde{X}^{(j,j')}\) .

  For \( u_t=0\) , this is simple: the last row is absent, so any linear combination of \( (\log(\mu)+\Gamma)^{(j,j')}\) for various indices \( (j,j')\) would yield a purely deterministic equation, involving only the second row of (62), as soon as the coefficients of the linear combination sum to zero. This also shows that the \( \mu_t^{j,j}\) are not stochastically independent either in this case, as we have already seen for the single oscillator model. For \( u_t \neq 0\) , we have to make this slightly more complicated.

  First, consider the real part of (62). We would define

  \[ \begin{equation} Y^{(j,j')}_{R,t} = \tilde{c}^{(j,j',1)} \mathfrak{R}(\log(\mu)+\Gamma)^{(j,j)}+\tilde{c}^{(j,j',2)} \mathfrak{R}(\log(\mu)+\Gamma)^{(j',j')}+\tilde{c}^{(j,j',3)} \mathfrak{R}(\log(\mu)+\Gamma)^{(j,j')} \end{equation} \]

  (63)

  with scalars \( \tilde{c}^{(j,j',k)} \) satisfying \( \sum_k \, \tilde{c}^{(j,j',k)} = 0\) , in order for the first row of (62) to cancel when writing \( dY^{(j,j')}_{R,t}\) . We further note that all variables associated to \( j'=j\) are real, so from (58) we can write for the last row of (62):

  \[ \begin{eqnarray*} \tfrac{-1}{2}\, \mathfrak{R}\left( 2\, (L_t^{(j,j)})^\dagger \, \tilde{J} \, u_t \right) &=& \mathfrak{R}\left( 2\, (C_0^{(j,j',1)}\, L_t^{(j,j)})^\dagger \, \tilde{J} \, u_t \right)\; ,\\ \tfrac{-1}{2}\, \mathfrak{R}\left( 2\, (L_t^{(j',j')})^\dagger \, \tilde{J} \, u_t \right) &=& \mathfrak{R}\left( 2\, (C_0^{(j,j',2)}\, L_t^{(j',j')})^\dagger \, \tilde{J} \, u_t \right) \; , \end{eqnarray*} \]

  while \( C_0^{(j,j',3)}= \mathbf{I}_{2n}\) trivially implies:

  \[ \mathfrak{R}\left( 2\, (L_t^{(j,j')})^\dagger \, \tilde{J} \, u_t \right) = \mathfrak{R}\left( 2\, (C_0^{(j,j',3)}\, L_t^{(j,j')})^\dagger \, \tilde{J} \, u_t \right) \; . \]

  Therefore, taking \( \tilde{c}^{(j,j',1)}= \tilde{c}^{(j,j',2)}=-1/2\) and \( \tilde{c}^{(j,j',3)}=1\) in (63) is compatible with a fully deterministic evolution:

  \[ \begin{eqnarray*} \tfrac{d}{dt} Y^{(j,j')}_{R,t} &=& \frac{-1}{2}\, \text{Tr}\left( B\, (\mathbf{I}_{2n}+\tfrac{Z^{(j,j)}_t}{2})\, (Z_t^{(j,j)})^{-1}\, \right) - \frac{1}{2}\, \text{Tr}\left( B\, (\mathbf{I}_{2n}+\tfrac{Z^{(j',j')}_t}{2})\, (Z_t^{(j',j')})^{-1}\, \right) \\ && + \mathfrak{R}\left(\, \text{Tr}\left( \,(B-\tilde{\omega}^{(j,j')})\, (\mathbf{I}_{2n}+\tfrac{Z^{(j,j')}_t}{2})\, (Z_t^{(j,j')})^{-1}\, \right) \,\right)\\ && + 2\, \mathfrak{R}(\tilde{X}_t^{(j,j')})^\dagger \, \tilde{J} \, u_t \, . \end{eqnarray*} \]

  Since \( \mathfrak{R}\log(\mu) = \log(|\mu|)\) , this shows how the weights \( |\mu_t^{(j\neq j')}|\) are deterministically linked to the \( \mu_t^{(j,j)}\) and \( \gamma_t^{(j,j)}\) .

  There remains to find deterministic variables among the \( \mathfrak{I}\log(\mu^{(j,j')}) = \text{phase}(\mu^{(j,j')})\) . Taking the imaginary part of (62), the first row drops rightaway. Furthermore, from (58), we note that for any triplet \( j,j',j''\) one has:

  \[ \mathfrak{I}(L^{(j,j')}_t + L^{(j',j'')}_t + L^{(j'',j)}_t ) = \mathfrak{I}(\tilde{X}^{(j,j')}_t + \tilde{X}^{(j',j'')}_t + \tilde{X}^{(j'',j)}_t )\; . \]

  Therefore, we can define

  \[ \begin{equation} Y^{(j,j',j'')}_{I,t} = \mathfrak{I}(\log(\mu)+\Gamma)^{(j,j')}+ \mathfrak{I}(\log(\mu)+\Gamma)^{(j',j'')}+ \mathfrak{I}(\log(\mu)+\Gamma)^{(j'',j)} \end{equation} \]

  (64)

  and obtain its deterministic evolution:

  \[ \begin{eqnarray*} \tfrac{d}{dt}\, Y^{(j,j',j'')}_{I,t} &=& \mathfrak{I}\left(\, \text{Tr}\left( \,(B-\tilde{\omega}^{(j,j')})\, (\mathbf{I}_{2n}+\tfrac{Z^{(j,j')}_t}{2})\, (Z_t^{(j,j')})^{-1}\, \right) \,\right)\\ && + \mathfrak{I}\left(\, \text{Tr}\left( \,(B-\tilde{\omega}^{(j',j'')})\, (\mathbf{I}_{2n}+\tfrac{Z^{(j',j'')}_t}{2})\, (Z_t^{(j',j'')})^{-1}\, \right) \,\right)\\ && + \mathfrak{I}\left(\, \text{Tr}\left( \,(B-\tilde{\omega}^{(j'',j)})\, (\mathbf{I}_{2n}+\tfrac{Z^{(j'',j)}_t}{2})\, (Z_t^{(j'',j)})^{-1}\, \right) \,\right)\\ && + 2 \, \mathfrak{I}(\tilde{X}^{(j,j')}_t + \tilde{X}^{(j',j'')}_t + \tilde{X}^{(j'',j)}_t )^\dagger \, \tilde{J} \, u_t \, . \end{eqnarray*} \]

  This expresses how the phases of the \( \mu^{(j,j')}_t\) are deterministically linked. One can check that for instance the phases of \( \mu^{(1,2)}_t,\,\mu^{(1,3)}_t,...,\mu^{(1,d)}_t\) can be chosen independently, then constraining all the other phases through the \( Y^{(j,j',j'')}_{I,t}\) , in agreement with the fact that there remained to find \( (d-1)\) stochastic variables and all other deterministic.

This concludes the computations. There remains to define the initial conditions, by imposing the Gaussian kernel to coincide with the Dirac-peak limit at \( t=0\) . We have thus established the following result.

Specializing these to a single cavity with \( B = \eta\, \mathbf{I}_{2n}\) and \( \tilde{R}=2\, \mathbf{I}_{2n}\) (and \( \Delta_k=0\) ) yields back the differential equations behind the solutions of Proposition 11. Some of the above expressions are undefined when \( S_t - \mathbf{I}_{2n}\, /2\) becomes singular; the reader is invited to go back through the derivations above in order to smoothen out these points.

Note that each \( R_l^{(re)}\) above is positive semidefinite, while \( q^\dagger (R_l^{(im)}\otimes J) \, q = \overline{q^\dagger (R_l^{(im)\dagger}\otimes J^\dagger) \, q } = \overline{q^\dagger (R_l^{(im)}\otimes (-J)) q}= -q^\dagger (R_l^{(im)}\otimes J) \, q = 0\) , proving that \( \tilde{R}\) is positive semi-definite. It governs the rate at which \( S_t\) converges towards \( \mathbf{I}_{2n}\, /2\) and \( M_t\) converges towards \( \mathbf{0}_{2n}\) .

Assuming that this has happened (e.g. \( \tilde{R}\) positive definite), the Gaussian kernel factorizes into a function of \( q_0\) and a Gaussian in \( q\) . The former integrates out with the initial state as a constant, and the former expresses the product of coherent states towards which the cavity has converged, indexed by \( j,j'\) . In this limit, from (53), the \( \gamma^{(j',j)}\) undergo no stochastic motion anymore: they follow deterministic dynamics which depend on \( (j,j')\) but not on the measurement results (nor, consistently, on the \( \eta_l\) ). These correspond to the coherent state amplitudes. The \( D^{(j,j')}\) and the \( \lambda^{(j,j')}\) (see (61) to resolve the undefiniteness) don't vary anymore either and the relative weight \( \text{Tr}(\rho_{t,(j,j)})\) of qudit level \( j\) is governed by \( \mu_{t}^{(j,j)}\) . Like for the single oscillator case, we can note \( \tilde{\mu}_{t}^{(j,j)} = \text{Tr}(\rho_{t,(j,j)})\) with \( \tfrac{d\tilde{\mu}}{\tilde{\mu}} = \tfrac{d\mu}{\mu}\) , such that \( Q_t = \sum_{j'} \tilde{\mu}_t^{(j',j')} \, \gamma_t^{(j',j')}\) while \( \gamma_t^{(j,j)} = \gamma_t^{(j,j)} \, \sum_{j'} \tilde{\mu}_t^{(j',j')}\) . From these observations, we obtain:

This generalization of (27) again matches the form (9). It expresses how the setting has become equivalent, in the limit, to a direct QND measurement of the \( \tilde{\mu}_t^{(j,j)}\) , with measurement operators \( L_l(t) = \sum_{j=1}^d \, |j\rangle\langle j|\,(\beta_{l,A}^\dagger \gamma_t^{(j,j)})\) .