Non-Markovian Quantum Dynamics in Strongly Coupled Multimode Cavities Conditioned on Continuous Measurement

An important challenge in non-Markovian open quantum systems is to understand what information we gain from continuous measurement of an output field. For example, atoms in multimode cavity QED systems provide an exciting platform to study many-body phenomena in regimes where the atoms are strongly coupled amongst themselves and with the cavity, but the strong coupling makes it complicated to infer the conditioned state of the atoms from the output light. In this work we address this problem, describing the reduced atomic state via a conditioned hierarchy of equations of motion, which provides an exact conditioned reduced description under monitoring (and continuous feedback). We utilise this formalism to study how different monitoring for modes of a multimode cavity affects our information gain for an atomic state, and to improve spin squeezing via measurement and feedback in a strong coupling regime. This work opens opportunities to understand continuous monitoring of non-Markovian open quantum systems, both on a practical and fundamental level.


I. INTRODUCTION
In quantum mechanics, a continuous measurement intrinsically influences the dynamics of a system such that its state depends on the particular measurement record [1][2][3][4][5][6][7][8][9]. Theoretical techniques to describe this conditioned state of the system have been well established [10][11][12], and play a central role in quantum control theory. In particular, the measurement record can be used as a feedback signal to drive the system in a desired state [13,14], counter decoherence [15] or increase entanglement [16,17]. In recent years, interest has shifted from few particle systems to the control of many-body quantum systems. For example, there are opportunities to explore particularly interesting many-body physics [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] with many atoms in optical cavities in regimes where the atoms interact strongly amongst themselves and also with the cavity mode(s), as depicted in Fig. 1. In principle, we would like to understand what information we can infer about the state of the atoms from light leaking out of the cavity. However, especially in the case of multimode cavities [18-20, 23, 25-27, 33], the problem becomes rapidly intractable as the system size and number of cavity modes grows.
In this work we present a new method to describe the information we obtain about the state of the atoms, which also contains a framework for understanding and utilising conditioned dynamics with continuous feedback. Our exact theory for the reduced state of the atoms (dashed blue box in Fig. 1) takes the form of a conditioned hierarchy of equations of motion, with the hierarchy capturing the information collectively present in the cavity mode(s). This theory thus accounts for a nontrivial interaction of the atoms with the cavity modes and vice versa, which makes the dynamics of the atoms non- gain of information about a many-body system of atoms in a multimode cavity, through detection of light leaking out of the cavity. While the combined cavity-atom system (red box) can be Markovian, it is a highly complex system with many modes. In the case of strong coupling between the atoms and the cavity, the reduced system of atoms (blue box) exhibits non-Markovian dynamics. Our formalism allows us to determine the information gained about this reduced system from the measurements, determining the measurement-conditioned dynamics, and allowing for continuous feedback.
Markovian. The efficiency of the hierarchical approach allows us to go beyond previous reduced descriptions for the state of the atoms in cases where the cavity modes can be adiabatically eliminated, addressing directly the strong coupling limit in a form that is also numerically tractable.
Although for concreteness we discuss the example of cavity QED, our theory can be applied to arbitrary quantum systems coupled to a non-Markovian bath which can be split into two parts with a larger Markovian open system in front of the detector. This includes networks of standard single-mode cavity QED, coupled resonator arrays and circuit QED [34][35][36], systems of electrons coupled to damped phonons (i.e., a dissipative Hubbard-Holstein model [37]) relevant for the solid-state, systems of cold atoms immersed in a BEC [38] or coupled coherently to an untrapped level [39,40], or quantum emitters coupled to plasmonic cavities [41] which can be modeled in terms of leaky quasinormal modes [42].
Furthermore, the formalism is itself an interesting extension to the general theory of open quantum systems. While established non-Markovian open system methods make it possible to compute the reduced state of the system (the atoms in Fig. 1) in the unmonitored case, making a connection to continuous measurement was possible only in some very special cases [43]. In contrast, here the exact non-Markovian dynamics of the atoms directly inherits a measurement interpretation from the full Markovian system of atoms and cavity modes. The embedding of a non-Markovian system in a larger Markovian system is a well known strategy [44][45][46][47][48][49][50][51]. However, our results provide a systematic exact theory to determine the conditioned states for the reduced system for different measurement schemes of the outgoing cavity field, connecting directly to various experimental setups.
Below in Section II, we introduce details on the cavity QED systems that we consider, and revise the corresponding continuous measurement theory. We also illustrate why in the strong coupling regime, standard adiabatic elimination of the cavity mode is not sufficient to accurately predict the dynamics of the atoms. In Section III, we present the derivation of our new exact description of the conditioned reduced atomic state, which takes the form of a conditoned hierarchy of mixed-state quantum trajectories for atom-only density matrices. In Section IV, the application of our method to multimode cavities is discussed, and we study how the information we gain about the atomic system and the resulting correlations depend on the way in which the cavity modes are monitored. In Section V, we extend the theory to include quantum feedback based on the continuous measurement, which has direct applications in quantum control. We show how squeezing of a collective spin in cavity QED can be improved via feedback, beyond the results known in an adiabatic regime [13]. Finally, in Section VI, we present our conclusions and discuss further future perspectives arising from this work.

II. CONDITIONED ATOM DYNAMICS
While the physics that we address in this work is far more general, we will use the language of cavity QED. First, we choose as a concrete example a simple and well known model: a collection of atoms A, coupled to a single mode of an optical cavity C via the Hamiltonian where H A and L are respectively the Hamiltonian and an arbitrary coupling operator for the atoms, a is the cavity mode annihilation operator, and g is the atomcavity mode coupling strength.
Multimode generalizations of this model constitute good candidates for example to explore many-body spin models from glassiness to associative memories [28] or to solve specific NP hard problems [29], as they are able to generate tunable-range interactions between atoms placed inside the cavities [25,52,53].
By way of an introduction think, however, of the simplest possible example, the Jaynes-Cummings model [54] consisting of a single two-level atom only. This is described by the choices where σ α (α = x, y, z) are the standard Pauli operators σ x = |0 1| + |1 0|, σ y = i(|0 1| − |1 0|) and σ z = |1 1| − |0 0| with |0 and |1 the two relevant atomic states, separated by the atomic transition frequency ω. This model is particularly simple because the Hamiltonian conserves the total number of atomic excitations and photons. In order to gain information on the state of the atoms via the cavity field, the outgoing cavity field can be continuously monitored, for instance with homodyne detection. The joint atom-cavity time evolution is then conditioned on the measured homodyne current J hom . For perfect detection efficiency, the detected signal can be written in the form [1] J hom dt = √ 2κ a + a † dt + dW.
Throughout the paper, the brackets ... denote the quantum expectation value with respect to the current conditioned state. The signal thus contains information on a particular quadrature of the mode. The second term that contributes to the homodyne current is white noise which arises due to the randomness of quantum measurement outcomes, written here in terms of increments dW of the Wiener process with zero mean E [dW ] = 0, obeying Ito's rule dW 2 = dt. The measurement strength κ is determined by the rate of light leaking out from the cavity mirrors. Through continuous measurement, information on the joint system of cavity field and atoms is continuously acquired. Thus, the state of atoms and field ρ AC is influenced by, or conditioned on, the measurement outcome and obeys the stochastic evolution equation [1] dρ Note that here we may in fact have pure state solutions ρ AC = |ψ AC ψ AC | (quantum trajectories, see later). To obtain the evolution of the average state, that is the average with respect to all measurement realizations, we simply have to omit the terms involving the stochastic increment dW . We are left with the standard master equation of cavity decay in Gorini-Kossakowski-Sudarshan-Lindblad (GKSL), or Lindblad, form [55,56]. The main goal of this work is to develop a theory describing the continuously monitored state of the atoms only, ρ A = tr C ρ AC , by tracing out the cavity field in (4), without further approximation. It is clear that while the monitored ρ AC may well be pure, ρ A will almost always be mixed. In the next section we will derive our exact hierarchical scheme for non-Markovian quantum trajectories of the monitored, mixed ρ A . Before doing so, we discuss two known approaches, where the desired ρ A is obtained through approximations: the simplest way is via adiabatic eliminationformally a second order perturbation theory in the coupling strength g, leading to an effective theory for the reduced state of the atoms described by the conditioned Redfield master equation Here, the operatorL is time-dependent and given as This equation is explicitly derived later in this work, and also in Ref. [49]. When the average over all measurement results is taken, i.e. the second line is neglected, equation (5) reduces to the deterministic Redfield master equation of atoms coupled to a non-Markovian reservoir with Lorentzian spectral density, reflecting the leaky cavity mode. Also, within the approximation, the measured homodyne current relates directly to the state of the atoms J hom dt ≈ −i √ 2κ g L −L † dt + dW , signaling that the cavity field is assumed to effectively follow the state of the atoms. The Redfield equation cannot capture strong memory effects which may occur in the interaction of the atoms with the cavity field. Recent work shows that for instance in the U (1)-symmetric Dicke model, a higher order perturbative evolution equation is required to correctly capture the state of the atoms and the phase transition [57], while in the standard dissipative Dicke model Redfield theory does yield accurate results [58]. Hence, even in relatively simple models, it is not obvious how to eliminate the cavity modes appropriately. We will see later that Eq. (5) will drop out naturally as only the first order approximation in our general hierarchical approach.
On top of the weak-coupling limit leading to the Redfield theory, one can make further simplifications based on assumptions on the timescales of the cavity and atom processes. A popular approximation is the 'bad cavity' limit, where the cavity decay is assumed to be the fastest timescale κ ω, ∆. Then we can approximate the Redfield operator asL ≈ g 2 L/κ which leads to an effective stochastic master equation of the same form as (4), In this approximation all memory effects of the cavity field vanish and the atoms obey GKS-Lindblad dynamics.
To give an example we consider a single realization of the Jaynes-Cummings model conditioned on a given homodyne detection signal, shown in Fig. 2 as a trajectory on or inside the Bloch sphere. The dynamics in the Jaynes-Cummings model is simply a relaxation of the atomic excitation to the ground state, which is a stationary state. The different plots in the figure compare an exact solution of the Jaynes Cummings model via master equation (4) with the Redfield approximation (5) and the bad cavity limit (7). Clearly, the bad cavity limit is not applicable in this parameter regime. Note that, in contrast to the bad cavity limit, the true conditioned atomic state does not remain pure, as indicated by a purity trρ 2 A (t) of less than one. This is because the atom-field interaction leads to finite entanglement between atom and the cavity mode. The Redfield approximation is close to the exact solution, but does not match perfectly even though, here, the cavity field can have at most a single photon occupation. This example highlights that a more systematic method is needed in order to compute the conditioned state of the atoms within an atom-only description.

III. ATOM-ONLY CONDITIONED HIERARCHICAL EQUATIONS OF MOTION
To overcome the clear limitation of an adiabtic treatment of the environmental modes, we introduce a fully non-Markovian theory for the atoms which is designed to capture the information we obtain from continuous monitoring of the environment. If the measurement outcome is averaged out, our theory reduces to the well known hierarchical equations of motion (HEOM) method, which has been used successfully for a numerical treatment of non-Markovian open quantum dynamics [59][60][61][62][63]. Our theory provides a general non-Markovian analogue to equations for the conditioned state in continuous measurement theory for Markovian systems [1]. In this analogy, the HEOM alone is the counterpart to a master equation that averages over the measurement results. In the limit of weak coupling, as well as for special integrable models, our approach reproduces the results of [49], and towards Markovian measurements with additional non-Markovian baths connects to the results of Refs. [64,65]. Our method does not rely on any assumptions on the atom-cavity coupling or the cavity quality, and can include measured and unmeasured environmental modes alike.
To sketch the derivation, we will focus on the simplest case of homodyne detection of the cavity output field from a single cavity mode only. The derivation of the various generalizations, which are presented later on, follows similar lines. For homodyne detection the joint atom and cavity state obeys the following stochastic Schrödinger equation [1,66,67], which is the pure state version of Eq. (4) 2κ a + a † dW is a factor ensuring normalization. The pure conditioned states |ψ AC are standard Markovian quantum trajectories, and have a clear physical interpretation in terms of continuous measurement.
Non-Markovian generalizations of quantum trajectories are well known in the literature, for instance the non-Markovian version of quantum jumps [68,69] and, more closely related to the results presented here, non-Markovian quantum state diffusion [70][71][72]. They have been used to compute the unmonitored dynamics of the atoms by propagating stochastic pure states. However, following [73][74][75][76][77][78], while a pure state solution of a general non-Markovian SSE can be interpreted as a conditioned state at any particular time [73,74], joining up these solutions to form a continuously monitored quantum trajectory is, up to special exceptions [43], generally impossible [76]. In our case, the joint atom-cavity quantum trajectory is Markovian, and, if pure, follows (8). We now aim to find an atom-only description of the non-Markovian dynamics of those atoms, leading to non-Markovian atomic quantum trajectories which are mixed states, and have a clear interpretation in terms of continuous measurement by construction.
In order to derive an equation for the state of the atoms, the cavity field in equation (8) must be traced out. For this we first project the equation onto a Bargmann coherent state of the cavity |y = exp(ya † )|0 , as in non-Markovian quantum state diffusion [70][71][72] or in Ref. [49], which yields To handle the derivative terms we define n-th order auxiliary states as |ψ (n) (y * , t) = (gi∂ y * ) n y|ψ AC (t) which themselves obey the coupled evolution equations reminiscent of the hierarchy of pure states (HOPS) in non-Markvoian quantum state diffusion [79,80]. Here, however, we determine conditioned atomic states under continuous (homodyne) measurement of the cavity modes. As with HOPS, the reduced state of the atoms is recovered by taking a Gaussian average over the y * variable upon acknowledging the completeness of coherent states with respect to a Gaussian measure Note that ρ A is the state of the atoms conditioned on the homodyne measurement record of the output field. Similar to [81], we can replace the hierarchy of conditioned pure states (10) by matrix hierarchichal equations of motion for the y-averaged, conditioned auxiliary matrices To find a closed evolution equation for these objects one has to employ partial integration under the Gaussian y-mean, which allows a main result of our work. Here, arising from the measured homodyne current, a term X = a + a † = tr ρ appears, that can be determined from the auxiliary matrices of first order. The zeroth order auxiliary state ρ (0,0) A is the exact physical state of the atoms conditioned on the measurement record, as can be seen from Eq. (11). As the above equation is expressed in Ito formalism, it is written in terms of the stochastic increment dW , related to the physical measurement current via Eq. (3). For completeness, this main result is presented in the Stratonovich formulation of stochastic calculus in appendix A, where the explicit dependence on the actual measurement current J hom becomes obvious.
To recover the unobserved average reduced atomic state, one additionally has to take the average over the homodyne current J hom , which amounts to taking the average with respect to the increments dW in (13): one simply has to omit all terms proportional to dW . As could be expected, this yields the standard HEOM for a quantum system in a bath with Lorentzian spectral density [81]. Here, remarkably, we obtain HEOM from an entirely different, Markovian continuous measurementbased approach.
Clearly, the derivation can be straightforwardly generalized to other measurement schemes such as direct photodetection or heterodyne detection. This is shown in appendix B, where also the corresponding hierarchical equations are provided.
While the full hierarchy (13) is in principle exact, the main practical advantage of the cHEOM arises from the fact that it can be truncated at finite order, and the consistency of that truncation can be checked: depending on the excitation of the modes, only a few auxiliary states need to be taken into account. In fact, the trace of auxiliary states gives the moments of cavity operators via a n (a † ) m (t) = trρ so that a neglect of high order auxiliary states amounts to neglecting corresponding higher order moments. As used earlier, Eq. (14) makes it possible to compute the homodyne current (3) from the first level auxiliary states. Note that in the following we consider an initial condition with no photons in the cavity, so that all higher order hierarchy states are initially zero. Their norm then increases over time, as the modes become occupied in the dynamics.
Nicely, a second order perturbation theory can be derived with ease from the full hierarchy by formally integrating the equations for the first level auxiliary states and neglecting all contributions of higher order. In this approximation, the first auxiliary states can be expressed as whereL(t) is the Redfield operator (6). Inserting this in the zeroth order equation of the hierarchy (13) gives the closed stochastic master equation for the conditioned state of the atoms (5). As expected, taking the average with respect to the increments dW results in the standard deterministic Redfield equation for a bath with Lorentzian spectral density. Equation (5) can thus be seen as a generalization of the Redfield theory to continuous measurement.
To showcase the developed cHEOM method, we go beyond the Jaynes-Cummings model (2) and include a driving term in the Hamiltonian ε 2 σ x . For a nonzero driving the number of excitations is no longer preserved and the hierarchy does not truncate at the second order. In Fig. 3 we show the convergence of the solutions with respect to the hierarchy depth, that is we simply truncate the hierarchy by setting ρ (n,m) A = 0 for n + m larger than the maximal depth k max . The figure clearly shows how a truncated cHEOM gives a systematic expansion beyond Redfield theory which converges to the exact result as the hierarchy depth is increased.
In fact, for this simple example with only a single cavity mode, solving the full Markovian stochastic master equation (4) is possible, and our hierarchical atom-only formulation is not required to numerically determine the atomic time evolution. In the following section we introduce a generalization of the hierarchy to multiple bath modes. In this case the cHEOM can open new possibilities to tackle the challenging description of cavity QED systems in multimode cavities.

IV. MULTIPLE CAVITY MODES
Exploiting the opportunities of multimode cavities or ensembles of coupled single-mode cavities experimentally has offered new fascinating possibilities for cavity QED physics simulations of many-body phenomena. The exponential size of the combined atom and cavity modes Hilbert space poses problems for a straightforward numerical treatment. Our cHEOM formalism can be generalized easily to multiple cavity modes which may or may not be continuously monitored by any of the measurement schemes discussed above, offering numerical advantages, as we shall see. The general Hamiltonian we like to consider reads Here, M cavity modes couple to the atoms with different strengths g k and possibly different coupling operators L k . For simplicity, we assume homodyne detection on all modes. In the multimode case the auxiliary density operators in the cHEOM acquire an index for each mode so that it is useful to define a vector notation where n = (n k ) and m = (m k ) are vectors of indices and w = (κ k + i∆ k ) is a vector storing the cavity mode detunings and decay rates. Then the cHEOM for the atom state conditioned on all homodyne currents J hom,k (t) reads one of the central results of our work. Here, we use the unit vectors e k = (δ kk ) and scalar prod- /(ig k ) arises from the homodyne current. In addition to taking into account multiple cavity modes with non-Markovian response, further Markovian dissipation channels of the atoms can be included simply by adding them to the hierarchy at each level. On the other hand also further non-Markovian baths could be accounted for by additionally employing any other suitable HEOM scheme, like the eHEOM method [82].
To demonstrate the applicability of our method to a nontrivial problem, we consider in the following a threemode cavity where three localized 'clusters' of atoms are trapped. Each of the identical atoms in a particular cluster interacts with the cavity field in the same way, as sketched in Fig. 4(a). Then the interaction is collective and the cluster can be described by a large spin J i of size j = N atoms /2. The atoms are assumed to interact with the individual cavity modes according to a generalized Dicke Hamiltonian which can be realized by a double Raman pumping scheme [28,83], and Markovian leakage of cavity photons with rate κ are assumed to be the only source of dissipation. Even for few cavity modes and moderate coupling, this system is numerically very challenging beyond an adiabatic regime where the cavity modes are either largely detuned or strongly damped. Thus it is a perfect setup to test the conditioned hierarchy (17).
As an interesting application, we study how monitoring of different cavity modes increases our knowledge about the state of the atoms in a single experimental run. This knowledge is quantified by the von-Neumann entropy of the state of the atoms S[ρ A ] = −trρ A ln ρ A . In case a mode of the cavity field is continuously monitored, this entropy decreases. The expected information gain from the monitoring is given by the difference of the entropy of the average state S[E [ρ A ]] and the mean entropy of the conditioned states E [S[ρ A ]] [10,84]. The latter is a nontrivial nonlinear average of conditioned states which  (19), where three clusters of atoms are trapped in a plane transverse to the cavity axis and interact with cavity modes localized around the different clusters. Such spatial-dependent intensity profiles could be obtained in a confocal cavity QED [23,25]. We assume that each cavity mode can be monitored independently. The panels (b), (c) and (d) show respectively the mean entropy, the mutual information and the negativity of the combined state of atomic cluster 1 and 3 depending on which modes of the cavity are continuously monitored via homodyning. The results are obtained from integrating 3000 samples of the conditioned HEOM (17) with parameters (19)- (21) and Natoms = 3. cannot be expressed by the reduced state, i.e. it is a property of the full ensemble and not just its mean. Thus it can only be evaluated by actually solving the full stochastic master equation for many different noise realizations. For an example calculation we choose three cavity modes and the following parameters which are far from an adiabatic regime. The coupling is chosen such that each cavity mode is localized around a single cluster in a way that the coupling strength decays towards neighboring clusters. This induces an effective interaction between the clusters, mediated by a common coupling to the cavity field. In particular, we focus on the combined state of the first and last cluster given by ρ 13 = tr 2 (ρ A ), where tr 2 denotes the trace over the second cluster. Fig. 4(b) shows the average entropy 13 ]] of the state ρ 13 depending on which modes of the cavity field are continuously monitored by homodyning. As should be, the more modes are monitored the more information is gained on the state and the lower the mean entropy. Because the atom clusters become entangled with the cavity modes and the cluster 2, the entropy is nonzero even if all modes are monitored. Note that, interestingly, monitoring mode 2, which is localized around atom cluster 2, does not give a significant information gain on the state of cluster 1 and 3, compared to measuring modes localized around the latter. As a second question of interest, we can study how the monitoring affects our knowledge on correlations between the clusters, characterized by the mutual information I 13 = S[ρ 1 ] + S[ρ 3 ] − S[ρ 13 ] between the first and the third cluster [10]. As seen in Fig. 4(c), showing the average mutual information E[I 13 ], barely any correlations between cluster 1 and 3 are present in case none of the modes is monitored. In contrast, in case both mode 1 and 3 are monitored with homodyning, significant correlations between the clusters are expected in the individual experimental realizations. If in addition also mode 2 is monitored, these correlations decrease again. Note however, that quantum correlations cannot increase from neglecting the measurement results of mode 2. This means that the increase in mutual information in the case where mode 2 is not monitored, versus the case where all modes are monitored, is exclusively due to classical correlations.
The third question we like to address is how quantum correlations between clusters 1 and 3 are affected by the monitoring. Specifically, we quantify this by the negativity N 13 ≡ N [ρ 13 ] = (||ρ Γ1 13 || 1 − 1)/2 [85], where ρ Γ1 13 denotes the partial transpose of ρ 13 with respect to the first cluster and ||O|| 1 = tr √ O † O the trace norm of an operator O. In Fig. 4(d), showing the mean negativity E[N 13 ], we observe how the completely monitored state contains the most negativity and then how it decreases when modes 1 and 3 and mode 1 alone are monitored. This is explained by the fact that quantum correlations cannot increase by classical averaging. When mode 2 alone or none of the modes are monitored, the state contains no negativity.
The effectiveness of our approach can be assessed when contrasted to the problem size when the cavity mode Hilbertspace would be directly truncated. The size of the space of the density matrices of the atomic system is r = (2j + 1) 2N clusters . Using our approach and truncating the hierarchy according to a triangular condition N modes k=1 (m k + n k ) ≤ k max where m k , n k are the hierarchy indices results to total of K auxiliary density matrices Thus, we need to solve a total of d = rK equations. A direct truncation of the Fock space of the modes at k max leads to a full state dimension of D = r(k max + 1) 2N modes . For k max = 3 we have K = 84 and this already gives two orders of magnitude reduction in the problem size d/D = 84/4096 ≈ 0.02. This reduction becomes even more dramatic when the number of modes is increased or the modes are coupled more strongly. Further improvements could be achieved with advanced truncation procedures as in Refs. [80,86].

V. FEEDBACK
From an experimental perspective, measuring single quantum trajectories is in general not feasible, as a particular trajectory cannot be prepared twice and any averaging required to reconstruct the state will reduce to the unmeasured case. This is different in case the measurement record is used as a feedback signal, which adjusts for instance the strength of a classical driving of the system. Then, average measurement results are predicted by a feedback master equation which deviates from the dynamics without feedback. Our formalism can be straightforwardly generalized to include feedback on the system based on continuous measurement. Here, we present such a theory in the simplest case of instantaneous feedback applied on the atomic system, but other generalizations such as delayed feedback can be derived similarly, following e.g. Ref. [1]. We want to consider in the following instantaneous feedback based on the homodyne signal of mode k. This corresponds to the feedback Hamiltonian where F k is an Hermitian operator which may be timedependent. We restrict ourselves to the case where the feedback is applied to the atoms, so that F k is an operator in the atom Hilbert space. It could, for instance, describe an external driving controlled by the experimenter. In the Stratonovich formulation of the conditioned HEOM provided in supplement A, the feedback can be trivially included simply by modifying the atom Hamiltonian accordingly. This can be converted to the Ito formalism which then allows to take the ensemble average. As shown in appendix C, the following contribution to the conditioned HEOM (17) arises due to the feedback on mode k: Most relevant is the feedback master equation, i.e. the equation for the averaged state. Again, in the Ito formalism, it is obtained by simply omitting all stochastic terms.
As an example, we apply the formalism to achieve unconditioned spin squeezing of an atomic system in a single-mode cavity via feedback of the homodyne current, following [13,87]. The results in Ref. [13], however, have been obtained for the bad cavity limit, where an adiabatic elimination of the cavity mode can be performed. By contrast, our approach allows us to study the case of a 'good' cavity and we demonstrate the possibility of achieving squeezing for longer times in this regime. This example shows that our formalism opens possibilities to investigate quantum state preparation in strong coupling regimes without any restrictions on cavity parameters. A detailed discussion of this example follows.
As in [13,88], we consider a system of N atoms two-level atoms (spin-1/2) collectively coupled to a single mode of a lossy cavity described by the Hamiltonian similar to the generalized Dicke Hamiltonian introduced in Eq. (18). In the following we consider N atoms = 10 and g = Ω/2. We assume that initially, the state of the atoms is a coherent spin state with all spins aligned along the direction x. The output field is measured via homodyningin a second step, the current is then fed back according to (23).
Since the coherent state is a minimum uncertainty state, the initial variances along the directions z and x are equal to j/2 = N atoms /4. As we will see below, the continuous measurement reduces the uncertainty along z, hence generating spin squeezing, which we quantify using the spin squeezing parameter [89,90] However, the spin squeezing of the conditioned states disappears after carrying out the ensemble average over all possible measurement records, because of the presence of a stochastic shift of J z which returns the unconditioned initial variance. Therefore, as in Ref. [13], one introduces a coherent feedback based on the measurement to counteract this stochastic shift, in order to maintain the spin squeezing. Our goal here is to apply this idea to achieve such spin squeezing in the good cavity limit, where the cavity mode cannot be eliminated. We will see that the feedback conditioned on the measurement will then not depend on quantities of the atoms only, but on quantities of both the atoms and the cavity, and more precisely on their correlations, as we show in the following. From the general form of the Hamiltonian as in Eq. (1) we can identify the correspondence A continuous measurement produces a stochastic shift of the z component J z . Using the conditioned HEOM equation (13), as well as (14), we find that this shift is given by where we used the quadrature X = a + a † . From the homodyne current in Eq. (3) we can extract the form of dW and use this in Eq. (28), obtaining which, using the fact that X = 0 (as shown in Figure 5), simplifies to In the absence of feedback, this stochastic shift, averaged over all the trajectories, will return a final non squeezed state, recovering the unconditioned variance for J z with the same value j/2 of the initial coherent state. This shift induced by the measurement is hence detrimental for the spin squeezing along z and must be counteracted with an additional feedback term that continuously acts on the system and induces an opposite shift of J z . Following the same protocol as in [13], we add a feedback that generates a rotation of the collective spin around the y axis: where we defined the feedback operator and where λ is the feedback strength. As the feedback Hamiltonian adds the extra terms (24) to the conditioned HEOM, the feedback induces another shift on J z For the total stochastic shift to vanish for a single trajectory, we have to impose d J z = 0. From Eq. (30) and Eq. (33) we therefore obtain the condition on the feedback strength It is worth highlighting that the feedback strength in Eq. (34) depends dynamically on conditioned expectation values of both a quantity of the atomic system alone, J x , and of a quantity that is related to the correlations between atoms and cavity, J z X, as shown in Fig. 5. Note that we can directly extract the relevant expectation value from the hierarchy (17) as In the bad cavity limit, where the cavity mode can be adiabatically eliminated, one can set optimised values for the feedback strength under the assumptions that perfect measurements on the system keep it in a pure state, so that the conditioned expectation values can be approximated with the unconditioned ones (indicated with the subscript u): J x J x u and J 2 z J 2 z u [88]. In the good cavity limit, however, the state becomes mixed, and we determine the optimal values of the feedback strength numerically. While we could choose a timedependent feedback, determining this would not be practical in an experimental setup. We therefore consider the simpler situation of constant feedback applied throughout the dynamics, using different values of the feedback strength λ within a given range, to see which values optimize the spin squeezing along z, Eq. (26). This is shown in Fig. 6 for both the good (red line) and bad (blue line) cavity limits.
First, we can observe that the good cavity case, in this regime of parameters, produces better spin squeezing. Second, while for the bad cavity regime, considering a feedback strength λ * set by the maximum of the correlations J z X is a good approximation, this is not true in the good cavity limit. In fact, if we use this approximation to set the value of the feedback strength, we do not obtain any spin squeezing, and this is due to two factors: first, the correlations are stronger and grow more slowly (as shown in Figure 5), introducing a limitation of the approximation to a constant value; second, in the good cavity limit the state of the system gets mixed quickly and the conditioned values obtained from the measurement differ from the unconditioned ones that we use to define the constant feedback strength. Fig. 6 also shows that in both the good and bad cavity limits, when we consider the full model, we have two local minima corresponding to values of λ having different signs. It is worth pointing out that if we adiabatically eliminate the cavity mode, for the same values of the parameters used here, we would obtain only one of the two minima and the physics of the second minimum would not be captured in this approximation.
The constant feedback applied with strengths set by the different values of the minima in Fig. 6 generates spin squeezing at different points in time, as can be observed from the dashed lines in Fig. 7. We can exploit this to contrast the increase of the the spin squeezing parameter after a transient [13] and decrease it further for longer times, by using a sequence of the optimal constant feedback strengths with different signs, switching from one to another at given times. The result is shown by the solid lines in Fig. 7.
Note that while our approach provides a simple physical picture of the mechanism leading to a better spin squeezing, one could potentially improve further these results by optimizing a continuous time-dependent feedback strength λ(t) via optimal control schemes. However, from an experimental point of view, this would be much harder to implement.

VI. CONCLUSIONS AND OUTLOOK
In this article we have developed a general theory describing the time evolution of a non-Markovian subsystem coupled to damped bosonic modes which are monitored continuously. This appears in many scenarios, including advanced cavity QED where multiple cavity modes couple to an interacting many-body quantum system (atoms) embedded in the cavity. We address the problem which information about the atoms can be inferred from continuous measurement of the leaking output light of the various cavity modes, and how that knowledge can be used to manipulate the many-body state using feedback. We refrain from using any of the usual approximations but succeed in formulating a theory that allows us to determine directly the exact, non-Markovian atomic (mixed) state dynamics, conditioned on the measurement record. We can thus study in detail the gradual information gain arising from monitoring more and more cavity modes. We show how the continuous observation of spatially selective cavity modes reveals information about (quantum) correlations between groups of atoms at different locations inside the cavity. Moreover, we can determine non-Markovian feedback hierarchical equations of motion as a starting point for control theory, allowing us to drive the many-body quantum system in desired states. Specifically, we consider a feedback scenario which aims to generate spin squeezing in a collective ensemble of atoms. Here, the mixedness of the conditioned atom state invalidates the previously known optimal feedback protocols of the adiabatic regime. In fact, with a more coherent cavity it is possible to generate stronger spin squeezing for significantly longer times.
Crucially, in our approach, no restrictions on coupling strengths or cavity qualities are required. Instead, our result takes the form of conditioned hierarchical equations of motion (cHEOM, Eq. (17)), proving to be an efficient scheme for tackling pressing issues in the highly complex quantum dynamics of strong atom-cavity coupling.
Our result paves the way for further exciting routes of study: our exact cHEOM could be naturally coupled to approximated schemes such as mean-field or cumulant expansion, to easily go towards larger many-body atomic systems while keeping the system-reservoir coupling exact. It could also be advantageous to formulate cHEOM with matrix product operators. In connection with the hierarchy of pure states method, recent works have demonstrated that such an ansatz can be used both to describe many-body dynamics in non-Markovian environments [37] and very strong system-bath coupling [91].
Our theory constitutes an ideal framework for the exploration of a wide range of many-body phenomena on the level of the reduced atomic non-Markovian dynamics in the strong-coupling regime, such as phase transitions [92], measurement-induced phase transitions [93,94], neural network like behaviors such as associative memories [28,95], cavity-enhanced transport [96][97][98][99][100] and superconductivity [101], continuous measurement of transport [102], or cavity cooling with higher capture range [103] and its monitoring [104]. New schemes for quantum information processing and production of entanglement [32] could also be investigated, exploiting the higher coherence achievable in the strong-coupling regime, the use of different cavity modes to realize quantum gates (or conversely the use of the atoms to realize quantum gates between photonic qubits [105,106]), or the potential of using the feedback formalism to implement error correction protocols. Finally, it is important to stress again that while we use the language of optical cavity QED, the underlying model is universal and can equally be applied to plasmonic cavities [41], cold atoms reservoirs [38][39][40], electron-phonon systems [37], or circuit QED [34,36]. Given the high cooperativity achievable in this latter platform, we expect our formalism to be indeed particularly useful for exploring control and readout of superconducting qubits. To derive the Stratonovich from of equation (13) we start from the conditioned atom and cavity evolution in Stratonovich convention [1,67] where N (t) = κ( a † a + 1/2 a 2 + (a † ) 2 ) − κ/2J hom (t) a + a † is a factor which ensures normalization. The measured homodyne current is given as where ζ(t) is a Gaussian white noise process in the Stratonovich sense with statistics E [ζ(t)] = 0 and E [ζ(t)ζ(s)] = δ(t − s). In order to obtain the hierarchy we follow the same derivation as in section III, to arrive at the Stratonovich version of Eq. (13) The advantage of the Stratonovich formulation is that it becomes obvious that the state depends on the homodyne current (A2) only. Similarly, the second order perturbation theory equation (5) is equivalent to the following Stratonovich equation which reduces to a Redfield master equation on average.

Appendix B: Different Detection Schemes
The derivations in section III can be straightforwardly applied to different measurement schemes. For example the formally very similar continuous heterodyne detection of the cavity output field yields the stochastic Schrödinger equation [1,66] This describes the evolution of atom and cavity state conditioned on the heterodyne current J het (t)dt = √ 2κ a dt + dW c , with the now complex Gaussian increment dW c , whose Another well known unraveling of master equation (16) are quantum jumps which are related to continuous direct photodetection [1,66]. In this case the conditioned evolution is piecewise deterministic until at a random time a photon is detected and the state jumps discontinuously. An evolution equation describing this can be formulated as follows (B6) One can check immediately that the last two lines drop out after taking the average with relations (B5) and the standard HEOM is recovered again. If a jump occurs then the hierarchy has to be modified according to i.e. the entire hierarchy is moving down one level.
Appendix C: Homodyne cHEOM with feedback To derive a HEOM that includes feedback we start from the SSE in Stratonovich convention as described in appendix A. There the rules for stochastic integration allow to simply add the feedback Hamiltonian (31), which describes linear feedback based on the measurement of a single mode k. We obtain As we consider only feedback on the atoms F K is an operator in the system Hilbert space. From Eq. (A2) we see that the above SSE depends on the real noise ζ(t). To convert it into Ito formalism one can then take an arbitrary basis, define ψ j = j|ψ , L j,k = j|L|k and make use of the conversion formula [67] for a Stratonovich equation of the form ∂ t ψ j = α j + β j ζ(t).
Applied to (C1) we arrive at the following SSE with feedback in Ito form: (C2) Starting from the equation above we now follow the same lines as in Sec. III, i.e. we project the equation onto a Bargmann coherent state |y , define the auxiliary states |ψ (n) = (ig∂ y * ) n y|ψ and obtain the HEOM from ρ , which is given in Eq. (24).