Thermalization with detailed-balanced two-site Lindblad dissipators

TL;DR

This paper introduces a detailed-balanced two-site Lindblad dissipator method for generating thermal states and analyzing heat transport, demonstrating high accuracy in thermalization at large temperatures.

Contribution

It proposes a novel variant of the Lindblad dissipator based on detailed balance, improving thermal state generation and characterization.

Findings

Steady state approaches thermal equilibrium with less than 1% error at high temperatures.

Effective temperatures from various observables are consistent.

Decay of two-time correlations is analyzed.

Abstract

The use of two-site Lindblad dissipator to generate thermal states and study heat transport raised to prominence since [J. Stat. Mech. (2009) P02035] by Prosen and \v{Z}nidari\v{c}. Here we propose a variant of this method based on detailed balance of internal levels of the two site Hamiltonian and characterize its performance. We study the thermalization profile in the chain, the effective temperatures achieved by different single and two-site observables, and we also investigate the decay of two-time correlations. We find that at a large enough temperature the steady state approaches closely a thermal state, with a relative error below 1% for the inverse temperature estimated from different observables.

Figures

Equations

$$H=\sum _{l=1}^L(h_x{\sigma }_l^x+h_z{\sigma }_l^z)+\sum _{l=1}^{L-1}J({\sigma }_l^x{\sigma }_{l+1}^x+{\sigma }_l^y{\sigma }_{l+1}^y)+\Delta {\sigma }_l^z{\sigma }_{l+1}^z,$$

$$\dot{\rho }$$

$$=-\frac{i}{\hslash }[H,\rho ]+\sum {\gamma }_{p,q}\left({\Gamma }_p\rho {\Gamma }_q^{†}-\frac{1}{2}\left\{{\Gamma }_p^{†}{\Gamma }_q,\rho \right\}\right),$$

$${\rho }_T=\sum _{j=0}^3{P}_j|j⟩⟨j|,$$

$$P_j=\frac{e^{-{\beta }_B{E}_j}}{{\sum }_i{e}^{-{\beta }_B{E}_i}},$$

$$\dot{\rho }$$

$$+\gamma \sum _{n,m>n}[{\lambda }_{(n,m)}^{+}\left({\Gamma }_{(n,m)}^{†}\rho {\Gamma }_{(n,m)}-\frac{1}{2}\left\{{\Gamma }_{(n,m)}^{†}{\Gamma }_{(n,m)},\rho \right\}\right)$$

$$+{\lambda }_{(n,m)}^{-}\left({\Gamma }_{(n,m)}\rho {\Gamma }_{(n,m)}^{†}-\frac{1}{2}\left\{{\Gamma }_{(n,m)}{\Gamma }_{(n,m)}^{†},\rho \right\}\right)]$$

$$0=$$

$$+P_3{\lambda }_{(0,3)}^{-}-P_0{\lambda }_{(0,3)}^{+}$$

$$0=$$

$$+P_3{\lambda }_{(1,3)}^{-}-P_1{\lambda }_{(1,3)}^{+}$$

$$0=$$

$$-P_2{\lambda }_{(0,2)}^{-}+P_0{\lambda }_{(0,2)}^{+}$$

$$0=$$

$$-P_3{\lambda }_{(1,3)}^{-}+P_1{\lambda }_{(1,3)}^{+}.$$

$${\lambda }_{(0,2)}^{+}=$$

$${\lambda }_{(0,3)}^{+}=$$

$$+\frac{P_1+P_2}{P_0+P_3}{\lambda }_{(1,2)}^{+}-\frac{P_2+P_3}{P_0+P_3}{\lambda }_{(2,3)}^{+},$$

$${\lambda }_{(1,3)}^{+}=$$

$$0$$

$${\lambda }_{(0,1)}^{+}$$

$${\lambda }_{(1,2)}^{+}$$

$${\lambda }_{(2,3)}^{+}$$

$${\lambda }_{(0,2)}^{+}=\frac{P_2}{P_0+P_2},$$

$${\lambda }_{(0,3)}^{+}=\frac{P_3}{P_0+P_3},$$

$${\lambda }_{(1,3)}^{+}=\frac{P_3}{P_1+P_3}.$$

$$H_{2S}$$

$$⟨{\sigma }_i^{+}(t){\sigma }_i^{-}(0){⟩}_s=\mathrm{tr}\left[{\sigma }_i^{+}{e}^{\mathcal{L}t}\left({\sigma }_i^{-}{\rho }_s\right)\right]$$

$${\rho }_T$$

$$=V\left(\sum _{j=0}^3{P}_j|j⟩⟨j|\right)V^{†}$$

$$V=\sum _{s_1,s_2,j}{\mu }_j^{s_1,s_2}|s_1,s_2⟩⟨j|$$

$$V{\Gamma }_{(n,m)}{V}^{†}=\sum _{s_1,s_2,t_1,t_2}{\mu }_p^{s_1,s_2}({\mu }_q^{t_1,t_2}{)}^{\ast }|s_1,s_2⟩⟨t_1,t_2|.$$

$$\displaystyle\left(\begin{array}[]{c}|\downarrow,\downarrow\rangle\\ |\downarrow,\uparrow\rangle\\ |\uparrow,\downarrow\rangle\\ |\uparrow,\uparrow\rangle\end{array}\right)=\left(\begin{array}[]{cccc}\sigma^{d}_{1}\sigma^{d}_{2}&\sigma^{d}_{1}\sigma^{-}_{2}&\sigma^{-}_{1}\sigma^{d}_{2}&\sigma^{-}_{1}\sigma^{-}_{2}\\ \sigma^{d}_{1}\sigma^{+}_{2}&\sigma^{d}_{1}\sigma^{u}_{2}&\sigma^{-}_{1}\sigma^{+}_{2}&\sigma^{-}_{1}\sigma^{u}_{2}\\ \sigma^{+}_{1}\sigma^{d}_{2}&\sigma^{+}_{1}\sigma^{-}_{2}&\sigma^{u}_{1}\sigma^{d}_{2}&\sigma^{u}_{1}\sigma^{-}_{2}\\ \sigma^{+}_{1}\sigma^{+}_{2}&\sigma^{+}_{1}\sigma^{u}_{2}&\sigma^{u}_{1}\sigma^{+}_{2}&\sigma^{u}_{1}\sigma^{u}_{2}\end{array}\right)\left(\begin{array}[]{c}|\downarrow,\downarrow\rangle\\ |\downarrow,\uparrow\rangle\\ |\uparrow,\downarrow\rangle\\ |\uparrow,\uparrow\rangle\end{array}\right),$$

Full text

I Introduction

The ability to control heat flow at the nanoscale can offer important technological opportunities. In many-body quantum systems the interactions can significantly change the transport properties, e.g., turning it from ballistic to diffusive Prosen (2011). Recently it was also shown that interactions can be instrumental in realizing highly performing diodes Balachandran et al. (2018, ). Of particular importance is the current between heat baths at different temperatures, i.e., heat current. However, the study of transport in the presence of heat baths requires an appropriate modelling of the baths which can be extended to intermediate-size and large many-body quantum systems.

Different approaches have been used to model the coupling of a many-body system to a thermal bath, some global and some local. In global approaches the model considers the Hamiltonian of the entire system, and it usually requires to diagonalize it, which is something that can be done only for limited system sizes. One of the global approaches uses master equations in Gorini-Kossakovski-Sudarshan-Lindblad (GKSL) form Gorini et al. (1976); Lindblad (1976) (see also Breuer and Petruccione (2007)). The use of this method, however, can result in an unphysical description of the current within a system. In fact, it can predict that a heat current flows from one bath to the system, and from the system to the other bath, while no heat current flows within the system Wichterich et al. (2007). Another global approach that gives more accurate results for the heat current within the chain is the Redfield master equation Redfield (1957). Recently it was even shown how to implement a Redfield master equation on a large system with matrix product states and operators Xu et al. (2019), although only for short-time evolutions. Other approaches to study many-body quantum systems in contact with a finite temperature bath of harmonic oscillators are based on the star-to-chain transformation, which converts the bath to a bosonic chain Chin et al. (2010); Prior et al. (2010); de Vega and Bañuls (2015); Guo et al. (2018); Cascio et al. . This strategy is particularly effective if preceded by a thermofield transformation which maps a finite temperature bath to two zero-temperature baths de Vega and Bañuls (2015); Guo et al. (2018); Cascio et al. . As of now, this approach is limited to relatively short-time evolutions and to noninteracting baths. For the cases in which it is possible to assume translational invariance in the transfer tensors describing the time evolution, it is possible to reach longer times Cerrillo and Cao (2014); Rosenbach et al. (2016).

A few approaches have been proposed for the long-time, or even steady state, study of many-body quantum systems in contact with a heat bath Gelman et al. (2004); Katz et al. (2008); Torrontegui and Kosloff (2016); Reichental et al. (2018); Prosen and Žnidarič (2009). Further investigation is needed to clarify which of these methods performs better, depending on the properties of both the systems and the baths studied. Here we focus on the two-site thermalization approach proposed in Prosen and Žnidarič (2009), which is readily implemented with matrix product states algorithms, and which has already been used to investigate transport physics Žnidarič et al. (2010); Žnidarič (2011a, b); Mendoza-Arenas et al. (2015). In this case, the dissipative part of the master equation in GKSL form acts simultaneously on the two sites closest to each edge trying to thermalize them. This approach improves significantly the results for local master equations whose dissipator acts on a single site. Here, we propose a modification of the scheme presented in Prosen and Žnidarič (2009) to obtain a modified two-site dissipator and study how it thermalizes the system of interest.

This paper is divided in the following sections: In Sec. II we describe the spin chain we aim to thermalize, and in Sec. III we describe our approach for two-site thermalizing GKSL master equations. In Sec. IV we characterize the performance of our method both with static and dynamical observables, and in Sec. V we draw our conclusions.

II Spin chain Hamiltonian

For the two-site dissipator to be able to thermalize the system, as we will see in more detail later, we need to consider a non-integrable Hamiltonian. We thus consider the following, one-dimensional spin Hamiltonian

[TABLE]

for a chain of length $L$, where $h_{x}$ is the magnitude of the transverse field, $h_{z}$ the longitudinal field, $J$ the tunneling amplitude and $\Delta$ the nearest neighbour interaction. The $\sigma_{l}^{a}$ with $a=x,\;y,\;z$ are operators associated to the Pauli matrices. This Hamiltonian is such that when $h_{x}=0$, or $h_{x}=J=0$ or $h_{z}=J=0$ the system is integrable. For $J=0$ we recover the Ising chain with transverse and longitudinal field, whose chaotic properties were analyzed in Ref. Mejía-Monasterio et al. (2005). For $J\neq 0$ we get an $XXZ$ chain with a transverse and a longitudinal field.

When the Hamiltonian is integrable, the statistics of the energy level spacing follows a Poisson distribution. On the contrary, non-integrable systems can have a level spacing statistics which follows a Wigner-Dyson distribution. Following Ref. Oganesyan and Huse (2007), it is possible to characterize how closely to a Poisson or a Wigner-Dyson distribution the system is. In Ref. Oganesyan and Huse (2007) the authors proposed to compute the average of the ratio between all the consecutive energy differences. This means that if $\delta_{n}$ is the energy difference between level $n$ and $n+1$, they compute $\bar{r}$, which is the arithmetic average of the $r_{n}$ where $r_{n}=\min(\delta_{n}/\delta_{n+1},\delta_{n+1}/\delta_{n})$ is the ratio between two consecutive energy difference. For a Hamiltonian whose level spacing follows a Poisson distribution, $\bar{r}=2\ln 2-1\approx 0.386$, while for the Wigner-Dyson $\bar{r}\approx 0.5295$ Oganesyan and Huse (2007). Note that for $J_{x}=0$ the Hamiltonian conserves the total magnetization, while for $J_{x}\neq 0$ the Hamiltonian has a reflection symmetry $j\rightarrow L-j$. We thus have to restrict ourselves to distinct number conserving sectors of the Hamiltonians to characterize the statistics of the energy levels.

III Two-site thermalizing baths

In order to thermalize the spin chain described in Sec. II, we apply a dissipator $\mathcal{D}$ in Lindblad form at one or two edges of the chain. The master equation is thus written as

[TABLE]

where $\rho$ is the density operator, $\Gamma_{p}$ are possible jump operators which act locally on one or the other edge, $\gamma_{p,q}$ the corresponding weight of each combination of jump operators and $\left\{\cdot,\cdot\right\}$ stands for the anti-commutator. With $\mathcal{L}$ we indicate the superoperator acting on the density operator $\rho$. Note that hereafter we set $\hbar=k_{B}=1$, where $k_{B}$ is the Boltzmann constant, and the unit of energy is set such that $\Delta=-1$. The matrix formed by the elements $\gamma_{p,q}$ has to be non-negative to guarantee that the density matrix $\rho$ is physical at all times Gorini et al. (1976); Lindblad (1976). In the following we use jump operators which consists of all the possible combinations of two-site operators composed by the identity matrix and the Pauli matrices. The task is then to figure out the elements $\gamma_{p,q}$ such that the overall chain approaches a thermal state in the steady state. In Ref. Prosen and Žnidarič (2009) the authors chose the $\gamma_{p,q}$ such that if the system comprised of only two sites, the steady state would be exactly the thermal state of an appropriately chosen two-site Hamiltonian. However this condition does not uniquely define the $\gamma_{p,q}$. In fact, this condition only ensures that the steady state is correct, but it does not fix the decay rate of all the other eigenmodes of Eq. (III), the so-called rapidities. Hence there is significant freedom in the choice of the jump operators $\Gamma_{p}$ and dissipation rate matrix $\gamma_{p,q}$. The choice to build the coefficients $\gamma_{p,q}$, and the operators $\Gamma_{p}$, done in Ref. Prosen and Žnidarič (2009) is such that all rapidities, non-corresponding to the steady state, decay quickly at the same rate. We will thus refer to this approach as the fast relaxation approach.

In the following we detail a different approach which is insipired by Ref. Prosen and Žnidarič (2009), but that uses a different condition to determine uniquely the master equation. This approach is based on maintaining detailed balance between the four energy levels of the two-site Hamiltonian which the bath is coupled to. We describe this method, which we refer to as detailed balance approach, in the following subsection.

III.1 Detailed balance approach

The dissipator is designed to generate a thermal state on the two sites it acts upon. The two sites form a four-level system (levels [math] to $3$ in increasing energies $E_{0}$ to $E_{3}$) and hence there are $6$ possible transitions and $12$ possible jump operators if we consider raising and lowering rates between all possible levels. In this rotated basis we can write the target density operator $\rho_{T}$ as

[TABLE]

where the $P_{j}$ are given by

[TABLE]

where the $E_{i}$ are the energies of the target two-site (i.e., four-level) Hamiltonian $H_{2S}$ and $\beta_{B}$ is the target inverse temperature that the bath will try to set. The master equation can then be written as

[TABLE]

where $\Gamma_{(n,m)}=|n\rangle\langle m|$ (or $\Gamma^{\dagger}_{(n,m)}=|m\rangle\langle n|$) is the operator that lowers (raises) the four-level system from level $m$ to $n$ ($n$ to $m$), $\gamma$ is a global relaxation rate, and $\lambda^{+}_{(n,m)}$ and $\lambda^{-}_{(n,m)}$ are respectively the relative rate of raising state $n$ to $m$ or lowering $m$ to $n$. Using Eq. (3) in Eq. (III.1), and imposing the steady-state condition (i.e., $d\rho/dt=0$), we obtain the following set of equations

[TABLE]

We choose $\lambda_{(n,m)}^{-}=1-\lambda_{(n,m)}^{+}$, and solve the above set of equations to get

[TABLE]

where the coefficients for close energy levels $\lambda^{+}_{(0,1)},\;\lambda^{+}_{(1,2)},\;\lambda^{+}_{(2,3)}$ are taken as free parameters. We now consider the same problem as if $\lambda^{+}_{(0,2)}=\lambda^{+}_{(0,3)}=\lambda^{+}_{(1,3)}=0$ (only nearest energy level transitions) and we get the simpler set of equations

[TABLE]

which have the unique solution (by considering again $\lambda_{(n,m)}^{+}=1-\lambda_{(n,m)}^{-}$)

[TABLE]

Using the results from Eq. (III.1) in Eq. (III.1) we get the remaining parameters

[TABLE]

We thus have described a method to obtain the desired target density operator $\rho_{T}$. It is simply necessary to rotate the operators from the two-site Hamiltonian diagonal basis to the computational basis to obtain the jump operators used in our calculations. We have given some details on how to implement this in Appendix A.

It is useful now to discuss about the target density operator $\rho_{T}$ as various approaches can be followed. Here we use the two-site Hamiltonian

[TABLE]

which is the Hamiltonian for each internal bond of the chain of Eq.(1). We then compute the thermal state $e^{-\beta H_{2S}}/\operatorname{tr}(e^{-\beta H_{2S}})$ and diagonalize it.

IV Results

In this section we study the steady state that is reached when applying the dissipator discussed in Sec. III.1. If the steady state is an exact thermal state, then the expectation value of any observable is determined once the temperature and the Hamiltonian are known. Consequently, by evaluating different observables, it is possible to check whether their expectation values are consistent with the steady state being thermal. In the following we consider different observables and we compute their expectation values in the steady state. Then, for each of them we check what would be the temperature of a thermal state which has the same expectation value for that observable. To do so, however, it is first necessary to have a reference temperature for the expectation values that the various observables can take. We discuss this in Sec. IV.1. We will then consider both static and dynamical observables, respectively in Sec. IV.2 and Sec. IV.3. Hereafter, for simulations with $L\geq 13$ we use a matrix product states algorithm for open systems Schollwöck (2011); Verstraete et al. (2004); bon .

IV.1 Reference temperature

From the Hamiltonian in Eq. (1), we compute thermal states for the reference inverse temperature $\beta_{R}\in[0,1]$ for chains of lengths from $L=5$ to $13$. The typical behavior of observables as a function of $\beta_{R}$ is depicted in Fig. 1$(a)$. In particular it shows $\langle\sigma^{z}_{l}\rangle$, $\langle\sigma^{x}_{l}\rangle$ and $\langle\sigma^{z}_{l}\sigma^{z}_{l+1}\rangle$ with $l=6$, the middle site of a chain with length $L=11$. At low $\beta_{R}$ (high temperatures) the observable changes sizeably with $\beta_{R}$, while at larger $\beta_{R}$ (lower temperatures), the observable is almost constant and it would be very hard to clearly identify the temperature corresponding to a given value of the observable. In Figs. 1$(b$-$d)$ we consider five different chain lengths $L(k)=2k+3$ with $k\in[1,5]$, and in each panel we show the absolute value of the difference between the expectation value of an observable in the middle of the chain for two different lengths $L(k+1)$ and $L(k)$, plotted versus the reference inverse temperature $\beta_{R}$. For example, in Fig. 1$(b)$ we show $|\Delta\langle\sigma^{x}_{i}\rangle|=|\langle\sigma^{x}_{k+3}\rangle_{k+1,\beta_{R}}-\langle\sigma^{x}_{k+2,\beta_{R}}\rangle_{k,\beta_{R}}|$ where with $\langle\cdot\rangle_{k,\beta_{R}}$ we mean the expectation value for a chain of length $L(k)$ and at an inverse temperature $\beta_{R}$. In Fig. 1$(c)$ we show $|\Delta\langle\sigma^{z}_{i}\rangle|=|\langle\sigma^{z}_{k+3}\rangle_{k+1,\beta_{R}}-\langle{\sigma^{z}_{k+2}}\rangle_{k,\beta_{R}}|$ and in Fig. 1$(d)$ we show $|\Delta\langle\sigma^{z}_{i}\sigma^{z}_{i+1}\rangle|=|\langle\sigma^{z}_{k+3}\sigma^{z}_{k+4}\rangle_{k+1,\beta_{R}}-\langle\sigma^{z}_{k+2}\sigma^{z}_{k+3}\rangle_{k,\beta_{R}}|$. Figs. 1$(b$-$d)$ show that for large enough temperatures and for $L=11$ the value of the observable in the middle of this non-integrable chain has converged. It is thus possible to use a chain of size $L=11$ to associate a temperature to a certain value of an observable.

In the following we use Eq.(III) to derive a steady state. Then, given and observable at a site $i$, we associate an effective inverse temperature which is given by the reference inverse temperature $\beta_{R}$ of a thermal state for which the expectation value of the same observable is identical. For example, when studying a chain of length $L=17$, we compute $\langle\sigma^{x}_{9}\rangle$ and check for which inverse temperature $\beta_{R}$ a thermal chain would have the same value for this observable. We then refer to this inverse temperature as $\beta_{x}$. Analogously we compute $\beta_{z}$ and $\beta_{zz}$. We shall note here that in Mendoza-Arenas et al. (2015) the authors suggest to compute the effective temperature by the minimum distance between a local reduced density matrix and the one from a reference thermal state.

IV.2 Static observables

We first consider the effective inverse temperatures $\beta_{z}$, computed from the middle site of the chain, versus the target inverse temperature $\beta_{B}$ that the dissipator $\mathcal{D}$ is trying to impose at the edge. In Fig. 2$(a)$ we show $\beta_{z}$ versus $\beta_{B}$. For $\beta_{B}\lesssim 1.5$ we find that $\beta_{z}$ grows linearly, whereas for larger values of $\beta_{B}$ it reaches a plateau. We have observed a similar behavior for all observables examined and for different types of chain Hamiltonians (here we show the XXZ chain with $h_{x}=3.375$, $h_{z}=2$, $J=-0.5$, $\Delta=-1$ and $L=17$). We then analyze the effective temperature reached in the system at different sites for low bath inverse temperature, $\beta_{B}=0.1$ in panel $(b)$, and larger bath inverse temperature, $\beta_{B}=1.5$ in panel $(c)$. As discussed previously, the effective temperature associated at each site to an observable is given by the temperature of a thermal state for which that observable has the same expectation value, e.g., $\beta_{x}$ is associated to $\langle\sigma^{x}_{l}\rangle$, $\beta_{z}$ to $\langle\sigma^{z}_{l}\rangle$ and $\beta_{z}$ to $\langle\sigma^{z}_{l}\sigma^{z}_{l+1}\rangle$. The difficulty here is that the size of the systems taken as reference is limited (here we consider $L=11$ for the reference thermal state). Nevertheless, the value of the observables in a non-integrable chain in a thermal state converges to a certain value within a few sites from the edges of the chain. Hence we can use the values at the edge of the reference thermal chains to obtain the effective temperature at the edge of the chain from the master equation, and for the other sites we use the bulk values. Note that in order to have accurate reference values for the inverse temperature, we use $10^{5}$ values of $\beta_{R}$.

In Figs. 2$(b,c)$ we plot the effective temperature for $\beta_{z}$ (red dot-dashed line with $\circ$), $\beta_{x}$ (blue dashed line with $\times$) and $\beta_{zz}$ (green solid line with $\square$). Note that here the dissipator $\mathcal{D}$ is acting only on sites $l=1,2$. At lower inverse temperature, panel $(b)$, the inverse temperatures $\beta_{z}$ and $\beta_{x}$ are practically indistinguishable for $l\geq 3$, and the temperature is fairly constant throughout the bulk of the chain. For an effective inverse temperature related to the correlation $\beta_{zz}$, we find that in the bulk there is an error of only a few percents, while at the edge the difference is much greater. At larger inverse bath temperature $\beta_{B}=1.5$, see Fig. 2$(c)$, the effective local inverse temperatures are much more clearly different from each other, highlighting the difficulty in reaching good thermal states at larger inverse temperatures.

As mentioned briefly earlier, it is expected that the chain reaches a steady state which resembles more a thermal state when it is non-integrable. Fig. 3 confirms this understanding. In both panel $(a)$ and panel $(b)$ we plot the inverse temperature for different observables $\beta_{M}$ against the Hamiltonian parameter $h_{x}$. In both panels we use $h_{z}=2$, $\Delta=-1$, $L=17$ and $\beta_{B}=0.85$, while in panel $(a)$ $J=-0.5$ and in panel $(b)$ $J=0$. We observe that in some regions the inverse temperatures of the different observables are close to one another, or even indistinguishable in these plots, while in other regions the inverse temperatures are very different. With Fig. 3$(c)$ we can compare the behavior of the effective inverse temperatures with the variation of the mean $\bar{r}$ as a function of $h_{x}$. In particular, the blue solid line shows $\bar{r}$ for $L=15$, $h_{z}=2$, $\Delta=-1$ and $J=-0.5$, while $L=16$ and $J=0$ for the red dashed line. The black dot-dashed horizontal lines show the values for $\bar{r}=0.386$ which would be obtained from a random Hamiltonian with the energy level difference following a Poisson distribution, and $\bar{r}=0.5295$ which is computed from large Gaussian orthogonal ensemble matrices. Although we do not observe a clear quantitative match between $\bar{r}$ and the proximity of the effective inverse temperatures, comparing Fig. 3$(c)$ with Fig. 3$(a,b)$ we observe that when $\bar{r}$ is largest, the inverse temperatures match better between each other, compared to when $\bar{r}$ is closer to $0.386$. In Fig. 3$(a)$ we see that for $h_{x}\lesssim 4$ the inverse temperature are closer than for $h_{x}\gtrsim 4$, which is consistent with the change in $\bar{r}$ shown in Fig. 3$(c)$ by the continuous blue line. In Fig. 3$(b)$ we zoom-in on the $y-$axis to see more clearly that the region of parameters of the transverse field $h_{x}$ in which the inverse temperatures match the most is approximately for $h_{x}\in(3,5)$, which is the region in which $\bar{r}$ is largest (red dashed line in Fig. 3$(c)$).

In order to find the best scenarios in which the steady state is closer to a thermal state, it is also important to check the role of the dissipation rate $\gamma$. In Fig. 4 we show intensity plot of $\beta_{z}$, panel $(a)$, $\beta_{x}/\beta_{z}$, panel $(b)$, and $\beta_{zz}/\beta_{z}$, panel $(c)$, versus inverse bath temperature $\beta_{B}$, and dissipation ratio $\gamma$. Fig. 4$(a)$ shows that when $\gamma$ is large enough (here e.g., for $\gamma\gtrsim 0.3$), the resulting inverse temperature $\beta_{z}$ shows little dependence on $\gamma$ itself. The behavior of the inverse temperature from other observables, $\beta_{x}$ and $\beta_{zz}$, with $\gamma$ has a similar behavior as $\beta_{z}$, as can be inferred from Figs. 4$(b,c)$. Going back to Fig. 4$(a)$, the dependence on $\beta_{B}$, for a given $\gamma$, seems to reach a plateau for $\beta_{B}\gtrsim 1.5$. In Fig. 4$(b,c)$ we note that the ratio between the inverse temperatures is closer to $1$ for larger temperatures, however the dependence of $\beta_{zz}/\beta_{z}$ with $\beta_{B}$ is, in general, non-monotonous.

To evaluate the performance of the detailed balance approach, we now compare it to the fast relaxation one. In Fig. 5, we plot the ratio of the inverse temperatures for different observables for both methods and for different non-integrable chains. More specifically we use blue squares for the detailed balance approach and red triangles for the results from fast relaxation. In Fig. 5 we depict $\beta_{x}/\beta_{z}$ with filled symbols and $\beta_{zz}/\beta_{z}$ with empty ones. In Fig. 5$(a,b)$, we consider chains with $h_{x}=3.375$, $h_{z}=2$, $\Delta=-1$, $L=17$ and $\gamma=0.35$ while in Fig. 5$(a)$ $J=-0.5$ and in Fig. 5$(b)$ $J=0$. We have chosen to represent these two cases as they are fairly representative of what we have observed. In particular, the case in Fig. 5$(b)$ is one of the best we have found for the fast relaxation approach, while in Fig. 5 we have chosen a set of parameters which gives the best performance for the detailed balance approach. In both cases, for the detailed balance approach both inverse temperature ratios are within few percents of $1$ for the inverse temperature of the bath $\beta_{B}\lesssim 1$. The ratios of inverse temperatures are in general further away from $1$ for the fast relaxation approach, and in particular the ratio of the inverse temperatures $\beta_{zz}/\beta_{z}$ is further away from $\beta_{x}/\beta_{z}$ compared to the detailed balance approach. It should be pointed out that the regime of parameters $\gamma$ and $\beta_{B}$ for which the relative error between the inverse tempratures is within $1\%$ is much broader for the detailed balance approach compared to the fast relaxation one. Moreover, the detailed balance approach can get simultaneously well bellow the $1\%$ threshold for both ratios. We should add that while reaching the steady state is faster, for all cases we checked, with the fast relaxation method, the time required by detailed balance approach to reach the steady state approaches that of the fast relaxation method when we increase the size of the system. Hence the improved accuracy of the detailed balance approach does not come at a particularly high cost in simulation time. It should also be mentioned that while we consider here chains of intermediate size, the ratio between effective temperatures from different observables does not improve significantly when increasing the system size.

IV.3 Dynamical correlations

While it is possible to produce steady states with targeted static properties, their dynamical properties may differ significantly from those of the target state Sciolla et al. (2015). For this reason we investigate the two-time correlation

[TABLE]

where $\langle\dots\rangle_{s}$ stresses that the two-time correlation is computed from the steady state $\rho_{s}$. In order to evaluate such a quantity, we first compute the steady state of Eq. (III), we then apply $\sigma^{-}_{i}$ to it, followed by evolving the system with Eq. (III). Last we apply $\sigma^{+}_{i}$ and take the trace.

We now consider the spin chain Hamiltonian of Eq. (1) and we couple the two-site bath with the parameters from the detailed balance approach only on one edge of the chain (sites $1$ and $2$). We then study the two-time correlation $\langle\sigma^{+}_{L}(t)\sigma^{-}_{L}(0)\rangle_{s}$ on the last site for different system lengths $L$, bath inverse temperatures $\beta_{B}$ and stronger coupling strengths $\gamma$. In Fig. 6 we consider different system size, in particular $L=13$ for panel $(a)$, $L=17$ for panel $(b)$ and $L=25$ for panel $(c)$. In each panel we show two curves: the red solid line represents the two-time correlations while keeping the system coupled to the Lindblad bath (i.e., $\gamma\neq 0$), while for the blue dashed line we take $\gamma=0$, i.e., the chain which was thermalized by the bath is now no longer coupled to it fn (1). We notice that at a certain time, indicated in the figures by a dot-dashed vertical line, the effects of the boundaries, in the form of beating-like oscillations, are felt in the correlation. In fact, the time at which the two-time correlation starts to have beatings increases with the system size. In absence of coupling to the bath the finite size effects are much larger, especially for shorter chains.

In Fig. 7 we study the evolution of the absolute values of the two-time correlations for different values of the inverse temperatures of the bath $\beta_{B}$. Similarly to Fig. 6 the red solid lines correspond to an evolution in which the chain and the bath are always connected, while the blue dashed lines (dotted line in panel $(d)$) are for evolutions with $\gamma=0$. The inverse temperature increases in the panel from $\beta_{B}=0.01$ in panel $(a)$, $\beta_{B}=0.05$ in panel $(b)$, $\beta_{B}=0.2$ in panel $(c)$ and $\beta_{B}=0.5$ in panel $(d)$. At larger $\beta_{B}$, the correlations are expected to decay more slowly, however, in this finite-size chain, the correlations seem to decay only for a finite time, and then stabilize to a non-zero value for long times (as expected the oscillations around the mean are smaller when $\gamma\neq 0$). At lower inverse temperatures, the correlations are expected to decay faster. We see that the correlations at the beginning follow a decay well approximated by a power law. At longer time the finite-size effects come into play. For $\gamma=0$ the correlations decay more slowly, while if the chain is still coupled to the bath, i.e., $\gamma=0.35$ (red solid lines), the correlations can decay even faster than the initial power-law regime. Using an intermediate strength of coupling to the bath, see yellows lines for $\gamma=0.1$ in Fig. 7$(a,b)$, results in a decay that resembles the power-law decay for longer times.

V Conclusions

We have presented an alternative method to produce thermal states with a two-site GKSL master equation for non-integrable systems. We have characterized the performance of the method by analysing the expectation values of single and two-site observables in the steady state, and we found that different observables can be consistent with those from a thermal state as long as we accept a relatively small error on its temperature.

We have also studied dynamical correlations, which give a clear indication, depending on the temperature, of how long a chain should be for its dynamical properties to look thermal enough. As expected, for larger temperature the chain can be shorter, while for smaller temperatures a much longer chain is required. The details depend on the particular parameters of the Hamiltonian.

Since it is possible to define many GKSL master equations to obtain the same target density operator $\rho_{T}$, further research could be focused on optimization of the parameters of the Hamiltonian, and/or of the dissipator to reach better thermal states. Here we have used, as target density operator, the thermal state of the Hamiltonian $H_{2S}$ in Eq. (11), however, for the fast relaxation method other target $\rho_{T}$ have been used as, for example, the thermal state for the first two sites of the Hamiltonian in Eq. (1) (Mendoza-Arenas et al., 2015), or the reduced density operator of the middle two sites of a long enough thermal chain with Hamiltonian (1) Žnidarič (2011a). When different target states are chosen, the relation between the effective inverse temperature obtained in the system (e.g., $\beta_{x}$) and the one of the bath $\beta_{B}$ will be quantitatively different from the one represented in Fig. 2$(a)$. We have run various tests using also the $\rho_{T}$ described in Mendoza-Arenas et al. (2015); Žnidarič (2011a), and we found that using Eq. 11 was consistently resulting in ratios of effective inverse temperatures closest to $1$. However a thorough investigation is required to shed light on which would be the best target state $\rho_{T}$.

Acknowledgments: We acknowledge fruitful discussions with M. Žnidarič. D.P. and M.P. acknowledge support from the Singapore Ministry of Education, Singapore Academic Research Fund Tier-II (project MOE2016-T2-1-065).

Appendix A Representation of the jump operators in the computational basis

We describe here how we write the jump operators in a way that it is conveniently implemented in our code. The operators which are implemented are the raising operator $\sigma^{+}=|\uparrow\rangle\langle\downarrow|$, the lowering operator $\sigma^{-}=|\downarrow\rangle\langle\uparrow|$, the projector on the $|\uparrow\rangle$ state $\sigma^{u}=|\uparrow\rangle\langle\uparrow|$, and the projector on the $|\downarrow\rangle$ state $\sigma^{d}=|\downarrow\rangle\langle\downarrow|$. We now write the density operator as

[TABLE]

where

[TABLE]

and where $s_{1}$, $s_{2}$, $t_{1}$, $t_{2}$ can take the values $\uparrow$ or $\downarrow$. The $\mu^{s_{1},s_{2}}_{j}$ are given by diagonalizing $\rho^{s_{1},s_{2}}_{t_{1},t_{2}}$ (e.g., after having combined the two indices $s_{1}$ and $s_{2}$ in a single one). The jump operators in the computational basis are then given by

[TABLE]

In order to get the jump operators for the code, we need to translate $|s_{1},s_{2}\rangle_{l}\langle t_{1},t_{2}|$ to lowering/raising/projection operators, which can be done, for example, from

[TABLE]

where $\sigma^{a}_{l}$ (with $a=u,d,+,-$) acts on site $l$. Since the master equation Eq. (III) uses two jump operators, we will need to consider the product of four $\sigma^{a}_{l}$ operators, for which there are $4^{4}=256$ combinations. The coefficients for these operators can be readily derived from the $16\times 16$ combinations of the coefficients $\mu^{s_{1},s_{2}}_{j}$.