Thermalization in open many-body systems based on eigenstate thermalization hypothesis

TL;DR

This paper explores how the eigenstate thermalization hypothesis (ETH) can predict the steady states of open quantum systems under weak dissipation, highlighting differences between bulk and boundary dissipation through theoretical and numerical analysis.

Contribution

It establishes a criterion linking ETH, dissipation strength, and the validity of perturbation theory in describing steady states of open quantum systems.

Findings

Gibbs state at an effective temperature describes steady states under ETH.

Perturbation theory validity depends on system size and dissipation type.

Numerical results confirm theoretical predictions about steady state behavior.

Abstract

We investigate steady states of macroscopic quantum systems under dissipation not obeying the detailed balance condition. We argue that the Gibbs state at an effective temperature gives a good description of the steady state provided that the system Hamiltonian obeys the eigenstate thermalization hypothesis (ETH) and the perturbation theory in the weak system-environment coupling is valid in the thermodynamic limit. We derive a criterion to guarantee the validity of the perturbation theory, which is satisfied in the thermodynamic limit for sufficiently weak dissipation when the Liouvillian is gapped for bulk-dissipated systems, while the perturbation theory breaks down in boundary-dissipated chaotic systems due to the presence of diffusive transports. We numerically confirm these theoretical predictions. This work suggests a connection between steady states of macroscopic open quantum systems and the ETH.

Equations

$$\{\begin{array}{rl}\frac{d\rho }{dt}& =\mathcal{L}\rho =-i[H,\rho ]+\gamma \mathcal{D}\rho ,\\ \mathcal{D}\rho & =\sum _{a=1}^{V_{\mathrm{D}}}\left(L_a\rho L_a^{†}-\frac{1}{2}\{L_a^{†}{L}_a,\rho \}\right),\end{array}$$

$$O_{nm}={⟨O⟩}_{{\beta }_n}{\delta }_{nm}+\frac{r_{nm}}{\sqrt{D_{{\bar{E}}_{nm}}}}{f}_O({\bar{E}}_{nm},{\omega }_{nm}),$$

$${\rho }_{\mathrm{s}}=\sum _{n=0}^{\infty }{\gamma }^n{\rho }_{\mathrm{s}}^{(n)}.$$

$$[H,{\rho }_{\mathrm{s}}^{(0)}]=0$$

$$-i[H,{\rho }_{\mathrm{s}}^{(1)}]+\sum _a\left[L_a{\rho }_{\mathrm{s}}^{(0)}{L}_a^{†}-\frac{1}{2}\{L_a^{†}{L}_a,{\rho }_{\mathrm{s}}^{(0)}\}\right]=0$$

$${\rho }_{\mathrm{s}}^{(0)}=\sum _n{P}_n|n⟩⟨n|.$$

$$\sum _m\left(W_{nm}{P}_m-W_{mn}{P}_n\right)=0,$$

$$\frac{W_{nm}}{W_{mn}}=e^{-\beta (E_n-E_m)}$$

$$\mathrm{Tr}O{\rho }_{\mathrm{s}}^{(0)}\approx {⟨O⟩}_{\beta }.$$

$$C(\beta )\equiv \sum _a{⟨[L_a^{†},H]L_a⟩}_{\beta }=0.$$

$$\underset{\gamma \to 0}{\lim }\underset{V\to \infty }{\lim }\mathrm{Tr}O{\rho }_{\mathrm{s}}=\underset{V\to \infty }{\lim }\underset{\gamma \to 0}{\lim }\mathrm{Tr}O{\rho }_{\mathrm{s}}=\underset{V\to \infty }{\lim }\mathrm{Tr}O{\rho }_{\mathrm{s}}^{(0)}.$$

$$s\equiv \frac{1}{V}S({\rho }_{\mathrm{s}}\Vert {\rho }_{\mathrm{s}}^{(0)})=\frac{1}{V}\mathrm{Tr}\left[{\rho }_{\mathrm{s}}(\ln {\rho }_{\mathrm{s}}-\ln {\rho }_{\mathrm{s}}^{(0)})\right].$$

$$\delta \rho =-{\tilde{\mathcal{L}}}^{-1}\mathcal{L}{\rho }_{\mathrm{s}}^{(0)}=-\gamma {\tilde{\mathcal{L}}}^{-1}\mathcal{D}{\rho }_{\mathrm{s}}^{(0)},$$

$$\mathcal{D}{\rho }_{\mathrm{s}}^{(0)}=$$

$$-\frac{1}{2}⟨n⟩L_a^{†}{L}_{a}m(P_n+P_m)]|n⟩⟨m|.$$

$$\mathcal{D}{\rho }_{\mathrm{s}}^{(0)}\sim \sum _{n\ne m,|E_n-E_m|\lesssim \Lambda }\frac{V_{\mathrm{D}}^{1/2}{R}_{nm}}{D^{3/2}}|n⟩⟨m|,$$

$$\Vert \delta \rho {\Vert }_F$$

$$\sim \frac{V_{\mathrm{D}}^{1/2}}{\Delta D^{3/2}}\sqrt{\sum _{n,m:|E_n-E_m|\lesssim \gamma \Delta }1}\sim \sqrt{\frac{V_{\mathrm{D}}}{D}\frac{\beta \gamma }{\Delta }}.$$

$$\left|\delta M\right|\le \Vert M{\Vert }_F\Vert \delta \rho {\Vert }_F\sim \sqrt{\frac{V_{\mathrm{D}}}{V}\frac{\beta \gamma }{\Delta }}.$$

$$s\sim \frac{V_{\mathrm{D}}}{V}\frac{\beta \gamma }{\Delta }.$$

$$\frac{V_{\mathrm{D}}}{V}\frac{\beta \gamma }{\Delta }\ll 1.$$

$$\{\begin{array}{rl}& H=\sum _{i=1}^V(h^z{S}_i^z+h^x{S}_i^x+g{S}_i^z{S}_{i+1}^z),\\ & L_i=S_i^{-}(i=1,\cdots ,V),\end{array}$$

$$H=$$

$$+\sum _{i=1}^{V-2}\left[-h^{\prime }(b_i{b}_{i+2}^{†}+b_i^{†}{b}_{i+2})+J^{\prime }\left(n_i-\frac{1}{2}\right)\left(n_{i+2}-\frac{1}{2}\right)\right],$$

$$L_1=\sqrt{1+\mu }{b}_1^{†},L_2=\sqrt{1-\mu }{b}_1,$$

$$L_3=\sqrt{1-\mu }{b}_V^{†},L_4=\sqrt{1+\mu }{b}_V,$$

$$H_{\mathrm{S}}(t)=\sum _{i=1}^V\left((h^z+\omega )S_i^z+\frac{h^x}{2}(S_i^{+}{e}^{-i\omega t}+S_i^{-}{e}^{i\omega t})+g{S}_i^z{S}_{i+1}^z\right),$$

$$H_{\mathrm{T}}(t)=H_{\mathrm{S}}(t)+\lambda H_{\mathrm{I}}+H_{\mathrm{B}},$$

$$\{\begin{array}{rl}{H}_{\mathrm{B}}=& \sum _{i=1}^V\sum _{\alpha }{\omega }_{\alpha }{a}_{i,\alpha }^{†}{a}_{i,\alpha },\\ H_{\mathrm{I}}=& \sum _{i=1}^V\sum _{\alpha }(S_i^{+}{a}_{i,\alpha }+S_i^{-}{a}_{i,\alpha }^{†}),\end{array}$$

$$U(t)=e^{-i\omega t{\sum }_{i=1}^V(S_i^z+{\sum }_{\alpha }{a}_{i,\alpha }^{†}{a}_{i,\alpha })}.$$

$$H_{\mathrm{T}}^{\mathrm{R}}=U^{†}(t)\left(H_{\mathrm{T}}(t)-i\frac{\partial }{\partial t}\right)U(t)=H+\lambda H_{\mathrm{I}}+H_{\mathrm{B}}^{\mathrm{R}},$$

$$H_{\mathrm{B}}^{\mathrm{R}}=\sum _{i=1}^V\sum _{\alpha }({\omega }_{\alpha }-\omega )a_{i,\alpha }^{†}{a}_{i,\alpha },$$

Full text

I Introduction

A quantum system that is weakly coupled to a large environment usually relaxes to a steady state due to dissipation Breuer and Petruccione (2002); Weiss (2012); Spohn (1980). When the environment is in thermal equilibrium, the steady state is universally described by the Gibbs state by virtue of the detailed balance condition Spohn and Lebowitz (1978). In contrast, when the environment is out of equilibrium, the detailed balance condition is violated and there is no simple criterion to determine the steady state. It is a challenge in statistical physics to predict the steady state in such a nonequilibrium situation Mitra et al. (2006); Prosen and Pižorn (2008); Prosen (2011a, b); Dhar et al. (2012); Sieberer et al. (2013); Torre et al. (2013); Sieberer et al. (2016); Foss-Feig et al. (2017). Recent experimental progress using ultracold atoms and trapped ions has enabled us to introduce controlled dissipation Barreiro et al. (2011); Barontini et al. (2013); Rauer et al. (2016); Tomita et al. (2017), which leads us to the possibility of designing dissipation so that the steady state has desired properties Diehl et al. (2008); Verstraete et al. (2009). This experimental background also motivates us to theoretically study steady states of quantum systems under dissipation not obeying the detailed balance condition.

In the present paper, the steady state of a quantum many-body system in a weak contact with an out-of-equilibrium environment is investigated. It turns out that the Gibbs state at a certain effective temperature well describes the steady state in some open quantum systems despite the violation of the detailed balance condition. We theoretically argue that there are two ingredients in the realization of a Gibbs steady state, i.e., the validity of the perturbation theory in weak dissipation and the eigenstate thermalization hypothesis (ETH), which is recognized as an important property of the Hamiltonian in explaining the approach to thermal equilibrium in isolated quantum systems Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Mori et al. (2018).

The ETH consists of two parts; the diagonal ETH and the off-diagonal ETH. The diagonal ETH states that the diagonal elements of a local operator in the energy basis are smooth functions of the energy. The off-diagonal ETH states that off-diagonal elements are irregularly fluctuating and exponentially small with respect to the system size. The ETH is expected to hold in a wide class of nonintegrable Hamiltonians.

On the other hand, the validity of the weak-dissipation perturbation theory in macroscopic systems is highly nontrivial because it is known that its convergence radius quickly shrinks in the thermodynamic limit Žnidarič (2015); Lemos and Prosen (2017). Here, we derive a theoretical criterion [Eq. (19) below] of the validity of the perturbation theory by using the off-diagonal ETH, and then the diagonal ETH ensures that the steady state is well described by a Gibbs state if this criterion is satisfied.

In this way, this paper suggests a connection between steady states of macroscopic open quantum systems and the ETH. In a recent work Ashida et al. (2018), it was shown that a Gibbs state with a time-dependent temperature emerges in the transient dynamics of an open quantum system in which the system of interest is finite and obeys the ETH. Our result should be distinguished from this recent result since we here focus on the steady state (i.e., the long-time limit) in a macroscopic open quantum system (i.e., the thermodynamic limit).

This paper is organized as follows. Section II gives a brief review of the weak-dissipation perturbation theory of the nonequilibrium steady state and the ETH. In Sec III, we derive a theoretical criterion of the validity of the perturbation theory, which clarifies under what condition the ETH is relevant in characterizing the nonequilibrium steady state. In Sec IV, we numerically confirm our theoretical prediction in a bulk-dissipated system and a boundary-dissipated system. In Sec V, this paper is summarized with some future directions.

II Perturbative expansion and the ETH

We consider a macroscopic system of the volume $V$ on $d$-dimensional cubic lattice that is in contact with an environmental system. We denote by $H$ and $\rho$ the Hamiltonian and the reduced density matrix of the system of interest, respectively. The dynamics of $\rho$ is assumed to be described by the Lindblad equation Lindblad (1976); Breuer and Petruccione (2002)

[TABLE]

where $[\cdot,\cdot]$ and $\{\cdot,\cdot\}$ denote the commutator and the anti-commutator, respectively, and we put $\hbar=1$. The dissipator $\mathcal{D}$ is characterized by the $V_{\rm D}$ Lindblad operators $\{L_{a}\}_{a=1,2,\dots,V_{\rm D}}$. When all the sites are subject to dissipation, $V_{\rm D}\propto V$, while dissipation acts only at the boundaries, $V_{\rm D}\propto V^{(d-1)/d}$. The superoperator $\mathcal{L}$ is referred to as the Liouvillian. Here, we consider the weak dissipation regime, i.e., small $\gamma$.

We assume that $H=\sum_{n}E_{n}\ket{n}\bra{n}$ obeys the ETH. The ETH states that matrix elements $O_{nm}=\braket{n}{O}{m}$ of any local operator $O$ takes the following form with $\bar{E}_{nm}=(E_{n}+E_{m})/2$ and $\omega_{nm}=E_{n}-E_{m}$:

[TABLE]

where $f_{O}(\bar{E}_{nm},\omega_{nm})$ is a smooth function of the argument and decays exponentially fast for $|\omega_{nm}|\gtrsim\Lambda$ with a cutoff frequency $\Lambda$, which is independent of $V$. The equilibrium expectation value of $O$ at the inverse temperature $\beta_{n}$ is denoted by $\braket{O}_{\beta_{n}}$, where $\beta_{n}$ is determined by the condition $E_{n}=\braket{H}_{\beta_{n}}$. The quantity $D_{E}$ is the number of energy eigenstates with eigenvalues between $E-\Delta E$ and $E$ with some width $\Delta E$ 111Since off-diagonal elements decay for $|E_{n}-E_{m}|\gtrsim\Lambda$, we should choose $\Delta E\sim\Lambda$. Otherwise, the function $f_{O}$ will diverge or vanish in the thermodynamic limit., and it scales as $D_{E}\sim e^{O(V)}$ with volume $V$. The last term of Eq. (2) expresses a small fluctuating part, and $\{r_{nm}\}$ behave as if their real and imaginary parts were random variables of mean zero and variance unity.

The steady state $\rho_{\rm s}$ is defined by $\mathcal{L}\rho_{\rm s}=0$. Since $\gamma$ is assumed to be small, we perform the perturbative expansion of $\rho_{\rm s}$ in $\gamma$:

[TABLE]

By substituting this expression into $\mathcal{L}\rho_{\rm s}$ and requiring that it vanishes in each order in $\gamma$, we obtain

[TABLE]

for $O(\gamma^{0})$, and

[TABLE]

for $O(\gamma^{1})$. Equation (4) implies that $\rho_{\rm s}^{(0)}$ is diagonal in the energy basis, i.e.

[TABLE]

By looking at the $n$th diagonal element of Eq. (5) in the energy basis, we obtain

[TABLE]

where $W_{nm}=\sum_{a}|\braket{n}{L_{a}}{m}|^{2}$ Thingna et al. (2013). Equation (7) determines the diagonal elements of $\rho_{\rm s}^{(0)}$, i.e., $\{P_{n}\}$. It is noted that $W_{nm}$ can be interpreted as the transition rate from the state $m$ to $n$. The transition rates satisfy the detailed balance condition

[TABLE]

when the environment is in thermal equilibrium at the inverse temperature $\beta$ Davies (1976, 1974); Spohn and Lebowitz (1978). As a result, the steady state is given by the Gibbs state $P_{n}=e^{-\beta E_{n}}/Z(\beta)$ with the partition function $Z(\beta)=\sum_{n}e^{-\beta E_{n}}$.

When the environment is out of equilibrium, the detailed balance condition is violated, and hence $\{P_{n}\}$ is not necessarily of the Gibbs form. Nevertheless, $\rho_{\rm s}^{(0)}$ is indistinguishable from the Gibbs state if the system Hamiltonian $H$ obeys the ETH. Since $\beta_{n}$ in Eq. (2) is almost constant, $\beta_{n}\approx\beta$, as long as the energy fluctuation in $\rho_{\rm s}^{(0)}$ is subextensive, we have

[TABLE]

In Appendix C, we show that the inverse effective temperature $\beta$ can be determined by numerically solving the following equation,

[TABLE]

Since $C(\beta)$ only depends on equilibrium values of local operators $[L_{a}^{\dagger},H]L_{a}$, numerical methods for thermal equilibrium states, e.g. quantum Monte Carlo method, can be used to calculate $C(\beta)$. In this way, the steady state is well described by the Gibbs state for small $\gamma$ as long as the naive perturbation theory is valid 222The effective temperature introduced here is operator independent in contrast to previous studies Sieberer et al. (2013, 2014).

III Validity of perturbation theory

In Ref. Lemos and Prosen (2017), it is numerically shown that the convergence radius of the perturbative expansion of Eq. (3) shrinks to zero in the thermodynamic limit, $V\to\infty$. This means that it is a nontrivial issue whether the thermodynamic limit commutes with the weak-dissipation limit. If they are commutable in evaluating the expectation value of an operator $O$, we have

[TABLE]

For macroscopic systems, the thermodynamic limit should be taken before the weak-dissipation limit, and hence, the left-hand side of Eq. (11) is the quantity we want. On the other hand, the right-hand side of Eq. (11) corresponds to the solution in the leading-order perturbation theory.

In the present paper, we discuss whether Eq. (11) holds by investigating the relative entropy density

[TABLE]

If $\lim_{\gamma\to 0}\lim_{V\to\infty}s=0$, we can conclude the macrostate equivalence between $\rho_{\rm s}$ and $\rho_{\rm s}^{(0)}$ Touchette (2015), i.e., Eq. (11) holds for intensive macroscopic observables $O$ that obey the large-deviation principle in the steady state Mori (2016).

Now we derive a criterion for the validity of the perturbation theory, i.e. Eq. (11). We assume the open boundary condition for simplicity, but we can also derive the identical result for the periodic boundary condition. The exact steady state $\rho_{\rm s}$ is decomposed as $\rho_{\rm s}=\rho_{\rm s}^{(0)}+\delta\rho$. Then, the equality $\mathcal{L}\rho_{\rm s}=0$ is rewritten as

[TABLE]

where $\tilde{\mathcal{L}}^{-1}$ is the pseudo-inverse of $\mathcal{L}$. From the definition of $\mathcal{D}$ and Eq. (6), we obtain

[TABLE]

Since the Lindblad operators are assumed to be local, they obey the ETH [see Eq. (2)]. We consider an energy shell that consists of the energy eigenstates $\{\ket{n}\}$ with $E_{n}\in[\bar{E}-\delta E/2,\bar{E}+\delta E/2]$. Here, $\bar{E}=\mathrm{Tr}\,H\rho_{\rm s}^{(0)}$ and $\delta E$ is chosen so that it is macroscopically small $\delta E=o(V)$ but large enough to ensure $\sum_{n:E_{n}\in[\bar{E}-\delta E/2,\bar{E}+\delta E/2]}P_{n}\approx 1$ (typically, $\delta E\sim V^{1/2}$) 333The difference between $\delta E\sim V^{1/2}$ here and $\Delta E\sim\Lambda$ in the ETH (2) does not matter unless $\beta^{-1}\gg\Lambda$. For simplicity, below we do not take care of the difference between $\Delta E$ and $\delta E$. In more careful analysis, it turns out that $\beta$ in Eqs. (17), (18), and (19) should be replaced by $\max\{\beta,1/\Lambda\}$.. Then, each $P_{n}$ is roughly equal to $1/D$, where $D$ is the number of eigenstates within the energy shell. We can then evaluate the order of magnitude of $\mathcal{D}\rho_{\rm s}^{(0)}$ as

[TABLE]

where $\{R_{nm}\}$ behave as random variables of mean zero and variance of $O(1)$, and $\Lambda$ is a cutoff frequency that is independent of $V$.

Next, we multiply $\tilde{\mathcal{L}}^{-1}$. Let us define the gap $\gamma\Delta$ of the Liouvillian as the nonzero smallest absolute value of the eigenvalues of $\mathcal{L}$. Then, the dominant contribution in multiplying $\tilde{\mathcal{L}}^{-1}$ comes from eigenmodes with eigenvalues close to $\gamma\Delta$, so we only consider such slow eigenmodes. It is expected that matrix elements corresponding to fast oscillations, i.e. $|E_{n}-E_{m}|\gtrsim\gamma\Delta$, do not contribute to slow eigenmodes near the gap $\gamma\Delta$, and hence we can roughly evaluate the Frobenius norm of $\delta\rho$ by using Eqs. (13) and (15) as follows:

[TABLE]

Here, let us consider a macroscopic quantity $M=(1/V)\sum_{i=1}^{V}O_{i}$, where $i$ is an index of the lattice sites and $O_{i}$ is a local operator acting to sites near $i$. In an energy shell, the diagonal elements of $M$ are roughly constant due to the ETH, so without loss of generality we put $\braket{M}_{\beta}=0$. Then, the Frobenius norm of $M$ within the energy shell is evaluated as $\|M\|_{F}\approx\sqrt{D/V}$. As a result, $\delta M:=\mathrm{Tr}\,M\delta\rho$ is evaluated as follows:

[TABLE]

Since the relative entropy density $s=S(\rho_{\rm s}\|\rho_{\rm s}^{(0)})/V$ is related to $\delta M$ by $s\gtrsim|\delta M|^{2}$ Mori (2016), Eq. (17) implies

[TABLE]

Equation (17) or Eq. (18) gives the following criterion for the validity of the perturbation theory:

[TABLE]

This is a main result of our paper. In the case of bulk dissipation, i.e. $V_{\rm D}\propto V$, as long as the Liouvillian is gapped $\Delta>0$ in the thermodynamic limit, the criterion (19) is satisfied for a small but finite $\gamma$ in the thermodynamic limit. The Liouvillian is expected to be gapped for a wide class of nonintegrable systems under bulk dissipation with no conserved quantity, and hence Eq. (19) implies that the steady state is described by a Gibbs state for an equally wide class of open systems. On the other hand, in the case of boundary dissipation, $V_{\rm D}\propto V^{(d-1)/d}$ and the criterion reads $\beta\gamma/(V^{1/d}\Delta)\ll 1$. The Liouvillian gap in a boundary-dissipated chaotic system typically behaves as $\Delta\sim V^{-\theta/d}$ with an exponent $\theta\geq 1$ Žnidarič (2015). In many cases $\theta>1$, and then our theory predicts that the perturbation theory may break down for an arbitrarily small $\gamma$ in the thermodynamic limit 444One might expect $\theta=2$ for chaotic systems because of the presence of diffusive transport, but it is not always the case Žnidarič (2015).. Below, theoretical predictions discussed here will be numerically confirmed.

IV Numerical result

IV.1 Bulk-dissipated system

We first study a nonequilibrium system under bulk dissipation with no conserved quantity. We consider the following dissipative Ising chain under the periodic boundary condition:

[TABLE]

where $\{\vec{S}_{i}\}_{i=1}^{V}$ are spin-1/2 operators and $S_{i}^{\pm}\equiv S_{i}^{x}\pm iS_{i}^{y}$. The parameters of the Hamiltonian are set as $(h^{z},h^{x},g)=(1.809,1.618,4)$, with which the ETH has been numerically shown to hold Kim et al. (2014). This open quantum system has been implemented using Rydberg atoms Carr et al. (2013); Letscher et al. (2017) and nonequilibrium phase transitions have been theoretically discussed Lee et al. (2011). In this system the up and down spin states correspond to the Rydberg state and the ground state of an atom, respectively, and $\{L_{i}\}$ describes the spontaneous emission in each atom. It is noted that the detailed balance condition is not satisfied in this model. In Appendix, the microscopic derivation of the Lindblad equation is given (Appendix B) and the violation of the detailed balance condition is demonstrated (Appendix D).

The system-size dependences of $s$ at $\gamma=0.01,0.03$ and $0.05$ are shown in Fig. 1 (a). We find the linear dependence of the distances on $1/V$ for large $V$ ($V\geq 9$). By using this linear dependence, we extrapolate the data to the thermodynamic limit, and obtain $\lim_{V\to\infty}s$ for several small values of $\gamma$ [Fig. 1 (b)]. We find $\lim_{V\to\infty}s\propto\gamma$ for $\gamma\lesssim 0.02$, which is consistent with Eq. (18) with $V_{\rm D}=V$. This conclusion is independent of the choice of the parameters $(h^{z},h^{x},g)$ as far as we have calculated. In Appendix E, we provide another example showing the same $\gamma$-dependence of $s$.

IV.2 Boundary-dissipated system

Next, we discuss a nonequilibrium system with boundary dissipation. As we have already argued, our criterion (19) tells us that the perturbation theory would break down and the steady state is not described by a Gibbs state in a boundary-dissipated system. In order to understand this result more intuitively, suppose a one-dimensional system in contact with two particle reservoirs with different chemical potentials at each end. The chemical potential difference drives the system, and particles will flow diffusively in the bulk. Such diffusive transports result in a gradient in the particle density profile. On the other hand, if the system Hamiltonian possesses translation invariance in the bulk, an individual energy eigenstate shows a uniform density profile, and hence its mixture like $\rho_{\rm s}^{(0)}$ cannot reproduce the expected gradient of the density profile in the steady state. This argument can be generalized to other conserved currents (e.g., an energy current between two thermal reservoirs at different temperatures). The perturbation theory fails in such a situation.

The failure of the perturbation theory is demonstrated for the hard-core Bose-Hubbard model driven by two environments with different chemical potentials:

[TABLE]

and

[TABLE]

where $b_{i}$ and $b_{i}^{\dagger}$ are annihilation and creation operators of a boson at site $i$, and $n_{i}=b_{i}^{\dagger}b_{i}=0$ or $1$. The parameters of the Hamiltonian are given by $(h,h^{\prime},J,J^{\prime})=(0.9167,0.2449,4,0.9045)$. The Lindblad operators $\{L_{a}\}_{a=1}^{4}$ act on the boundaries of the lattice and $\mu$ effectively controls the chemical potential of the environments. We set $\mu=0.1$. In the steady state, we have a nonuniform particle density profile.

The system-size dependences of $s$ at $\gamma=0.01,0.03$ and $0.05$ are shown in Fig. 2 (a). Again, we find linear dependence on $1/V$ for $V\geq 9$, so the value in the thermodynamic limit is estimated by extrapolating the data. In this way, we obtain the $\gamma$-dependence of $\lim_{V\to\infty}s$ [Fig. 2 (b)], showing that the distance is finite in the limit of $\gamma\to 0$: $\lim_{\gamma\to 0}\lim_{V\to\infty}s\neq 0$. This result clearly shows the failure of the perturbation theory in a boundary-dissipated system.

V Summary

In the present paper, we have investigated steady states of macroscopic quantum systems under dissipation not obeying the detailed balance condition. We have theoretically argued that even in such nonequilibrium situations, the Gibbs state at effective temperature is a good description of the steady states. There are two ingredients in emergence of the Gibbs state: the validity of the weak-dissipation perturbation theory and the ETH.

We have derived a criterion for the validity of the perturbation theory beyond the convergence radius, which shrinks to zero in the thermodynamic limit Žnidarič (2015); Lemos and Prosen (2017). It tells us that the perturbation theory works well for sufficiently weak bulk dissipation as long as the Liouvillian is gapped. On the other hand, the perturbation theory breaks down for an arbitrarily weak boundary dissipation because of the vanishing gap of the Liouvillian due to the presence of diffusive transport. Our numerical calculations have confirmed those theoretical predictions.

There remain some issues to be studied. The effect of an extensive number of conserved quantities in integrable models on the steady states should be studied. Our theoretical criterion (19) has been derived by using the off-diagonal ETH, which is not valid in integrable systems. It is expected that the steady state is well described by a generalized Gibbs ensemble under a certain condition Lenarčič et al. (2018). The extension of our theory to systems with finite dissipation strength are also important open problems.

Acknowledgements.

The authors would like to thank Yuto Ashida, Ryusuke Hamazaki, Tomotaka Kuwahara and Keiji Saito for useful discussion. The numerical calculations have been done mainly on the supercomputer system at Institute for Solid State Physics, University of Tokyo. This work was supported by the Japan Society for the Promotion of Science KAKENHI Grant No. JP18K13466 and No. JP19K14622.

Appendix A Convergence radius of perturbative series

Lemos and Prosen Lemos and Prosen (2017) developed a numerical method to estimate the convergence radius of the perturbative expansion, $\gamma_{\rm c}$, and they argued that $\gamma_{\rm c}$ shrinks to zero with the system size for generic open quantum systems. In the main text, we have performed the linear fitting of the numerical data for large system size at different values of $\gamma$ to obtain $\lim_{V\to\infty}s$. In the argument, we have assumed that $\gamma=0.005$ is greater than $\gamma_{\rm c}$ for $V=9$. Here, we show that it is true for our model [Eq. (20)].

In Fig. 3(a), the system-size dependence of $\gamma_{\rm c}$ is presented up to $V=5$. The convergence radius $\gamma_{\rm c}$ shows the exponential decay with the system size, and it suggests that $\gamma_{\rm c}$ shrinks to zero in the thermodynamic limit. However, as in Ref. Lemos and Prosen (2017) we could not obtain $\gamma_{\rm c}$ for larger system size.

In order to estimate $\gamma_{\rm c}$ for larger system size $(V\geq 6)$, we compare two perturbative solutions, which are obtained by truncation of the perturbation series [Eq. (3)] up to fourth order and sixth order. In Fig. 3(b), we plot the expectation values of $\sum_{i=1}^{V}S_{i}^{x}/V$ over each perturbative solution by solid curve (fourth order) and dotted curve (sixth order), respectively. We found that two curves start to deviate at a certain value of $\gamma$, and the value is close to $\gamma_{\rm c}$ for $V=5$. We use this relation to estimate the approximated values of $\gamma_{\rm c}$ for larger system size. In Fig. 3(b), we found the approximated value of $\gamma_{\rm c}$ decreases with the system size and it is smaller than $0.005$ for $V=9$.

Appendix B Derivation of the Lindblad equation

We give a microscopic derivation of the Lindblad equation, Eq. (20). The following derivation is essentially the same as the one found in Refs. Verso and Ankerhold (2010); Torre et al. (2013) although the model is different. Our microscopic model consists of a chain of laser-driven Rydberg atoms in contact with a bath of harmonic oscillators Lee et al. (2011); Carr et al. (2013); Letscher et al. (2017). Each atom is regarded as a spin-1/2 (a two-level system). The up and down spin states correspond to an excited Rydberg state and the ground state, respectively. Let us denote by $h^{x}$ and $\omega$ the amplitude and the frequency of the laser driving. Then, the Hamiltonian of the Rydberg atoms alone is given by

[TABLE]

where $h^{z}$ is the detuning frequency and $g$ is the strength of interactions between neighboring atoms.

The Hamiltonian of the total system, including the bath, can be written as

[TABLE]

where $H_{\rm B}$ and $H_{\rm I}$ are the Hamiltonians of the bath and the interaction between the Rydberg atoms and the bath, respectively. Here, $\lambda$ represents the interaction strength, which is assumed to be small. The explicit forms of the Hamiltonians $H_{\rm B}$ and $H_{\rm I}$ are given by

[TABLE]

where $a_{i,\alpha}$ and $a_{i,\alpha}^{\dagger}$ are bosonic annihilation and creation operators of the bath modes, respectively. It is assumed that each spin is coupled to its own thermal bath independently.

The time dependence of the Hamiltonian can be eliminated by moving to a rotating frame. For this purpose, we introduce a unitary operator,

[TABLE]

The total Hamilton in the rotating frame reads

[TABLE]

where $H$ and $H_{\rm B}^{\rm R}$ are given by Eq. (20) and

[TABLE]

respectively. Note that the energy of each bath mode in $H_{\rm B}^{\rm R}$ is shifted by $\omega$ with respect to that in $H_{\rm B}$. Then, the bath in the rotating frame does not satisfy the Kubo-Martin-Schwinger (KMS) relation [see Eq. (34) below], which results in the violation of the detailed balance condition [see Appendix D]. Therefore, in the rotating frame, the problem is equivalent to that of an open quantum system in contact with an out-of-equilibrium environment.

The relaxation dynamics of the system of interest with a weak system-bath coupling is described by a Markovian quantum master equation. By applying the Born-Markov approximation, we obtain the following master equation Kubo et al. (1991); Breuer and Petruccione (2002),

[TABLE]

where $\rho$ is the density matrix of the system of interest in the rotating frame and $S_{i}^{\pm}(u)=e^{iHu}S_{i}^{\pm}e^{-iHu}$. Here, $G_{1}^{\rm R}(\epsilon)$ and $G_{2}^{\rm R}(\epsilon)$ denote the bath correlation functions in the rotating frame,

[TABLE]

where $\rho_{\rm B}$ is a Gibbs state at inverse bath temperature $\beta_{\rm B}$: $\rho_{\rm B}=e^{-\beta_{\rm B}H_{\rm B}}/{\rm Tr}e^{-\beta_{\rm B}H_{\rm B}}$. It is noted that the Gibbs state is invariant under the unitary transformation, i.e. $U(t)\rho_{\rm B}U^{\dagger}(t)=\rho_{\rm B}$. The bath correlation functions in the rotating frame are related to those in the original frame $G_{1}(\epsilon)$ and $G_{2}(\epsilon)$ as

[TABLE]

where

[TABLE]

The correlation functions in the original frame satisfy the KMS relation Kubo (1957),

[TABLE]

while those in the rotating frame do not satisfy this relation;

[TABLE]

As is mentioned before, this implies that the bath is in thermal equilibrium in the original frame but not in the rotating frame.

The integration over $u$ in Eq. (29) gives

[TABLE]

where $H\ket{n}=E_{n}\ket{n}$ and ${\rm P}$ denotes the Cauchy principal value. Let us denote by $D[\rho]$ the dissipation part of Eq. (35). Then, Eq. (35) is written as $d\rho/dt=-i[H,\rho]+D[\rho]$.

The master equation, Eq. (35), is simplified under the assumptions that $\omega\gg\max\{h^{z},h^{x},g\}$. It is proved that matrix elements $\braket{n}{S_{i}^{\pm}}{m}$ decay exponentially in $|E_{n}-E_{m}|$, and the dominant contribution in the integral over $\epsilon$ in Eq. (35) comes from $|\epsilon|\lesssim\max\{h^{z},h^{x},g\}\ll\omega$. Therefore, we can neglect the dependences of correlation functions on $\epsilon$,

[TABLE]

In this approximation, the Cauchy principal integrals become zero, and thus the dissipation part of Eq. (35) reads

[TABLE]

where we have used the KMS condition, Eq. (33). If we further assume that $\omega\gg 1/\beta_{\rm B}$, we obtain

[TABLE]

where $\gamma\equiv 2\pi\lambda^{2}G_{1}(-\omega)$. This is identical to the Lindblad form in main text, Eq. (20). The Lindblad operators $\{S_{i}^{-}\}$ describe the transition from the up spin state to the down spin state due to dissipation.

Appendix C Effective temperature

We have shown in the main text that the steady state is locally indistinguishable from the Gibbs state $\rho_{\beta}=e^{-\beta H}/\mathrm{Tr}\,e^{-\beta H}$ for some open quantum systems that obey the Lindblad equation

[TABLE]

Here, $\beta$ is an inverse effective temperature that characterizes the steady state.

We provide a method to estimate $\beta$. The steady state of the Lindblad equation $\rho_{\rm s}$ is obtained by setting the right-hand side of Eq. (39) to zero. By multiplying it by $H$ from the left and then taking its trace, we obtain the following equation:

[TABLE]

Since $[L_{a}^{\dagger},H]L_{a}$ is a local operator, $\rho_{\rm s}$ in Eq. (40) can be replaced by $\rho_{\beta}$. Then, we obtain Eq. (10) in the main text, i.e.,

[TABLE]

where $\mathrm{Re}$ has been omitted since $\mathrm{Tr}\,[L_{a}^{\dagger},H]L_{a}\rho_{\beta}$ is always real. We can estimate the value of $\beta$ by using Eq. (41).

There exists at least one solution because $C(\beta)$ is a continuous function of $\beta$ and satisfies $C(+\infty)\geq 0$ and $C(-\infty)\leq 0$. We can prove $C(+\infty)\geq 0$ as follows. At $\beta=+\infty$, $\rho_{\beta}=\ket{\Psi_{\rm min}}\bra{\Psi_{\rm min}}$ where $\ket{\Psi_{\rm min}}$ is the ground state of $H$ with energy $E_{\rm min}$. Then,

[TABLE]

where

[TABLE]

We can also prove $C(-\infty)\leq 0$ in the similar way.

We apply this method to the Lindblad equation in main text [see Eq. (20)]. In Fig. 4 (a), we plot $C(\beta)$ for the system size $V=4$. As it is mentioned above, $C(\beta)$ is negative at $\beta=-\infty$ $[C(-\infty)\simeq-22.6]$ and positive at $\beta=+\infty$ $[C(+\infty)\simeq 3.83]$. Between them, $C(\beta)$ has a maximum around $\beta\simeq 2.46$, and thus it is not a monotonic function of $\beta$. In the present case, there is only one solution for $C(\beta)=0$, that determines the inverse effective temperature $\beta$. In Fig. 4 (b), the system-size dependence of $\beta$ is depicted by red circles. The estimated values of $\beta$ show a weak system-size dependence, and they are almost converged to $\beta\simeq 0.256$ for $V\geq 7$.

The inverse effective temperature $\beta$ is also evaluated by comparing the expectation values of the energy between $\rho_{\rm s}^{(0)}$ and $\rho_{\beta}$:

[TABLE]

In Fig. 4 (b), we also plot the system-size dependence of $\beta$ obtained in this method by black squares. The estimated values of $\beta$ in two methods approach each other with increasing system size. It seems that $\beta$ obtained by $C(\beta)=0$ quickly converges to the value in the thermodynamic limit. It should be emphasized that the estimation of $\beta$ using Eq. (44) is a numerically hard task since it requires the evaluation of the transition rates $W_{ij}$ for all the pairs of energy eigenstates (to do so, we have to diagonalize the Hamiltonian), while solving Eq. (41) is much easier since it only requires the calculation of the equilibrium expectation values of $[L_{a}^{\dagger},H]L_{a}$.

Appendix D Violation of the detailed balance condition

In this appendix, we demonstrate the violation of the detailed balance condition in the model given by Eq. (20) in the main text. The “transition rate” calculated by applying the perturbation theory is given by

[TABLE]

and the detailed balance condition with respect to the Gibbs state is expressed by $\frac{W_{nm}}{W_{mn}}=e^{-\beta(E_{n}-E_{m})}$ or equivalently,

[TABLE]

for any pair of eigenstates $n$ and $m$.

In order to judge whether the detailed balance condition holds in the model given by Eq. (20) in the main text, we make a scaled histogram for all the pairs of $n$ and $m$,

[TABLE]

where $\ell\in\mathbb{Z}$ and $\Delta_{\rm d}=0.1$ is the bin size of the histogram and $\beta=0.256$ is the inverse effective temperature [see Appendix C]. The detailed balance condition holds when $h_{\rm all}(\ell)=\delta_{\ell,0}$. In Fig. 5(a), we plot the scaled histogram for different system sizes, $V=10$ (red solid line) and $V=16$ (blue dotted line). They are almost overlapped with each other, and have a peak at $\ell=0$ with finite width. Thus, the detailed balance condition is violated irrespective of the system size.

In the histogram $h_{\rm all}(\ell)$, all the transition rates between the energy eigenstates are taken into account, but transition rates between the states with macroscopically different energies will be irrelevant to determine the steady state. Thus, we produce another scaled histogram that omit such irrelevant contributions:

[TABLE]

where $\ell\in\mathbb{Z}$, $\bar{E}={\rm Tr}H\rho_{\rm s}^{(0)}$, and $\Delta E=4\simeq 1/\beta$. Normalization constant ${\cal N}$ is determined from the condition that $\sum_{\ell}h_{\rm part}(\ell)=1$. In Fig. 5(b), we plot the scaled histogram $h_{\rm part}(\ell)$ for different system sizes. As in Fig. 5(a), it implies the violation of the detailed balance condition irrespective of the system size.

We also study the detailed balance condition with respect to $\{P_{n}\}$ in Eq. (6). We make the corresponding scaled histograms,

[TABLE]

where $\ell\in\mathbb{Z}$ and ${\cal N}^{(0)}$ is the normalization constant determined by $\sum_{\ell}h_{\rm part}^{(0)}(\ell)=1$. In Figs. 6, we plot the scaled histograms $h_{\rm all}^{(0)}(\ell)$ and $h_{\rm part}^{(0)}(\ell)$ for different system sizes, which also imply the violation of the detailed balance condition with respect to $\{P_{n}\}$.

Appendix E Emergence of the Gibbs state in another dissipative spin system

We provide another model with bulk dissipation that shows the same $\gamma$-dependence of $\lim_{V\to\infty}s$ as the model in the main text. The Hamiltonian and the Lindblad operators are given by

[TABLE]

where $(h^{z},h^{x},g)=(1.809,1.618,4)$. We assume that $V$ is even. There are two types of Lindblad operators: $L_{i}^{(1)}$ acting on every site and $L_{i}^{(2)}$ acting on even sites. The system does not have any conserved current in the steady state, which implies that the steady state is described by a Gibbs state.

In Fig. 7 (a), we show the system-size dependences of $s$ at $\gamma=0.01$, $0.03$, and $0.05$. We find the linear dependence of the distance on $1/V$ for large $V$ ($V\geq 8$), and again we extrapolate the data for each $\gamma$ to the thermodynamic limit. In Fig. 7 (b), we present $\lim_{V\to\infty}s$ as a function of $\gamma$. The figure shows $\lim_{V\to\infty}s\propto\gamma$, which is the same dependence as the model of the main text.