Interacting lattice systems with quantum dissipation: a quantum Monte Carlo study

TL;DR

This paper investigates the effects of quantum dissipation on one-dimensional interacting lattice systems using Quantum Monte Carlo simulations, revealing novel phases and phase transitions induced by dissipation.

Contribution

It introduces a simple thermodynamical approach to model dissipation in quantum lattice systems and explores its impact through numerical simulations.

Findings

Discovery of a compressible zero-temperature state with symmetry breaking

Identification of a thermal phase transition in a 1D dissipative quantum system

Counterintuitive effects of dissipation on quantum many-body states

Abstract

Quantum dissipation arises when a large system can be split in a quantum system and an environment where the energy of the former flows to. Understanding the effect of dissipation on quantum many-body systems is of particular importance due to its potential relations with quantum information processing. We propose a conceptually simple approach to introduce the dissipation into interacting quantum systems in a thermodynamical context, in which every site of a 1d lattice is coupled off-diagonally to its own bath. The interplay between quantum dissipation and interactions gives rise to counterintuitive interpretations such as a compressible zero-temperature state with spontaneous discrete symmetry breaking and a thermal phase transition in a one-dimensional dissipative quantum many-body system as revealed by Quantum Monte Carlo path integral simulations.

Equations

$$H_s=\sum _i-J(a_i^{†}{a}_{i+1}+h.c)+V(n_i-\frac{1}{2})(n_{i+1}-\frac{1}{2}),$$

$$H_b=-\sum _{i,j}{J}_b(b_{i,j}^{†}{b}_{i,j+1}+h.c)$$

$$H_{sb}=-\sum _i{J}^{\prime }(a_i^{†}{b}_{i,0}+h.c).$$

$$U_2=\frac{3}{2}(1-\frac{1}{3}\frac{⟨m^4⟩}{⟨m^2{⟩}^2})$$

$${\xi }_a=\frac{1}{q_1}\sqrt{\frac{S(\pi )}{S(\pi +q_1)}-1}$$

$$Z$$

$$W=⟨i_1|\hat{T}|i_2⟩e^{-({\tau }_2-{\tau }_1)E_{i_2}}⟨i_2|\hat{T}|i_3⟩e^{-({\tau }_3-{\tau }_2)E_{i_3}}\cdots$$

$$\times e^{-({\tau }_n-{\tau }_{n-1})E_{i_n}}⟨i_n|\hat{T}|i_1⟩e^{-(\beta +{\tau }_1-{\tau }_n)E_{i_1}}$$

$${\rho }_s=\frac{L}{\beta }⟨W^2⟩,$$

$$\kappa =\frac{1}{L}\frac{d}{d\mu }N(\mu ){|}_{\mu =0}$$

$$N(\mu )=\frac{1}{Z(\mu )}Tr(\sum _i{\hat{n}}_i)e^{-\beta {\hat{H}}_{tot}+\beta \mu {\sum }_i{\hat{n}}_i}$$

$$\frac{dN(\mu )}{d\mu }$$

$$\chi =\frac{d}{d{\mu }_s}m({\mu }_s){|}_{{\mu }_s=0}=\frac{L}{T}(⟨m^2⟩-⟨m{⟩}^2)$$

$${\chi }^{\prime }=\frac{L}{T}(⟨m^2⟩-⟨|m|{⟩}^2)$$

$${\rho }_s^c(L)={\rho }_s^c(L\to \infty )(1+\frac{1}{2\ln (L)+C}),$$

Full text

I Introduction

Almost all the quantum systems are inevitably coupled to their surroundings. Understanding a quantum system immersed in an environment is not only of immense practical significance in quantum simulation and information processing Schlosshauer (2007), but also in understanding fundamental questions such as the quantum measurement process, quantum-to-classic crossover/transitions, and the mechanism of decoherence Hänggi and Ingold (2006); Hänggi et al. (2008); Caldeira and Leggett (1981); Caldeira and Leggett (1983a, b); Leggett et al. (1987); Weiss (1999); Breuer and Petruccione (2002); Gardiner and Zoller (1999); Leggett et al. (2004). The scenario is further complicated if the open quantum system itself is a (strongly) interacting many-body system, which is indeed the case for the majority of current quantum computing systems including trapped ions, Rydberg atoms, and solid-state quantum computers/simulators. In such systems, the interplay between the quantum many-body effects and the dissipation gives rise to a plethora of novel phenomena Chakravarty et al. (1986); Fisher (1987); Cazalilla et al. (2006); Dalla Torre et al. (2010); Zhu et al. (2015); Werner et al. (2005); Werner and Troyer (2005); Sperstad et al. (2011); Stiansen et al. (2012); Cai et al. (2014); Diehl et al. (2008); F. Verstraete and Cirac (2009); Prosen and Pižorn (2008); Daley et al. (2009); Wilming et al. (2017); Rançon et al. (2013); Schindler et al. (2013); Johnson et al. (2011) in fields as diverse as solid state physics, quantum information and atomic physics.

A variety of theoretical and numerical methods for open quantum systems has been employed. Most of them follow the spirit of taking the environment into account in an exact or approximate way through deriving an effective action by integrating out the bath degrees of freedom Caldeira and Leggett (1981); Caldeira and Leggett (1983a, b); Leggett et al. (1987); Weiss (1999); Breuer and Petruccione (2002); Gardiner and Zoller (1999). Despite its great success, this treatment is difficult – or at least impractical in dealing with the complex situations where the system itself is a quantum many-body system. Except for some special cases Exc (a), often in the field of open quantum systems, tracing out the bath degrees of freedom produces an effective action with unequal-time (retarded) interactions in imaginary or real time. This complicates the quantum many-body system usually to such a degree that numerical methods currently used in strongly correlated physics can no longer be applied Exc (b). A different route to study open quantum many-body systems follows the strategy of quantum optics by generalizing the Born$-$Markov master equation to the many-body case, which is restricted to those open quantum systems with a weak system-bath (SB) coupling to a Markovian environment, neglecting the time delay in the interactions.

In this paper, we propose a conceptually simple approach to study the dissipative interacting quantum system by treating the bath degrees of freedom on the same footing as the system variables, even though we are only interested in the properties of the quantum system. An essential ingredient of this approach is the separation of a global system into a quantum subsystem and the bath, where the choice of the bath Hamiltonian needs to be simple enough to be tractable by conventional many-body numerical methods, yet complicated enough to capture the essential physics of such environments occurring in nature. Motivated by a recent intriguing proposal of modeling the environment by an engineered spin chain Ramos et al. (2016); Vermersch et al. (2016), we perform a numerically exact Quantum Monte Carlo (QMC) path integral simulation with the worm algorithm Prokof’ev et al. (1998) (here in the implementation of Ref. Pollet et al. (2007)) and the stochastic series expansion algorithm Syljuåsen and Sandvik (2002) to study a composed quantum many-body system with a special lattice geometry, which is interpreted as an interacting lattice system with quantum dissipation.

The main points of our results are highlighted as follows. The phase diagram of an interacting lattice system with quantum dissipation is investigated by a numerically exact method. At zero temperature, it is shown that the dissipation can fundamentally alter the properties of the 1D system, and induce a novel state absent as a ground state in typical 1D closed quantum systems. At finite temperature, we find that even though both the dissipation and the thermal fluctuations are individually detrimental to long-range order, their conspiracy can facilitate it and give rise to a finite-T phase transition in this 1D dissipative quantum many-body system. Our results also show that spontaneous symmetry breaking can only occur in a subsystem spatially embedded in a larger system with an inhomogeneous Hamiltonian.

II Model and method

The Hamiltonian describing a dissipative system contains three parts: $H_{tot}=H_{s}+H_{b}+H_{sb}$. For the system Hamiltonian $H_{s}$ we choose a prototypical example of an interacting quantum model: a 1D hard-core boson chain with nearest-neighboring density-density interactions (or, equivalently, an XXZ model in the spin language), which reads:

[TABLE]

where $a_{i}^{\dagger}$ ($a_{i}$) is the creation(annihilation) operator of a hard-core boson on site $i$, $J$ is the tunneling amplitude, and $V>0$ denotes the repulsive interaction strength.

Apart from the intrinsic difficulties in solving strongly correlated physics, a proper modeling and dealing with a realistic environment is also a theoretical challenge. For many realistic quantum systems we do not have a proper understanding of the microscopic origin of dissipation. A desired bath model consists of a quantum system with gapless excitation spectrum and its number of degrees of freedom should be much larger than that of the system it is coupled to. In addition, we devise the baths surrounding different system sites in such a way that they do not influence each other, reflecting the situation in scalable quantum computing with solid-state devices. Finally, the bath model needs to be as simple as possible to be tractable by numerical methods. A minimal bath model satisfying the above requirements is a set of independent chains of hard-core bosons, each coupled to a system site as shown in Fig. 1, with Hamiltonian:

[TABLE]

where $b^{\dagger}_{i,j}$ ($b_{i,j}$) denotes the creation(annihilation) operator of a hard-core boson at site $j$ of the bath linked to the system on site $i$. However, we should emphasize that the bosons in the bath need not have the meaning of physical particles, nor the 1D structure that of the real geometry of the baths in nature. $H_{b}$ in Eq.(2) plays the role of a quantum reservoir with continuum spectrum that can absorb extra energy from the system. In passing we note that the idea of independent baths for composite systems was first introduced in laser theory and later in modeling heat conduction Davies (1978). In mesoscopic physics, Buttiker has proposed the idea of using a local reservoir to study the environment-induced inelastic scattering in quantum transport Büttiker (1985, 1986). However, in none of the above cases have interactions been considered.

Regarding the SB coupling, a system can either interact diagonally with the bath preserving its particle number, or off-diagonally, exchanging both energy and particles simultaneously. The former case has been investigated previously by two of us using a different algorithm Cai et al. (2014). We therefore focus here on the off-diagonal, particle exchange SB coupling with the Hamiltonian:

[TABLE]

We assume that each system site couples only to the central site $j=0$ of the bath chain with coupling strength $J^{\prime}$. In spite of their simplicity, Eq.(2) and Eq.(3) capture some of the essential physics and represent at least part of the physical reality in both a quantum optics and solid state context; e.g., the spectrum of the bath resembles the one of a photonic band gap material, whereas the SB coupling, after a proper transformation, is reminiscent of the the matter-light interaction in quantum optics. In addition, both the SB coupling strength and the bath properties (e.g. the memory time) are highly tunable, which enables us to investigate both Markovian and non-Markovian physics ranging from weak to strong SE coupling regimes in a unified picture.

The most obvious advantage of using Eq.(2) and Eq.(3) to mimic dissipation is that it enables us to treat the bath on the same footing as the system. In the following, we use Quantum Monte Carlo (QMC) simulations to study the composed quantum system in the comb-like geometry, in which the 1D subsystem is regarded as “system”and the remainder as “bath(s)”. Even though we solve the global system, we focus only on the properties of the system. In our simulations, we assume that the composed system (system + bath) is weakly coupled to a superbath with working temperature $T$ and is thus in thermal equilibrium even though the system is entangled with the bath. In all simulations, we focus on half-filling for the composed system. Without loss of generality, we assume that the system and bath chain have the same length $L$. Periodic boundary conditions are used. Both the ground state and the finite-temperature properties of the composed system are investigated.

III Phase diagram at zero temperature: an infinitely compressible insulating phase with density modulations

We first focus on the ground state ($T=0$) of the composed system. In the QMC simulations, we scale the inverse temperature as $\beta=L$, as is relevant for Luttinger liquids, thus setting the dynamical critical exponent $z=1$. The thermodynamic limit is approached in the limit $L\rightarrow\infty$. Without the SB coupling ($J^{\prime}=0$) it is known that 1D hard-core bosons with NN interaction at half-filling undergo a Kosterlitz-Thouless (KT)-type quantum phase transition at $V=2J$ from a Luttinger liquid to a gapped Mott insulator with charge-density-wave (CDW-MI) order. We investigate how the SB coupling changes the above picture. To gain insight, we first discuss qualitatively the effect of the SB coupling on the respective phases. Since a particle which escapes into the $y$-direction must come back to the same site, the bath can be seen as enlarging the unit cell from a single site to an (infinitely long) chain. Hence, we do not expect the 1d nature of the Luttinger liquid to change substantially, but since the system in the $y$-direction is gapless and compressible the charge per unit cell does not need to be quantized.

The first quantity we analyze is the superfluid density, defined as $\rho_{s}=L/\beta\langle W^{2}\rangle$ with $W$ the winding number defined along the $x$-direction Pollock and Ceperley (1987). In Fig.2 (a), we fix the SB coupling $J^{\prime}=0.2J$ and plot $\rho_{s}$ as a function of $V$. A Weber-Minnhagen fit Weber and Minnhagen (1988) shows that $\rho_{s}$ exhibits a sudden drop from a finite value $\rho_{s}^{c}=0.228$ to zero at $V=2.85J$(see the Suppl. Mat. Sup ), indicating that there is still a KT-type quantum phase transition. The non-Luttinger liquid phase on the other side of the transitions is also characterized by spontaneous translational symmetry breaking through the emergence of density-wave (DW) order characterized by the static DW structure factor, $S(\pi)=1/L^{2}\sum_{ij}(-1)^{i-j}\langle(n_{i}-\frac{1}{2})(n_{j}-\frac{1}{2})\rangle$, as shown in Fig 2 (b). Up to now, it seems that the dissipation doesn’t bring any qualitative change to the system except shifting the position of the phase transition point. A striking difference can however be found in the compressibility of the DW phase, defined as $\kappa=\beta/L(\langle N^{2}\rangle-\langle N\rangle^{2})$ with $N$ the total particle number within the system chain (not the total system). As shown in Inset (ii) of Fig. 2 (c), the compressibility $\kappa$ increases linearly with $L$ in the dissipative density wave (DDW) phase, indicating an infinitely compressible state in the thermodynamic limit which makes it fundamentally different from the CDW Mott insulating state at $J^{\prime}=0$. The divergence of the compressibility is due to the fact that every system site is coupled to an infinite number of degrees of freedom; thus, every unit cell can be doped with no energy cost. In passing, we note that anomalously large isochoric compressibilities have in the past been found in totally different contexts such as the superclimb of edge dislocations in solid 4He Ray and Hallock (2008); Söyler et al. (2009). Furthermore, the vanishing of superfluid density in the direction of the DW modulation rules out a supersolid. Such a compressible insulating phase is absent as a ground state in typical closed quantum spin-systems. The SB coupling allows particles to delocalize over a larger system and is expected to enhance the superfluid phase. Indeed, when increasing $J^{\prime}$ at fixed $V>2J$, the system goes over from the CDW-MI to the DDW phase and finally to a superfluid, as can be seen in Fig. 2, panel (c). However, further increasing $J^{\prime}\gg J$ (which no longer corresponds to a physically motivated SB coupling) enhances singlet formation and leads to a strong suppression of $\rho_{s}$ , as can be seen in inset (i), but we expect no phase transition in the thermodynamic limit.

IV Finite temperature phase diagram: a thermodynamic phase transition in a 1D dissipative interacting quantum system

The most striking effect of the SB coupling on the quantum many-body system can be found at finite temperature. On the one hand, it is well known that there is no spontaneous symmetry breaking at any finite temperature for a closed 1D system with local Hamiltonian (while long-range interactions may change this scenario Thouless (1969); Dyson (1969); Kosterlitz (1976); Imbrie and Newman (1988); Spohn (1989); Luijten and Meßingfeld (2001)). The argument goes as follows, for instance for the CDW phase: Thermal fluctuations induce pairs of kink-antikink domain walls which cost only a finite energy. At any finite temperature, the entropy gain by deconfining the excitations overwhelms the energy cost. The leads to a proliferation of domain walls, which destroys the long-range CDW order. Furthermore, quantum fluctuations induced by the SB coupling are also detrimental to the CDW order. On the other hand, the infinite compressibility of the DW phase may alter the domain-wall proliferation argument.

In order to better understand the latter, we switch off the quantum fluctuations in the system ($J=0$), and calculate the DW (staggered) susceptibility of the system, defined as $\chi=\frac{\beta}{L}(\langle m^{2}\rangle-\langle|m|\rangle^{2})$Kotze (2008), with $m=\frac{1}{L}\sum_{i}(-1)^{i}n_{i}$. We plot $\chi$ as a function of $T$ in Fig. 3 (a). Other than the upturn associated with the ground state previously discussed, we observe a peak which keeps growing with system size and whose peak position extrapolates to a finite value of temperature in the thermodynamic limit (see the inset of Fig. 3 (a) for the finite size scaling of the peak position), which is a signature of a finite temperature phase transition. To verify this point, we also calculate the Binder cumulant $U_{2}$ and the correlation length $\xi_{a}$ defined as

[TABLE]

where $S(Q)=1/L^{2}\sum_{ij}e^{iQ(i-j)}\langle(n_{i}-\frac{1}{2})(n_{j}-\frac{1}{2})\rangle$ and $q_{1}=2\pi/L$. The $U_{2}$ and $\xi_{a}/L$ as a functions of $T$ with different system sizes are plotted in Fig. 3 (b) and (c), both of which have exhibited the signature of phase transition.

However, it is difficult to infer the universality class of the transition from the numerical data: the lowest system sizes in Fig. 3 (a) are apparently not in the scaling regime and the values of $\chi$ attained are rather low. Due to the presence of the bath we expect the finite size effects to be more important than for the usual spin or bosonic models with one particle. The finite size effects can also been seen in Fig. 3 (b) and (c), where we don’t find an strictly size independent (common crossing) point. To study the universal class of this phase transition, we need to simulate significantly larger system within the scaling regime. This is difficult in our current QMC simulations with worm-type update due to the special comb-like geometry of the system lattice, especially in the case we studied with $J=0$, where the update may be inefficient even though the model is free from sign problem. Also, due to the limited system sizes we can study and the strong finite size effects, we can’t preclude the possibility of an crossover, where the position of peak in Fig. 3 (a) may shift extremely slowly to zero with increasing system size.

If we start from a system whose ground state has no DW order (e.g. $V=2J^{\prime}$) and increase the temperature, there is no finite temperature phase transition, as shown in Fig. 4 (a). Hence, even though thermal fluctuations and the SB coupling are both detrimental to the DW long-range order, they can conspire to facilitate it. Similar phenomena have observed before in the finite temperature phase transition of a quantum compass model in a square lattice with $\sigma^{z}_{i}\sigma^{z}_{i+e_{x}}$ coupling along the horizontal bonds and $\sigma^{x}_{i}\sigma^{x}_{i+e_{y}}$ coupling along the vertical bonds playing a similar role of quantum fluctuationsWenzel and Janke (2008). The reason behind this is that the total system we studied has a comb-like geometry in 2D. The domain wall excitations are no longer deconfined: due to the SB coupling, separating a kink-antikink pair over some distance $d$ will inevitably perturb the bath chains between them, and cost an energy depending on $d$. In the case of models with long range interactions Thouless (1969); Dyson (1969); Kosterlitz (1976); Imbrie and Newman (1988); Spohn (1989); Luijten and Meßingfeld (2001)), the interactions between a kink-antikink pair are also confined, which originates from the long-range nature of the interactions, instead of the system-bath couplings as in our case. In the previous discussion in this section, we set $J=0$ thus all the quantum fluctuations come from the off-diagonal SB coupling terms. To complete our discussion, we switch on the single particle hopping in the system chain, and study the finite temperature phase diagram of the model in the $T-J$ plane with a fixed $V=4J^{\prime}$, the result is shown in the Fig. 4 (b).

V Discussion

While it is widely believed that dissipation leads to dephasing which drives a quantum system towards a classical one, here we show an example that the quantum dissipation induced by an off-diagonal SB coupling to a quantum environment, can enhance the quantum fluctuations in the system and give rise to counterintuitive phenomena. One of the most important questions is to what extent real dissipation can be modeled by the specific choice of the bath and SB coupling used in this paper, which provides the simplest example of gapless quantum modes coupling locally and off-diagonally to the system. Usually, the randomness feature of the environment can be absorbed into the spectrum function of the bath, while a proper revision of the bath Hamiltonian in our model can mimic more general dissipation with various spectral functions. However, our model doesn’t apply for those open quantum many-body systems with non-local dissipation where the entire system shares the same environment (as e.g. in cavity QEDCormick et al. (2013)). Also, the open quantum systems coupled to a classical environment e.g. a classical noise involves the real time evolutionsPoletti et al. (2012); Cai and Barthel (2013); Sieberer et al. (2013); Marino and Diehl (2016), thus are beyond the scope of our current scheme.

VI Conclusion and Outlook

We propose a conceptually simple approach to deal with quantum dissipation in open interacting quantum systems, which allows us to investigate both Markovian and non-Markovian physics ranging from weak to strong SB coupling regimes in a unified picture. We explored the properties of an embedded 1D quantum many-body system by numerically solving a special comb-like lattice geometry. Counterintuitive zero-temperature and thermal behavior has been discovered, indicating that dissipation can fundamentally alter the way we look at the properties of quantum many-body systems: spontaneous symmetry breaking can occur in a subsystem of reduced dimensionality spatially embedded in a larger system with an inhomogeneous Hamiltonian. Some avenues for further work immediately suggest themselves. First, a generalization to higher dimensions is straightforward. Notice that both the bath and the SB coupling term lead to positive definite expansions in the path integral formulation. Therefore, as long as the system Hamiltonian is free from the sign problem, the total Hamiltonian can be solved in a numerically exact way by QMC. This allows one in principle to investigate the effect of quantum dissipation on various systems with long-range order (symmetry breaking) and even topological order, but also on how quantum coherence and entanglement properties may change. Even though in this paper we only study the equilibrium state (for the total system), it would be very interesting to explore the non-equilibrium dynamics of the dissipative quantum many-body system in the current scheme, which is beyond the scope of QMC simulations but may be accessible by other methods (e.g. an exactly solvable quadratic fermionic model defined on a similar comb-like lattice). The far-from-equilibrium dynamical and steady state properties of a driven-dissipative quantum many-body system remain an open and elusive questionHartmann et al. (2016).

Appendix A A Brief introduction of the quantum Monte Carlo method with worm algorithm

We consider a Hamiltonian $\hat{H}$ which can be decomposed as $\hat{H}=-\hat{T}+\hat{V}$ , where the $\hat{T}$ represents the off-diagonal terms in the Hamiltonian and $V$ is the diagonal ones in the Fock basis $|i\rangle=|n_{1}n_{2}\cdots n_{L}\rangle$. The partition function $Z=\mathbf{Tr}e^{-\beta\hat{H}}$ can be expanded in terms of the probability functions of space-imaginary time configurations as

[TABLE]

where $\tau_{n}$ is the time of the $n$th tunneling event in the world-line of the particles, and $|i_{m}\rangle$ denotes the particle configurations between the imaginary time $\tau_{m-1}$ and $\tau_{m}$ (as shown in Fig.5), and for a given space-imaginary time configuration: $\{\tau_{1},\cdots,\tau_{n},|i_{1}\rangle,\cdots,|i_{n}\rangle\}$, the corresponding probability $W(\tau_{1},\cdots,\tau_{n},|i_{1}\rangle,\cdots,|i_{n}\rangle)$ can be written as:

[TABLE]

where $E_{i_{m}}=\langle i_{m}|\hat{V}|i_{m}\rangle$ is the interaction (diagonal) energy for particle configuration between $\tau_{m}$ and $\tau_{m-1}$. As long as $\hat{T}$ is a positive definite operator, we can prove that $W(\tau_{1},\cdots,\tau_{n},|i_{1}\rangle,\cdots,|i_{n}\rangle)$ is always positive, which enables us to perform the importance sampling and evaluate the average value of physical quantities over a limited number of space-(imaginary) time configurations. A worm algorithm is a update algorithm where the partition function in Eq.(6) is sampled indirectly in the extended configuration space of open world-line configurations by performing local movementProkof’ev et al. (1998).

Appendix B Physical quantities obtained by quantum Monte Carlo (QMC)

In this section, we will discuss the derivations of several physical quantities from the QMC simulations. We should emphasize that all these quantities are only defined for the system instead of the system+bath.

B.1 Superfluid density

The superfluid density $\rho_{s}$ is defined as

[TABLE]

where $W$ is the winding number along the system chain, which can be obtained in the quantum Monte Carlo (QMC) simulationsPollock and Ceperley (1987). Considering a world line of a hard-core boson in the imaginary time-space configurations, even though a system boson can escape into the bath from certain system site, it will finally return to the system via the same system site, since each bath chain is independent and a world line should be a closed curve.

B.2 Compressibility

The system compressibility $\kappa$ is defined as the system particle number $N=\sum_{i}\langle n_{i}\rangle$ in response to the perturbation $H^{\prime}=-\sum_{i}\mu n_{i}$, which only operates in the system sites:

[TABLE]

We assume that the bosons in the bath chain don’t feel the chemical potential. As we demonstrated in the main text, the bosons in the bath chains don’t necessarily mean the real bosonic particles as in the system, they can be any bosonic energy reservoir, e.g. phonon, photon and so on, which is insensitive to the external fields operating on the system bosons. We can prove that $\kappa=\frac{\beta}{L}\Delta N$, where $\beta=1/T$ is the inverse temperature, $L$ the system size, and $\Delta N=\langle\hat{N}^{2}\rangle-\langle\hat{N}\rangle^{2}$ the variance of the system particle number. Notice that

[TABLE]

with $Z(\mu)=Tre^{-\beta\hat{H}_{tot}+\beta\mu\sum_{i}\hat{n}_{i}}$, thus

[TABLE]

In the numerical simulation, the particle number fluctuations can be derived from the histogram of the system particle number distributions in the QMC simulations, which is assumed to be a Gaussian distribution for sufficient large system. We fit the distribution of the particle number by a Gaussian distribution as $P(N)=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(N-\bar{N})^{2}}{2\sigma}}$ and the system particle number fluctuation is the width of the continuous Gaussian distribution (in the thermodynamic limit), $\Delta N=\sigma$.

B.3 Charge-density-wave susceptibility

The charge-density-wave (CDW) susceptibility of the system is defined as the static linear susceptibility of the CDW order parameter $m=\frac{1}{L}\langle\sum(-1)^{i}n_{i}\rangle$ in response to the perturbation $H^{\prime}=-\mu_{s}\sum_{i}(-1)^{i}n_{i}$ only operating on the system sites :

[TABLE]

For any finite system size $\langle m\rangle=0$, thus in the QMC simulation, we use a revised CDW susceptibility $\chi^{\prime}$ defined as:

[TABLE]

which agrees with $\chi$ the in thermodynamic limit.

Appendix C The Koserlitz-Thouless transition

At the Kosterlitz-Thouless (KT) phase transition point, the superfluid density flows as a function of $L$ as Weber and Minnhagen (1988):

[TABLE]

where $\rho_{s}^{c}(L\rightarrow\infty)$ is the critical value of the superfluid density at the transition point in the thermodynamic limit, and $C$ is a non-universal constant. To determine the $\rho_{s}^{c}$ and the critical point $V_{c}$, we use the Weber-Minnhagen fit of Eq(14) to obtain the superfluid density calculated in our QMC simulations, $\rho_{s}^{\rm QMC}(L)$, and calculate the average error of the fit, $\delta=\frac{1}{4}\sum_{L}|\rho_{s}^{\rm fit}(L)-\rho_{s}^{\rm QMC}(L)|$ with $L=32,48,64,96$, for different interaction strengths $V$. Since the flow equation Eq.(14) is only valid at the phase transition point, we expect that the the average error will reach it minimum at the critical point $V=V_{c}$, as shown Fig.6.

*Acknowledgements – * ZC wishes to thank P. Zoller for fruitful discussions. ZC acknowledges the support from the National Key Research and Development Program of China (grant No. 2016YFA0302001), the National Natural Science Foundation of China under Grant No.11674221 and the Shanghai Rising-Star Program. LP is supported by FP7/ERC starting grant No. 306897 and acknowledges funding by DFG through NIM-2. X. Wang was supported by the National Program on Key Research Project (Grants No. 2016YFA0300501), and by the National Natural Science Foundation of China (Grant No. 11574200). ZY, JL and YC are supported by the National Natural Science Foundation of China (Grant No. 11625416 and 11474064) and the State Key Programs of China (Grant No. 2016YFA0300504). We acknowledge the support from the Center for High Performance Computing of Shanghai Jiao Tong University.