Experimental observation of a dissipative phase transition in a multi-mode many-body quantum system

We characterize the dissipative phase transition in a driven dissipative Bose-Einstein condensate of neutral atoms. Our results generalize the work on dissipative nonlinear Kerr resonators towards many modes and stronger interactions. We measure the effective Liouvillian gap and analyze the microscopic system dynamics, where we identify a non-equilibrium condensation process.

One of the paradigmatic examples exhibiting a dissipative phase transition is a single mode quantum system, realized, e.g., with a nonlinear optical medium in a cavity, which is coherently pumped by a laser and subject to cavity losses [6,[15][16][17]21]. The Hamiltonian part of such a Kerr type system readŝ H = ∆â †â + Uâ †â †ââ + F * â + Fâ † , where ∆ is the detuning between the mode energy and the photon energy of the pump laser,â is the annihilation operator of the mode, F is the amplitude of the pumping laser and U denotes the interaction. The cavity losses are included via a master equation in Lindblad form, where γ diss is the loss rate andρ is the density operator. This model exhibits several prominent features: (i) There are two classes of states (with low and high occupation number) which are connected by a first order dissipative phase transition. (ii) While a mean-field treatment predicts optical bistability, the full quantum description (Eqs. (1), (2)) allows for tunneling between the steady states, thus erasing bistability and rendering the states metastable. (iii) As a consequence of (i) and (ii), single trajectories show a pronounced switching behavior between the two classes of states [6].
For multi-mode systems instead, much less is known. There exist only a few theoretical works [3,4,[22][23][24][25], mostly on the driven-dissipative 2D Bose-Hubbard model. The extension to many modes in the interacting many-body regime poses many challenges to its theoretical description (see methods for an explicit form of the master equation). In practice, already a few tens of modes and particles make exact diagonalization intractable. Experimental studies of a dissipative phase transition in a multi-mode system therefore not only generalizes the results of the single mode Kerr oscillator, but also provide benchmarks for advanced theoretical methods to describe driven-dissipative systems.
In this work, we investigate the multi-mode version of the nonlinear Kerr resonator in an atomic quantum gas experiment. First results for the system under investigation have been reported previously by one of the authors [26,27]. We present experimental indications that this system undergoes a first order dissipative phase transition when tuning the ratio between the coherent drive and the dissipation. We observe the same characteristic switching behavior as in the single mode system and find signatures of a closing of the Liouvillian gap. On intermediate timescales, the system exhibits hysteresis and metastability. In particular, we identify and analyze a non-equilibrium condensation process. The system has triggered recent theoretical work using c-field methods [24] and a single-mode variational truncated Wigner approach [25].
The experimental system (Fig. 1a) consists of a Bose-Einstein condensate loaded into a 1D optical lattice. This can be considered as an array of coupled two-dimensional Bose-Einstein condensates of neutral atoms. Each condensate occupies about 80 transverse harmonic oscillator modes. With a total of about 800 atoms in each condensate, the system is deep in the multi-mode regime with high mode occupation numbers. The central lattice site is subject to a homogeneous loss process with rate γ diss . The loss mechanism is realized by a focused electron beam, which is scanned across the central site and removes the atoms. A fraction of the atoms is converted into ions [28], which form the basis of our measurement signal. This way, we can measure the time-resolved atom  1. (a) A Bose-Einstein condensate is loaded into a onedimensional optical lattice. The individual sites of the condensate array are coupled via the coherent tunneling coupling J. The central site is subjected to particle loss, realized with an electron beam. The nonlinear interaction leads to an effective tunneling coupling J eff to the neighboring sites (see text). The produced ions are guided towards a detector and provide a time-resolved measurement signal. (b) Single experimental trajectory (amber line, moving average is shown as dotted line) for an initial state with low atom number in the central site. After 200 ms evolution time, the trajectory switches to a high atom number, which is comparable to that of the neighboring sites. This jump is fitted with a sigmoid function (indigo dashed line). (c) Same as (b) for an initially high atom number in the central site.
number and the transverse extent of the atomic cloud in the central site. The lattice sites are connected via quantum mechanical tunneling with rate J, which constitutes the drive into the central site. The experimental sequence starts by initially preparing the system either with an "empty" site (low occupation number in the central site) or with a "full" site (same occupation number in all sites). We then record the temporal evolution of the site occupancy under the simultaneous influence of drive and dissipation.
In Fig. 1, we show typical individual trajectories observed in the experiment for an initially empty (Fig. 1b) and full central site (Fig. 1c). The trajectories are characterized by a switching behavior, indicated by pronounced jumps in the site occupancy between the two states. The absolute magnitude of these jumps is large: about 700 atoms tunnel within a short time into or are removed from the central site. This switching behavior is reminiscent of that of the corresponding single mode system [6] and hints at a dominant role of fluctuations. It has also been seen in a c-field treatment of this system [24].
The nature of the appearing states between which the system switches has been identified previously [27,29]. For high occupation number, when the coherent drive dominates over the dissipation, there is a condensate in the central site, which is phase locked to the neighboring condensates. The whole system is in a superfluid state, all sites have similar atom numbers and the system exhibits steady-state transport known as coherent perfect absorption [29]. In the opposite limit of large loss rate, the central site features only a small occupation in each mode. The atomic ensemble is in a normal state and macroscopic phase coherence is absent. Hence, superfluid transport into that site is inhibited [27]. In the vicinity of the phase transition, where dissipation and drive are of similar strength, jumps are induced between the states.
To quantify the switching dynamics, we measure the time stamps of the jump for each trajectory. We then determine the time scales τ normal→superfluid = τ ns and τ superfluid→normal = τ sn for switching between the two configurations as well as the corresponding rates γ ns = τ −1 ns and γ sn = τ −1 sn . The physical meaning of the switching rates becomes clear if one looks at the spectral decomposition of the Liouvillian Here, λ i denotes the complex eigenvalue, and Re(λ i ) ≤ 0 is the damping rate ofρ i . Eigenvalues with Re(λ i ) = 0 denote steady states. Note that for any non-zero Re(λ i ), theρ i are not proper density matrices, as the trace is not preserved in time. A proper density matrix therefore always contains a contribution from at least one steady state. If two eigenvalues of Eq. (1) are zero, the system exhibits true bistability. If only one eigenvalue is zero, the spectrum of Eq. (3) shows a gap, the Liouvillian gap, which is defined by the smallest nonzero eigenvalue, min(− Re(λ i )). The Liouvillian gap therefore also defines the slowest relaxation rate of the system. As the full many-body problem is intractable, we adapt a simplified toy model, which was developed in [30] for the driven dissipative Bose-Hubbard model and applied in [16]. The toy model considers only two states between which the system can switch: the normal state and the superfluid state, represented by the density matrices ρ s and ρ n . The corresponding rate equation model features two steady states. One is a stable steady-state with eigenvalue λ = 0: ρ ss = (γ ns ρ s +γ sn ρ n ). The second one is an exponentially decaying state, ρ dec = (γ ns ρ s − γ sn ρ n ), with eigenvalue λ gap = −(γ ns + γ sn ). Hence, the toy model features an effective Liouvillian gap, given by λ gap . Because the steady-state ρ ss is a mixture of the normal and the superfluid state, a time resolved measurement of the site occupancy is characterized by jumps between the two states.
In the experiment, we can observe the system only for a finite time because of atom loss. As a consequence, we typically observe only one jump in a single trajectory (in very rare cases, we observe two). Therefore, we chose two different initial conditions to determine the rates γ ns and γ sn individually, see also Fig. 1. In Fig. 2, we show the resulting effective Liouvillian gap, γ ns + γ sn , as a function of the tunneling coupling for fixed dissipation strength. At the minimum, the gap is smaller than one Hertz. This is two to three orders of magnitude smaller than the dissipation rate γ diss and the tunneling coupling J/ , which are the intrinsic time scales of the system. The large difference in time scales is a strong indication of a dissipative phase transition and we take the minimum of the relaxation rate to be the position of the phase transition. Around the critical point we find metastability of the two initial states, which manifests itself in the emergence of intermediate hysteresis for finite measurement times (Fig. 3). Note that only in the thermodynamic limit, which corresponds in our case to an infinite transverse extent of all lattice site, the Liouvillian gap is expected to close completely [1,31], thus providing bistability at the critical point. The toy model introduced above neglects all density matrices except of ρ ss and ρ ns . Therefore, the effective Liouvillian gap is only an upper limit for the true Liouvillian gap. Strictly speaking, the large number of allowed density matrices makes it practically impossible to experimentally measure the true Liouvillian gap. Interestingly, experiment and theory encounter similar constraints in this respect.
So far, our findings are qualitatively similar to those for a single mode nonlinear Kerr resonator. We now turn to the question, what genuine characteristics emerge from the multi-mode nature of our system and make it distinct from the single-mode system. We start our discussion by analyzing the jumps in more detail. Because we are using interacting bosons, we expect to pass a condensation process during a jump from the normal to the superfluid state. Such a condensation process happens under non-equilibrium conditions, characterized by the competition between the coherent drive, atomic collisions/interactions and dissipation. Information on this non-equilibrium condensation process can be obtained from the time resolved site occupancy during the change and the transverse extent of the atom cloud. Figure 4 shows the averaged time-resolved detection signal as well as the transverse extent of the atomic cloud during state changes from the normal to superfluid state. We can identify three different stages. After an initial slow exponential growth of the atom number a sudden rise in the growth rate is apparent. An exponential fit to the initial part of the averaged trajectory is shown to highlight the speed-up in the filling rate above a critical atom number. In a final step, the filling of the central site eventually equals that of the neighboring sites and remains constant. It is instructive to look at the timescale of this process by analyzing the duration of the state change in each trajectory individually. We find that the shortest duration for a state change from the normal to the superfluid state is given by the tunneling time, τ ≈ /J (Fig. 5a). This can be explained as follows. For a small atom number in the central site we only expect an incoherent transport of particles into that site. However, when the density in the central site exceeds a critical value, high mode occupancy builds up. Fluctuations then trigger a condensation process, thus allowing for a coherent Josephson type of transport into the central site. For fully established phase coherence across the central site, the maximum possible coherent transport is expected. For scrambled and fluctuating phases, much longer durations are possible. We observe that the average duration of the state change is about four tunneling times and the longest duration for a state change can be even 20 times larger. In the other direction, from the superfluid to the normal state, the dissipation rate is the determining factor (Fig. 5b). No state change can be faster than the inverse of the dissipation rate and this is also our experimental observation. The average state change is again about two to five times longer and few trajectories show even longer durations. For both types of jumps, such large variations hint at the dominant role of fluctuations for the system dynamics.
Another important observation is that history appears to be relevant in our system. If the state changes were fully independent of the elapsed evolution time, the temporal distribution of the state changes would follow an exponential decay [32]. In our experiment, however, the temporal distribution of state changes is non-trivial (Fig. 6). We find small probabilities of a state change for early times, increasing probabilities as time elapses further followed by a final slow decay of probability. This supports our previous observation in Fig. 4, that a state change has preceding internal dynamics.
To explain the above reported observations for state changes from the normal to the superfluid state, we have developed a simple stochastic rate model, which describes the microscopic system dynamics in the normal state. As indicated in Fig. 1, the mean-field energy (chemical potential) in the neighboring sites renders the direct tunneling process into the ground state of the central site offresonant. In a single-mode Josephson array, this would lead to macroscopic quantum self-trapping, thus inhibiting any transport [33]. The multi-mode nature of our system changes this picture. Atoms are allowed to tunnel into higher orbital states of the central site, where they are subject to collisions between particles in different orbital states. Due to a Franck-Condon overlap between the condensate wave function in the neighboring site and the single particle states in the central site, the effective tunneling coupling is reduced and depends on the atom number in the central site [26]. Overall, we get an effective tunneling rate Γ eff , which depends on the dissipation rate γ diss and the effective tunneling coupling J eff (N ), where N is the number of atoms in the central site. We thereby assume that the effective tunneling coupling depends linearly on the atom number and reaches the value of J for equal filling with the neighboring sites [26]. We find (see methods for further details): where γ coll is the collision rate in the lossy site in the normal state. The stochastic differential equation describing the population N in the lossy site then reads where N F is the atom number in the neighboring sites. The term σdW t with σ = √ γ diss describes the fluctuations induced by the dissipation. The dynamics resemble those of geometric Brownian motions (for more details on the model see methods section). Note that due to the N -dependence of the effective tunneling rate, the rate equation is nonlinear. To simulate the system we time evolve many individual trajectories by Monte-Carlo simulations. When a trajectory reaches a threshold filling for the first time macroscopic phase coherence builds up and the jump occurs. In Fig. 7a, we show a series of simulated trajectories. We take the time of the jump to be the time the threshold is reached. The introduction of a threshold filling followed by a condensation process gets backing from the measurement of the width of the atomic cloud within the central site during the jump (Fig. 4): All atoms tunneling from a neighboring reservoir into the central site have total energy equal to the chemical potential µ of the neighboring sites. In the beginning of the refilling process, when the central site contains barely any atoms, the interaction energy is converted into kinetic energy, leading to an effective finite temperature. From an approximation of the density distribution of the thermal part inside the central side by a Gaussian function the width of the cloud in the central site is calculated to be σ th = k B T / (mω 2 ).
Here ω denotes the radial trap frequency, k B the Boltzmann constant, m the mass of rubidium atoms, and T the temperature in the site. The condensate's radial extent in the central site is approximated by the Thomas-Fermi radius r TF = µ/ (2mω 2 ). The calculated widths for our system are 2.1(1) µm for the normal/thermal gas and 1.7(1) µm for the condensate. Thus, the transverse extent of the cloud shrinks with increasing atom number. This is a consequence of the condensation process and the reduction of thermal fluctuations. The experimental results (indigo line in Fig. 4) are in good agreement with this prediction.
In Fig. 7b we compare the temporal distribution of the simulated jumps with our measurement and in Fig. 7c we show the comparison for the fraction of trajectories exhibiting a jump. The only fitting parameter of our model is the threshold filling and we find a critical value of 300 atoms, which is in good agreement with the value extracted from the measurement (Fig. 4).
This has to be compared to the critical atom number for Bose-Einstein condensation in thermal equilibrium. 2D BEC theory [34] gives a value of 100 atoms. This confirms that the fluctuations increase the critical atom number for condensation in the central site and underlines the qualitative difference between a dissipative phase transition and a phase transition in thermal equilibrium.
In summary, we have experimentally characterized a first-order dissipative phase transition in a multi-mode many-body quantum system. Our results extent the experimental studies on nonlinear optical Kerr resonators. While some features are qualitatively similar, the microscopic dynamics are much richer due to the multi-mode nature. For the future, we expect to see more experiments working on related systems, such as the driven dissipative 2D Bose-Hubbard model or corresponding fermionic counterparts. For the system at hand, a detailed study of the scaling of the Liouvillian gap with the system size would allow for a finite size scaling towards the thermodynamic limit. The measurement of critical exponents in the vicinity of the phase transition would allow for a deeper understanding of the classification of the phase transition [11].

The master equation
The master equation of our system is constructed in the following way. We first note that the atoms occupy the lowest band along the lattice direction. We are well within the tight binding approximation along this direction and introduce the discrete site index i. The central site corresponds to i = 0. The kinetic energy in the lattice direction results in the tunneling coupling J. In the two transverse directions (x, y) each site experiences a harmonic confinement with radial frequency ω ⊥ . The Hamiltonian of the central site in terms of bosonic field operators (Ψ 0 (x, y) =Ψ 0 ) readŝ The effective interaction strength in the lattice site is given by g 2D . The coherent drive is given by a mean-field wave function of the condensate in the two neighboring sites given by Φ(x, y)e −iµt/ , where µ is the chemical potential. The large number of additional adjacent reservoir sites justifies to set Φ(x, y) constant in time.
The losses are included via a master equation, As outlined in the main text, an exact theoretical treatment is not possible.
A stochastic mean-field description of this master equation has been developed and analyzed in detail in Ref. [24].

Experimental setup
We use a Bose-Einstein Condensate (BEC) of 87 Rb containing about 150 × 10 3 atoms. The atoms are trapped in a focused CO 2 -laser beam with a beam waist of 34 µm. The condensate is cigar-shaped with dimensions (115 × 6 × 6)µm in an approximately harmonic trap with frequencies 2π · (9.8, 197, 183)Hz. After its preparation the BEC is loaded into a one-dimensional blue detuned optical lattice created by two beams (λ = 774 nm, beam waist 500 µm) intersecting at 90 • . The resulting lattice has a period of d = 547 nm. Each site contains a small 2D BEC with about 800 atoms at the center of the trap. The total number of lattice sites is about 200. In addition, an electron column implemented in the vacuum chamber provides a focused electron beam which is used to introduce a well-defined local particle loss as a dissipative process at one site of the lattice. In order to ensure a homogeneous loss process over the whole extent of this lattice site we rapidly scan the electron beam in the transverse direction (3 kHz scan frequency) with a sawtooth pattern. To adjust the dissipation strength γ, we vary the transverse extent of the scan pattern. In the experiments considered in this paper we start out with a BEC, load it into the 1D lattice at a depth of s = 30 (we give the values of the lattice depth V 0 = s E rec in units of the recoil energy E rec = π 2 2 / 2md 2 ) and -depending on the initial conditions -deplete one site. We then reduce the lattice depth to a value between s = 6 and s = 12 while keeping the electron beam scanning over the depleted site. Our measurement signal is the time-resolved ion signal from the scan at the lower lattice depth. Only in the case of the dynamical hysteresis measurement, the lattice depth is varied in time (for the measurement presented in Fig. 3a the lattice depth was varied at a constant rate from s = 4 to s = 20 and back).

Microscopic model
In order to develop a rate model for the dynamics in the central site in the normal state (low filling), we start out from a two-level system. The neighboring site corresponds to the ground state. Coherent coupling to a radially excited state in the central site is provided by an effective tunneling coupling J eff , which takes into account the finite Franck-Condon overlap between the neighboring reservoir site and the radially excited state [27]. Atom loss from this radially excited state is decribed by the dissipation rate γ diss and dephasing caused by atomic collissions is included via γ coll . The rate equation for the central site is then given by where N is the number of atoms in the central site and N F is that of the neighboring sites. Due to the large number of reservoir sites to which the neighboring site is coupled we set N F as constant. The effective tunneling rate Γ eff in this effective two-level system is then given by In our system atoms can tunnel between the full site (|Φ F ) and the resonant mode (|Φ E ) of the empty site with the effective tunnel rate Γ eff . Introducing the overlap integral η = | Φ F |Φ E |, and assuming a linear dependence on the atom number for finite filling, we rewrite the effective tunnel as follows [26]: For low filling in the central site, as it is the case in the normal state, we linearize the above expression: To model the fluctuations induced by the dissipation we introduce a stochastic term and obtain the stochastic differential equation dN = (λN + a) dt + σdW t (12) with a Wiener process W t and σ = √ γ diss .
To solve Eq. (12) we use the Itô formula, a variant of the chain rule for stochastic multi-variable functions [35][36][37]. The solution to this SDE is then given by with a standard normal variable x and c = a/λ.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.