Feedback-stabilized dynamical steady states in the Bose-Hubbard model

The implementation of a combination of continuous weak measurement and classical feedback provides a powerful tool for controlling the evolution of quantum systems. In this work, we investigate the potential of this approach from three perspectives. First, we consider a double-well system in the classical large-atom-number limit, deriving the exact equations of motion in the presence of feedback. Second, we consider the same system in the limit of small atom number, revealing the effect that quantum fluctuations have on the feedback scheme. Finally, we explore the behavior of modest sized Hubbard chains using exact numerics, demonstrating the near-deterministic preparation of number states, a tradeoff between local and non-local feedback for state preparation, and evidence of a feedback-driven symmetry-breaking phase transition.

The implementation of a combination of continuous weak measurement and classical feedback provides a powerful tool for controlling the evolution of quantum systems. In this work, we investigate the potential of this approach from three perspectives. First, we consider a double-well system in the classical large-atom-number limit, deriving the exact equations of motion in the presence of feedback. Second, we consider the same system in the limit of small atom number, revealing the effect that quantum fluctuations have on the feedback scheme. Finally, we explore the behavior of modest sized Hubbard chains using exact numerics, demonstrating the near-deterministic preparation of number states, a tradeoff between local and non-local feedback for state preparation, and evidence of a feedback-driven symmetry-breaking phase transition.
Statistical mechanical equilibrium is based on the straightforward observation that a system's eigenstates are probabilistically occupied according to the thermal Boltzmann distribution, subject to any relevant constraints. The resulting density operator affords a timeindependent description of the system as a thermalequilibrium steady state. Conversely, any system described by a time-independent density operator may be associated with thermal equilibrium, for an effective Hamiltonian proportional to the logarithm of that density operator, although detailed balance may not necessarily be satisfied. In this sense then, is it possible to identify and explore systems maintained in low-entropy steady states associated with exotic effective Hamiltonians with highly non-local or N -body interactions? Optical pumping [50] and laser cooling [51]-both described by the physics of open quantum systems-are iconic examples where large ensembles of atoms enter low-entropy steady states well described by single-atom physics with little correlation between atoms. * jeremy.young@colorado.edu † ian.spielman@nist.gov In contrast, quantum error correction (QEC) codes are highly specialized and sophisticated examples of measurement-feedback systems that can dynamically stabilize strongly entangled arrays of qubits [52,53]. A digital quantum computation device consists of an interconnected collection of qubits, and QEC utilizes multiqubit syndrome measurements, often realized via "ancilla" qubits that are first coupled to the error-prone qubits and then measured. The coupling and measurement are designed to be sensitive to select errors, but not to the quantum state involved in computation. In some forms of QEC [54], the classical information resulting from measuring the ancilla can inform subsequent error-correcting feedback stages, thereby maintaining the quantum computation device in a type of dynamical steady state.
Quantum state preparation is closely related to QEC: both attempt to drive a system towards a particular state. In the case of QEC, this is the state prior to an error; in the case of state preparation, it is a particular state of interest. State preparation is particularly important in the field of metrology [55,56]. By preparing highly entangled states, such as squeezed states [35,36,[57][58][59][60][61], it becomes possible to make measurements that are far more accurate than in an unentangled classical system. While a variety of sophisticated techniques have resulted in experimental measurements with incredible accuracy, incoherent effects often provide a fundamental limitation. Like with typical QEC for quantum computation, measurement and feedback may provide a means of reducing the problems introduced by incoherent processes [62][63][64][65][66][67][68][69][70][71][72].
Here, we study this problem from three perspectives progressively increasing in complexity. (I) we begin by considering a minimal two-site system in the large-atomnumber classical limit-a type of non-linear top-and derive the exact equations of motion including feedback. This model predicts what types of feedback can drive the system into different fixed points, providing guidance for what type of quantum states will appear. (II) We then consider a stochastic Schrödinger equation description of the same double well. By investigating both the classical large-atom limit and the quantum smallatom limit, we connect the classical fixed points to their associated quantum dynamical steady states. (III) We conclude with the numerically exact time dynamics for modest sized Hubbard chains in the quantum low-density limit (see Refs. [40] and [41] for investigations in the more classical high-density limit) and show that nearly arbitrary distributions of number states can be neardeterministically prepared, and that non-local feedback algorithms can both improve fidelity and accelerate state preparation. Moreover, we show that the corresponding feedback can be modified in a simple way to realize a symmetry breaking transition.
From a broader perspective, this manuscripts shows that straightforward feedback derived from highly idealized continuous quantum measurements can stabilize low-entropy dynamical steady states. It is the task of future work to expand upon the impact of increasingly realistic experimental parameters.

I. MODEL
Here we focus on the idealized 1D Bose-Hubbard model subjected to continuous weak measurement [73,74] to illustrate the possible dynamical steady states in manybody quantum systems. As depicted in Fig. 1 this model consists of: a 1D Bose-Hubbard chain (the system) dispersedly coupled to a transverse laser field (reservoir); measured via phase contrast imaging (optical homodyne detection); and with local tunneling strengths and on-site potentials that are dynamically adjusted based upon the measurement outcome (feedback). At time t, the Bose-Hubbard chain is governed by the system Hamiltonian expressed in terms of the bosonic field operatorsâ † j describing the creation of an atom at site j andn j = a † j a j the number operator at site j; with positive real-valued tunneling strength J j (t) between sites j and j +1; on-site energy V j (t); and interaction strength U . The chemical potential is absent in this microcanonical ensemble study.
Measurement model. The system is coupled to the reservoir by the system-reservoir Hamiltonian that describes a dispersive off-resonance interaction of light with two-level atoms, detected by homodyne measurement [75] or equivalently phase contrast imaging, shown in Fig. 1. Here, the bosonic field operatorb † j describes the addition of a photon into mode j associated with lattice site j (here we focus on an idealized case with mode functions associated one-to-one with lattice sites); g captures the strength of the system-reservoir coupling; and the reservoir modesb j are each in the coherent state |α prior to interacting with the system. Physically this operator describes rotations in theX j -P j quadrature plane for each reservoir mode, in proportion to the local number of atoms. After interacting with the system for a time t m , the reservoir modes are projectively measured by optical homodyne measurement (in practice implemented by phase contrast imaging) giving access to theX j = (b † j +b j )/2 quadrature observables. A strong measurement on the reservoir in effect affords a single weak measurement of the system observablesn j with strength governed by g, the reservoir field amplitude α, and the measurement time t m . A measurement outcome at time t contains a contribution from the system observable's instantaneous expectation value n j (t) , along with a noise contribution δn j = m j /ϕ from projectively measuring the reservoir coherent state |α j . Here m j is a classical random variable with zero mean m j = 0 and variance m j m j = δ j,j /2. We consider an optical mode function consisting of a traveling pulse of light with extent ct m , where c is the speed of light. Converting from the continuum mode functions to this compact function introduces a factor of 1/(ct m ) into the definition of the generalized measurement strength parameter. This leads to the generalized measurement strength parameter ϕ 2 = g 2 |α| 2 t m /c, defining the measurement noise δn 2 j = ∆n 2 = 1/2ϕ 2 and the back-action of this measurement on the system as described by the Kraus operator [76] K(n j ) = exp conditioned on the measurement outcomes n j . Given an initial state |ψ and a measurement outcome n j , the postmeasurement state of the system is given byK(n j )|ψ . As one might intuitively expect from a simple classical measurement, this gives a distribution of measurement outcomes with a width ∝ 1/ϕ that decreases linearly with the system-reservoir coupling g but as the square root of the measurement time t m and intensity of the probe field |α| 2 . In the following analysis, we consider the continuous limit of many such weak measurements, taking the time between subsequent measurements to be t m and applying measurements with a strength associated with a measurement time t m . Feedback model. For sufficiently weak measurements, the quantum projection noise δn j present in any individual measurement can greatly exceed the contribution of the instantaneous expectation value n j (t) . As a result, we borrow from classical control theory [77], and derive FIG. 1. System and measurement schematic. An ensemble of locally interacting (with strength U ) neutral atoms described by the bosonic field operatorsâj are confined in a lattice potential with tunneling strength Jj between sites j and j + 1. These atoms interact with independent modes of an optical field described by bosonic field operatorsbj. After interacting, the optical modes are imaged via phase contrast imaging, implementing a spatially resolved optical homodyne detection providing a measurement sensitive to the quadrature field op-eratorXj = (bj +b † j )/2. This idealized measurement assumes that the incident laser field is only forward scattered, essentially making the paraxial approximation.
an error signal (t) from a temporal low-pass filter applied to the direct measurement outcomes. This consists of an integrator with a low-frequency gain limit, with time constant τ . This filter retains the low frequency 1/τ components of n j (t) (containing contributions from both system dynamics and projection noise), but rejects the high frequency components (dominated by projection noise). Here we assume that the measurement time t m τ , implicitly making Eq. (3) a stochastic differential equation, see Ref. [41] for a more detailed discussion. The resulting error signal j (t) thus approximates the continuous measurement limit of a sequence of weak measurements of n j (t) with measurement times τ regardless of the physical measurement time t m τ . Note that because this filter is applied in real time, j (t) will lag behind n j (t) by τ , and we take j (0) = 0 in our simulations. Because the goal is neither to perform quantum state estimation nor to drive the system into a predefined state, we do not employ Kalman filter tech-niques [78].
Here we focus on possible dynamical steady states when this classical information is then fed back into the Hamiltonian in one of two ways. We shift either the onsite energy (4) in proportion to the local feedback signal (with gain G V ), or the tunneling based on the difference along that link (with gain G J ). These two forms of feedback are the most simple physically realistic options: the potential V j (t) simply shifts in proportion to the error signal; however, because negative tunneling is difficult to achieve, quadratic feedback is the most simple realistic feedback to the tunneling term. Both of these can be realized using quantum gas microscopes [79,80] where a digital mirror device can locally control both the intensity of the lattice light (increasing or decreasing tunneling [81]) as well as changing the intensity on the lattice sites (altering the on-site potential [30]).

II. TWO-SITE LATTICE
The two-site Bose-Hubbard model with N atoms can be recast as an f = N/2 angular momentum system with mappingsn 1 −n 0 → 2F z ,â † 1â 0 +â † 0â 1 → 2F x , and a † 1â 0 −â † 0â 1 → 2iF y , that obey a dimensionless angular momentum commutation relation [F i ,F j ] = i ijkFk . This gives the Hamiltonian where we have absorbed constant terms into an overall shift in the zero of energy and defined ∆ = V 1 − V 0 . This also results in the single filter equation for (t) ≡ 1 (t) − 0 (t) and feedback equations which are slightly adjusted with an additional factor of 2 on the right-hand side of the filter equation, resulting from the details of the angular momentum mapping. The on-site measurement outcomes at time t can be re-expressed as with contribution [82] from the instantaneous expectation value F z (t) . As before, the noise is defined by a classical random variable m with zero average m = 0, and variance m 2 = 1/2. Back-action is described by the Kraus operator Taken together, these expressions provide an equivalent formulation of the two-site system in the language of angular momentum. Classical limit. The Heisenberg equations of motion in the large-N limit (i.e., ignoring quantum fluctuations and thus measurement uncertainty) reduce to the classical equations of motioṅ where F ≈ ( F x , F y , F z ) denotes the classical angular momentum vector, ignoring quantum fluctuations. Note that because the interaction U enters via a quartic term in the Hamiltonian, it enters the classical equations of motion nonlinearly. Thus in order to properly compare systems with different f , we keep U f fixed to a constant value. Similarly, since the potential and tunneling feedback terms behave like effective fourth-and sixth-order terms in the Hamiltonian, respectively, we must also keep f G V and f 2 G J fixed to constant values. Tunneling causes the macroscopic spin vector to precess around e x ; a potential imbalance leads to precession around e z ; and the non-linear interaction term effectively drives precession around e z with angular frequency in proportion to F z . Figure 2(a) plots the energy [from Eq. (6)] and trajectories [from Eq. (11)] for ∆ = 0, as a function of the polar angles φ and θ with respect to an e x aligned spherical coordinate system (φ = 0 in the e ydirection). The displayed trajectories illustrate the two classes of orbits. (1) Near F x = f (i.e., θ = 0), the phase φ exhibits running solutions [83,84] orbiting around e x : Josephson oscillations (JO). (2) For larger θ, closed orbits around F x = −J/U , F y = 0 and F 2 z = f 2 − F 2 x correspond to oscillations about a density imbalanced excited state [85]: macroscopic self trapping (MST), present only for |J/U | < f . For |J/U | > f the two MST solutions merge to form a second stable JO fixed point at Linearized dynamics. The impact of feedback is most clearly understood by first linearizing the dynamics around each of the stable fixed points of the JO and MST orbits, in general giving elliptical orbits. The displacements δF x,y,z each obey second-order differential equations such asδ F z + ω 2 δF z = 0 (12) describing motion in an effective harmonic potential, with angular frequencies F z and shifts the location of the fixed points (potentially even eliminating them entirely), as well as altering the frequency and ellipticity of the orbits. Still, the form of the differential equations is unchanged: no damping. In the above example, conventional damping would arise from an additional friction termδF z /τ , and more generally damping (or anti-damping) will occur only when odd and even derivatives are mixed. In our case, (t) mixes these derivatives when τ = 0, changing the relationship between (t) and f z (t) according to Eq. (7), thus changing the effect of δF z on Eqs. (8,11). The effect of feedback can be quantified in linear response about both the JO and the MST fixed points, allowing us to derive a damping rate Γ in each case, in the limit of weak damping ωΓ 1, where ω is the angular frequency of the relevant fixed point. In this limit, the damping rates are ωΓ = γωτ /(1 + τ 2 ω 2 ), with the most effective damping when ωτ = 1, and where the strength γ governs the system dependence of the feedback. For feedback in the potential, the resulting strengths are showing that, depending on the sign of G V , the system can be damped effectively. Focusing on the case of positive interactions U > 0 and positive base tunneling J 0 > 0, either only the JO fixed point is stable (0 < −G V < U ), only the MST fixed points are stable (0 < G V < U ), or both types of fixed points are stable (−G V > U ). In contrast, for feedback in the tunneling channel, the quadratic response δJ = G J (t) 2 leads to very slow damping as the density-balanced JO fixed point is approached, i.e., at the JO fixed point, (t) = 0, and the linear contribution to δJ is zero. Thus, only the MST state has damping or anti-damping in linear response, and for the physical sign of the gain coefficient G J > 0, this results in anti-damping.
Fluctuations. Although the classical limit omits quantum fluctuations (by assuming a spin-coherent state), it need not omit measurement back action. A prescription similar to that of Ref. [40] for coherent states gives a final state conditioned on the random variable m for measurements ofF z described by Eq. (9), where θ z is defined with respect to e z . This prescription takes an initial spin coherent state, then applies the Kraus operator; because the resulting state is not necessarily a spin-coherent state, we find the coherent state with the largest overlap as the updated state. The core message of this expression is that F z is updated as informed by the measurement outcome, but for initial states nearly polarized along e z , the state is nearly unaltered. Physically, this is expected because we are working at fixed total angular momentum f , and a state already polarized along e z cannot become more polarized. Lastly, equating Eqs. (9) and (13) suggests an optimal measurement strength for which the true value of F z following the measurement [Eq. (9)] is equal to the measurement outcome [Eq. (13)], i.e., the measurement outcome accurately reflects the new state of the system. This occurs only for an optimal measurement strength ϕ 2 opt = (2f sin 2 θ z ) −1 , that would be selected to yield optimal performance near the desired fixed point. We also note that with Eqs. (9) and (13), this implies that for the optimal measurement strength, the measurement noise and back-action both scale like f −1/2 , fractionally going to zero in the f → ∞ limit as one would expect.
The low-pass filter from which we derive the error signal has time-constant τ , which plays the role of an effective measurement time. Assuming no delay between measurements, this implies an individual-measurement coupling ϕ 2 opt,0 (2f sin 2 θ z ) −1 × (t m /τ ), so that the τ /t m measurements which take place in the time interval τ give a single average outcome with the optimal strength ϕ opt .
Quantum double well. With these basic insights from the classical model in mind, we compare to a quantumtrajectories stochastic Schrödinger equation simulation of this double-well problem, including all effects of backaction resulting from measurement to the classical spin model, including measurement noise and back-action. Figure 3(a) shows the coherent evolution of a spin coherent state with N = 64 or f = 32, according to the classical spin model (red) and the quantum simulation (blue), with nearly perfect overlap. The agreement continues to improve with increasing N , as neglecting quantum fluctuations becomes an increasingly good approximation. This overlap persists in Fig. 3(b) which includes the effects of measurement, back-action, and feedback, again showing good correspondence between the classical and quantum models. Lastly, as the number N is reduced from 64 to 2 (and U is correspondingly increased by a factor of 64/2 to keep U f = 2), neglecting fluctuations becomes a very poor approximation, and the classical and quantum trajectories deviate.
The classical-limit modeling provided valuable insight in guiding parameter selection even for quantum systems: (1) feedback can change the dynamical steady state by introducing effective potentials; and (2) the optimal filter time τ is inversely proportional to the dynamical time scale of excitations in the measurement channel. A basic intuition for the latter point comes from a continuously monitored classical harmonic oscillator. The measured position converts into knowledge of momentum one-quarter period later, allowing a time-delayed feedback force to reduce that known velocity. In the present case, τ in the lowpass filter gives the same outcome by introducing a phase shift for signals with angular frequency 1/τ . We certainly expect that more complex filters, proportional-integral-differential (PID) or Kalman [78] filters for example, would give improved performance.

III. EXTENDED LATTICE
In this section, we extend our analysis to the case of small one-dimensional Bose-Hubbard chains. Unlike the case of two sites, there is no angular momentum map, and the appropriate mean-field description is a discrete Gross-Pitaevskii equation, similar to that discussed in Refs. [40,41]. Here, we focus on the quantum region by considering states with mean densities of just a few particles per-site and studying the system's behavior numerically using a quantum trajectories approach [86][87][88]. Correspondingly, the fluctuations in all figures are a result of sampling error. We will also consider a range of target density distributions, in which the error signal derives from the difference between the observation and the target. Throughout, we will express everything in units where the interaction strength U = 1. Recent closely-related work investigates preparing similar Mottlike states via global measurements rather than the local measurements considered here [39].
Feedback models. We consider two possible tunneling feedback schemes. The first is the nearest-neighbor approach introduced in Eq. (5). The second, the "imbalance approach", uses highly non-local information, and relies instead on imbalances between the left and right sides of the system. For the double-well system considered earlier, these two schemes are identical. The tunneling strengths are, respectively, where and N j indicates the target occupation number of site j.
Since the tunneling strength without feedback is J 0 = 0 and measurement occurs in the number basis, this means that the target state |N 1 N 2 · · · will be an evolutionfree fixed point in the absence of feedback induced by noise from the measurements. Given the structure of the imbalance approach, which entails separating the system into left and right sides, we will consider open boundary conditions when discussing state preparation. The purpose of the imbalance approach is to take advantage of the global knowledge of the target state, knowledge that is not used in the nearest-neighbor approach. For example, consider the scenario in which the target state is |43210 and the system is in the state |04321 . In this case, the nearest-neighbor approach will only initially turn on the tunneling between the first two sites and gradually turn on the other tunneling terms from left to right as the bosons move to the left. In contrast, the imbalance approach will turn on all tunneling terms and begin to transfer atoms from the right side of the system to the left. However, the use of this more global knowledge comes at a cost to the feedback uncertainty. Since j,L/R involves a sum of N different measurement records, the resulting uncertainty will be enhanced by a factor of √ N . As a result, there is a trade-off between using global knowledge and the uncertainty in the feedback. Finally, we will discuss these two feedback approaches in terms of a continuous measurement rate κ, where ϕ ≡ √ 2κt m and we use t m = 0.01 in our numerical simulations.
State preparation. We compare the performance of these two approaches for three different target states: ψ 1 = |11111 , ψ 2 = |20202 , and ψ 3 = |30003 . We fix the filter time τ = 0.3 and the evolution time T = 10 + τ [89] in units of U −1 and examine the behavior of the final target state overlap | ψ i |ψ(T ) | 2 as a function of κ and G. For each trajectory, we initialize the system in a Haar random state. While we could consider the steady-state behavior of this system, from the perspective of creating a target state in a physical system, it is more useful to consider the finite-time behavior. Fig. 4 illustrates the target state overlap for ψ 2 at different values of κ and G using the nearest-neighbor approach. The behavior for other target states and the imbalance approach are qualitatively similar. For all values of the measurement strength κ, we see that there is an optimal choice of the gain G for which the overlap is maximized. This optimal value of the gain represents a balance between the need to direct the system to the desired target state quickly and the fact that measurement uncertainty can lead to "errors" in the application of the feedback. We identify the optimal value of the gain and the corresponding target state overlap for the different target states and measurement strengths in Fig. 5. Several trends can be identified in the optimal feedback behavior.
The first trend is that the optimal gain grows linearly with the measurement strength. This reflects the fact that stronger measurements correspond to more accurate knowledge of the system's state, so the system can be more strongly driven to the target state without making errors due to the measurement uncertainty. Additionally, the gain must increase with the measurement strength to avoid quantum Zeno effects, i.e., the dynamics due to the feedback must become sufficiently fast so that the repeated weak measurements do not effectively become a strong measurement before the system can evolve coherently.
The second trend is that optimal gain is larger for the nearest-neighbor approach than the imbalance approach. This is most likely a result of the fact that the uncertainty is larger for the imbalance approach. Since the imbalance approach uses 5 measurements for each tunneling link while the nearest-neighbor approach uses 2, the relative uncertainties in the tunneling strength differ by roughly a factor of 2.5 [90], which is consistent with the observed behavior.
The third trend is that for the parameters and target states considered, the imbalance approach performs much better than the nearest-neighbor approach. The nearest-neighbor approach performs better the more homogeneous the target state is. However, this begins to change at lower values of the measurement strength, where the target state overlap sharply begins to fall. This is a consequence of the fact that, at sufficiently small measurement strengths, the uncertainty in the measurement-and therefore the feedback signaloverwhelms the useful information. Since the imbalance approach involves more measurements, this behavior occurs earlier than for the nearest-neighbor approach.
Nearby state preparation. For larger systems, the effect of the increased uncertainty that comes with the imbalance approach becomes far more important. A natural approach to avoid this issue is to use a more hierarchical approach to the application of feedback. When the system starts initially far from the desired target state, a more global approach to feedback should be used. As the system gets closer to the target state according to the measurement record, then increasingly local feedback approaches can be used.
While the systems we can consider numerically are too small to apply such a hierarchical feedback approach, we can explore how the performance of the two feedback approaches can change if the system is already closer to the target state. To do so, we consider initial states whose only difference from the target state is that a single bo- son has been moved one site away. Additionally, we allow the feedback to be applied for only one unit of time after the initialization of the measurement record (which requires time τ ), reflecting the fact that the local feedback would be applied for a shorter time in this hierarchical approach. The results of our numerics are shown in Fig. 6 for κ = U , demonstrating that the nearest-neighbor approach performs better when the system is already close to the desired state due to the reduced uncertainty in the feedback. Note that this is in spite of the fact that the nearest-neighbor approach will initially turn on three coupling terms while the imbalance approach will turn on only one in the absence of measurement uncertainty.
The emergence of phase transitions via feedback has been investigated in a variety of diverse contexts [41][42][43][44][45][46][47][48][49], and here we explore a symmetry-breaking transition in an extended Hubbard chain. To observe this, we will consider a six-site lattice with periodic boundaries whose tunneling feedback drives the system towards the states |202020 or |020202 . Due to the periodic boundary conditions and type of symmetry breaking, the left-right imbalance approach is not applicable, so we consider a generalization of the nearestneighbor feedback. In particular, the tunneling feedback is with 7 ≡ 1 . This is essentially a product of the tunneling feedback in Eq. (15) for N j − N j+1 = ±2. The Gaussian term is included in order to increase the tunneling for j − j+1 ≈ 0 compared to just the quartic potential so that it is comparable to the corresponding tunneling when only one of the two aforementioned states is a target state. This will serve as a barrier which hinders the system from moving from |202020 to |020202 and vice-versa (see Fig. 7). An additional important feature to note for this choice of feedback is that these are not the only two states which will lead to small tunneling. States of the form |024200 can also be relatively stable. While there is one pair of sites with large tunneling, this pair of sites has no bosons, so the large hopping is not sufficient on its own to change the state. However, there are two factors that can make such states unfavorable. The first is that the creation of these states requires four bosons on a single site, which is energetically unfavorable. The second is that once measurement-induced fluctuations cause one of the bosons at sites with two bosons to move into an empty site, the boson will quickly move to the next site.
In order to investigate the symmetry-breaking behavior of this feedback, we consider the behavior of the steady-state density matrix ρ ss for fixed κ as a function of G, which we obtain via a combination of ensemble and ergodic averaging. Additionally, we define an effective Hamiltonian H ss according to e −Hss ≡ ρ ss . In analogy to equilibrium systems, a symmetry breaking phase will occur when the two lowest energy levels of H ss become gapped from the rest of the spectrum. In terms of ρ ss , the two corresponding states will have a much larger probability than the rest of the eigenstates, and the gap in H ss describes the exponential suppression in probability of these other eigenstates. Hence in the absence of a gap, no particular configuration is favored, and the system appears disordered. Fig. 8 plots the eigenvalues of H ss as a function of G for κ = 5U using the initial state ψ(0) = |111111 , which will prefer both symmetry-breaking states equally, although we expect the same steady state for Haar random initial states. We see that an effective gap opens for a range of G, with two effective energies much smaller than the rest. As expected, these two eigenvectors correspond approximately to |202020 and |020202 , while the other possible stable states (e.g., |024200 ) are not strongly occupied. Additionally, as this gap increases, the auto-correlation times increase, necessitating more sampling for larger gaps. This corresponds to the fact Hss as a function of G. Eigenvalues are calculated after symmetrizing ρss with respect to translations and reflections, leading to degeneracy for states related through these symmetries; quantitatively similar behavior emerges without symmetrization, although the eigenvalues are no longer as degenerate. The orange dots are doubly degenerate and correspond approximately to |202020 , |020202 , with the strongest overlap occurring for the larger gaps. For G < .5U , ρss is averaged over 128 trajectories from t = 100U −1 to t = 400U −1 . For G ≥ .5U , ρss is averaged over 32 trajectories from t = 20U −1 to t = 320U −1 . The initial time delay for averaging ensures that the system has relaxed to the steady state. that when the system is in either of the two symmetrybroken states, the feedback makes it difficult to leave this state, so the system spends a long time in the same state, thus slowing down the ergodic sampling of ρ ss . In the limit of large G, even small fluctuations will drastically modify J j (t), preventing the system from preferring any particular state and causing the gap to close.
As G is increased, this gap eventually closes and does not reopen. Physically, this is because in the limit of large G, the hopping becomes very large, so small fluctuations in the measurements lead to large fluctuations in the hopping, preventing the system from preferring any given state. This is qualitatively similar to an infinite temperature limit in which no particular state is energetically preferred. In the limit of small G, the gap similarly closes because the measurement strength overwhelms the feedback, causing the system to relax to a random number state. The start of this closure can be observed in Fig. 8, although due to the quantum Zeno effect, identifying the steady state for smaller values of G becomes numerically prohibitive as it requires increasingly long relaxation times.

IV. OUTLOOK
This study combines concepts from many-body physics with the techniques of quantum control and feedback to create and characterize low entropy many-body dynamical steady states. Creating such dynamical steady state in the laboratory environment hinges on exploring the implication of experimental realities, for example: What are the consequences of imperfect detection? Of limited spatial resolution? What is the impact of limited feedback bandwidth? Real implementation must face these questions head-on.
From a foundational perspective, it remains to be seen if there exist dynamical steady states which are forbidden in associated thermal equilibrium systems. For example, nonreciprocal couplings can be realized using feedback, which can lead to rich non-equilibrium phenomena [91][92][93][94]. Alternatively, feedback that is non-local in time can be used to engineer non-Markovian baths even though the measurements themselves are Markovian. Similarly, even with local measurements of density, non-local feedback can emulate aspects of long-range potentials [37]. Stochastic descriptions such as ours can be reframed in terms of the so-called feedback master equation [95], in which the time delay introduced by the low pass filter introduces non-Markovian terms that cannot be described in a Lindblad form. This at least provides the opportunity for creating dynamics and steady states that are not achievable using realistic Lindblad terms. Moreover, even if the steady state may look thermal, its linear response may still violate the fluctuation-dissipation theorem [96,97], indicating the persistence of non-equilibrium features. Additionally, new types of dynamics are possible when the bosons include a spin degree of freedom, such as quantum Zeno-like effects which can emerge if one spin state is subjected to stronger measurements than the other.
While we have demonstrated how the use of non-local feedback can be employed to prepare desired target states more efficiently than local feedback, this was done for the relatively simple case of Mott-like number states. An interesting next direction, then, is to consider preparing more complex states, such as superfluid-like states with long-range coherences. By adding density modulations through the use of non-local feedback as we did for the number states, this opens the possibility of preparing supersolid and supersolid-like states [98][99][100][101]. More broadly, this would provide further insight into how the use of measurement and feedback affects the ability to realize intrinsically quantum dynamics.
Although we have illustrated how a symmetrybreaking phase transition may in principal emerge due to measurement and feedback in a many-body quantum system, several open questions remain. Here, our analysis has been restricted to small 1D chains, so an important question is whether this corresponds to a phase transition in the thermodynamic limit. Even if it does not exist in 1D, an analogous transition may exist in higher-dimensional systems. Additionally, understanding the behavior of critical phenomena in these systems is another important direction, such as whether they are equivalent to a quantum phase transition, a classical phase transition, or something entirely different. In Ref. [46], it was shown that in zero-dimensional systems, modifying the form of non-Markovian feedback can lead to changes in the critical exponents of the phase transition, so similar phenomena may arise in a many-body context as well.
A natural direction to explore in order to investigate these questions are measurement-induced phase transitions, which have been the subject of intense theoretical research in recent years and are defined by the scaling behavior of the entanglement entropy [102][103][104][105][106][107][108][109][110][111][112][113][114][115][116][117][118], with recent initial experimental evidence [119]. Aside from the key difference that these systems are not subject to feedback, these often involve coherent evolution defined by random unitary circuits and strong measurements, although similar behavior has been shown to emerge for weak measurements as well [114][115][116][117][118]. A natural question is how these phase transitions can be modified through different applications of feedback; will the phase transition only shift, or can qualitatively new physics emerge?
Another promising approach to exploring the above questions is to investigate the similarities that systems subject to measurement and feedback have with drivendissipative systems, which are systems that are subjected to coherent drive and incoherent dissipation that have been studied extensively in a variety of contexts [10-21, 93, 120-136]. From an abstract standpoint, these two types of systems are very similar, with dissipation playing a role analogous to continuous measurements and drive playing a role analogous to the feedback. As a result, insights or techniques from studying one type of system can lead to insights in the other. For example, phase transitions in driven-dissipative systems have been studied extensively using a Keldysh-Schwinger and functional integral formalism [93,[122][123][124][125][126][127][128][129][130][131], so it is important to explore how these same techniques can be applied to systems subject to measurement and feed-back. Similarly, effective classical equilibrium criticality is observed generically in many driven-dissipative phase transitions [127][128][129], with some important exceptions [93,[130][131][132], so measurement-feedback systems may provide new avenues to realize novel forms of nonequilibrium criticality and quantum criticality. Moreover, recent research indicates that dissipative phase transitions and measurement-induced phase transitions need not coincide [136], which means that measurementfeedback phase transitions may similarly lead to distinctive forms of criticality.