Feedback Induced Magnetic Phases in Binary Bose-Einstein Condensates

Weak measurement in tandem with real-time feedback control is a new route toward engineering novel non-equilibrium quantum matter. Here we develop a theoretical toolbox for quantum feedback control of multicomponent Bose-Einstein condensates (BECs) using backaction-limited weak measurements in conjunction with spatially resolved feedback. Feedback in the form of a single-particle potential can introduce effective interactions that enter into the stochastic equation governing system dynamics. The effective interactions are tunable and can be made analogous to Feshbach resonances -- spin-independent and spin-dependent -- but without changing atomic scattering parameters. Feedback cooling prevents runaway heating due to measurement backaction and we present an analytical model to explain its effectiveness. We showcase our toolbox by studying a two-component BEC using a stochastic mean-field theory, where feedback induces a phase transition between easy-axis ferromagnet and spin-disordered paramagnet phases. We present the steady-state phase diagram as a function of intrinsic and effective spin-dependent interaction strengths. Our result demonstrates that closed-loop quantum control of Bose-Einstein condensates is a powerful new tool for quantum engineering in cold-atom systems.


I. INTRODUCTION
Quantum gas experiments have exquisite control over the low-energy Hamiltonian governing system dynamics, providing demonstrated opportunities to study interacting many-body quantum systems with great precision. As a result, ultracold atoms have emerged as a leading platform in 'analog quantum simulation' [1][2][3][4][5][6], where experiments have successfully explored condensedmatter phenomena such as the superfluid-Mott insulator transition [7], the BEC-BCS crossover [8,9], and spinorbit coupling [10]. Cutting-edge experiments now realize systems with long-range interactions [11] or novel non-equilibrium dynamics [12,13]. In contrast, quantum simulation of open systems remains relatively unexplored [14], and careful application of feedback control to many-body quantum systems is a new approach toward this goal.
Feedback control of many-body systems could enable observation of a wide range of new phenomena in dynamical steady state, where a potentially larger class of states are possible than in thermal equilibrium [15,16]. Existing proposals include preparation of many-body pure states via reservoir engineering [17][18][19][20], nonthermal steady states [21,22], stable non-Abelian vortices [23], or time crystals [24]. Here, we showcase the flexibility of weak measurements coupled with spatially resolved feedback for quantum simulation of time-dependent effective Hamiltonians using a two-component Bose-Einstein condensate (BEC) as a model spinor system [25][26][27].
We develop a theory of weak measurement and classical feedback in weakly interacting quantum systems framed in the context of quantum control theory [28]. Using our general formalism we investigate the steadystate phases of a two-component BEC subject to weak measurement and classical feedback via a spin dependent applied potential, enabling both density and spin dependent feedback protocols.
Spatially local feedback can result in spin-dependent effective interaction terms in the stochastic equation governing condensate dynamics. Depending on the interplay of intrinsic and effective (i.e. feedbackinduced) spin-dependent interactions, the condensate steady-state phase is either an easy-axis ferromagnet or spin-disordered paramagnet. The effective interaction is tunable via the gain of the feedback signal, enabling a reversible, feedback-induced phase transition. The transition is reminiscent of what is achieved by tuning intrinsic interactions via a spin-dependent Feshbach resonance [29], however here the atomic scattering lengths remain unchanged. We develop a signal filtering and cooling scheme to minimize heating and show that the condensate remains intact under feedback and measurement backaction. Our result opens the door to engineering new dynamical and/or spatially dependent effective interactions in quantum gases via closed-loop feedback control.
The paper is structured as follows: In Sec. II we present our main formal results, including the stochastic equation describing condensate dynamics, and introduce a toy model illustrating the salient features of the control protocol. We show that locally applied feedback induces a phase transition between easy-axis ferromagnetic and disordered paramagnetic phases in a two-component condensate.
In Sec. III we elaborate on our feedback cooling protocol and characterize the resulting steady state via condensate fraction, Von Neumann entropy, and energy. We show that heating due to measurement backaction can be effectively mitigated by feedback cooling. In Sec. IV we discuss the feedback-induced steady-state phases in more detail and elucidate the nature of the phase transition in our system. We conclude in Sec. V.

A. General Formalism
We model dispersive imaging of a quasi-onedimensional (1D) multicomponent Bose-Einstein condensate of length L via spin resolved phase-contrast imaging [39] and we label individual components by an index s. We consider time and space resolved measurements of atomic densityn s (x, t) in each component using the Gaussian measurement model developed in detail in Ref. [40]. Stroboscopic weak measurements with strength ϕ result in the measurement signal where m s (x) describes spatiotemporal quantum projection noise associated with the measurement. The measurement is characterized by Fourier domain Gaussian statisticsm s,k = 0 andm s,kms,k = LΘ(|k| − k c )dW s,k dW s ,k /2dt 2 , where dW s,k is a Wiener increment with dW s,k = 0 and dW s,k dW s ,k = dtδ ss δ kk for a time increment dt [41]. The Heaviside function Θ enforces a momentum cutoff at k c = 2π/λ, accounting for the fact that the physical measurement process can only resolve information with length scales larger than λ/2π. The observer does not directly obtain information about the condensate phase using this protocol. We use the aggregate measurement result M, a function of x and s, to generate feedback signals in the form of a single-particle potentialV [M], where· indicates an operator in component space. In this work we consider a potential which is local in space.
We describe the condensate in the mean-field approximation using a complex spinor order parameter Ψ(x) = (ψ 1 (x), ψ 2 (x), . . .) T , where ψ s (x) is a classical field describing the dynamics of component s. The total density is n(x) = Ψ † (x)1Ψ(x) and the order parameter is normalized to the number of particles, N = dx n(x). From Eq. (1) the measurement results at the meanfield level therefore depend on the field amplitude via n s (x) → |ψ s (x)| 2 . Measurement backaction leads to stochastic evolution of the order parameter, which results in condensate heating [40,42] in the absence of a cooling protocol, which we describe in Sec. III.
The combined measurement and quantum control process is described by a stochastic equation of motion for the condensate order parameter Ψ(x). Here denote contributions from unitary (i.e. closed system) evolution, measurement backaction, and feedback, respectively and µ is the chemical potential. We adopt the implied summation convention over repeated indices and set = 1.
Using this general formalism we study a condensate of 87 Rb atoms from which we select two hyperfine states, yeilding a two-component condensate [27,43] with components denoted by s = ↑, ↓. The Hamiltonian in Eq. (3) is the usual Gross-Pitaevskii equation (GPE) describing closed system dynamics, which takes the explicit form for two component condensates, with (x, t) indices suppressed for clarity. Here, S z (x) = Ψ † (x)σ z Ψ(x) indicates the spin density andσ = (σ x ,σ y ,σ z ) is a vector of the Pauli operators. The single particle Hamiltonian isĤ 0 =p 2 /2m a for atoms of mass m a . The intrinsic spin-independent u 0 and spin-dependent u 2 interaction strengths serve to define ξ = 1/ √ 2m a µ and ξ s = ξ u 0 /2|u 2 |, the healing length and spin-healing length respectively. Equation (4) describes measurement backaction. Separate measurements of each condensate component result in independent backaction noise m s (x). Equation (5) describes feedback, applied via the potential termV [M]. The feedback potential combines a deterministic part containing information about the condensate dynamics with a stochastic part due quantum projection noise. Therefore, both dΨ| F and dΨ| M contribute to stochastic condensate dynamics. When each individual measurement is very weak, the density of noncondensed particles remains low. Therefore we assume Ψ(x) to be well described by a lowest order Hartree-Fock theory throughout it's evolution. This assumption is validated in Sec. III B and III C.

B. Key Feedback Concepts
Our aim is to develop feedback schemes which add new effective interaction terms to the Hamiltonian while min-imizing quantum projection noise. We illustrate the core concept of feedback using a toy model. The toy model is a simplified version of the feedback protocols developed in later sections, that nonetheless illustrates a key result: weak measurements combined with feedback can be used to engineer new effective Hamiltonians.

Toy Model
Here we construct a minimal model of measurement and feedback for single component systems, and therefore suppress the component index s. We weakly measure the density, then apply a proportional feedback potential where the gain parameter g 0 denotes the feedback strength. Inserting Eq. (1) into Eq. (7) gives a feedback potential with two contributions. The first is an effective mean-field interaction and the second is a stochastic contribution The effective HamiltonianĤ eff (x) has the same form as the spin-independent term in Eq. (6), but with u 0 replaced by an effective interaction constant u eff 0 = u 0 + g 0 . Likewise, the noise in the stochastic evolution is modified due to the contribution of V fluct (x, t). This simple model illustrates how feedback can be used to create new effective Hamiltonians with modified interaction terms.
Returning to the two-component case, we consider the spin-dependent feedback potentiaľ describing separate contributions to the density and spin sectors controlled by independent gain parameters g 0 and g 2 , respectively. Measurement signals M s are used to calculate total density and spin density, given by M n = M ↑ + M ↓ and M z = M ↑ − M ↓ , respectively. Following the same algebraic arguments, the feedback potential (12) leads to effective interaction strengths u eff 0 = u 0 +g 0 , u eff 2 = u 2 +g 2 , along with modified stochastic noise on each component ψ s .
In the following, we use this guiding prinicple to develop a measurement and feedback scheme which controls the magnetic properties of a two-component condensate without changing the internal interaction parameters. The simplified protocol presented in this section is impractical due to runaway heating [40], from the repeated and uncompensated application of the stochastic potential in Eq. (11). In Sec. III we introduce a feedback cooling protocol that prevents runaway heating and thus completes our toolbox for quantum feedback control.

Signal Filtering
In the toy model above, the feedback potential is governed only by local in time measurement results. Because Eqs. (3)-(5) describe continuous time evolution, the effect of V fluct (x, t) in Eq. (9) would seem to diverge as dt → 0. However, any measurement signal M i (x, t) can be filtered in time to provide a running best estimate of the measured observable i (where i = n, z, etc.).
The resulting estimator ε i is derived from M i via the low-pass filter i.e.
where τ i is the filter time constant and M i (x, t) indicates the unfiltered measurement signal. This process filters the contribution of projection noise present at timescales below τ i , making τ i the effective measurement time associated with the estimator ε i . We derive all of our feedback potentials using estimators ε i instead of measurement signals M i , thereby controlling the noise applied to the system via feedback. In our feedback scheme we use separate estimators of the total density, spin density, or density in component s, denoted ε n , ε z , ε s , respectively, which can have different filter time constants τ n , τ z , and τ s .

C. Feedback Induced Magnetic Phases
We now focus on feedback-tuned spin-dependent interactions with g 2 = 0 and g 0 = 0. Guided by our toy model, we expect the steady-state phase diagram of a two-component BEC to resemble the ground state phase diagram for u 2 . The ground state density n(x) and spin density S z (x) are shown in in Fig. 1 (a). For u 2 > 0, the ground state is an easy-plane ferromagnet with S z (x) = 0, while for u 2 < 0 the ground state is an easy-axis ferromagnet, consisting of spin-polarized domains [25,[43][44][45], separated by a domain wall.
Using the measurement and feedback procedure outlined in Sec. II B 1, we apply a forcing potentiaľ along with a cooling potentialV c , to be described in Sec. III. Equation (15) changes the effective spin dependent interaction strength via the gain g 2 , based the estimator of the spin density ε z . The effective Hamiltonian for this protocol iš The phase diagram is now a function of two variables: spin-dependent interaction strength u 2 and signal gain g 2 , which give an effective interaction strength u eff 2 ≈ u 2 + g 2 . Examples of the two steady-state phases are shown in Fig. 1 (b). Both phases have uniform density, but with very different spin character. For u eff 0, the system is an easy-axis ferromagnet with well defined, spin polarized domains. For u eff 2 0 the system enters a spin-disordered paramagnetic phase, with large spin fluctuations. Fig. 1 (b) shows the spin density averaged over 100 ms (darker solid curve) and ten individual time traces (semi-transparent curves). The individual time traces show that the spin is essentially static in the ferromagnetic phase but has large spatiotemporal fluctuations in the paramagnetic phase. Figure 1 (c) shows the steady-state phase diagram as a function of u 2 /u 0 and g 2 /u 0 . As expected, the phase diagram is divided into two regimes delineated by u eff 2 = 0 (black dashed curve). We quantify the steady-state phase using a time-separated correlation function of magnetization, where A is an overall normalization factor. A condensate with well defined domains gives η 0.5; for the ground state with a single domain wall η ≈ 1. The disordered paramagnet phase with fluctuating magnetization has η ≈ 0, because the local magnetization at any point x fluctuates strongly in time.
Like many magnetic systems, this system exhibits hysteretic behavior. When g 2 < 0, the easy-axis phase is robust to the initial condition of the system and over many different repetitions of the simulation with different noise realizations. The phase in the region where u eff 0 with u 2 < 0 and g 2 > 0 is sensitive to the initial state, denoted by the hatched region in Fig. 1 (c). In this region, the steady state of the system is an easy-axis ferromagnet only if it was initially in the ferromagnetic ground state with u 2 < 0, as in 1 (a.i). For the easy-plane ground state, as in Fig. 1 (a.ii), domains do not form. We discuss this steady-state behavior for the easy-plane initial condition in Appendix B.
In the following sections we examine the robustness of the feedback-induced magnetic phases and feedback cooling. We show that despite repeated weak measurements and feedback the condensate remains largely intact over the ∼ 4 s time period of the simulation. Furthermore, by changing the effective interaction via feedback we demonstrate tunability between different steady-  state phases. Spatially-resolved, time dependent feedback therefore provides a tool to dynamically change effective interactions in cold atom systems.

III. FEEDBACK COOLING
Measurement backaction adds excitations to the condensate. The aim of feedback cooling is to apply feedback using information from the measurement signal to suppress the excitations, thereby stabilizing the condensate and preventing runaway heating. In this section we develop a feedback cooling protocol for single and multicomponent condensates which ensures the stability of the condensate during measurement and feedback. We connect the continuous measurement limit presented in Sec. II A to the experimental reality of discrete measurements. We then develop a feedback cooling protocol using a single discrete measurement as a building block. Finally, we show that during this protocol the condensate fraction and entropy reaches a steady-state, but the GPE energy functional continues to slowly increase.

A. Single Measurement Protocol
The continuous measurement limit is typically assumed a priori by taking dt → 0. Since the variance of the measurement signal in Eq. (1) is ∝ 1/dt, the variance in the measurement record diverges in this limit. However, no physical measurement is infinitely fast. Integrating Eq. (1) over a small time window therefore yields a 'single measurement'. By considering this type of measurement, we can quantify a measurement protocol which extracts maximal information from the condensate while minimizing the negative effects of backaction. As in Sec. II B 1, here we consider measurements of a single component condensate and drop the s index. It is straightforward to generalize this procedure to multicomponent condensates.
Consider a time-integrated version of Eq. (1) over an interval ∆t, giving a single measurement of density. The measurement result is M(x) = n(x) +m(x)/κ, where the measurement strength κ = √ ∆tϕ. The spatial quantum projection noise ism(x) wherem k has the same Fourier space statistics previously discussed, withm k = 0 andm kmk = Lδ kk Θ(|k| − k c )/2. Directly after measurement, the updated wavefunction is such that the measurement outcome matches the postmeasurement density n |M exactly, i.e. M(x) = n |M (x). In principle, the optimal measurement strength depends on the local density, however as this is difficult to implement experimentally we instead approximate κ * to be constant. We then use this coupling value for feedback cooling.
If we could find a potential V c|M (x) for which the post measurement state is the ground state, ψ |M (x) would sat-isfy the stationary GPE In our feedback cooling protocol, we first apply the potential V c|M (x) for which the post-measurement state would be the ground state (assuming a uniform phase). Then we approach the initial state by slowly -adiabaticallyramping off the applied cooling potential. We approximate V c|M using the Thomas-Fermi (TF) approximation of Eq. (19), giving V c|M (x) = µ − u 0 n |M (x). We then make the substitution u 0 n |M (x) → g c M(x) where g c is the cooling gain, an externally adjustable parameter (for which the expected value of u 0 is found to be optimal). This gives the feedback cooling potential function where t m is the time of the measurement and f (t) is a ramp off function where f (0) = 1 and f (t → ∞) = 0. In where γ is the ramp-off rate.

B. Bogoliubov Theory for Single Measurement Protocol
Here we provide an analytical solution of the singlemeasurement-feedback protocol described above using Bogoliubov theory [46], with periodic boundary conditions. After making the Bogoliubov transformation, small excitations above the ground state of a weakly interacting spinless BEC with density n are described by the HamiltonianĤ whereb † k describes the creation of a Bogoliubov phonon with momentum k and energy k = µξ|k| ξ 2 k 2 + 2. To facilitate our analytic treatment, we focus on the weak measurement regime, in which at most one phonon mode is occupied, leading to wavefunctions of the form |ψ = α|vac + k β k |k , where |k =b † k |vac , and |vac is the phonon vacuum.
Measurement backaction is described by the Kraus operator with the density difference operator δn(x) ≡n(x) − n.
In the phonon basis δn(x) can be expressed as a sum δn(x) = n/L k (c k e −ikxb k + h.c.) of phonon creation and annihilation operators, with c k = [1 + 2/(ξk) 2 ] −1/4 . In this representation, the feedback cooling operator derived from (20) iŝ Assuming adiabatic evolution, with ramp-off rate γ → 0, and using first order perturbation theory, the operator describing the cooling protocol iŝ This expression is valid for g c c k √ n κ k √ L. The probability of finding a phonon in state |k after a measurement-feedback cycle is (25) We draw two conclusions from this result: (1) Setting g c = 0 gives the probability nκ 2 c 2 k /2 that the measurement created a phonon in state |k ; and (2) the phonon mode with energy k,opt = g c κ −2 can be perfectly cooled with this protocol. Figure 2 (a) compares Eq. (25) with our stochastic GPE simulation with a linear ramp-off function f (t). The analytic calculation exactly reproduces the numerically predicted phonon distribution immediately following a single measurement (red curve), while the results with cooling have additional periodic features resulting from the finite ramp-off rates in the simulations. The shaded region denotes the parameters for which our perturbation theory is inapplicable.
In the thermodynamic limit L ξ, the per-particle energy after one measurement-feedback cycle is parabolic. With ξ 1/k c , the minimal per-particle energy increase ∆E * /µ = κ 2 φ 2 c (πφ c −6 √ 2)/(6π 2 ξ) occurs for a gain where φ c = k c ξ/ √ 2 parameterizes the cutoff and A = (4 √ 2κ 2 µξ) −1 . Figure 2 (b) compares the optimal energy increase predicted by Eq. (26), with that obtained from numerical simulations of the stochastic GPE (horizontal black dashed line and black squares, respectively), and the corresponding optimal gains are denoted by the red circles. The GPE simulation exhibits three regimes: (1) For very rapid ramps γ → ∞, the adiabatic assumption is invalid, and the GPE optimal gain is larger than anticipated from analytic model. (2) In the adiabatic ramping regime where γ → 0, we find both g c * and ∆E * converge, with ∆E * greater than our predicted value. This results from phonon-phonon scattering processes redistributing phonons between modes, which is not included in our Bogoliubov theory. And, (3) in the intermediate regime (γ between 3 ms −1 and 10 ms −1 ) our theory performs optimally and ∆E * coincides with the analytic prediction, albeit with much higher gain. We note that the optimal gain g c = u 0 obtained in Sect. III A is close to that predicted by Eq. (27), where for the parameters in Fig. 2, g c * ≈ 2.8u 0 .

C. Continuous Feedback Cooling Protocol
The single measurement procedure described in Sec. III A is a building block for continuous feedback cooling. We periodically measure the condensate with measurement strength κ = κ * ∆t/τ where κ * is the ideal single measurement strength in Eq. (18) and τ is the filtering time constant for the measurement signal. The cooling potential is derived from the density estimator ε(x, t) [47] and is decreased between measurements, as described by Eq. (20).
The effect of the cooling potential is to drive ψ(x) toward it's ground state between measurements. This procedure leverages the optimal single measurement strength and signal filtering to measure the condensate more weakly. We implement this protocol numerically and simulate condensate evolution under measurement and feedback using Eqs. (3)-(5).
Here we simulate an elongated condensate with N = 10 5 particles, healing length ξ = 0.8 µm and total system size L = 80 µm, computed for k c = 2π/λ with λ = 780 nm. The interval between measurements is set to dt = 200 µs to match typical image acquisition times in experiment, and the estimator time constant and cooling ramp-off rate were set to τ = 1/γ = 4.6 ms. We characterize the quasi-steady state by three metrics: conden-sate fraction, Von Neumann entropy, and energy, and find that the condensate remains remarkably coherent throughout the feedback cooling protocol. Upon implementing continuous feedback cooling, the condensate fraction and Von Neumann entropy reach a steady state while the GPE energy functional slowly increases, as shown in Fig. 3.
We calculate the condensate fraction using the Penrose-Onsager criteria [48]. Per this criteria, upon diagonalizing the one body density matrixρ asρ|n = N n |n , a condensate is present in mode |n if it's eigenvalue is N n ∼ O(N ) where N is the total number of particles. We obtainρ from an ensemble of stochastic trajectories of pure states [49], starting from the GPE ground state. In Fig. 3 (a) we show the four largest eigenvalues ofρ, normalized by N , giving a measure of the fractional occupation in each mode. The condensate fraction is the largest eigenvalue, which stabilizes at ≈ 0.99, with a secondary mode having an occupation fraction of ≈ 0.01. The remaining eigenvalues are orders of magnitude smaller than the leading two; therefore those modes have negligible occupation.
The second metric we use to characterize the steady state is the Von Neumann entropy, defined as S = Tr [ρ logρ]. As shown in Fig. 3 (b), S saturates at ≈ 0.01 of it's maximum possible value log(D), where D is the Hilbert space dimension. This is consistent with the final condensate fraction of ≈ 0.99. We extract an equilibration time τ eq ≈ 200 ms by fitting S to the function The third metric, energy, does not reach a constant value, rather it slowly increases even after the condensate fraction and entropy saturate, as shown in Fig. 3 (b). Here we define energy in terms of the per-particle GPE energy without any feedback terms present. The final energy after 4 s of evolution is ∼ 0.15 µ, indicating a 15% increase from the ground state value throughout the protocol. We determined that this energy increase is due to the gradual population of modes above the momentum cutoff which cannot be directly addressed by feedback cooling. However, this increase is slow enough to provide ample time (on the order of seconds) for additional experiments while the condensate is being measured.
Cooling for the two-component case proceeds similarly, but with cooling applied in the spin and density channels separately. Weak measurements add magnons (spin waves) in addition to phonons [27]. For the easy-axis ground state with u 2 < 0, the results are qualitatively the same as as the single component case, with the final condensate fraction reduced to ≈ 0.85, indicating cooling is not quite as efficient for the two-component system. However, in the easy-plane case (i.e. u 2 > 0), cooling is not as effective at long times and the condensate enters a spin-disordered phase with large spin fluctuations and a lower condensate fraction of ≈ 0.35. The cooling protocol for two-component condensates is discussed in Appendix C.  Fig. 1 (c) phase diagram with u eff 2 < 0. The colored markers indicate the calculated spin healing length averaged over 1.6 s window. The black markers indicate the spin healing length for a ground state system (i.e., no feedback) with u2 equal to the marked value of u eff 2 . The dashed curve indicates the predicted spin healing length ξs = ξ/ 2|u eff 2 /u0| with no fitting parameters.

IV. FEEDBACK INDUCED MAGNETIC PHASES
In this section, we elaborate on the steady state magnetic phases and their measurement signatures. The phase diagram in Fig. 1 (c) was computed for a gas of N = 10 5 87 Rb atoms with healing length ξ = 0.8 µm and total system length L = 80 µm, with feedback both to control the effective interactions and cool the system. In all of our simulations, feedback cooling is continuously applied. We add the forcing feedbackV f (x, t) = g 2 ε z (x, t)σ z in the time window from 1 s to 3 s and allow the simulations to continue until the total run time reaches 4 s. Figure 1 (c) shows that the magnetic phase of the system reaches a steady-state governed by the effective spindependent interaction strength u eff 2 = g 2 + u 2 while the forcing potential is on, leading to the easy-axis ferromagnet and spin-disordered paramagnetic phases discussed in Sec. II C. The spin-dependent interaction strength u 2 and gain g 2 serve as tunable parameters.
The easy axis ferromagnetic phase for u eff 2 < 0 exhibits well defined, spin-polarized domains. The order parameter η for this phase is the time-separated correlation function of the magnetization, given in Eq. (17). We find that η 0.5 indicates the existence of persistent domains. We can identify an effective spin healing length ξ s ∝ 1/ |u eff 2 | in this phase, similar to the spin healing length in closed two-component systems [43]. Changing u eff 2 via the feedback strength thus alters the spin healing length in the steady state. Figure 4 shows the effective spin healing length, obtained by fitting the spin density S z (x) to a function with N d domains, where Here, x n are the positions of each domain wall, S is the overall amplitude of domains, and ξ s is the spin healing length. The ± sign in front accounts for the polarity of the domain signal (i.e. which domain is at the edge), as the measurement and feedback process spontaneously breaks a Z 2 symmetry to determine the domain orientations [40,50]. The spin healing length diverges upon approaching the transition at u eff 2 = 0, indicated system behavior that is analogous to the expected phase transition from changing the interaction parameters. The markers in Fig. 4 are color-coded based on the value of the η, where we can see that for lower values there is more variability in the data. This is because lower values of η generally correspond to a spin texture with multiple domains, where there is movement of the domain boundaries over time due to fluctuations parameterized by the nonzero enropy [40]. The black diamonds in Fig. 4 show the spin healing length obtained for the corresponding closed system ground state, and the dashed curve is the computed functional dependence ξ s = ξ[u 0 /2|u eff 2 |] 1/2 for u eff 2 < 0, which shows excellent agreement with the simulations.
The disordered paramagnetic phase is characterized by a spatially and temporally fluctuating spin structure. An example of these fluctuations in real space is shown in Fig. 5 (a). In the disordered paramagnetic phase, a spin healing length is not well defined. The power spectral density (PSD) of the spin, provides a measure of how much the spin fluctuates [43]. HereS z (x) is the time-averaged value of the spin density andS z (k, t) is the Fourier transform of S z (x, t). Figure 5 (a) shows PSD z (k) in the steady-state magnetic phase averaged over 1 s. At low momenta the signature for the disordered phase is significantly higher than for the easy-axis ferromagnetic phase. The large fluctuations in spin are thus a signature of the paramagnetic phase which can be deduced from the measurement signals. Above the cutoff k c λ = 200π indicated by the black, dashed line, we see additional spectral features at multiples of k c , indicating higher-order resonances due to the measurement process. Population of modes above the cutoff leads to a gradual increase in energy and affects cooling, as discussed in Sec. III C.

V. OUTLOOK
Hamiltonian engineering for multicomponent Bose gases has been achieved at the level of the single-particle Hamiltonian via synthetic gauge fields [51,52], spinorbit coupling [10,53,54], and spin-dependent potentials [55,56]. The ability to tune the character and strength of interactions beyond those already present in the system has heretofore been limited to using Feshbach resonances [29], which typically change only one interaction constant at a time, or via coupling to an external cavity field [57][58][59]. In contrast, our feedback technique can simultaneously change all the spin-dependent effective interaction strengths in situ: not possible with Feshbach resonances or cavity mediated interactions.
Our result shows that spatially local feedback control based on a record of weak measurements is a viable route toward engineering effecting interactions in quantum gases. We demonstrated that a dynamical steady state can be engineered in a two-component Bose-Einstein condensate where the magnetic phase is determined by the interplay of the intrinsic and feedbackinduced interaction strengths.
Going beyond previous works [34,40], we implemented a cooling scheme which avoids runaway heating of the condensate during the feedback process. Further optimization of the cooling protocol will be important for experimental implementation. For example, Eq. (25) suggests that the k dependent gain g c (k) = nκ 2 k would lead to near-perfect cooling for all momentum states. Actual measurements have limited resolution, detector inefficiencies, and technical noise, which could possibly be addressed by further optimizing the cooling protocol.
The feedback control method of engineering effective Hamiltonians is flexible and allows for the introduction of tailored, spatially dependent effective interaction terms. Future work could implement nonlocal or timedependent interactions which have no analogue in closed systems. Our protocols can be generalized to higher dimensions, and could stabilize topological defects such as non-Abelian vortex anyons which are unstable in closed systems [23]. Finally, our methods enable real-time feedback control, so over the course of one experiment we can study both quasi-steady-state behavior and dynamics.

ACKNOWLEDGMENTS
This work was partially supported by NIST and NSF through the Physics Frontier Center at the JQI. HMH acknowledges the support of the NIST NRC postdoctoral program.

Appendix A: Simulation Parameters
Here we briefly review the simulation method for Eqs. (3)-(5) and the parameters we use in this work. All simulations have N = 10 5 atoms and we consider a quasi-1D system of length L = 80 µm with hard wall boundary conditions such that Ψ(x = −L/2) = Ψ(x = L/2) = 0. Hard-wall boundaries can be implemented using flatbottomed traps instead of a harmonic one [60]. The momentum cutoff is k c = 2π/λ with λ = 780 nm being the wavelength of imaging light. We simulate a single component condensate in order to study steady state behavior under feedback cooling in Sec. III. Elsewhere, we simulate a two-component condensate with an easy-axis magnetic ground state, i.e. u 2 < 0, or easy-plane ground state with u 2 > 0. In the main text results are presented using the easy-axis ground state with u 2 = 0.01u 0 as the initial condition.
The system is initialized in it's ground state by solving the GPE in imaginary time. The natural units for this set up are the total system length L and the chemical potential µ = 2 /2mξ 2 as the unit of energy where ξ = 0.8 µm is the healing length. Upon re-scaling the variables to unitless quantities x → xL, t → t(2m a ξ 2 / ), ψ ↑(↓) → N/Lψ ↑(↓) , the Hamiltonian in Eq. (19) iš where dx n(x) = 1. Therefore, the spinless case has one free parameter ξ/L and the two-component case has the additional free parameter u 2 /u 0 . For our parameters  we have ξ/L = 0.01 and we consider different values of u 2 . We simulate the nonlinear dynamics using a secondorder symplectic integration method [61]. In these units it is natural to express u 2 and the gain strengths g 0 , g 2 , etc in units of u 0 .
In order to simulate a small measurement interval (approaching the continuous measurement limit), we consider a separation of timescales dt τ such that the measurement interval dt of the system is much shorter than the signal filtering timescale τ for any observable. This enables us to write the evolution Eq. (3)-(5) as continuous time stochastic differential equations. As indicated by the hatched region in Fig. 1 (c) , the steady state phase diagram has a region of bistability depending on the initial state of the system. In this Appendix we present the results for the phase diagram calculated using the easy-plane ground state as the initial condition, shown in Fig. 6. In the steady-state magnetic phase, the system forms domains for u eff 2 < 0 and g 2 < 0. An example of the density and spin density in this region is shown in Fig. 6 (a.i), where we see that there are multiple domains in the spin texture. This is in contrast to the case presented in the main text where there is only one domain, due to the single-domain being the ground state. The number of domains depends on many parameters including u 2 , g 2 , and the timescale over which feedback is turned on. We consider further investigation of these variables to be outside the scope of this work.
Unlike the easy-axis initial condition, the spindisordered phase occurs for a wider range of parameters, most notably in the hatched region where u eff 2 = 0 but g 2 > 0. The spin texture in this regime is shown in Fig. 6 (a.ii), which indicates relatively uniform density but a highly fluctuating spin texture. We suspect that the observed bistability could be due in part to the underlying cooling protocol for the two-component system, which can also affect the spin texture, as discussed in Appendix C.

Appendix C: Two-Component Feedback Cooling
The density is measured in each component s with strength κ = κ * ∆t/τ n where ∆t is the measurement duration and τ n is the low-pass filtering time constant for the total density. Measurements M ↑ and M ↓ are then combined to give a measurement of total density (M ↑ + M ↓ ) or spin density (M ↑ -M ↓ ), which is used in a low-pass filter to calculate the estimators ε n and ε z . Crucially, the filtering works best when ε n and ε z have different filtering time constants; we use τ n = 4.6 ms and τ z = 46 ms, respectively. This is due to the different types of excitations in the two-component case, which can be phonons or magnons. Phonons have faster time dynamics than magnons, which necessitates different time constants in each channel.
The spin-dependent cooling potential iš As in the spinless case, the potentials V c,n and V c,z are calculated after each measurement and then exponentially ramped off between measurements. Cooling in the density channel is done via the potential where g c is the gain. This potential drives the total density toward a uniform state based on estimator ε n with Fractional occupation of first four modes in the single-particle density matrix for (a) u2 < 0 and (b) u2 > 0. The eigenvalue of the four highest-occupied modes is pictured. The condensate fraction (solid curve) is ≈ 0.85 in the steady state for u2 < 0 and ≈ 0.35 for u2 > 0. (c) Average energy (black) for a condensate with u2 < 0 (solid curve) and u2 > 0 (dashed curve) calculated from 124 independent stochastic trajectories. As in the spinless case, energy computed from the GPE energy functional increases slowly (d) Von Neumann entropy for a condensate with u2 < 0 (solid curve) and u2 > 0 (dashed curve).
ramp-off rate γ n . Cooling for the spin sector is via the spin-dependent potential V c,z (x, t) = g c,z [ε z (x, t) − ε z (x, t)] e −γz(t−tm) , (C3) where γ z is the spin ramp-off rate, g c,z is the cooling gain for the spin sector, andε z indicates a running time average of ε z . This potential drives the spin density S z (x) toward it's time-averaged value, effectively cooling short wavelength (high momentum) spin fluctuations but allowing long-wavelength spin textures such as domain walls to remain intact. In practice we use γ −1 n = τ n and γ −1 z = τ z , with the other parameters the same as for the spinless case. We calculateε z by averaging the original signal over a 120 ms time window. Cooling is most effective when the gain parameters are g = u 0 and g c,z = u 2 .
As in the spinless case, feedback cooling drives the twocomponent condensate to a quasi-steady state. Condensate fraction and Von Neumann entropy stabilize around constant values and the energy per particle increases slowly over the course of the simulation. We compute the energy from the GPE energy functional without any feedback terms present. The steady-state properties for cooling a two-component condensate are presented in Fig. 7. The results are qualitatively different for the case with u 2 < 0 (easy-axis ground state) and u 2 > 0 (easy plane ground state).
The easy-axis case is similar to the spinless cooling results presented in the main text. In Fig. 7 (a) we present the condensate fraction for u 2 < 0, which can also be calculated for multicomponent condensates [62]. The condensate fraction is ≈ 0.85 in the steady state with one additional mode having occupation ≈ 0.15 and other modes having negligible occupation. The energy increase, shown in Fig. 7 (c) is ≈ 0.25µ. The Von Neumann entropy, shown in Fig. 7 (d) (solid curve) increases to about 10% of it's maximum value. These metrics indicate that the cooling protocol is effective for two-component condensates with u 2 < 0. Furthermore, we find that at the end of the cooling protocol the domain wall is still intact, showing that this spin dependent cooling protocol is effective both at maintaining a high level of condensation and preserving the spin structure. The equilibration time extracted from the entropy is τ eq ≈ 400 ms.
In the case of an easy-plane initial condition (i.e. u 2 > 0), the cooling protocol is not as effective. In Fig. 7 (b) we show the fractional occupation of the first four modes from the one-body density matrix. The condensate fraction (blue, solid curve) decreases to ≈ 0.35 while the other modes also have fractional occupations of O(0.1). This indicates that the Penrose-Onsager criterion for condensation is violated in this regime. Furthermore, we find that the entropy S increases considerably more than the easy-axis case, reaching a constant value of ≈ 0.4 log(D) after 2 s of time evolution. The entropy increase is likely being driven by an instability toward spin separation in the condensate. Under our current feedback protocol, the easy-plane ground state eventually enters a spin-disordered phase with large spin fluctuations, which accounts for the higher entropy and lower condensate fraction we observe. Future work could develop a feedback cooling protocol specifically for u 2 > 0 systems to combat this instability more effectively.