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

Dissipative phase transitions [1–6] are characteristic of open systems, where the coupling to an external reservoir provides a source of dissipation which competes with an internal or external pumping process. Depending on the relative strength of the two processes, different steady states emerge or even coexist. Most notably, the steady-states can show properties which cannot be observed in equilibrium [7]. The generic properties of first and second order dissipative phase transitions can be derived from a spectral analysis of the system's Liouvillian [1, 2]. Critical behavior at the phase transition has been studied both, theoretically [2, 4, 8–12] and experimentally [6, 12–20]. 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–17, 21]. The Hamiltonian part of such a Kerr type system reads

$\begin{equation}\hat{H}={\Delta}{\hat{a}}^{{\dagger}}\hat{a}+U{\hat{a}}^{{\dagger}}{\hat{a}}^{{\dagger}}\hat{a}\hat{a}+{F}^{\ast }\hat{a}+F{\hat{a}}^{{\dagger}},\end{equation} \tag{ 1 }$

where Δ is the detuning between the mode energy and the photon energy of the pump laser, $\hat{a}$ 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,

$\begin{equation}\begin{aligned}\hfill \dot{\hat{\rho }}& =\hat{\mathcal{L}}\hat{\rho }\hfill \\ \hfill & =-\frac{\text{i}}{\hslash }\left[\hat{H},\hat{\rho }\right]+\frac{{\gamma }_{\text{diss}}}{2}\left(2\hat{a}\hat{\rho }{\hat{a}}^{{\dagger}}-{\hat{a}}^{{\dagger}}\hat{a}\hat{\rho }-\hat{\rho }{\hat{a}}^{{\dagger}}\hat{a}\right),\hfill \end{aligned}\end{equation} \tag{ 2 }$

where γ_diss is the loss rate and $\hat{\rho }$ 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 (equations (1) and (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–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 generalize 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 (figure 1(a)) consists of a Bose–Einstein condensate (BEC) loaded into a 1D optical lattice. This can be considered as an array of coupled two-dimensional BECs 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 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.

Zoom In Zoom Out Reset image size

In figure 1, we show typical individual trajectories observed in the experiment for an initially empty (figure 1(b)) and full central site (figure 1(c)). 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 ${\gamma }_{\text{ns}}={\tau }_{\text{ns}}^{-1}$ and ${\gamma }_{\mathrm{sn}}={\tau }_{\mathrm{sn}}^{-1}$. The physical meaning of the switching rates becomes clear if one looks at the spectral decomposition of the Liouvillian

$\begin{equation}\hat{\mathcal{L}}\hat{{\rho }_{i}}={\lambda }_{i}{\hat{\rho }}_{i}.\end{equation} \tag{ 3 }$

Here, λ_i denotes the complex eigenvalue, and Re(λ_i ) ⩽ 0 is the damping rate of ${\hat{\rho }}_{i}$. Eigenvalues with Re(λ_i ) = 0 denote steady states. Note that for any non-zero Re(λ_i ), the ${\hat{\rho }}_{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 equation (1) are zero, the system exhibits true bistability. If only one eigenvalue is zero, the spectrum of equation (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. This characteristic behaviour of the Liouvillian gap is sketched in the inset of figure 2.

Zoom In Zoom Out Reset image size

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 ${\hat{\rho }}_{\text{s}}$ and ${\hat{\rho }}_{\text{n}}$. The corresponding rate equation model features two steady states. One is a stable steady-state with eigenvalue λ = 0: ${\hat{\rho }}_{\text{ss}}$ = (γ_ns ${\hat{\rho }}_{\text{s}}$ +  γ_sn ${\hat{\rho }}_{\text{n}}$). The second one is an exponentially decaying state, ρ_dec = (γ_ns ${\hat{\rho }}_{\text{s}}$ − γ_sn ${\hat{\rho }}_{\text{n}}$), with eigenvalue λ_gap = −(γ_ns + γ_sn). Hence, the toy model features an effective Liouvillian gap, given by λ_gap. Because the steady-state ${\hat{\rho }}_{\text{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 figure 1. In figure 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 1 Hz. 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 (see also the inset of figure 2, where the typical behaviour of a first order dissipative phase transition is sketched).

In the thermodynamic limit only at the critical point of the first order phase transition will the steady state be a mixture of the two phases. Right before and after the critical point the steady state is one of the two different phases of the system. For a finite size system a small region around the critical point is expected to be metastable, which may introduce hysterical behavior. This is typical of first-order phase transitions. In figure 3(a) we show that the normal (superfluid) state is the steady state of the system on the right (left) of the critical point. A region of metastability is found around this point. To explicitly show the emerging hysterical behavior, we measured the dynamic hysteresis by ramping the tunneling coupling up and back down (figure 3(b)).

Zoom In Zoom Out Reset image size

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 ρ_s and ρ_n. 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. The possibility of a change of the transverse extent of the cloud is a direct consequence of the multi-mode nature of the system and is therefore not expected to occur in single-mode systems.

Zoom In Zoom Out Reset image size

The average width shows an inverted behavior, however. 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: 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 site by a Gaussian function the width of the cloud in the central site is calculated to be ${\sigma }_{\text{th}}=\sqrt{{k}_{\text{B}}T/\left(m{\omega }^{2}\right)}$. 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}_{\text{TF}}=\sqrt{\mu /\left(2m{\omega }^{2}\right)}$. The calculated widths for our system are 2.1(1) μm for the normal/thermal gas and 1.7(1) μm for the condensate. The experimental results (indigo line in figure 4) are in good agreement with these predictions.

It is instructive to look at the timescale of this process by analyzing the duration of the state change in each trajectory individually. The duration of the individual state changes is taken to be the 2σ width of the fitted sigmoid function (see methods for details on the fitted function). We find that the shortest duration for a state change from the normal to the superfluid state is given by the tunneling time, τ ≈ ℏ/J (figure 5(a)). 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 [32]. For fully established phase coherence across the central site, the maximum possible coherent transport is expected [29]. 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 (figure 5(b)). 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 small tunneling couplings and thus little coherent transport into the lattice site the dissipation is strong enough that even small fluctuations lead to an emptying of the site. With increasing tunneling coupling the coherent transport into the site becomes stronger thereby on average decreasing the rate at which the atom number in the site can be decreased. This leads to an increase of the average duration of jumps. For both types of jumps, such large variations hint at the dominant role of fluctuations for the system dynamics.

Zoom In Zoom Out Reset image size

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 [33]. In our experiment, however, the temporal distribution of state changes is non-trivial (figure 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 figure 4, that a state change has preceding internal dynamics. To the best of our knowledge, such complex internal system dynamics has not been reported on single-mode systems.

Zoom In Zoom Out Reset image size

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 figure 1, the mean-field energy (chemical potential) in the neighboring sites renders the direct tunneling process into the ground state of the central site off-resonant. In a single-mode Josephson array, this would lead to macroscopic quantum self-trapping, thus inhibiting any transport [34]. 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):

$\begin{equation}{{\Gamma}}_{\text{eff}}(N)=\frac{{J}_{\text{eff}}^{2}(N)}{({\gamma }_{\text{diss}}+{\gamma }_{\text{coll}})+{J}_{\text{eff}}^{2}(N)/{\gamma }_{\text{diss}}},\end{equation} \tag{ 4 }$

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

$\begin{equation}\mathrm{d}N=\left({{\Gamma}}_{\text{eff}}(N){N}_{\text{F}}-{\gamma }_{\text{diss}}N\right)\mathrm{d}t+\sigma \,\mathrm{d}{W}_{t},\end{equation} \tag{ 5 }$

where N_F is the atom number in the neighboring sites. The term σ dW_t with $\sigma =\sqrt{{\gamma }_{\text{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 figure 7(a), 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 (figure 4), as discussed above.

Zoom In Zoom Out Reset image size

In figure 7(b) we compare the temporal distribution of the simulated jumps with our measurement and in figure 7(c) we show the comparison for the fraction of trajectories exhibiting a jump. The only fitting parameter of our model is the threshold filling. We use it to match the experimentally measured temporal distribution of state changes and the fraction of trajectories that exhibit a jump (see figures 7(b) and (c) for a direct comparison). We find a critical value of 300 atoms, which is in good agreement with the value extracted from the measurement (figure 4).

This has to be compared to the critical atom number for Bose–Einstein condensation in thermal equilibrium. 2D BEC theory [35] 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].

1. Methods

1.1. 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 $({\hat{{\Psi}}}_{0}(x,y)={\hat{{\Psi}}}_{0})$ reads

$\begin{align}\hfill \hat{H}& =\int \mathrm{d}x\,\mathrm{d}y{\hat{{\Psi}}}_{0}^{{\dagger}}\left[-\frac{{\hslash }^{2}}{2m}{\nabla }_{\perp }^{2}+\frac{1}{2}m{\omega }_{\perp }^{2}({x}^{2}+{y}^{2})\right]{\hat{{\Psi}}}_{0}\hfill \\ \hfill & \quad +\frac{{g}_{2\text{D}}}{2}\int \mathrm{d}x\,\mathrm{d}y{\hat{{\Psi}}}_{0}^{{\dagger}}{\hat{{\Psi}}}_{0}^{{\dagger}}{\hat{{\Psi}}}_{0}{\hat{{\Psi}}}_{0}\hfill \\ \hfill & \quad -J\int \mathrm{d}x\,\mathrm{d}y\left({{\Phi}}^{\ast }(x,y){\text{e}}^{\text{i}\mu t/\hslash }{\hat{{\Psi}}}_{0}+{\Phi}(x,y){\text{e}}^{-\text{i}\mu t/\hslash }{\hat{{\Psi}}}_{0}^{{\dagger}}\right).\hfill \end{align} \tag{ 6 }$

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,

$\begin{equation}\begin{aligned}\hfill \dot{\hat{\rho }}& =\hat{\mathcal{L}}\hat{\rho }\hfill \\ \hfill & =-\frac{\text{i}}{\hslash }\left[\hat{H},\hat{\rho }\right]+\frac{{\gamma }_{\text{diss}}}{2}\int \mathrm{d}x\,\mathrm{d}y\left(2{\hat{{\Psi}}}_{0}\hat{\rho }{\hat{{\Psi}}}_{0}^{{\dagger}}-{\hat{{\Psi}}}_{0}^{{\dagger}}{\hat{{\Psi}}}_{0}\hat{\rho }-\hat{\rho }{\hat{{\Psi}}}_{0}^{{\dagger}}{\hat{{\Psi}}}_{0}\right).\hfill \end{aligned}\end{equation} \tag{ 7 }$

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 reference [24].

1.2. Fitting model

We fit a sigmoid function

$\begin{equation}{y}_{\text{fit}}=A\ast \left(\pm \mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{x-{x}_{0}}{\sigma }\right)+1\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{x}{\tau }\right)+\mathrm{o}\mathrm{f}\mathrm{f}\mathrm{s}\mathrm{e}\mathrm{t}\end{equation} \tag{ 8 }$

to the measured trajectories y_meas. Here, x₀ gives the center of the error function (erf), which we take to be the time of the jump. The duration of the jump is defined as twice the width σ of the error function. The atom loss from the atomic cloud leads to a decline in signal strength over time, especially in the superfluid state. The exponential function takes this atom loss of the overall system into account. An offset is added in the model function to account for signal measured when probing an empty site.

1.3. Experimental setup

We use a BEC of ⁸⁷Rb containing about 150 × 10³ atoms. The atoms are trapped in a focused CO₂-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₀ = s E_rec in units of the recoil energy ${E}_{\text{rec}}={\pi }^{2}{\hslash }^{2}/\left(2m{d}^{2}\right)$) 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 figure 3(a) the lattice depth was varied at a constant rate from s = 4 to s = 20 and back).

1.4. 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 described by the dissipation rate γ_diss and dephasing caused by atomic collisions is included via γ_coll. The rate equation for the central site is then given by

$\begin{equation}\frac{\mathrm{d}N}{\mathrm{d}t}={{\Gamma}}_{\text{eff}}{N}_{\text{F}}-{\gamma }_{\text{diss}}N,\end{equation} \tag{ 9 }$

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

$\begin{equation}{{\Gamma}}_{\text{eff}}=\frac{{J}_{\text{eff}}^{2}}{({\gamma }_{\text{diss}}+{\gamma }_{\text{coll}})+\frac{{J}_{\text{eff}}^{2}}{{\gamma }_{\text{diss}}}}.\end{equation} \tag{ 10 }$

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 rate as follows [26]:

$\begin{equation}{{\Gamma}}_{\text{eff}}=\frac{{\left(\frac{N}{{N}_{\text{F}}}(J-\eta J)+\eta J\right)}^{2}}{\left({\gamma }_{\text{diss}}+{\gamma }_{\text{coll}}(N)\right)+\frac{{\left(\frac{N}{{N}_{\text{F}}}(J-\eta J)+\eta J\right)}^{2}}{{\gamma }_{\text{diss}}}}.\end{equation} \tag{ 11 }$

For low filling in the central site, as it is the case in the normal state, we linearize the above expression:

$\begin{equation}\frac{\mathrm{d}N}{\mathrm{d}t}=\lambda N+a\end{equation} \tag{ 12 }$

To model the fluctuations induced by the dissipation we introduce a stochastic term and obtain the stochastic differential equation

$\begin{equation}\mathrm{d}N=\left(\lambda N+a\right)\mathrm{d}t+\sigma \,\mathrm{d}{W}_{t}\end{equation} \tag{ 13 }$

with a Wiener process W_t and

$\begin{equation}\sigma =\sqrt{{\gamma }_{\text{diss}}}.\end{equation} \tag{ 14 }$

To solve equation (13) we use the Itô formula, a variant of the chain rule for stochastic multi-variable functions [36–38]. The solution to this SDE is then given by

$\begin{equation}N(t)=\left(N(0)+c\right){\text{e}}^{\left(\lambda -\frac{1}{2}{\sigma }^{2}\right)t+\sigma \sqrt{t}x}-c\end{equation} \tag{ 15 }$

with a standard normal variable x and c = a/λ.

Acknowledgments

Data availability statement