Observation of the Mott Insulator to Superfluid Crossover of a Driven-Dissipative Bose-Hubbard System

TL;DR

This paper demonstrates how engineered dissipation can control quantum phase transitions in a Bose-Hubbard system, revealing dissipation's role in delaying the Mott insulator to superfluid crossover in cold atom experiments.

Contribution

It introduces a method to manipulate quantum states using controllable dissipation in a Bose-Hubbard system, a novel approach in quantum many-body physics.

Findings

Strong on-site dissipation delays Mott insulator melting

Dissipation suppresses phase coherence growth

Quenching dissipation enables non-equilibrium studies

Abstract

Dissipation is ubiquitous in nature and plays a crucial role in quantum systems such as causing decoherence of quantum states. Recently, much attention has been paid to an intriguing possibility of dissipation as an efficient tool for preparation and manipulation of quantum states. Here we report the realization of successful demonstration of a novel role of dissipation in a quantum phase transition using cold atoms. We realize an engineered dissipative Bose-Hubbard system by introducing a controllable strength of two-body inelastic collision via photo-association for ultracold bosons in a three-dimensional optical lattice. In the dynamics subjected to a slow ramp-down of the optical lattice, we find that strong on-site dissipation favors the Mott insulating state: the melting of the Mott insulator is delayed and the growth of the phase coherence is suppressed. The controllability of the dissipation is highlighted by quenching the dissipation, providing a novel method for investigating a quantum many-body state and its non-equilibrium dynamics.

Equations

$$\hslash \frac{d}{dt}{\hat{\rho }}_{\mathrm{eff}}=-i\left[{\hat{H}}_{\mathrm{eff}},{\hat{\rho }}_{\mathrm{eff}}\right]+L_2({\hat{\rho }}_{\mathrm{eff}}),$$

$${\hat{H}}_{\mathrm{eff}}=\sum _j\frac{U}{2}{\hat{n}}_{\mathrm{A},j}({\hat{n}}_{\mathrm{A},j}-1)-\sum _{⟨j,k⟩}J\left({\hat{a}}_j^{†}{\hat{a}}_k+h.c.\right),$$

$$L_2({\hat{\rho }}_{\mathrm{eff}})=\frac{\hslash {\Gamma }_{\mathrm{PA}}}{4}\sum _j\left(-{\hat{a}}_j^{†}{\hat{a}}_j^{†}{\hat{a}}_j{\hat{a}}_j{\hat{\rho }}_{\mathrm{eff}}-{\hat{\rho }}_{\mathrm{eff}}{\hat{a}}_j^{†}{\hat{a}}_j^{†}{\hat{a}}_j{\hat{a}}_j+2{\hat{a}}_j{\hat{a}}_j{\hat{\rho }}_{\mathrm{eff}}{\hat{a}}_j^{†}{\hat{a}}_j^{†}\right),$$

$${\Gamma }_{\mathrm{PA}}=8\frac{g^2}{{\hslash }^2{\Gamma }_{\mathrm{M}}}.$$

$$\hslash \frac{d}{dt}\hat{\rho }=-i\left[\hat{H},\hat{\rho }\right]+L_{\mathrm{M}}(\hat{\rho }).$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\mathrm{am}}+{\hat{H}}_{\mathrm{hop}},$$

$${\hat{H}}_0$$

$${\hat{H}}_{\mathrm{am}}$$

$${\hat{H}}_{\mathrm{hop}}$$

$$L_{\mathrm{M}}(\hat{\rho })=\frac{\hslash {\Gamma }_{\mathrm{M}}}{2}\sum _j\left(-{\hat{n}}_{\mathrm{M},j}\hat{\rho }-\hat{\rho }{\hat{n}}_{\mathrm{M},j}+2{\hat{m}}_j\hat{\rho }{\hat{m}}_j^{†}\right).$$

$$\hslash \frac{d}{dt}{\hat{\rho }}_{\mathrm{eff}}=-i\left[{\hat{H}}_{\mathrm{eff}},{\hat{\rho }}_{\mathrm{eff}}\right]+L_2({\hat{\rho }}_{\mathrm{eff}}),$$

$${\hat{H}}_{\mathrm{eff}}=\sum _j\frac{U}{2}{\hat{n}}_{\mathrm{A},j}({\hat{n}}_{\mathrm{A},j}-1)-\sum _{⟨j,k⟩}J\left({\hat{a}}_j^{†}{\hat{a}}_k+h.c.\right),$$

$$L_2({\hat{\rho }}_{\mathrm{eff}})=\frac{\hslash {\Gamma }_{\mathrm{PA}}}{4}\sum _j\left(-{\hat{a}}_j^{†}{\hat{a}}_j^{†}{\hat{a}}_j{\hat{a}}_j{\hat{\rho }}_{\mathrm{eff}}-{\hat{\rho }}_{\mathrm{eff}}{\hat{a}}_j^{†}{\hat{a}}_j^{†}{\hat{a}}_j{\hat{a}}_j+2{\hat{a}}_j{\hat{a}}_j{\hat{\rho }}_{\mathrm{eff}}{\hat{a}}_j^{†}{\hat{a}}_j^{†}\right),$$

$${\Gamma }_{\mathrm{PA}}=8\frac{g^2}{{\hslash }^2{\Gamma }_{\mathrm{M}}}.$$

$$|n_a,n_m{⟩}_j=\frac{1}{\sqrt{n_a!n_m!}}({\hat{a}}_j^{†}{)}^{n_a}({\hat{m}}_j^{†}{)}^{n_m}|vac{⟩}_j$$

$$\displaystyle|l\rangle_{j}=\left\{\begin{array}[]{cc}|n_{a}=l-1,n_{m}=0\rangle_{j},&{\rm when}\,\,1\leq l\leq d,\\ |n_{a}=l-1-d,n_{m}=1\rangle_{j},&{\rm when}\,\,d+1\leq l\leq 2d.\end{array}\right.$$

$$\hat{\rho }=\sum _{l_1,l_2,\dots }\sum _{m_1,m_2,\dots }{\rho }_{l_1,l_2,\dots }^{m_1,m_2,\dots }\prod _j|l_j⟩⟨m_j{|}_j.$$

$$\hslash \frac{d}{dt}\mathbf{\rho }=\left({\hat{\mathcal{M}}}_0+{\hat{\mathcal{V}}}_{\mathrm{am}}+{\hat{\mathcal{V}}}_{\mathrm{hop}}\right)\mathbf{\rho },$$

$$|l_j,m_j{)}_j=|l_j⟩⟨m_j{|}_j.$$

$$(l_j^{\prime },m_j^{\prime }|l_j,m_j)={\delta }_{l_j^{\prime },l_j}{\delta }_{m_j^{\prime },m_j}.$$

$${\hat{\mathcal{M}}}_0=\sum _j{\hat{\mathcal{M}}}_{0,j}^{\mathrm{loc}},$$

$${\hat{\mathcal{M}}}_0^{\mathrm{loc}}=\sum _{l=1}^{2d}\sum _{m=1}^{2d}|l,m)(l,m|i\left({\left(E_l^{(0)}\right)}^{\ast }-E_m^{(0)}\right)+\sum _{l=1}^d\sum _{m=1}^d|l,m)(d+l,d+m|\hslash {\Gamma }_{\mathrm{M}},$$

$$\displaystyle E_{l}^{(0)}=\left\{\begin{array}[]{cc}\frac{U}{2}(l-1)(l-2),&{\rm when}\,\,1\leq l\leq d,\\ D+\frac{U}{2}(l-d-1)(l-d-2)+W(l-d-1)-i\frac{\hbar\Gamma_{\rm M}}{2},&{\rm when}\,\,d+1\leq l\leq 2d.\end{array}\right.$$

$${\hat{\mathcal{V}}}_{\mathrm{am}}=\sum _j{\hat{\mathcal{V}}}_{\mathrm{am},j}^{\mathrm{loc}}$$

$${\hat{\mathcal{V}}}_{\mathrm{am}}^{\mathrm{loc}}$$

$${\tilde{g}}_l=g\sqrt{(l-1)(l-2)}.$$

$${\hat{\mathcal{M}}}_0^{\mathrm{loc}}|v_{\alpha })$$

$$(w_{\alpha }|{\hat{\mathcal{M}}}_0^{\mathrm{loc}}$$

$$(w_{{\alpha }^{\prime }}|v_{\alpha })={\delta }_{\alpha ,{\alpha }^{\prime }}.$$

$${\hat{\mathcal{P}}}_{\alpha }^{\mathrm{loc}}=|v_{\alpha })(w_{\alpha }|.$$

$${\hat{\mathcal{P}}}_{\mathrm{eff}}^{\mathrm{loc}}$$

$${\hat{\mathcal{P}}}_{\mathrm{eff}}=\prod _{j=1}^M{\hat{\mathcal{P}}}_{\mathrm{eff},j}^{\mathrm{loc}},$$

Full text

I Introduction

Dissipation -coupling to the environment- plays an essential role in quantum systems. On the one hand it causes decoherence of quantum states thus it limits the coherent dynamics. Therefore, protection of the quantum states from the coupling to the environment has been a crucial issue in quantum engineering. On the other hand, the dissipation can be used as an efficient tool for preparation and manipulation of particular quantum states of interest Daley-14 ; Muller-12 . Understanding and controlling non-equilibrium dynamics of correlated quantum many-body systems with dissipation are indeed an imperative issue shared in common among experimental systems in diverse areas of physics, including ultracold gases Daley-14 ; Muller-12 ; Sieberer-16 , Bose-Einstein condensates placed in optical cavities Sieberer-16 ; Ritsch-13 , trapped ions Blatt-12 ; Bohnet-16 , exciton-polariton BEC Sieberer-16 ; Carusotto-13 , microcavity arrays coupled with superconducting qubits Houck-12 ; Fitzpatrick-17 .

Cold atoms, which attract much attention for investigation of many-body quantum systems owing to high controllability of various parameters, are often regarded as an ideal closed (or isolated) quantum system. However, this controllability allows also for creating open quantum systems by introducing dissipation processes. So far, various kinds of theoretical works on the effect of dissipation have predicted novel quantum states engineered by dissipation due to photon scattering and particle loss Diehl-08 ; Witthaut-08 ; Verstraete-09 ; Han-09 ; DallaTorre-10 ; Diehl-10 ; Tomadin-11 ; LeBoite-13 ; Bidanovic-14 ; Stannigel-14 ; Ashida-16 . Experimentally, a one-body dissipation has been introduced in a controlled manner with several methods. The utility of an electron beam has been demonstrated in Ref. Barontini-13 ; Labouvie-15 ; Labouvie-16 . With a well-designed photon scattering process, measurement backaction on the many-body state Patil-15 and the many-body localization in open quantum systems Luschen-17 have been investigated. In the case of three-body loss process, controlling the strength of three-body recombination by Feshbach resonance and realization of a novel metastable many-body state have been demonstrated Mark-12 .

Because the two-body interaction is fundamental and crucial for the emergence of the novel quantum states and many-body physics such as quantum phase transitions, it is important to investigate the effect of two-body dissipation on quantum many-body systems in the way that the strength of dissipation can be widely controlled. While the pioneering works were reported in which the two-body loss process was realized by using intrinsic nature of molecules, such as vibrational quenching Syassen-08 and chemical reaction Yan-13 , and the lifetime of the molecules was investigated, a systematic study of the effect of two-body dissipation on quantum many-body physics has not been reported.

Here we report an investigation of a Bose-Hubbard system using ultracold atoms in a three-dimensional (3D) optical lattice, in which we introduce engineered dissipation of the two-body particle losses. By exploiting the highly controllable nature of the dissipation that we introduce, we successfully reveal the effect of the dissipation on the quantum phase transition from a Mott insulator to a superfluid state in a systematic manner. In particular, we observe in the ramp-down dynamics across the crossover from the Mott insulator to the superfluid states that the melting of the Mott state is delayed and the growth of the phase coherence is suppressed for the strong dissipation. Note that the type of the dissipation introduced in this work is an on-site one. The highly controllable on-site dissipation allows us to study the quench dynamics as the novel method of the initial state preparation, providing a new way for investigating non-equilibrium quantum dynamics. The success in engineering the controllable dissipation of the Bose-Hubbard system offers new opportunities for exploring novel roles of the dissipation in quantum many-body systems.

II Results

Engineered two-body dissipation

A crucial part of the research of a driven-dissipative quantum many-body state is the design of engineered dissipation for a quantum many-body system. In our experiment, two-body inelastic atom loss with controllable strength is successfully implemented by introducing a single-photon photo-association (PA) process for ultracold ytterbium (174Yb) atoms in a 3D optical lattice (see Materials and Methods). The PA beam drives the intercombination transition of ${}^{1}S_{0}\leftrightarrow{}^{3}P_{1}$, by which two atoms in the doubly-occupied sites are photoassociated into the ${}^{1}S_{0}+{}^{3}P_{1}$ molecular state and immediately dissociated into the two ground state atoms (see Fig. 1 (A) and Materials and Methods). This process gives high kinetic energy to the dissociated atoms, which thus results in the escape from the lattice. In this way, the PA laser induces the two-body inelastic collision loss between the two atoms occupying the same site. We can realize controllable strength of inelastic collision coefficient $\beta_{\rm{PA}}$ by changing the intensity of the PA laser $I$, as shown in Fig. 1 (B). $\beta_{\rm{PA}}$ is determined through the relation $\Gamma_{\rm{PA}}=\beta_{\rm{PA}}\int|w({\bm{r}})|^{4}d{\bm{r}}$, where the inelastic collision rate $\Gamma_{\rm{PA}}$ is measured from the exponential fit of collisional loss dynamics as shown in the inset of Fig. 1 (B), and $w({\bm{r}})$ is the Wannier function of the lowest band. Note that the measurement is done in a deep lattice potential of $V_{0}=14~{}E_{R}$, so that atom tunneling is suppressed in the timescale of this measurement. Here $E_{R}=h^{2}/(2m\lambda_{L}^{2})$ is a recoil energy with $m$ the mass of the 174Yb atom, $h$ the Planck’s constant, and $\lambda_{L}=532$ nm the wavelength of the lattice laser. In this way, we can realize controllable strength of inelastic collision up to $\beta_{\rm{PA}}\sim 1.2\times 10^{-10}~{}{\rm cm}^{3}/{\rm s}$ corresponding to $\Gamma_{\rm{PA}}\sim$70 kHz in the lattice depth of $V_{0}=14~{}E_{R}$ and $\gamma\sim 5$, where $\gamma=\hbar\Gamma_{\rm PA}/U$ is the dimensionless dissipation strength which is independent of the lattice depth and $U$ denotes the on-site interaction.

Model

It is important to understand how this PA process is effectively described in a dissipative Bose-Hubbard model. Theoretically, the system of bosonic atoms in a sufficiently deep optical lattice coupled coherently to the ${}^{1}S_{0}+{}^{3}P_{1}$ molecular state via the PA laser is well described by the Markovian master equation for the coupled atom-molecule mixture model Rousseau-08 with a one-body molecular loss term. By adiabatically eliminating the molecular degrees of freedom on the basis of a second-order perturbation theory for the master equation Brennen-05 ; GarciaRipoll-09 , we derive the effective master equation (see S1 in the Supplementary Materials for details),

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

$\hat{a}_{j}$ denotes the annihilation operator of atoms at site $j$ and $\hat{n}_{{\rm A},j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$. $\langle j,k\rangle$ represents nearest-neighboring pairs of lattice sites. This model is nothing but the single-component Bose-Hubbard model with a two-body loss term GarciaRipoll-09 , where $J$, and $\Gamma_{\rm PA}$ denote the hopping energy, and the strength of the two-body inelastic collision induced by PA. In Eq. (4), $g$ and $\Gamma_{\rm M}$ denote the strengths of the atom-molecule coupling and the one-body molecular loss, respectively. While $g$ is controllable by varying the intensity of the PA laser, $\Gamma_{\rm M}$ is fixed for a specific molecular state. Note that the effective master equation (1) is valid only when $\hbar\Gamma_{\rm M}\gg\max(|g|,|D|,|W|,|U|,J)$, where $D$ and $W$ denote the detuning of the PA coupling and the on-site interaction between an atom and a molecule, respectively. Since $\max(|g|,|D|,|W|,|U|,J)/\hbar\sim 100\,{\rm kHz}$ at most and $\Gamma_{\rm M}\sim 1\,{\rm MHz}$ in our experiments (see S2 in the Supplementary Materials), this condition is safely satisfied.

Stability of the atoms with a unit-filling initial state

Before studying the effect of the dissipation on the quantum phase transition, we investigate the stability of the atoms with a unit-filling initial state at a fixed lattice depth. Here, the strength of the dissipation is varied in a wide range from the weak region, in which the dissipation acts as perturbation, to the strong region, in which it exceeds any other energy scale. In contrast to the previous works in which the experiments were done only in the limited range of the dissipation strength, the wide range of our engineered dissipation enables us to observe a crossover between qualitatively different roles of the dissipation. This gives a clue to understand the many-body physics with the dissipation.

In contrast to the measurement of $\beta_{\rm{PA}}$ shown in Fig. 1 (B), this measurement is done at a shallow lattice depth of $V_{0}=8~{}E_{R}$ in which the tunneling rate $6J/\hbar$ is 4.7 kHz. Thus, in the absence of dissipation, atoms can tunnel to neighboring sites in a timescale of this measurement, so that we can investigate how the atom tunneling, which is the only mechanism causing doubly-occupied sites, is modified by the dissipation. In our experiment, we first adiabatically load a Bose-Einstein condensate (BEC) of $1.0\times 10^{4}$ atoms into a 3D optical lattice with $V_{0}=15~{}E_{R}$, in which the state is a singly-occupied Mott insulator. Subsequently, we ramp down the lattice to $V_{0}=8~{}E_{R}$ in 0.2 ms, and apply the PA laser at the same time. The initial atom number at each site is at most unity, confirmed by the absence of the atom loss by the PA laser and also the occupancy-sensitive high-resolution laser spectroscopy Kato-16 .

The measured two-body loss rate $\kappa$ is shown in Fig. 1 (C). For weak dissipation, the two-body loss rate $\kappa$ grows as $\Gamma_{\rm{PA}}$ increases, which reflects the increase of the detection rate of tunneling. For strong dissipation, however, the two-body loss rate decreases when $\Gamma_{\rm{PA}}$ increases. Namely, atom loss is suppressed by the strong on-site dissipation. This counter-intuitive behavior is a manifestation of the continuous quantum Zeno effect misra1977zeno , that is, the strong two-body inelastic collision plays a role of the strong measurement and suppresses the coherent process of tunneling. From the comparison between the theory and experiment in a wide range of dissipation strength, we confirm that the measured loss behavior correctly captures the theoretical prediction (see S3 in the Supplementary Materials). Note that we observe unexpectedly large atom loss for much higher intensity of PA laser, which prevents the suppression of the two-body loss rate for the strong dissipation region over $\gamma\sim 5$ from clear observation (see S4 in the Supplementary Materials). Therefore, in our experiment, we restrict the region of the dissipation strength under $\gamma\sim 5$.

Effect of the dissipation on the quantum phase transition

We next investigate the effect of the on-site dissipation on the quantum phase transition from the Mott insulator to the superfluid greiner2002quantum , which is the main topic of the present work. Specifically, starting with a singly-occupied Mott insulating state, we analyze the dynamics of the atoms subjected to the PA laser during a ramp-down of the lattice depth. The ramp-down speed is $-2\,E_{\rm R}/{\rm ms}$, which is much slower than the case of the two-body loss measurement discussed above. Before presenting the experimental observation, we theoretically analyze such dynamics by assuming a homogeneous system and solving the effective master equation (1) within the Gutzwiller mean-field approximation Diehl-10 ; Tomadin-11 in order to obtain some insights on the problem. Details of the theoretical analyses are shown in S3 in the Supplementary Materials.

An important effect of the dissipation on the quantum phase transition is that it explicitly breaks the conservation of the particle number of the system. Since the superfluid-Mott insulator transition at $\gamma=0$ is originated from the U(1) symmetry associated with the particle-number conservation, the introduction of finite $\gamma$ changes the transition to a crossover. Notice, however, that the two-body loss term does not explicitly break the U(1) symmetry. Indeed, the master equation (1) is invariant under the U(1) transformation, $\hat{a}_{j}\rightarrow\hat{a}_{j}e^{i\varphi}$, where $\varphi$ is an arbitrary constant.

This crossover can be theoretically characterized by the growth rate of the superfluid order-parameter amplitude; when the growth rate is smaller, the system is deeper in the Mott insulator region. According to this characterization, we find that in the strong on-site dissipation region, where $\gamma\gg 1$, the Mott insulating state is more favored for larger $\gamma$ (see Fig. S10 (A) of the Supplementary Materials). This effect originates from the quantum Zeno suppression of the tunneling, which is observed in the two-body loss rate measurement.

This interesting effect of the on-site dissipation on the crossover manifests as the delay in the melting of the singly-occupied Mott insulator in the ramp-down dynamics. In Figs. 2 (A), we show the atom number per site $\langle\hat{n}_{\rm A}\rangle$ and the condensate fraction $|\psi|^{2}/\langle\hat{n}_{\rm A}\rangle$ as functions of the instantaneous lattice depth during the ramp-down dynamics. We clearly see that in the strong dissipation region the onset of the atom loss or the order parameter growth shifts to the side of small lattice depth as $\gamma$ increases. This result suggests that one may experimentally observe the delay in the Mott-insulator melting by measuring the time evolution of the atom number and the momentum distribution during the ramp-down dynamics.

Having the above theoretical insights in mind, we perform the experiment for measuring ramp-down dynamics across the crossover from the Mott insulator to the superfluid. The atom number and the momentum distribution during ramp-down dynamics are obtained from the fluorescence detection and the density distribution of the time-of-flight (TOF) absorption image, respectively. Our experiment starts with ramping up the lattice to $V_{0}=20~{}E_{R}$ for preparation of the singly-occupied Mott insulator state. The atom number is tuned to be small enough that no doubly-occupied site exists. Subsequently we ramp down the lattice with applying the PA laser. The lattice ramp-down speed is $-2\,E_{R}$/ms. After ramping down the lattice to the final lattice depth, we perform the fluorescence detection for measuring the atom number, or we suddenly turn off all the trap and take the absorption image after 8 ms ballistic expansion for obtaining the density distribution.

We first focus on the atom loss measurement during ramp-down dynamics. Figure 2 (B) shows the atom number measured with various dissipation strengths. The experimental result well reproduces the overall features of the calculation shown in Fig. 2 (A) (left). Specifically, the significant atom loss starts around $V_{0}=10~{}E_{R}$ in the presence of weak dissipation (\footnotesize2⃝), while the atom number is conserved during ramping down the lattice without dissipation (\footnotesize1⃝). This onset shifts to the deep lattice side as $\gamma$ increases (\footnotesize3⃝) for weak dissipation ($\gamma<2$). However, when $\gamma$ increases further from $\gamma\sim 2$, the onset shifts to the shallow lattice side (\footnotesize4⃝). In order to identify the onset, we fit the double linear function to the data (Fig. 2 (C)), which are shown in Fig. 2 (B). In the presence of on-site dissipation, the atom loss is correlated with the melting of the Mott insulator in ramp-down dynamics because the melting creates the double occupation which is blasted out by the PA laser. Our result suggests that the melting of the Mott insulator is delayed for strong on-site dissipation. Quantitatively, the onset changes from $V_{0}=11.7(4)~{}E_{R}$ to $V_{0}=9.2(4)~{}E_{R}$ at the maximum as $\gamma$ increases. This corresponds to the increase of $zJ/U$ by a factor of 2.2. These behaviors capture the essence of the theoretical predictions mentioned in the above.

Figure 3 (A) shows a series of TOF absorption images obtained by changing the final lattice depth from Mott insulator regime to superfluid regime with various strength of the dissipation. Without dissipation, we clearly observe the transition from a Mott insulator state to a superfluid state as shown in Fig. 3 (A) of $\gamma$ = 0: in the deep lattice such as $V_{0}=20~{}E_{R}$, we obtain a broad distribution with no pattern, which indicates that the atoms have no phase coherence corresponding to the Mott insulator state. As ramping down the lattice across the critical depth of $V_{0}=11.3~{}E_{R}$, which is calculated from the scattering length of ${}^{1}S_{0}$ state of 174Yb PhysRevA.77.012719 , we obtain a clear interference pattern characterizing the presence of the phase coherence of the superfluid state. In the presence of dissipation, the observed transition is significantly modified, as shown in Fig. 3 (A). As the strength of the dissipation increases, the interference pattern becomes unclear in the shallow lattice regime. For strong dissipation such as $\gamma\sim 5$, any pattern cannot be observed. This result indicates that the growth of the phase coherence is suppressed by the strong dissipation.

In order to evaluate the phase coherence quantitatively, we introduce the visibility of the interference peaks as $v=(N_{\rm{max}}-N_{\rm{min}})/(N_{\rm{max}}+N_{\rm{min}})$ PhysRevA.72.053606 . Here, $N_{\rm{max}}$ is the sum of the number of atoms in the regions of first-order interference peaks and $N_{\rm{min}}$ is the sum of the number of atoms in the regions at the same distance from the central peak along the diagonals. While the visibility increases with the ramp-down of the lattice, this increase becomes more moderate in the stronger dissipation, as shown in Fig. 3 (B). Especially, clear dependence on the strength of the dissipation is observed below the depth of $V_{0}=11~{}E_{R}$, which is around the calculated critical depth at $\gamma=0$. As shown in Fig. 3 (B), the effect of the dissipation on the width of the crossover region is observed as more moderate growing of the visibility below the depth of $11~{}E_{R}$ associated with the increase of $\gamma$. In addition, as shown in Fig. 3 (C), the dissipation moderates the narrowing of the width (in $\mu$m) of the density distribution and the slope of the width (in $\mu$m) with respect to the lattice depth becomes less steep as $\gamma$ increases. Narrowing the width of the distribution indicates the localization of the state in the momentum space. All of these measurements support the delay in the melting of the singly-occupied Mott insulator in the ramp-down dynamics as an effect of the on-site dissipation as we see in the calculation of the condensate fraction shown in Fig. 2. Notice that the observation of the excitation gap, which is the direct evidence of the formation of the Mott insulator, is difficult in this dissipative system because the excitation spectrum should have a broad linewidth determined by the inelastic collision rate $\Gamma_{\rm PA}$ of a few tens of kHz.

Quenching the dissipation

It is important to experimentally check whether this behavior is attributed to some heating effect by the PA laser. For this purpose, we measure the phase coherence after turning off the PA laser. If the absence of the interference pattern is attributed to the heating, the phase coherence is no longer restored after the PA laser is turned off. In contrast, if the state after the ramp-down of the lattice is still a Mott insulator, the phase coherence can be restored. Similarly to the measurement of Fig. 3 (A), we ramp down the lattice to a final lattice depth in $-2\,E_{R}$/ms with the maximum strength of dissipation $\gamma$ = 4.6(4). Then, we suddenly turn off the PA laser and investigate the subsequent time evolution of the atoms in the lattice by observing the phase coherence through a TOF absorption image at some hold time.

The result for the case of the final lattice of $V_{0}=8~{}E_{R}$ is shown in Fig. 4 (B) for the observed TOF images and Fig. 4 (C) for the evolution of the visibility and width of the density distribution. After some hold time, an interference pattern grows. It serves as a direct signature of the restoration of the phase coherence, indicating that the absence of the interference pattern in Fig. 3 (A) is not completely attributed to the heating. We confirm that the total atom number is conserved in this dynamics as shown in Fig. 4 (E). This means that the evaporative cooling which could possibly explain the observed behavior does not occur during the dynamics. Importantly, the atoms after turning off the dissipation can be considered as an isolated (closed) system. Therefore, the observed dynamics in our experiment should not be considered as the usual thermal relaxation with the environment, but the relaxation in the isolated quantum system, which is a hot topic actively studied in recent experiments and theories Polkovnikov-11 ; Langen-16 ; Eisert-15 ; Nandkishore-15 . Here, we consider the tunneling time as a relevant time scale because the superfluid state is realized through the process of delocalization of the particles by the tunneling. As shown in Fig. 4 (C), the time constant of the increase of the visibility and the decrease of the width is comparable to the tunneling time $(6J/\hbar)^{-1}$ = 0.21 ms. Figure 4 (D) shows the visibility and the width of the density distribution, 0 ms and 4 ms after the ramp-down for various final lattice depths, similarly indicating the restoration of the coherence.

III Discussion

We have realized the engineered dissipative Bose-Hubbard system by introducing a controllable strength of two-body inelastic collision with use of a PA laser. By exploiting the highly controllable nature of the dissipation, we have investigated the effect of the dissipation on the quantum phase transition from the Mott insulator state to the superfluid state in the lattice ramp-down dynamics. We have observed that the melting of the Mott state is delayed and the growth of the phase coherence is suppressed for the strong on-site dissipation. The favored state depends on the type of the dissipation. For example, the stabilization of the superfulid state with use of a well-designed off-site dissipation is proposed Diehl-08 . Thanks to the dramatic change in the onset of the Mott-insulator melting, as shown in the fact of increase of $zJ/U$, we can access the interesting problem of quenching the dissipation across the crossover from the Mott insulator to the superfluid, where turning off the dissipation corresponds to a sudden parameter change of the Bose-Hubbard system Altman-02 . In this method, the required time for turning off the dissipation could be very short while the sudden change of the depth of optical lattice needs a certain time in order to stabilize the power of the lattice laser as well as to prevent non-adiabatic inter-band transition. Moreover, while we have used 174Yb, which is a bosonic isotope of an alkali-earth-like species, to demonstrate our method for controlling the dissipation, it is applicable generally to other atomic species that can be coupled to a state of lossy PA molecule. The crossover properties can also be caused by varying the on-site interaction (see Fig. S10 (B) of the Supplementary Materials). Controlling the on-site interaction with Feshbach resonance, for example using alkali atoms, enables us to investigate wider range of strength of dissipation including infinitely strong, since weakening the on-site interaction corresponds to strengthen the dissipation $\gamma$. Our work opens a new way to study the quantum many-body system by controlling the dissipation.

IV Materials and Methods

Preparation of ultracold 174Yb atoms in an optical lattice

After collecting atoms with a magneto-optical trap (MOT) with the intercombination transition of ${}^{1}S_{0}\leftrightarrow{}^{3}P_{1}$, we load the atoms into a crossed far-off-resonant trap (FORT). Subsequently, an evaporative cooling is performed, resulting in an almost pure BEC with no discernible thermal component. The atom number of BEC is controlled by changing the collecting time of atoms in MOT. The trap frequencies of the FORT at the final stage of the evaporative cooling are $(\omega_{x^{\prime}},\omega_{y^{\prime}},\omega_{z})/2\pi=(162,31,166)$ Hz. Here, the $x^{\prime}$- and $y^{\prime}$-axes are tilted from the lattice axes ($x$ and $y$) by 45*∘* in the same plane. Then, the BEC is loaded into a 3D optical lattice with the wavelength of the lattice laser $\lambda_{L}$ of 532 nm. A typical atom number of BEC is $6\times 10^{4}$ in a measurement of the photo-association rate of doubly-occupied sites, or $1\times 10^{4}$ in a measurement of two-body loss rate and observation of the quantum phase transition for atoms without multiply-occupied sites.

Detail of PA

By applying the PA laser to the atoms in the lattice, two atoms in the doubly-occupied sites are photoassociated into the ${}^{1}S_{0}+{}^{3}P_{1}$ molecular state whose vibrational quantum number $v^{\prime}_{e}=16$. These molecules are immediately dissociated into the two ground state atoms. We note that the photon scattering of the atoms due to the PA laser can be negligible because the selected PA line is 3.7 GHz below the atomic transition line, which is about $2\times 10^{4}$ times larger than its natural line width. $\beta_{\rm{PA}}$ for various intensities is determined through the loss dynamics of the atoms by measuring the remaining atom number with the fluorescence detection method after applying the PA laser to the atoms in the lattice with depth of $V_{0}=14~{}E_{R}$. In this lattice depth, the system is in a state with singly- and doubly-occupied sites. After the PA laser is applied, the remaining atom number $N(t)$ decreases as $N(t)=N_{1}+N_{2}\exp(-\Gamma_{\rm{PA}}t)$, where $N_{1}$ and $N_{2}$ are the initial atom number in the singly and doubly-occupied site respectively.

Acknowledgements

We thank Y. Ashida, M. Barbier and P. Naidon for fruitful discussions. This work was supported by MEXT/JSPS KAKENHI Grant Numbers JP25220711, JP26247064, JP16H00990, JP16H01053, JP16H00801, CREST, JST No. JPMJCR1673, and Matsuo Foundation. T.T. acknowledges support from the Japan Society for the Promotion of Science (KAKENHI Grant Number JP16J01590).

Author contributions

T. Tomita, S. Nakajima and Y. Takasu carried out experiments and the data analysis. I. Danshita carried out the theoretical calculation. Y. Takahashi conducted the whole experiment. All the authors contributed to the writing of the manuscript.

Competing interests

The authors declare no competing financial interests.

Data and materials availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from T.T. (tomita@scphys.kyoto-u.ac.jp).

S1 Derivation of the dissipative Bose-Hubbard model

In this section, we present a detailed derivation of the dissipative Bose-Hubbard model with a two-body loss term from the coupled atom-molecule mixture model with a one-body molecular loss term. The derivation is based on a kind of perturbative approach developed in Ref. garcia-09 .

S1.1 Models

We start with the Markovian master equation for the coupled atom-molecule mixture model with a one-body molecular loss term,

[TABLE]

The Hamiltonian is given by

[TABLE]

where

[TABLE]

The quantum phases of this atom-molecular Hamiltonian have been theoretically studied in Ref. rousseau-08 . We anticipate that the molecule consists of a ${}^{1}S_{0}$ atom and a ${}^{3}P_{1}$ atom as in the experiment such that its linewidth is on the order of 1 MHz due to the short lifetime of the latter state of atom. Hence, we have to include the one-body loss term of molecules in the master equation,

[TABLE]

Equation (S1.1) describes the dynamics of ultracold bosonic atoms in an optical lattice coupled with a molecular state via photo-association (PA) laser. $\hat{a}_{j}$ and $\hat{m}_{j}$ annihilate an atom and a molecule on site $j$ while $\hat{n}_{{\rm A},j}=\hat{a}^{\dagger}_{j}\hat{a}_{j}$ and $\hat{n}_{{\rm M},j}=\hat{m}^{\dagger}_{j}\hat{m}_{j}$ are the density operators of atoms and molecules. $\Gamma_{\rm M}$, $D$, $U$, $W$, $g$, and $J$ denote the one-body loss of molecules, the detuning of the PA coupling from the molecular state, the on-site interaction between two atoms, the on-site interaction between an atom and a molecule, the atom-molecule coupling, and the hopping of atoms. $\langle j,k\rangle$ represents nearest-neighboring pairs of lattice sites. Since $\hbar\Gamma_{\rm M}\gg\max(|D|,|g|,|J|)$, a molecule created on a lattice site via the PA laser decays much earlier than the creation of another molecule on the same site. In this sense, we can safely assume the hardcore constraint on the molecules, which forbids more than one molecules to occupy a single site. We aim to show that when $\hbar\Gamma_{\rm M}\gg\max\left(|D|,|U|,|W|,|g|,|J|\right)$, the molecular degrees of freedom can be properly projected out by means of a perturbation theory garcia-09 . such that the system is well approximated by the following effective master equation,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Equations (S1.7)-(S1.10) are written also in the main text as Eqs. (1)-(4).

In order to express the density matrix of the system more explicitly, we define the local Fock state on site $j$,

[TABLE]

We set the maximum number of atoms per site to be $d-1$. Hence, the dimension of the local Hilbert space is $2d$. While the maximum number of bosonic atoms per site is in principle the total number of atoms, at a finite filling factor the occupation probability of large-$n_{a}$ states decays exponentially. This means that setting the cutoff of the local Hilbert space at $n_{a}=d-1$ does not affect results of numerical calculations in practice as long as $U>0$ and $d$ is sufficiently large kuhner-98 .

For convenience, we introduce a simpler notation for the local state,

[TABLE]

Using the local states defined above, we express a general form of the density matrix as

[TABLE]

Regarding the density matrix $\hat{\rho}$, which is a $(2d)^{M}\times(2d)^{M}$ matrix, as a $(2d)^{2M}$-dimensional vector ${\bm{\rho}}$, we can rewrite the master equation in the following form,

[TABLE]

where the superoperators $\hat{\mathcal{M}}_{0}$, $\hat{\mathcal{V}}_{\rm am}$, and $\hat{\mathcal{V}}_{\rm hop}$ are $(2d)^{2M}\times(2d)^{2M}$ matrices originated from $\hat{H}_{0}$ and $L_{\rm M}(\hat{\rho})$, $\hat{H}_{\rm am}$, and $\hat{H}_{\rm hop}$, respectively. Associated with the change of the notation from $\hat{\rho}$ to ${\bm{\rho}}$, we also express the local matrix as the following vector,

[TABLE]

This vector satisfies the orthonormality condition,

[TABLE]

S1.2 Local projection

In order to derive the effective master equation (S1.7), we need to express some of the superoperators explicitly and introduce the projection superoperator. The non-perturbative superoperator $\hat{\mathcal{M}}_{0}$ can be expressed as a sum of commuting local superoperators $\hat{\mathcal{M}}^{\rm loc}_{0,j}$,

[TABLE]

where

[TABLE]

and

[TABLE]

The perturbative superoperator $\hat{\mathcal{V}}_{\rm am}$ originated from the atom-molecule coupling $\hat{H}_{\rm am}$ can be also expressed as a sum of commuting local superoperators,

[TABLE]

where

[TABLE]

and

[TABLE]

We omitted the site index $j$ of $\hat{\mathcal{M}}^{\rm loc}_{0,j}$ and $\hat{\mathcal{V}}^{\rm loc}_{{\rm am},j}$ in Eqs. (S1.20) and (S1.25) because they do not depend on $j$.

We do not write an explicit expression of the superoperator $\hat{\mathcal{V}}_{\rm hop}$ because it is unnecessary for our purpose. Nevertheless, it is worth noting that $\hat{\mathcal{V}}_{\rm hop}$ changes neither the number of atoms nor the number of molecules. This means that this superoperator does not have matrix elements connecting the effective Hilbert space with the truncated one.

In order to construct the projection superoperator, we need to solve the following eigenvalue problem,

[TABLE]

Notice that the left eigenvector $(w_{\alpha}|$ in general is not equal to the conjugate of the right eigenvector $|v_{\alpha})$ because $\hat{\mathcal{M}}_{0}^{\rm loc}$ is not Hermitian. The eigenvectors satisfy the following orthonormality condition,

[TABLE]

When we derive the effective model, we utilize the fact that $d^{2}$ eigenvalues have the property $|\lambda_{\alpha}|=O(\max(|D|,|U|,|W|))$ and the other eigenvalues have $|\lambda_{\alpha}|=O(\hbar\Gamma_{\rm M})$. The local subspaces that include states with the former and latter properties are denoted by $\mathcal{D}_{\rm eff}^{\rm loc}$ and $\mathcal{D}_{\rm trc}^{\rm loc}$, respectively.

We define the local projection superoperators as

[TABLE]

This superoperator projects a state or a superoperator on state $\alpha$. From these projectors, we construct the projection superoperator on the effective Hilbert space,

[TABLE]

Notice that from the first line to the second line of Eq. (S1.31), we neglected the terms on the order of $\max(|D|,|U|,|W|)/(\hbar\Gamma_{\rm M})$ on the basis of the assumption that $\hbar\Gamma_{\rm M}\gg\max(|D|,|U|,|W|)$.

S1.3 Second-order perturbation

From the local projection operator of Eq. (S1.31), we construct the projection operators for the entire system as

[TABLE]

where $\hat{\mathcal{I}}$ is the global identity matrix. Multiplying $\hat{\mathcal{P}}_{\rm eff}$ on Eq. (S1.16) from the left and using the facts that $\hat{\mathcal{P}}_{\rm eff}\hat{\mathcal{V}}_{\rm am}\hat{\mathcal{P}}_{\rm eff}=0$ and $\hat{\mathcal{P}}_{\rm eff}\hat{\mathcal{V}}_{\rm hop}\hat{\mathcal{P}}_{\rm trc}=0$, we obtain

[TABLE]

where

[TABLE]

In Eq. (S1.34), it is obvious that the first and second terms in the right hand side already agree with the term $-i\left[\hat{H}_{\rm eff},\hat{\rho}_{\rm eff}\right]$ in Eq. (S1.7).

We will next derive $L_{2}(\hat{\rho}_{\rm eff})$ in Eq. (S1.7) from the last term in the right hand side of Eq. (S1.34). The components in $\hat{\mathcal{P}}_{\rm trc}$ that can give finite contributions to $\hat{\mathcal{P}}_{\rm eff}\hat{\mathcal{V}}_{\rm am}\hat{\mathcal{P}}_{\rm trc}$ are the ones written as

[TABLE]

where

[TABLE]

and $\beta\in\mathcal{D}_{\rm trc}^{\rm loc}$. Hence, Eq. (S1.34) can be rewritten as

[TABLE]

In order for Eq. (S1.39) to be closed within the effective Hilbert space, we need to express ${\bm{\rho}}_{\beta}$ in terms of ${\bm{\rho}}_{\rm eff}$. For this purpose, we look into the equation for ${\bm{\rho}}_{\beta}$ given by

[TABLE]

where

[TABLE]

We neglect the last term in the right hand side of Eq. (S1.40) because it gives higher-order contributions with respect to $\hat{\mathcal{V}}$. Making a variable transformation,

[TABLE]

Eq. (S1.40) is simplified a little,

[TABLE]

Formally solving Eq. (S1.46), we obtain

[TABLE]

Performing a partial integral, this solution becomes

[TABLE]

The remaining integral can be neglected because it is of higher order in $\max(|D|,|U|,|W|,|g|,|J|)/(\hbar\Gamma_{\rm M})$. The second term in Eq. (S1.48), which includes $e^{\lambda_{\beta}t/\hbar}\tilde{\bm{\rho}}_{\rm eff}(0)$, decays very quickly on the order of $1/\Gamma_{\rm M}$ so that it can be also neglected as long as we are interested in much longer time scale than $1/\Gamma_{\rm M}$. Moreover, $\lambda_{\beta}=-\frac{\hbar\Gamma_{\rm M}}{2}$ in its leading order. Hence, ${\bm{\rho}}_{\beta}$ is well approximated as

[TABLE]

Substituting Eq. (S1.49) into Eq. (S1.39), we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Rewriting Eq. (S1.7) with use of superoperators and the vector form of the density matrix, and setting $\Gamma_{\rm PA}=8g^{2}/(\hbar^{2}\Gamma_{\rm M})$, we recognize that the derived effective master equation, which is Eq. (S1.50), is equivalent to that for the dissipative Bose-Hubbard model with the two-body loss term. Thus, we have successfully derived the dissipative Bose-Hubbard model from the coupled atom-molecule mixture model with the one-body molecular loss term.

S2 Loss dynamics from the Mott insulating state with double filling

In this section, we analyze the dynamics of the Mott insulator with two bosonic atoms per site in an optical lattice subjected to a sudden increase of the atom-molecule coupling $g$ from zero.

As an initial condition, we assume that the system is deep in a Mott insulating state of atoms with double filling and that the atom-molecule coupling $g$ is zero. In such a situation we can safely neglect the hopping term of atoms as long as we are interested in the atom-loss dynamics and its timescale $1/\Gamma_{\rm PA}$ is much shorter than that of the hopping, which is on the order of $\hbar U/J^{2}$ in the Mott insulator. Thus, the system can be described by the following single-site master equation,

[TABLE]

The Hamiltonian is given by

[TABLE]

where

[TABLE]

The one-body molecular loss term is given by

[TABLE]

Since the initial state is $|n_{\rm A}=2,n_{\rm M}=0\rangle$, the Hilbert space necessary for describing the dynamics of Eq. (S2.1) is spanned by only three states, namely, $|n_{\rm A}=0,n_{\rm M}=0\rangle$, $|n_{\rm A}=2,n_{\rm M}=0\rangle$, and $|n_{\rm A}=0,n_{\rm M}=1\rangle$. When $\hbar\Gamma_{\rm M}\gg|g|$, we can properly eliminate the state $|n_{\rm A}=0,n_{\rm M}=1\rangle$ by means of the perturbation theory used in the previous section, to derive the effective master equation,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

In contrast to the case in the previous section, the effective model is valid without the condition that $\hbar\Gamma_{\rm M}\gg\max(|D|,|U|,|W|)$ and we could obtain the analytical expression of $\tilde{\Gamma}_{\rm PA}$, which includes the explicit dependence on $U$ and $D$. It is obvious that when $\hbar\Gamma_{\rm M}\gg\max(|D|,|U|,|W|)$, $\tilde{\Gamma}_{\rm PA}$ coincides with $\Gamma_{\rm PA}$.

The dynamics of the effective master equation (S2.6) involves only the two states such that we can easily obtain its analytical solution,

[TABLE]

In order to check the validity of the effective master equation (S2.6), in Fig. S1 we compare Eq. (S2.10) with the numerical solution of the original master equation (S2.1) including the molecular degree of freedom. We see that the analytical and numerical results agree when $\hbar\Gamma_{\rm M}\gg g$.

The above results indicate that when the condition $\hbar\Gamma_{\rm M}\gg g$ is safely satisfied, one can determine $\tilde{\Gamma}_{\rm PA}$ in experiment by measuring the atom-loss dynamics in the double-filling Mott insulating state and extracting the exponent in the exponential decay of the atom number. Once $\tilde{\Gamma}_{\rm PA}$ is measured as a function of the detuning $D$, $\Gamma_{\rm M}$ can be also determined by fitting $\tilde{\Gamma}_{\rm PA}(D)$ to the Lorentzian function of Eq. (S2.9). In Fig. S2, we show measured $\tilde{\Gamma}_{\rm PA}$ as a function of $D-U$. Fitting the data to Eq. (S2.9), we determine $\Gamma_{\rm M}=2\pi\times$185(13) kHz = 1.16(8) MHz. On the other hand, $g$ can be estimated from the experimentally determined $\Gamma_{\rm M}$ and $\Gamma_{\rm PA}$, resulting in $g/\hbar\sim 100~{}{\rm kHz}$ at most. Therefore we can confirm that the condition $\hbar\Gamma_{\rm M}\gg g$ is satisfied, as well as the condition $\hbar\Gamma_{\rm M}\gg\max(|D|,|U|,|W|)$ since $D/\hbar$, $U/\hbar$ and $W/\hbar$ are at most a few 10 kHz in our experiment.

S3 Details of the theoretical analyses using the Gutzwiller variational approach

In this section, we present some details of theoretical calculations regarding dynamics starting from the Mott insulator with one bosonic atom per site in the presence of the atom-molecule coupling that leads to the loss of atoms. The corresponding dynamics is experimentally analyzed in the main text. For this purpose, we first review how to solve the master equations (S1.1) and (S1.7) within the Gutzwiller mean-field approximation  diehl-10 ; tomadin-11 . Using the introduced prescription, we next compute the atom-loss dynamics after a fast ramp-down of the lattice depth and that during a slow ramp-down of the lattice depth.

Here we explicitly explain the Gutzwiller mean-field theory applied to the atom-molecule mixture model of Eq. (S1.1) because it contains the molecular degrees of freedom, which has not been taken into account in previous studies. We note that one can easily apply the same prescription to the effective Bose-Hubbard model of Eq. (S1.7) in a very similar manner. In the Gutzwiller mean-field approximation, the many-body density matrix is assumed to be a single product of local density matrices,

[TABLE]

where

[TABLE]

From the Gutzwiller density matrix, we define the local superfluid order parameter as

[TABLE]

In the Gutzwiller mean-field approximation, we ignore the second order terms with respect to the fluctuation of $\hat{a}_{j}$ from its mean value $\psi_{j}$. In this way, the Hamiltonian is simplified as

[TABLE]

where

[TABLE]

[TABLE]

In Eq. (S3.6), $\langle k\rangle_{j}$ means sites nearest-neighboring to $j$. It is worth noting that $\hat{H}_{j}^{\rm GW}$ includes only local operators at site $j$ because one of field operators in the hopping term is replaced with its mean value. Thanks to this property, the master equation under the Gutzwiller mean-field approximation is closed within local site $j$,

[TABLE]

such that we can solve the master equation at a very low numerical cost.

While the Gutzwiller approximation is a simple mean-field theory, it has been extensively used to study various phenomena and properties of Bose gases in optical lattices, including the quantum phase transitions rokhsar-91 ; sheshadri-93 ; iskin-11 , the elementary excitations kovrizhin-05 ; krutitsky-11 , the superfluid critical momentum altman-05 ; saito-12 , and the non-equilibrium dynamics snoek-07 ; bissbort-11 ; Snoek-11 . Recently, it has been applied for solving the master equation of the Bose-Hubbard system with dissipation terms diehl-10 ; tomadin-11 ; leboite-13 ; vidanovic-14 . This approximation is more accurate in higher dimensions, where there are more mean fields to interact with. In the case of the Bose-Hubbard model on a cubic lattice, for instance, the Gutzwiller approximation gives the critical point for the superfluid-Mott insulator quantum phase transition at unit filling as $zJ/U=0.1716$ while that by the quantum Monte Carlo method is $zJ/U=0.2045$ capogrosso-07 , where $z$ is the coordination number. Since the experimental system considered here is three dimensional, the Gutzwiller approximation can give reliable results at least qualitatively.

Assuming that the system is homogeneous, we further simplify the master equation, i.e., we drop the dependence on the site index $j$. This assumption means that we neglect fluctuations other than the zero-momentum one and the effect of the trapping potential. Since the trapping potential is present in the actual experiment, our theoretical analyses within this simplification do not correspond to the experiment at a quantitative level. We emphasize that the main purpose of our theoretical analyses is to provide qualitative explanations for the interesting effects of the engineered dissipation observed in the experiment.

S3.1 Dynamics after a fast ramp-down of the lattice depth

We consider the atom-loss dynamics from the initial unit-filling Mott insulator induced by the PA laser after a fast ramp-down of the lattice depth. The corresponding dynamics is experimentally analyzed in the main text (see, e.g., Fig. 2). The optical lattice potential is given by

[TABLE]

where $k=\pi/d$ with the lattice constant $d=266\,{\rm nm}$. We set the time sequence of the lattice depth,

[TABLE]

where the ramp-up and ramp-down speeds are given by

[TABLE]

Once the lattice depth is given, we numerically compute the Wannier function $w_{j}({\bm{r}})$ localized at site $j$, from which we can determine the on-site interaction and the hopping energy as

[TABLE]

Sites $j$ and $k$ are assumed to be nearest neighboring. The $s$-wave scattering length of 174Yb used in the experiment is $a_{s}=5.55\,{\rm nm}$ kitagawa-08 . In Fig. S3, we plot the on-site interaction $U$ and the hopping energy $J$ as a function of the lattice depth.

Let us elaborate the time sequence of the dynamics computed here. We start with the superfluid ground-state at $V_{0,{\rm ini}}=5E_{\rm R}$ and $t_{0}=-100.2\,{\rm ms}$. Strictly speaking, in the experiment, there is an additional ramp-up process from $V_{0}=0E_{\rm R}$ to $5E_{\rm R}$, but we can not take into account this process because the Bose-Hubbard model is invalid for such a shallow optical lattice. Nevertheless, the ramp-up speed in the experiment is so slow that the superfluid state prepared at $V_{0,{\rm ini}}=5E_{\rm R}$ can be regarded as the ground state. We next prepare the Mott insulating state at $V_{0,{\rm max}}=15E_{\rm R}$ and $t_{1}=-0.2\,{\rm ms}$, which imply that the ramp-up speed is $v_{\rm up}=0.1\,E_{\rm R}/{\rm ms}$. Right after the preparation of the Mott insulating state, we ramp down the lattice depth to $V_{0,{\rm fin}}=8E_{\rm R}$ in $0.2\,{\rm ms}$, implying that $v_{\rm down}=-35\,E_{\rm R}/{\rm ms}$ and $t_{\rm 2}=0$. Finally, we turn on the dissipation term $\Gamma_{\rm PA}$ and keep $V_{0}(t)=V_{0,{\rm fin}}$ until $t_{\rm fin}=20\,{\rm ms}$. The time sequence for $V_{0}/E_{\rm R}$, $zJ/U$, and $\gamma=\hbar\Gamma_{\rm PA}/U$ is summarized in Fig. S4.

In Fig. S5, we plot the time evolution of the atom number per site $\langle\hat{n}_{\rm A}\rangle$ at $t>t_{2}=0$, where $\gamma=0.05$ (A) and $2$ (B). We obviously see that $\langle\hat{n}_{\rm A}\rangle$ decays as the time evolves. In order to extract the loss rate from the numerical data, we use the fitting function,

[TABLE]

where the loss rate $\kappa$ is treated as a free parameter. Notice that the same fitting function is used to extract the loss rate from the experimental data (see Fig. 1 (C) of the main text). In Fig. S6, the extracted loss rate is plotted as a function of $\gamma$. While the function of Eq. (S3.17) is well fitted to $\langle\hat{n}_{\rm A}\rangle(t)$ at small $\gamma$ as shown in Fig. S5(A), it is worse for relatively large $\gamma$ (see, e.g., Fig. S5(B)). This leads to small and large error bars for small and large $\gamma$, respectively.

We see from Fig. S6 that when the dissipation strength $\gamma$ increases, the loss rate $\kappa$ initially increases but starts to decrease around $\gamma=1$. A similar behavior is also seen in the experiment as shown in Fig. 1 (C) of the main text. The decreasing loss rate at $\gamma\gtrsim 1$ can be attributed to the suppression of the double occupancy due to the strong two-body loss, namely the continuous Zeno effect. In order to corroborate this interpretation, we show the time evolution of $\rho_{3,3}$, which is an element of the local density matrix corresponding to the probability of the double occupancy, in Fig. S7. We take the maximum value of $\rho_{3,3}$ during the time evolution at $t>t_{2}=0$ and plot it as a function of $\gamma$ in Fig. S8(A). We clearly see that the double occupancy is strongly suppressed due to the strong dissipation. Roughly speaking, the loss rate can be estimated by $\Gamma_{\rm PA}\times\rho_{3,3}^{\rm max}$, which is plotted as a function of $\gamma$ in Fig. S8(B). There we see the tendency similar to the actual loss rate shown in Fig. S6.

S3.2 Dynamics during a slow ramp-down of the lattice depth

We consider the dynamics of the initial unit-filling Mott insulator during a slow ramp-down of the lattice depth towards the supefluid state. By analyzing this type of dynamics, we aim to understand the effect of the two-body loss term $\Gamma_{\rm PA}$ on the quantum phase transition between the superfluid and Mott insulator. As mentioned in the main text, an important effect is that the two-body loss term explicitly breaks the conservation of the particle number of the system. Since the superfluid-Mott insulator transition at $\Gamma_{\rm PA}=0$ is originated from the U(1) symmetry associated with the particle-number conservation, the introduction of finite $\Gamma_{\rm PA}$ changes the transition to a crossover. Notice, however, that the explicit breaking of the particle-number conservation does not mean that of the U(1) symmetry in our dissipative system. This is in clear contrast to a closed system with no dissipation.

Another important effect is that the two-body loss term makes the “superfluid” state at unit filling so dissipative that it cannot carry dissipationless superflow. In this sense, even a small loss term immediately breaks the superfluidity. However, one can distinguish such a lossy gas with delocalized atoms and long-range coherence from the Mott insulating state and study the crossover from the latter state to the former in the dynamics subjected to a slow ramp-down of the lattice depth. Hereafter, for convenience we use the term “superfluid” to describe the former state.

Before recreating the experimental situation, we present some important properties of the crossover phenomenon that are independent of either the preparation procedure of the initial Mott insulating state or the ramp-down speed within the Gutzwiller approximation. As elaborated below, we specifically focus on the growth rate of the superfluid order parameter amplitude during a slow increase of the hopping energy or the on-site interaction. The time sequence of the hopping energy is given by

[TABLE]

while the other parameters are fixed to be time-independent, where $\tau$ denotes the total evolution time. In the case that we vary the on-site interaction, its time sequence is given by

[TABLE]

In our experiment shown in the main text, we ramp down the lattice depth to cause the crossover from the Mott insulator to the “superfluid”. Since the dominant effect of the ramp down of the lattice depth is the exponential increase of the hopping energy, the case of $J$ increase is closer to the experimental situation. Notice that our 174Yb atoms in their electronic ground state do not have usable Feshbach resonance such that we can not dynamically control $U$ to realize the case of $U$ decrease.

We set the initial and final values of the hopping energy as $zJ_{\rm ini}/U=0$ and $zJ_{\rm fin}/U=0.5$, and those of the on-site interaction energy as $U_{\rm ini}/(zJ)=50$ and $U_{\rm fin}/(zJ)=0.5$. At $zJ/U=zJ_{\rm ini}/U$ or $zJ/U_{\rm ini}$ at unit filling, the ground state is the Mott insulating state, i.e., $\rho_{l,m}^{\rm gs}=\delta_{l,2}\delta_{m,2}$. As an initial state of the dynamics, we add small random noise terms to $\rho_{l,m}^{\rm gs}$ as

[TABLE]

where $\epsilon_{l,m}^{\rm re}$ and $\epsilon_{l,m}^{\rm im}$ are assumed to be independent random variables with zero average and a box distribution from $-\varepsilon$ to $\varepsilon$. In the absence of the noise terms ($\varepsilon=0$), the system remains in the initial state and the evolution towards the “superfluid” state can not be captured because $\rho_{l,m}^{\rm gs}$ is a time-independent solution of the effective master equation (S1.7) within the Gutzwiller approximation.

In Fig. S9(A), we show the time evolution of the amplitude of the superfluid order parameter $|\psi|^{2}$ for several values of the noise strength $\varepsilon$. We see that $|\psi|^{2}$ significantly depends on $\varepsilon$. In contrast, as shown in Fig. S9(B), we find in the time evolution of the rate of the exponential growth in $|\psi|^{2}$, namely $G=\frac{d}{dt}\ln|\psi|^{2}$, that there is a time region where $G$ increases and is independent of $\varepsilon$. Notice that $G$ in such a time region is also independent of $\tau$ as long as $\tau$ is sufficiently large. We use the values of $G$ in the time region to characterize the time scales of the crossover from the Mott insulator to the “superfluid” that are independent of either $\varepsilon$ or $\tau$. As indicated in Fig. S9(C), for given $\gamma$ we determine the value of $zJ/U$ at which $G$ in the universal time region takes a certain value, e.g., $\hbar G/U=0.1$ (red dashed line) or $0.05$ (blue dotted line). In Fig. S10(A), we show a contour plot of $\hbar G/U$ in the $(\gamma,zJ/U)$-plane. There we see that when $\gamma$ increases from zero, the contour lines become more distant from one another, i.e., the transition is changed to a crossover. A similar behavior is also seen in Fig. S10(B), where a contour plot of $\hbar G/(zJ)$ in the $(\hbar\Gamma_{\rm PA}/(zJ),zJ/U)$-plane is shown. We also see that $zJ/U$ on each contour line exhibits a non-monotonic behavior as a function of $\gamma$; when $\gamma$ increases, it initially decreases but starts to increase above a certain $\gamma$. This result indicates that the strong two-body loss term, i.e., $\gamma\gg 1$, favors the Mott insulating state over the “superfluid”.

We argue that this interesting effect of the two-body loss term on the crossover from the Mott insulator to the “superfluid”, namely the suppression of the dynamical melting of the Mott insulating state, can be observed through the measurement of the atom number or the momentum distribution during the slow ramp-down of the lattice depth. In order to corroborate this argument, we compute the dynamics associated with the change of the lattice depth in time,

[TABLE]

where the ramp-up and ramp-down speeds $v_{\rm up}$ and $v_{\rm down}$ are given by Eqs. (S3.13) and (S3.14). We start with the superfluid ground state at $V_{0,{\rm ini}}=5\,E_{\rm R}$ and $t_{0}=-100\,{\rm ms}$ while setting $\gamma=0$. We slowly ramp up the optical lattice in $100\,{\rm ms}$ to $V_{0,{\rm max}}=20\,E_{\rm R}$, implying that $t_{1}=0\,{\rm ms}$ and $v_{\rm up}=0.15\,E_{\rm R}/{\rm ms}$, in order to prepare a Mott insulating state. Right after preparing the Mott insulating state, we turn on $\gamma$ to be a finite value and ramp down the optical lattice to $V_{0,{\rm fin}}=5\,E_{\rm R}$ in $7.5\,{\rm ms}$, implying that $t_{2}=7.5\,{\rm ms}$ and $v_{\rm down}=-2\,E_{\rm R}/{\rm ms}$. Notice that in contrast to the dynamics subjected to the hopping ramp-up analyzed above, we do not explicitly include small random noise terms in the initial condition. Instead, the finite-time ramp-up process creates small excitations in the prepared Mott insulating state at $t=t_{1}$, which practically take a role of small initial noise terms needed for dynamically melting the initial Mott state into the “superfluid” state.

In Fig. S12, we show the atom number per site during the ramp-down of the lattice depth, where $t_{1}\leq t<t_{2}$, as a function of the instantaneous value of $V_{0}/E_{\rm R}$. In Fig. S12(A), we see that the onset of the atom loss shifts to the side of large $V_{0}/E_{\rm R}$ when $\gamma$ increases up to $\gamma=0.5$. In contrast, as shown in Fig. S12(B), the onset significantly shifts to the side of small $V_{0}/E_{\rm R}$ when $\gamma$ increases further from $\gamma=0.5$. This means that the melting of the initial Mott insulating state is delayed due to the effect of the strong two-body loss term. As shown in Fig. S13, a similar tendency is also seen in the dynamics of the condensate fraction $|\psi|^{2}/\langle\hat{n}_{\rm A}\rangle$, which qualitatively corresponds to the strength of the coherence peak in the momentum distribution. When $\gamma$ increases from $\gamma=0.5$, the onset of the growth of $|\psi|^{2}/\langle\hat{n}_{\rm A}\rangle$ shifts significantly to the side of small $V_{0}/E_{\rm R}$. Note that the oscillation of the condensate fraction originates from non-adiabaticity of the ramp down of the lattice depth. Because the gap of the amplitude mode is small in the crossover region [41], a relatively fast ramp-down across the crossover excites the amplitude mode. In contrast, such an oscillation is not observed in the experiment likely because of the combined effect of quantum and thermal fluctuations, and the spatial inhomogeneity due to the trap potential. Specifically for the inhomogeneity, the frequency of the amplitude mode significantly depends on the chemical potential, which varies in space in the presence of a trap potential, and this leads to the dephasing of the oscillation.

S4 Unexpectedly large atom loss for strong intensity of PA laser

In this section, we discuss the possible origin of the unexpectedly large atom loss rate $\kappa$ observed for much higher intensity of PA laser corresponding to $\gamma>5$. This additional atom loss which is not taken into account in the present theory prevents the suppression of two-body loss rate from clear observation, and is also observed in the ramp-down dynamics in the deep lattice region. For $V_{0}=20~{}E_{R}$, the measured loss rate is about 30 Hz, while the loss rate expected from the theory is about 3 Hz. We confirm that this loss is not attributed to the photon scattering: the photon scattering rate we measure is about 3 Hz for intensity $I\sim 30$ W/cm2.

Our calculation shows that Raman-assisted tunneling aidelsburger-11 due to the PA laser can explain the observed additional loss. Our high intensity PA laser not only induces the molecular formation but also the coupling between ${}^{1}S_{0}$ and ${}^{3}P_{1}$ atomic states with the detuning of 3.7 GHz, which results in the enhanced atomic tunneling by Raman process even though the overlap of the Wannier functions between the nearest neighboring site is small. From our calculation, for $\gamma\sim 5$ in the lattice depth of $V_{0}=20~{}E_{R}$, this enhancement amounts to the two-body loss rate $\kappa\sim 50$ Hz, which is consistent with observation. To avoid the effect of this additional loss, we restrict the region of the dissipation strength under $\gamma\sim 5$ in this work.