Interacting bosons in two-dimensional lattices with localized dissipation

Although dissipation in quantum systems leads to decoherence of quantum states, recent years have witnessed an upsurge in allowing dissipation on purpose to study non-equilibrium dynamics in various physical systems such as optical cavities [1], trapped ions [2, 3], exciton–polariton BECs [4], and microcavity arrays coupled with superconducting qubits [5]. This is in part, because dissipation can be used as an efficient tool for preparing and manipulating quantum states [6], and in part because the interplay between unitary and dissipative dynamics leads to the emergence of non-equilibrium steady states [7]. Among various experimental platforms for studying dissipative dynamics, the most promising candidate turns out to be cold atoms due to its high degree of experimental controllability. This has led to several cold atom experiments where single [8–10] and two-body particle losses have been investigated with widely controllable dissipation strength, revealing the melting of quantum phases across the superfluid-Mott transition [11].

In this work, we consider a two-dimensional (2D) Bose–Hubbard (BH) model and examine the effect of dissipation on the quantum phase transition from a Mott insulator (MI) to a superfluid phase (SF) regime. In particular, we focus on the different parameter regimes of the interaction for a fixed chemical potential. It is well established that the superfluid phase can itself be classified into two regimes depending on the relative strength of the two momentum scales, namely the Ginzburg momentum ${k}_{G}\approx \sqrt{n{(U/J)}^{3}}$ and healing momentum ${k}_{h}\approx \sqrt{n(U/J)}$ with $n,U$ and J as the Bogoliubov approximated superfluid density, interaction and tunneling strength respectively. For k_G/k_h ∼ 1, a weakly-correlated superfluid (WCSF) regime is obtained, where superfluid stiffness or condensate density reaches to the maximum Bogoliubov-approximated values. In contrast, for k_G/k_h ≪ 1, a strongly-correlated superfluid regime (SCSF) emerges, where no Bogoliubov theory exists and the condensate density and superfluid stiffness are found to be highly suppressed [12]. Furthermore, the SCSF regime has particle or hole character of phonons, and a gapped mode exhibit therein as compared to WCSF [13, 14]. In view of that, we aim to address the following issues: Is the dynamics different in these two regimes of superfluid states? Can we probe these two regimes by studying dissipative dynamics? We note that although dissipative dynamics has been studied both in few-sites and extended BH model [15–18], the effect of dissipation on the different superfluid regimes is yet to be explored in detail.

To study dissipative dynamics of the BH model, we use the Gutzwiller mean-field (GMFT) theory for density matrices. Previous works have shown that the GMFT can validate correct dynamics of the BH model in the clean [19, 20] and disordered limit [21–23]. Moreover, the study of nonequilibrium phenomena such as expansion dynamics of Mott phase [24], dynamical generation of molecular condensates [25], dipole oscillations [26], dissipative Rydberg atoms [27], etc can well be described in the framework of Gutzwiller approximations. Also, it has been shown that GMFT can reasonably describe particle loss in the quantum Zeno regime of bosonic systems with long-range interaction [17]. Additionally, more recent experimental results of Mott-superfluid crossover in the presence of localized dissipation has been explained well within the GMFT [11].

In what follows, we use the GMFT for density matrices and Lindblad master equation to study the particle loss in the different correlated superfluid regimes. In the WCSF regime, we find that, for weak dissipation, the loss rate is proportional to the strength of the dissipation. However, for stronger dissipation, the systems protects itself from being dissipated, hence the loss rate decreases, corroborating the commonly known Zeno effect in the superfluid phase when the onsite repulsive interactions are weak [49, 50]. The non-monotonic nature of the particle loss is in stark contrast to the superfluid regime when interaction is relatively strong. In particular, we see Zeno-anti-Zeno crossover as a function of dissipation strength. We furthermore show that, in the SCSF, the long-time dynamics differs significantly from the WCSF. These facts also translate into the fidelity obtained in these two regimes. Thus the dissipative dynamics in the 2D BH model may indicate the effect of interaction in the coherent many-body phases of this system.

The rest of the paper is organized as follows. In section 2, we discuss the model and formalism in the presence of local dissipation. In particular, we demonstrate single-site particle loss using density matrix formalism. This is followed by section 3, where we discuss particle loss in the different regimes of the onsite interaction of the model Hamiltonian. We also present the long-time dynamics in those regimes. We furthermore show the evolution of fidelity for different strengths of the dissipation. Finally, we conclude with a discussion on the possible future directions in section 4.

We consider a homogeneous bosonic gas at zero temperature trapped in a 2D optical lattice potential. Employing the tight-binding approximation, the statics and dynamics of the bosonic atoms can be well described by the single lowest band BH Hamiltonian with nearest-neighbor (NN) hopping between the lattice sites, and local onsite repulsive interactions (U) between the atoms:

$$\begin{eqnarray}&&\hat{H}=-J\displaystyle \sum _{ \langle {ll}^{{\rm{{\prime} }}} \rangle }{\hat{a}}_{l}^{\dagger }{\hat{a}}_{{l}^{{\rm{{\prime} }}}}+\displaystyle \frac{U}{2}\displaystyle \sum _{l}{\hat{n}}_{l}({\hat{n}}_{l}-1)-\displaystyle \sum _{l}\mu {\hat{n}}_{l},\end{eqnarray} \tag{ 1 }$$

where $\langle {ll}^{\prime} \rangle$ refers to NNs l and $l^{\prime}$, J is the hopping strength between two NN sites, and μ is the chemical potential; ${\hat{a}}_{l}^{\dagger }({\hat{a}}_{l})$ are the bosonic creation (annihilation) operators with ${\hat{n}}_{l}={\hat{a}}_{l}^{\dagger }{\hat{a}}_{l}$ as the occupation number. Depending on the relative values of J/U, this model supports two distinct phases. For J/U ≪ 1 , the system exhibits insulating phase, namely Mott-insulating (MI) phase with commensurate integer fillings and vanishing order parameter. In contrast, J/U ≫ 1 leads to superfluid (SF) phase with non-vanishing compressibility [28–31]. As mentioned before, the superfluid medium can be further divided into two regimes. Near the MI–SF boundary (at the onset of superfluidity), the superfluid states are strongly correlated with highly suppressed superfluid stiffness, whereas J/U ≫ 1 corresponds to WCSF with stiffness close to unity. Note that such correlated phases can be distinguished by the non-perturbative renormalization group analysis as shown in [12]. Although distinction between these phases are beyond the scope of Gutzwiller mean-field theory, however, it can estimate the superfluid stiffness [20] which qualitatively agrees with the stiffness presented in [12].

To obtain the ground state of BH model, we employ the mean-field approximation by decomposing ${\hat{a}}_{l}={\phi }_{l}+\delta {\hat{a}}_{l}$, where the order-parameter ${\phi }_{l}=\langle {\hat{a}}_{l}\rangle$. With this, equation (1) can further be simplified as a sum of single-site mean-field Hamiltonian $\hat{H}={\sum }_{l,l^{\prime} }{\hat{h}}_{l,l^{\prime} }^{\mathrm{MF}}$ with

$$\begin{eqnarray}\begin{array}{rcl}{\hat{h}}_{l,l^{\prime} }^{\mathrm{MF}} & = & -J\left[({\phi }_{l+1,l^{\prime} }^{* }{\hat{b}}_{l,l^{\prime} }-{\phi }_{l+1,l^{\prime} }^{* }{\phi }_{l,l^{\prime} }+{\phi }_{l,l^{\prime} -1}^{* }{\hat{b}}_{l,l^{\prime} }\right.\\ & & -\left.\,{\phi }_{l,l^{\prime} -1}^{* }{\phi }_{l,l^{\prime} })+{\rm{h}}.{\rm{c}}.\right]+\frac{U}{2}{\hat{n}}_{l,l^{\prime} }({\hat{n}}_{l,l^{\prime} }-1)-\mu {\hat{n}}_{l,l^{\prime} },\end{array}\end{eqnarray} \tag{ 2 }$$

where l and $l^{\prime}$ are the lattice site indices along the x and y directions, respectively. In the Fock basis of occupation numbers, equation (2) can eventually be diagonalised independently using the Gutzwiller variational ansatz

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{GW}}\rangle =\displaystyle \prod _{\otimes l,l^{\prime} }| \psi {\rangle }_{l,l^{\prime} }=\displaystyle \prod _{\otimes l,l^{\prime} }\displaystyle \sum _{n=0}^{{N}_{{\rm{b}}}}{c}_{n}^{(l,l^{\prime} )}| n{\rangle }_{l,l^{\prime} },\end{eqnarray} \tag{ 3 }$$

where $| \psi {\rangle }_{l,l^{\prime} }$ is the ground state of the lattice site with index $(l,l^{\prime} )$, $| n{\rangle }_{l,l^{\prime} }$ is the corresponding occupation number basis, ${c}_{n}^{(l,l^{\prime} )}s$ are the associated complex coefficients of the expansion, satisfying ${\sum }_{n}| {c}_{n}^{(l,l^{\prime} )}{| }^{2}=1$, and N_b is the total number of basis states. Using these definitions and considerations, the order parameter at the $(l,l^{\prime} )$th lattice site is given by

$$\begin{eqnarray}&&{\phi }_{(l,l^{\prime} )}=\langle {{\rm{\Psi }}}_{\mathrm{GW}}| {\hat{b}}_{l,l^{\prime} }| {{\rm{\Psi }}}_{\mathrm{GW}}\rangle =\displaystyle \sum _{n=1}^{{N}_{{\rm{b}}}}\sqrt{n}{c}_{n-1}^{* (l,l^{\prime} )}{c}_{n}^{(l,l^{\prime} )}.\end{eqnarray} \tag{ 4 }$$

As mentioned earlier, ${\phi }_{(l,l^{\prime} )}=0$ in the MI phase, whereas it is finite in the SF phase. Note that, in our GMFT analysis, the critical values of hopping ${J}_{c}/U$, at which the SF order parameter starts to develop, is obtained to be 0.04 at μ/U = 0.5. Thus, within the framework of GMFT, we consider J/U ∼ 0.04–0.05 as strongly-correlated SF regime, whereas J/U ≳ 0.1 as weakly-correlated SF regime for μ/U = 0.5. It is worth mentioning that in addition to equilibrium GMFT, the time-dependent GMFT has been successful in capturing the essential physics related to amplitude modes on interacting bosons in a lattice model [32].

2.1. One-body loss

Having obtained the many-body ground states for the BH model through GMFT, we wish to look at the effects of dissipation on these phases. We first introduce a local loss term that acts on a single site of the lattice system to remove particles. The loss can either be one-body, two-body or higher orders. For the present work, we have considered only one-body loss. Invoking Born-Markov approximation and rotating-wave approximation [33], the lossy dynamics of the BH model can thus be represented through the following Lindblad equation with local dissipation rate γ.

$$\begin{eqnarray}&&{\partial }_{t}\hat{\rho }=-{\rm{i}}[\hat{H},\hat{\rho }]+\gamma \left[{\hat{a}}_{j}\hat{\rho }{\hat{a}}_{j}^{\dagger }-\displaystyle \frac{1}{2}\{{\hat{a}}_{j}^{\dagger }{\hat{a}}_{j},\hat{\rho }\}\right].\end{eqnarray} \tag{ 5 }$$

The first term on the right represents unitary dynamics of the BH model, whereas the second term describes the non-unitary evolution. Experimentally, the dissipation rate can be tuned by varying the intensity of the applied electron beam. The local atom-electron collisions leave the atoms ionized and by the use of an electrostatic field the atoms escape from the lattice. Thus a local dissipation is generated with the aid of anelectron beam [34]. Applying the Gutzwiller approximation to the density operator $\rho ={\prod }_{\otimes l,l^{\prime} }{\sum }_{n,m}{c}_{{nm}}^{(l,l^{\prime} )}| n{\rangle }_{(l,l^{\prime} )}\langle m{| }_{(l,l^{\prime} )}$, the evolution of the matrix elements of the density operator corresponding to a lattice site is obtained to be

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{c}_{{nm}}^{(l,l^{\prime} )} & = & {\rm{i}}\left\{\Space{0ex}{3.0ex}{0ex}J(\sqrt{n}{c}_{n-1,m}^{(l,l^{\prime} )}-\sqrt{m+1}{c}_{n,m+1}^{(l,l^{\prime} )})({\phi }_{(l+1,l^{\prime} )}\right.\\ & & +\,{\phi }_{(l-1,l^{\prime} )}+{\phi }_{(l,l^{\prime} +1)}+{\phi }_{(l,l^{\prime} -1)})+J(\sqrt{n+1}{c}_{n+1,m}^{(l,l^{\prime} )}\\ & & -\,\sqrt{m}{c}_{n,m-1}^{l,l^{\prime} })({\phi }_{(l+1,l^{\prime} )}^{* }+{\phi }_{(l-1,l^{\prime} )}^{* }+{\phi }_{(l,l^{\prime} +1)}^{* }\\ & & +\,{\phi }_{(l,l^{\prime} -1)}^{* })+\left(\mu (n-m)-\displaystyle \frac{U}{2}[n(n-1)-m(m-1)]\right.\\ & & +\,\left.\left.{\rm{i}}\displaystyle \frac{\gamma }{2}(n+m)\right){c}_{n,m}^{(l,l^{\prime} )}\right\}+\gamma \sqrt{(n+1)(m+1)}{c}_{n+1,m+1}^{(l,l^{\prime} )}.\end{array}\end{eqnarray} \tag{ 6 }$$

Clearly equation (6) is a coupled equation of multiple variables which entails analytical solutions at some special values or cases which will be evident shortly. Equation (6) also involves several parameters of the Hamiltonian. Thus, the dynamics is expected to depend on the interplay between those parameters, which may lead to interesting dissipative dynamics.

To solve equation (6) numerically, we use the fourth-order Runge-Kutta method with open boundary condition with the order parameter ϕ of the BH Hamiltonian being obtained by imaginary time propagation. Using the coefficients of expansion c_nm, one can obtain the total number of lost particles

$$\begin{eqnarray}&&{\rm{\Delta }}N=N(t=0)-N(t)=\langle n\rangle (t=0)-\langle n\rangle (t),\end{eqnarray} \tag{ 7 }$$

where $N(t)={\sum }_{l,l^{\prime} }{\sum }_{n}{c}_{{nn}}^{(l,l^{\prime} )}(t)| n{| }^{2}$ is total number of particles in the system.

We finally compute the time evolution of fidelity [35] $F=\mathrm{Tr}(\hat{\rho }(t)\hat{\rho }(0))$ as a measure of probability to recover the initial state at a given time, t. For the present work, we evaluate F(t) for different strengths of the dissipation γ and discuss its non-monotonic dependence on γ as a manifestation of unconventional particle loss in the correlated superfluid regime.

We consider a 2d lattice potential in which the central site is subjected to dissipation mediated through one-body loss as shown in figure 1. Theoretically, the dynamics of the composite system is governed by the Lindblad equation in the Gutzwiller approximation following equation (6). On the experimental front, such a setup is feasible in which the central site is coupled to a localized bath and the resulting dynamics is monitored [34]. For the current study, we choose a finite-sized 41 × 41 lattice system and N_b = 10. We checked that this system size is sufficient to reveal correct dynamics in the parameter regime considered here. We then vary J to obtain the initial homogeneous ground state in the three different regimes, namely ground states of (a) the deep SF, (b) the strongly-correlated SF, and (c) the Mott insulating regime.

Zoom In Zoom Out Reset image size

3.1. Particle loss

The coupling between the system and bath is turned on at t = 0, which drives away the atoms from the central site, and a defect is thus formed. The central site tends to lose particles, and eventually gets refilled by the nearest neighbor tunneling processes. For t > 0, the dynamics includes hopping processes propagating through the system into the lossy site which gives rise to an effective coupling between the NN sites. This leads to the whole system getting affected as a result of localized dissipation. The difference in the rate of loss of particles and the rate of refilling leads to a variation in the dissipative dynamics which forms the motivation of our present investigation. To quantify this behavior, we compute the total number of particles lost from the system, and demonstrate that the intermediate and long-time behavior for the three different phases is indeed different upon the application of a localized dissipation. Furthermore, we also compute the fidelity to substantiate our findings.

Mott phase—When $J/U\ll 1$, that is in the MI regime with $\langle n\rangle (t=0)=1$, the loss rate monotonically increases with γ, which occurs due to randomization of phases in this regime [36, 37]. The corresponding plot is shown in figures 2(a), (d). On varying γ, and choosing an intermediate time tU = 20, we find ΔN to be increasing and gets saturated. This can be understood from equation (6). Since the superfluid order parameter ϕ remains zero in the MI phase for all values of γ, equation (6) for n = m leads to

$$\begin{eqnarray}&&{c}_{{nn}}^{(l,l^{\prime} )}(t)={c}_{n,n}^{(l,l^{\prime} )}(0){{\rm{e}}}^{-\gamma t},\end{eqnarray} \tag{ 8 }$$

where we have used ${c}_{n+1,n+1}=0$ due to the Mott nature of the wavefunction. This gives ${\rm{\Delta }}N=1-{{\rm{e}}}^{-2\gamma t}$ and validates the result of figures 2 (a), (d). We note that results obtained using GMFT have reasonably good agreement with the analytical results discussed for boson in double well trap with large filling fraction as discussed in [15]. However, the method discussed here applies to any filling and all range of interaction and valid for any system size.

Zoom In Zoom Out Reset image size

Deep SF phase—In the deep SF regime, we can ignore the term involving U in equation (6). Then for fixed μ, the competition between J and γ leads to non-monotonic particle loss as a function γ. Figures 2(c), (f) reflect this behavior, where the loss rate increases linearly with γ for small γ. This can be qualitatively understood from equation (6) and using some perturbative argument. We first set γ = 0 to find c_nm(t). This corresponds to equilibrium dynamics of a closed BH model. In this limit, c_nn(t) turns out to be constant as also evident from equation (5). This is essentially because the order parameter ϕ is real, and ρ is symmetric. Albeit, $\gamma \ne 0$ does not render c_nn(t) to be constant. Using this solution as an initial solution for finite but very small γ in equation (6), we obtain ${c}_{{nn}}(t)\sim {\alpha }_{0}+\gamma t$, where α₀ is a constant function of J. Note that, the order parameter ϕ of the central site and its nearest sites have been considered to be constant for some accessible time domain as seen in figure 3. Note also, the approximate solution of c_nn(t) is valid only up to some intermediate time. Thus for fixed time we obtain ΔN ∼ γ for small γ and it qualitatively agrees with the perturbative analysis of one-dimensional hard-core bosons in optical lattices [38].

Zoom In Zoom Out Reset image size

On the other hand, strong dissipation γ ≫ J leads to quantum Zeno effect which describes the suppression of dynamics due to the application of a continuous projective measurement. In this limit, the system ceases to get affected by attaining a dark state [38–40]. Consequently, the loss of atoms will be suppressed as evident from figures 2(c), (f). To understand this, we neglect terms involving J since γ ≫ J. Then c_nn(t) can be approximated as ${c}_{{nn}}(t)\sim {{\rm{e}}}^{-\gamma t}$, considering minimal contribution from nearest sites. Note that for stronger dissipation, the order parameter reaches to steady states much faster than weak dissipation. However, the Zeno behavior is reflected only in the order parameter of the nearest sites. Figure 3(b) clearly shows the non-monotonic behavior of $| {\phi }_{l,l+1}(t){| }^{2}$ for strong dissipation.

Close to SF-Mott boundary—A dramatic change in the dynamical behavior occurs for the phase close to the SF–MI transition, that is when J/U = 0.05. This phase also referred to as the strongly-correlated superfluid phase where the condensate density is suppressed, exhibiting a different behavior upon the application of a localized dissipation. It is evident from figures 2(b), (e) that, for intermediate time, ΔN initially increases with the increase of γ as before. However, for moderate range of γ, we see decaying and increasing trend of ΔN. This is in contrast to the particle loss both in the MI and deep SF phase. In addition, for stronger dissipation, ΔN decays notably slow as compared to the deep SF regime. In a nutshell, we find Zeno-anti-Zeno behavior at the onset of superfluidity (near the MI–SF phase boundary). Indeed, here the dynamics is governed by the combination of all parameters J, U, γ. As a result, the evolution of order parameters in the central site and its nearest site differs significantly than their dynamics deep into the superfluid (see figures 3(c), (d)), hence the particle loss. We would like to point out that the resulting dynamics emanates from the interplay of J and U, but there has been no analytical or numerical evidence related to the actual trend in particle loss in this regime even for few-wells bosonic systems. We also point out that the particle loss in the deep superfluid regime is itself non-monotonic, but there are further modifications to this non-monotonicity in this domain.

To this end, we comment on the dynamics away from μ = 0.05. For fixed J/U = 0.05, the Zeno-anti-Zeno regime is retained with slight variation of μ. However, such regime turns out to be very narrow and asymmetric with respect to μ. Indeed, this is attributed to the asymmetric SF–MI phase boundary.

3.2. Long-time dynamics

In this section, we discuss the long-time dynamics of particle loss for different strengths of dissipation γ. In the WCSF regime, the long-time dynamics turns out to be similar to the dynamics of intermediate time. This is apparent in figure 4 where the particle loss initially increases with γ and acquires some maximum value at some particular γ, and finally it decreases as a consequence of Zeno effect. This resembles to the particle loss at tU = 80 as seen from figure 2(d). In contrast, the Zeno-anti Zeno kind of behavior is not prominent in the long-time dynamics of the strongly-correlated regime (see figure 2(e)). In fact, the loss is almost same for all range of γ. The reason for such behavior can again be attributed to the randomization of phases close to the phase boundary similar to the Mott phase. As the oscillation amplitudes of the order parameter increase in the long-time limit (see Figure 3(c), (d)), we can approximate the average values of the order parameter to be zero. This leads to ${\rm{\Delta }}N\sim { \mathcal O }(1)$ similar to the Mott phase. Thus this distinction of particle loss in the long-time limit is in conjunction with our earlier observations that the dissipative dynamics in bosonic systems strongly depends on the local interaction.

Zoom In Zoom Out Reset image size

3.3. Fidelity

To substantiate the results obtained in the preceding sections, we discuss the evolution of fidelity both in the WCSF and SCSF regimes. Figures 5(a), (b) illustrate the fidelity in both the superfluid regimes for various γ. Clearly, F(t) decays exponentially as a function of t and finally saturates after an accessible time, corroborating the non-equilibrium steady state behavior of the system[41]. However, the degree of decay depends on the strength of the dissipation. For stronger dissipation, F(t) saturates even before relatively weaker dissipation in the SCSF regime (see figure 5(b)). This is in contrast to the case of WCSF where no such behavior is observed (see figure 5(a)). This is again attributed to the different dynamical behavior of the particle loss in different regime of the SF–MI transition. Since it is experimentally feasible to measure fidelity in open quantum dynamics, it can be used as a probe to identify different interacting regimes of the BH model [42–46].

Zoom In Zoom Out Reset image size

In conclusion, we present the dynamics of a 2D BH model in the presence of a local dissipation. We find that the particle loss in this setting strongly depends on the on-site bosonic interactions. While the particle loss in the Mott phase has monotonic behavior as a function of dissipation, it significantly differs in the superfluid regimes with non-monotonic dependence on the strength of the dissipation. In addition, we find notably different dissipative dynamics in the strongly and weakly correlated superfluid regimes. Such differences are reflected in the fidelity and long-time dynamics.

Although the results reported here are explained qualitatively within the framework of Gutzwiller mean-field theory for density matrices and Lindblad master equation, it would be desirable in future to verify it using more sophisticated tools like truncated Wigner approximations [47], quantum trajectories method [48] or by constructing density matrix using exact diagonalization techniques, alleviating the limitations of GMFT. A natural extension of the study presented here would be to explore dissipative dynamics in disordered systems, and in particular, the Bose glass phase.

We thank J Marino, A Lazarides and S Diehl for useful discussions. AR thanks R Bai, S Bandopadhyay, K Suthar, S Pal, and D Angom for insightful discussions on GMFT.