Origin of Temperature Gradient in Nonequilibrium Steady States in Weakly Coupled Quantum Spin Systems

TL;DR

This paper develops an analytical perturbation theory to understand the origin of temperature gradients in nonequilibrium steady states of weakly coupled quantum spin chains, revealing conditions for nonzero temperature gradients independent of system integrability.

Contribution

It introduces a perturbation approach for weakly coupled quantum spin systems to analytically determine the conditions for temperature gradients in NESSs.

Findings

Derived a simple, integrability-independent condition for nonzero temperature gradients.

Established an analytical framework for NESSs in weakly coupled quantum spin chains.

Demonstrated the applicability of the theory to various quantum spin models.

Abstract

We study nonequilibrium steady states (NESSs) in quantum spin-1/2 chains in contact with two heat baths at different temperatures. We consider the weak-coupling limit both for spin-spin coupling in the system and for system-bath coupling. This setting allows us to treat NESSs with a nonzero temperature gradient analytically. We develop a perturbation theory for this weak-coupling situation and show a simple condition for the existence of nonzero temperature gradient. This condition is independent of the integrability of the system.

Equations

$$H$$

$${\alpha }_{zz}\ne 0$$

$$H$$

$$H_{\mathrm{SB}}$$

$$\frac{d}{dt}\rho (t)$$

$${\overline{H}}_{SB}$$

$$\mathcal{D}\rho$$

$$S_{jl}$$

$${\Phi }_{jl}(t)$$

$$H_{\mathrm{S}}$$

$$H_0$$

$$S_{jl}$$

$${\Xi }_{jl}(\omega )$$

$$\frac{d}{dt}\rho$$

$$\mathcal{L}$$

$${\mathcal{L}}_0\rho$$

$${\mathcal{L}}_1\rho$$

$${\mathcal{L}}_2\rho$$

$$\frac{d}{dt}\rho =\mathcal{L}\rho =0.$$

$$\rho$$

$${\mathcal{L}}_0{\rho }_0$$

$${\mathcal{L}}_0{\rho }_1+{\mathcal{L}}_1{\rho }_0$$

$${\mathcal{L}}_0{\rho }_2+{\mathcal{L}}_1{\rho }_1+{\mathcal{L}}_2{\rho }_0$$

$$E_{n,k}$$

$$H_0|E_{n,k}⟩$$

$$⟨E_{n,k}|H_1|E_{n,l}⟩$$

$${\mathcal{L}}_0|E_{n,k}⟩⟨E_{m,l}|$$

$$\mathtt{Ker}{\mathcal{L}}_0$$

$${\mathcal{P}}_0{\mathcal{L}}_0={\mathcal{L}}_0{\mathcal{P}}_0=0$$

$${\mathcal{P}}_0{\rho }_0={\rho }_0$$

$${\mathcal{P}}_0{\mathcal{L}}_1{\rho }_0={\mathcal{P}}_0{\mathcal{L}}_1{\mathcal{P}}_0{\rho }_0=0.$$

Full text

I Introduction

The microscopic description and understanding of nonequilibrium steady states (NESSs) is one of the most challenging problems in physics. For heat-conducting systems, Fourier’s law states that the heat current is proportional to the temperature gradient. Although it is a very universal empirical law, there has been no microscopic derivation of it.

In a normal heat-conducting NESS, we observe a uniform nonzero temperature gradient and heat current. However, when considering a simple theoretical model, we often observe a flat temperature profile in the system, which represents ballistic transport. It is commonly believed that normal heat-conducting states are realized in non-integrable systems.lepri However, there are some exceptions to this rule. prosen ; michel ; wu ; mendoza Therefore, the actual condition for realizing normal NESSs is still an open problem.

In this paper, we consider heat-conducting NESSs in one-dimensional quantum spin chains in contact with two heat baths. In the study of such NESSs, the coupling between the system and the baths is often assumed to be weak. yuge ; sugita However, if the system-bath coupling is weak and the spin-spin interaction is not weak, the system will have a flat temperature profile because the thermal resistance at the boundary is large compared with that in the system (Fig. 1). To study NESSs with nonzero temperature gradients, we consider the weak-coupling limit both for spin-bath coupling and for spin-spin coupling. In this limit, the thermal resistance uniformly becomes infinitely large both in the system and at the boundary. Therefore, we expect a “frozen” NESS to appear, which has a nonzero temperature gradient and no current. We will show that the “frozen” solution is not just a tensor product of local equilibrium states.

Heat-conducting states with nonzero temperature gradients are difficult to study analytically, because they are usually realized in non-integrable systems. However, our setting allows some analytical treatment of such states. In this paper, we consider one-dimensional uniform quantum spin-1/2 systems with nearest-neighbor interactions, whose Hamiltonians are of the form

[TABLE]

Here, $\sigma_{l}^{j}\;(j=x,y,z)$ represents a Pauli matrix at the $l$th site. We assume $h\neq 0$ and $\alpha_{ij}=\alpha_{ji}$. Then we show that

[TABLE]

is necessary to form a nonzero temperature gradient. Note that this condition is independent of the integrability of the system. For example, it is consistent with the observation that the $XXZ$ model has normal conducting states in spite of its integrabilityprosen ; michel ; wu ; mendoza . We can understand intuitively that Eq. (2) must be satisfied so that the local temperature at a site can affect the temperatures of neighboring sites, because in the weak-coupling limit the local temperature is determined only from $\langle\sigma_{z}\rangle$. We also show that some three-body correlations play important roles in forming the temperature gradient. The forms of the correlation terms are universal, i.e., they do not depend on the form of the interaction.

This paper is organized as follows. In Sect. II, we derive the quantum master equation (QME) in the weak-coupling limit. We show that the heat bath superoperator becomes local when the spin-spin coupling is weak. In Sect. III, we derive a perturbative expansion of the stationary solution of the QME. In Sect. IV we study quantum spin-1/2 chains. First we analyze two-spin cases and show an explicit example of a frozen NESS with a nonzero temperature gradient. Then we study $N$-spin systems and derive Eq. (2). We also show some numerical results to verify our theoretical considerations. Section V is devoted to a summary.

II Quantum Master Equation

In this section, we derive the QME for weakly coupled spin systems. For simplicity, we first derive the QME for a system coupled to a single heat bath. Generalization to multiple heat baths is straightforward. We denote the Hilbert space of the system by $\mathsf{H}$ and assume that its dimension is finite. The set of all linear operators on $\mathsf{H}$ is denoted by $\mathsf{L}$. It becomes a Hilbert space by introducing the Hilbert–Schmidt inner product as ${\rm Tr}\left(X_{1}X_{2}\right)$ for any $X_{1},X_{2}\in\mathsf{L}$. We refer to a linear operator on $\mathsf{L}$ as a superoperator.

We start with the Hamiltonian

[TABLE]

where $\lambda$ is a small coupling parameter, $H_{\rm S}$ is the Hamiltonian for the system, and $H_{\rm B}$ is that for the heat bath. $H_{\rm SB}$ represents the interaction between the system and the heat bath, which can be written as

[TABLE]

Here, $X_{j}$ and $Y_{j}$ are Hermitian operators on the system and on the heat bah, respectively. By using the Born-Markov approximationkubo ; breuer , we obtain the standard Redfield-type QME

[TABLE]

Here, $\overline{H}_{\rm SB}\in\mathsf{L}$ represents the interaction averaged over the heat bath, whose explicit form is

[TABLE]

where $\overline{Y}_{j}={\rm Tr}\left\{\rho_{\rm B}Y_{j}\right\}$ and $\rho_{\rm B}$ is the equilibrium state of the heat bath $\rho_{\rm B}=e^{-\beta H_{\rm B}}/Z_{\rm B}$. $\mathcal{D}$ is the heat bath superoperator

[TABLE]

where

[TABLE]

$\Phi_{jl}$ is a correlation function for the heat bath,

[TABLE]

where $\Delta Y_{k}(t)=Y_{k}(t)-\overline{Y}_{k}$. Note that Eq. (5) is correct up to $O(\lambda^{2})$. We also use the Markovian approximation in deriving Eq. (5), which is correct when the density matrix in the interacting picture $e^{iH_{\rm S}t/\hbar}\rho(t)e^{-iH_{\rm S}t/\hbar}$ is slowly varying compared with the correlation time of the heat bath kubo . This assumption holds for NESSs in a weakly nonequilibrium regime. In our perturbation theory, the zeroth order solution satisfies this condition exactly even when it is far from equilibrium.

Now we consider a weakly coupled one-dimensional spin system with nearest-neighbor interactions. We assume that the spin-spin interaction is of the same order as the system-bath interaction. Then the system Hamiltonian can be written as

[TABLE]

where $H_{0}$ represents the external field applied to the spins and $H_{\rm I}$ is the interaction among the spins. We also assume that the $X_{j}$ in $H_{\rm SB}$ act only on the edge spin (Fig. 2) and write

[TABLE]

where $H_{\rm e}$ acts on the edge spin and $H_{0}^{\prime}$ acts on the other spins. Since $[H_{0}^{\prime},X_{j}]=0$, we have

[TABLE]

The $O(\lambda)$ term is negligible in the QME in Eq. (5), and now $S_{jl}$ is a local operator acting on the edge spin. It can also be written as

[TABLE]

where $\Xi_{jl}$ is the Fourier–Laplace transform of the correlation function

[TABLE]

$E_{r}$ is an eigenvalue of $H_{\rm e}$ and $|E_{r}\rangle$ is the corresponding eigenvector.

If we place heat baths at both ends of the spin chain, the QME for this weak-coupling limit is

[TABLE]

[TABLE]

[TABLE]

where $H_{1}=H_{\rm I}+\overline{H}_{\rm SB}$, and $\mathcal{D}_{\rm L}$ and $\mathcal{D}_{\rm R}$ are heat bath superoperators for the left and right heat baths, respectively. They are of the form given in Eqs. (7) and (13).

III Perturbative Expansion

III.1 Zeroth-order solution

The steady state is determined by the equation

[TABLE]

We calculate the steady state by expanding $\rho$ with respect to the coupling parameter $\lambda$:

[TABLE]

By substituting it into Eq. (20), we obtain the following equations up to the second order:

[TABLE]

In usual perturbation theories, the zeroth-order solution is determined from the zeroth-order equation. In this case, however, we cannot determine $\rho_{0}$ simply from Eq. (22), and the first order-equation Eq. (23) also fails to determine $\rho_{0}$. Hence we need all three equations Eqs. (22)-(24) to determine the zeroth-order solution. Physically speaking, this is because the Liouvillian $\mathcal{L}_{0}+\lambda\mathcal{L}_{1}$ represents Hamiltonian dynamics governed by $H_{0}+\lambda H_{1}$, and all energy eigenstates of $H_{0}+\lambda H_{1}$ are stationary in this dynamics. We require the heat baths, which are represented by $\mathcal{L}_{2}$, to make the stationary state unique.

Let us consider the eigenvalue problem of the Hamiltonian $H_{0}+\lambda H_{1}$ up to the first order of $\lambda$. Since $H_{0}$ is a non-interacting Hamiltonian and highly degenerate, we must use degenerate perturbation theory. Then we obtain the eigenvalues

[TABLE]

where $E_{n}$ is an eigenvalue of $H_{0}$ and $e^{(n)}_{k}$ is obtained by diagonalizing $H_{1}$ in the eigenspace belonging to $E_{n}$. Therefore, the corresponding eigenvector $|E_{n,k}\rangle$ satisfies

[TABLE]

and

[TABLE]

We assume the $E_{n,k}$ are nondegenerate.

Since $|E_{n,k}\rangle\langle E_{m,l}|\in\mathsf{L}$ satisfies the equation

[TABLE]

it is an eigenoperator of $\mathcal{L}_{0}$ with the eigenvalue $\frac{E_{n}-E_{m}}{i\hbar}$, and the operators of this type form a complete set in $\mathsf{L}$. The zero eigenspace is spanned by the operators with $n=m$:

[TABLE]

We denote the projection superoperator to $\mathtt{Ker\,}\mathcal{L}_{0}$ by $\mathcal{P}_{0}$. Then we have

[TABLE]

and

[TABLE]

because $\rho_{0}\in\mathtt{Ker\,}\mathcal{L}_{0}$. By applying $\mathcal{P}_{0}$ to Eq. (23), we obtain

[TABLE]

Hence, $\rho_{0}$ is in $\mathtt{Ker\,}\left(\mathcal{P}_{0}\mathcal{L}_{1}\mathcal{P}_{0}\right)$, which is spanned by the diagonal elements:

[TABLE]

We denote the projection superoperator to $\mathtt{Ker\,}\left(\mathcal{P}_{0}\mathcal{L}_{1}\mathcal{P}_{0}\right)$ by $\mathcal{P}_{1}$. Then

[TABLE]

We also define projection superoperators

[TABLE]

Note that

[TABLE]

and $\mathtt{Im\,}\mathcal{P}_{1}\subset\mathtt{Im\,}\mathcal{P}_{0}$, $\mathtt{Im\,}\mathcal{Q}_{1}\subset\mathtt{Im\,}\mathcal{P}_{0}$.

Applying $\mathcal{Q}_{0}$ to Eq. (23), we obtain

[TABLE]

Hence,

[TABLE]

Note that the superoperator $\left(Q_{0}\mathcal{L}_{0}Q_{0}\right)^{-1}$ is defined on ${\tt Im}Q_{0}$. We implicitly assumed in Eq. (39) that $\left(Q_{0}\mathcal{L}_{0}Q_{0}\right)^{-1}\rho=0$ if $\rho\notin{\tt Im}Q_{0}$. Hence,

[TABLE]

By applying $\mathcal{P}_{1}$ to Eq. (24), we obtain

[TABLE]

Since $\mathcal{P}_{1}\mathcal{L}_{1}\mathcal{P}_{0}=\mathcal{P}_{1}\mathcal{P}_{0}\mathcal{L}_{1}\mathcal{P}_{0}=0$, we have

[TABLE]

We can show that

[TABLE]

by a straightforward algebra (see AppendixA). Hence, we obtain the following three equations to determine the zeroth-order solution:

[TABLE]

III.2 Higher-order solutions

In a similar way, we can obtain higher-order solutions. If the Liouvillian is expanded with respect to $\lambda$ as

[TABLE]

the $m$th-order equation is

[TABLE]

We decompose the $m$th-order solution as

[TABLE]

and represent the three terms by lower-order solutions in the following.

First we apply $\mathcal{Q}_{0}$ to Eq. (49). Then we have

[TABLE]

Hence,

[TABLE]

Then we apply $\mathcal{Q}_{1}$ to the $(m+1)$th-order equation. Since $\mathcal{Q}_{1}\mathcal{L}_{0}=0$, we have

[TABLE]

Because

[TABLE]

we obtain

[TABLE]

Next we apply $\mathcal{P}_{1}$ to the $(m+2)$th-order equation. Then we have

[TABLE]

The first term is

[TABLE]

where we used Eq. (52) with $m$ replaced by $m+1$ in the second line and Eq. (44) in the last line. Therefore,

[TABLE]

The superoperator $\mathcal{P}_{1}\mathcal{L}_{2}\mathcal{P}_{1}$ is not invertible since it has a one-dimensional kernel spanned by $\rho_{0}$. Therefore, we define $\left(\mathcal{P}_{1}\mathcal{L}_{2}\mathcal{P}_{1}\right)^{-1}$ as the Moore-Penrose pseudoinverse of $\mathcal{P}_{1}\mathcal{L}_{2}\mathcal{P}_{1}$. Then we have

[TABLE]

where the constant $c$ is determined by the condition

[TABLE]

Combining Eqs. (52), (55), and (62), we obtain the $m$th-order solution $\rho_{m}$. For example, the first-order solution is given by

[TABLE]

[TABLE]

Here, we give some comments about the range of applicability of our perturbation theory. Our method has two parameters, the system size and the coupling constant, which compete with each other. If the system size is large, $\mathcal{Q}_{1}\rho_{m}$ may be very large because we have $(\mathcal{Q}_{1}\mathcal{L}_{1}\mathcal{Q}_{1})^{-1}$ in Eq. (55), whose eigenvalues may be very large. Then our method could work only for a very small coupling constant $\lambda$ and might be practically useless for large systems. However, it is also possible that even when the higher-order terms of the density matrix are large, they give only small contributions to local physical quantities such as local temperature and current. (Note that very different states, such as a pure state and the Gibbs state, can give the same expectation values for local observables in macroscopic systems sugita2 ; sugiura .) Then our method may be useful even for large systems.

To consider this point, we will give some numerical results for spin systems in the next section. The results show that, for local temperatures, the zeroth-order approximation is equally good for the system sizes $N=2$ and $N=6$. Therefore, we hope that our method is useful for larger systems.

Another important point to note is that a NESS in a rather small system may be similar to a part of a macroscopic NESS, at least in some aspects. Thus, we will be able to obtain some insights on macroscopic NESSs by studying those in small systems, where our method will be helpful.

IV Spin-1/2 Chain

In this section, we consider a quantum $N$-spin chain with the Hamiltonian

[TABLE]

[TABLE]

In this section we set $\hbar=1$.

IV.1 Basis of $\mathtt{Ker\,}\mathcal{L}_{0}$

We represent spin up and down states by $|1\rangle$ and $|-1\rangle$, respectively. Then $H_{0}$ has eigenvectors of the form

[TABLE]

where $s_{l}=\pm 1$, and its eigenvalue is $h\sum_{l}s_{l}$. $\mathtt{Ker\,}\mathcal{L}_{0}$ is spanned by elements of the form

[TABLE]

with

[TABLE]

Note that the number of sites with $s_{l}\neq s^{\prime}_{l}$ must be even to satisfy this condition.

Let us consider some simple cases. If $N=1$, $\mathtt{Ker\,}\mathcal{L}_{0}$ is spanned by $|1\rangle\langle 1|$ and $|-1\rangle\langle-1|$. In this case, only diagonal terms appear. Therefore, in this case, $\mathtt{Ker\,}\mathcal{L}_{0}$ is spanned by $\sigma^{0}$ and $\sigma^{z}$, where $\sigma^{0}$ is the $2\times 2$ unit matrix.

If $N=2$, $\mathtt{Ker\,}\mathcal{L}_{0}$ is spanned by six elements, in which there are four diagonal elements,

[TABLE]

and two off-diagonal elements,

[TABLE]

The diagonal terms can be represented by linear combinations of $\sigma_{\alpha}\otimes\sigma_{\beta}$, where $\alpha,\beta$ are [math] or $z$. The off-diagonal terms are represented as

[TABLE]

where

[TABLE]

Therefore, the off-diagonal terms can be represented by the two terms $\sigma_{1}^{x}\otimes\sigma_{2}^{x}+\sigma_{1}^{y}\otimes\sigma_{2}^{y}$ and $\sigma_{1}^{x}\otimes\sigma_{2}^{y}-\sigma_{1}^{y}\otimes\sigma_{2}^{y}$.

For the $N$-spin case, a basis element in Eq. (72) has an even number of off-diagonal sites where $s_{l}\neq s^{\prime}_{l}$. The diagonal parts can be represented by combining $\sigma_{0}$ and $\sigma_{z}$, and the off-diagonal parts can be obtained by combining the following two types of correlation terms:

[TABLE]

IV.2 Zeroth-order solution for the two-spin case

Let us consider the two-spin case. In this case $\mathtt{Ker\,}\mathcal{L}_{0}$ is six dimensional, and we take the basis as

[TABLE]

Here we used the notation

[TABLE]

$\rho^{\rm A}$ and $\rho^{\rm B}$ are correlation terms in the form of Eqs. (79) and (80), respectively, with $l=1$ and $m=2$. $\rho^{\rm C}$ is defined as

[TABLE]

Then we consider Eq. (46). $\mathcal{L}_{1}$ is represented by two terms,

[TABLE]

where

[TABLE]

By straightforward algebra, we obtain

[TABLE]

[TABLE]

where $\alpha_{+}\coloneqq\alpha_{xx}+\alpha_{yy}$.

We write

[TABLE]

which represents the first-order effect of the system-bath coupling from the left and right heat baths. Then we have

[TABLE]

[TABLE]

where $\Delta\alpha^{z}\coloneqq\alpha_{\rm L}^{z}-\alpha_{\rm R}^{z}$. Therefore, we have

[TABLE]

where

[TABLE]

[TABLE]

Let us consider the local equilibrium form

[TABLE]

Here, $\rho_{\rm LE}$ represents a single-spin equilibrium state, whose explicit form is

[TABLE]

where $t_{l}\coloneqq\tanh(\beta_{l}h)$. Then we apply $\mathcal{P}_{0}\mathcal{L}_{1}$ to Eq. (99). Since

[TABLE]

we have

[TABLE]

Therefore, we see that the local equilibrium form in Eq. (99) is not the zeroth-order solution unless $t_{1}=t_{2}$. The general form of the solution of the equation $\mathcal{P}_{0}\mathcal{L}_{1}\rho_{0}=0$ is

[TABLE]

where the three parameters $t_{1},t_{2}$, and $g$ are to be determined from $\mathcal{P}_{1}\mathcal{L}_{2}\rho_{0}=0$. Note that the correlation term $\rho^{\rm A}$ appears when there is a temperature gradient in the system.

Let us consider the tilted Ising model as an example:

[TABLE]

This model is known to be integrable at $\theta=\pi/2$ (transverse Ising model), and the NESS has a flat temperature profile in this case saito1 .

By taking the $z$-axis in the direction of the external field, the Hamiltonian can be transformed as

[TABLE]

where

[TABLE]

The interaction term can also be written as

[TABLE]

Therefore, we have

[TABLE]

We take the spin-bath coupling as having the same form as the spin-spin coupling (Fig. 3). This situation may be regarded as similar to that for the two spins in a NESS of a long spin chain. We also assume that the temperature difference between the two heat baths is small. Then we obtain

[TABLE]

where $t_{l}=\tanh(\beta_{l}h)\;\;(l={\rm L,R})$, and $\beta_{\rm L}$ ($\beta_{\rm R}$) is the inverse temperature for the left (right) heat bath. $K$ is defined as

[TABLE]

where $c$ and $d$ are non-negative constants. (See AppendixB for details of the calculation.) We see that a nonzero temperature gradient is formed unless $\alpha_{zz}=\cos^{2}\theta=0$, which corresponds to the transverse Ising model.

Figure 4 shows the temperature profiles for the tilted Ising model with $N=2$. The numerical results are obtained by diagonalizing the Liouvillian in Eq. (16). We set $h=1$ and take the heat bath parameters as $c_{L}=c_{R}=c_{L}=c_{R}=1$ in Eqs. (191) and (192). Then Eqs. (111) and (112) should hold exactly in the limit $\lambda\rightarrow 0$. In the case with $\theta=\pi/4$, the numerical results are in good agreement with the analytical results for $\lambda=0.01$ with relative errors of order $10^{-2}$. For the case with $\theta=\pi/2$ (transverse Ising model), where the temperature gradient vanishes, the convergence to the weak-coupling limit is much faster.

IV.3 $N$-spin case

Next we consider $N$-spin systems. Let us first consider the product of local equilibrium states:

[TABLE]

We write $\mathcal{L}_{1}$ as

[TABLE]

where $\mathcal{L}_{{\rm I}(l,l+1)}$ represents a local interaction,

[TABLE]

[TABLE]

and $\mathcal{L}_{{\rm E}(l)}$ represents the first-order contribution from a heat bath,

[TABLE]

[TABLE]

Then we have

[TABLE]

and

[TABLE]

Since $\mathcal{P}_{0}\mathcal{L}_{1}\rho_{\rm PLE}\neq 0$, $\rho_{\rm PLE}$ is not the zeroth-order solution. Therefore we have to add some correction terms to $\rho_{\rm PLE}$ to obtain $\rho_{0}$. To cancel the two-body correlated term in Eq. (121), we have to add three-body correlated terms, but if we apply $\mathcal{P}_{0}\mathcal{L}_{1}$ to them, slightly different residual three-body terms appear. Then we have to add four-body correlated terms to cancel the three-body residual terms, which give new four-body residual terms, and so on. Thus, we obtain a BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon)-like hierarchy of equations, which is difficult to solve.

In this paper, we do not try to obtain $\rho_{0}$ explicitly. Instead, we consider local projections of the equation $\mathcal{P}_{0}\mathcal{L}_{1}\rho_{0}=0$. A one-site projection $\mathcal{P}_{(l)}$ traces out all spins except for that on the $l$th site. It gives no information, however, since $\rho_{\rm PLE}$ already satisfies

[TABLE]

Next we consider a two-site projection $\mathcal{P}_{(l,l+1)}$, which traces out all spins except for those at the $l$th and $(l+1)$th sites. Since

[TABLE]

if $l+1<m$ or $m+1<l$, only three interaction terms, $\mathcal{L}_{{\rm I}(l-1,l)}$, $\mathcal{L}_{{\rm I}(l,l+1)}$, and $\mathcal{L}_{{\rm I}(l+1,l+2)}$ matter when we consider $\mathcal{P}_{(l,l+1)}$. Hence, the equation to solve is

[TABLE]

We write

[TABLE]

Without loss of generality, we can assume that the additional term $\rho_{a}$ does not change the local state:

[TABLE]

Then the general form of $\rho_{a}$ that is relevant to Eq. (125) is

[TABLE]

In the coefficients $c^{\alpha}_{p,q,r}$, the overline attached to a suffix shows the site where $\sigma^{z}$ is. We assume that, at a site whose state is unspecified, the state is $\sigma^{0}/2$, whose trace is unity. For example,

[TABLE]

Substituting Eq. (128) into Eqs. (125) and (126), we obtain

[TABLE]

where $\rho^{\rm C}_{l,m}\coloneqq\sigma^{0}_{l}\sigma^{z}_{m}-\sigma^{z}_{l}\sigma^{0}_{m}$. Hence,

[TABLE]

At the edges of the system, we have slightly different equations. For example, at the left edge ($l=1$), we have

[TABLE]

The general form of the relevant additional term is

[TABLE]

Substituting Eq. (135) into Eq. (134), we obtain

[TABLE]

Hence,

[TABLE]

In the same way, we obtain the following three equations at the right edge:

[TABLE]

Thus, we have obtained the following set of equations related to the temperature gradient:

[TABLE]

Summing these equations, all terms with $\alpha_{+}$ vanish and we obtain

[TABLE]

If we assume that the system-bath coupling has the same form as the spin-spin coupling as in Sect. IV.2, we have $\alpha_{1}^{z}=-\alpha_{zz}t_{\rm L}$ and $\alpha_{N}^{z}=-\alpha_{zz}t_{\rm R}$. Replacing the coefficients with the expectation values of the correlations, we obtain

[TABLE]

This is our main result. The left-hand side (LHS) is related to the temperature difference between the two ends of the chain. Here, we see that $\alpha_{zz}\neq 0$ is necessary to form a nonzero temperature gradient. We can also see that the three-body correlations

[TABLE]

and

[TABLE]

play essential roles in forming the temperature gradient. Note that the two-body correlation terms related to the heat baths also have similar forms since $t_{l}=-\langle\sigma_{l}\rangle$.

In Fig. 5 we show numerically obtained temperature profiles for the tilted Ising model with $N=6$. We see that the temperature gradient is nonzero for $\theta=\pi/4$, where $\alpha_{zz}\neq 0$, while it vanishes for $\theta=\pi/2$, where $\alpha_{zz}=0$, as expected from Eq. (150). In the case with $\theta=\pi/4$, the numerical results agree with the weak-coupling limit for $\lambda=0.01$, with relative errors of order $10^{-2}$. The convergence is as fast as for the two-spin case. For $\theta=\pi/2$, the convergence is much faster, which is also the same as in the two-spin case.

In Fig. 6 we show numerical results for the LHS and the right-hand side (RHS) of Eq. (150). We see that both sides converge to the same value in the weak-coupling limit, although the convergence is slightly slower than that for the temperature.

V Summary

We considered NESSs in weakly coupled one-dimensional quantum spin-1/2 chains. The spin-spin coupling and system-bath coupling were assumed to be of the same order. We derived the QME for this case, in which the heat bath superoperator becomes local. We developed a perturbative expansion for the stationary solution of the QME. For the two-spin case, we gave an explicit example in which the zeroth-order solution has a nonzero temperature gradient. For the general $N$-spin case, we showed that $\alpha_{zz}\neq 0$ is necessary to obtain a nonzero temperature gradient. We also showed that some three-body correlations are essential for the formation of a temperature gradient.

Acknowledgements.

The authors thank Akira Terai, Katsuhiro Nakamura, and Shogo Tanimura for helpful discussions and comments.

Appendix A Proof of Eq. (44)

Let us denote a matrix element of $H_{1}$ as

[TABLE]

It satisfies

[TABLE]

We apply $\mathcal{L}_{1}\left(Q_{0}\mathcal{L}_{0}Q_{0}\right)^{-1}\mathcal{L}_{1}\mathcal{P}_{1}$ to a basis element $|E_{n,k}\rangle\langle E_{n,k}|$ in $\mathtt{Im\,}\mathcal{P}_{1}$. Then we obtain

[TABLE]

Applying $\mathcal{P}_{1}$ to the above equation, we have

[TABLE]

Hence, Eq. (44) is proven.

Appendix B Two-spin tilted Ising model

In this appendix, we show the details of the calculation of the zeroth-order solution for the two-spin tilted Ising model.

First we consider a single bath and assume spin-bath coupling of the form $\sigma^{z^{\prime}}\otimes\sigma^{z^{\prime}}$, where

[TABLE]

Then the heat bath superoperator becomes

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Here, $\langle\cdot\rangle_{\beta}$ represents the equilibrium average in the heat bath at the inverse temperature $\beta$. We removed suffices $j,l$ from $\Xi$ and $\Phi$ since the system-bath coupling is represented by a single term, and use the subscript $\beta$ to show the inverse temperature of the heat bath explicitly. Note that we take $\hbar=1$ here.

$\Xi_{\beta}$ can be written as

[TABLE]

where the real part is represented by the Fourier transform of $\Phi_{\beta}(t)$:

[TABLE]

It satisfies the KMS (Kubo-Martin-Schwinger) condition kubo

[TABLE]

and $\tilde{\Phi}(\omega)\geq 0$. We ignore the imaginary part $\Psi_{\beta}$ because usually it does not give a significant contribution.saito2

Then

[TABLE]

where

[TABLE]

Thus, we have

[TABLE]

Next we apply $\mathcal{P}_{1}\mathcal{L}_{2}=\mathcal{P}_{1}(\mathcal{D}_{\rm L}+\mathcal{D}_{\rm R})$ to Eq. (103). Then we obtain

[TABLE]

with obvious notations. Since in this case

[TABLE]

we have

[TABLE]

and

[TABLE]

We can determine $t_{1},t_{2}$, and $g$ by solving the following three equations:

[TABLE]

[TABLE]

When $t_{\rm L}-t_{\rm R}\ll 1$, $a_{\rm L}\simeq a_{\rm R}$ and $b_{\rm L}\simeq b_{\rm R}$. Hence, we obtain

[TABLE]

ignoring higher-order terms with respect to $t_{\rm L}-t_{\rm R}$. Here,

[TABLE]

and

[TABLE]

where

[TABLE]

and $\alpha={\rm L},{\rm R}$. Note that $a_{\alpha}=c_{\alpha}\sin\theta$ and $b_{\alpha}=d_{\alpha}\cos\theta$. Since $c,d\geq 0$, $K$ satisfies $0\leq K\leq 1$. $K=0$ holds for $\theta=\pi/2$ (horizontal magnetic field), and $K=1$ for $\theta=0$ (vertical magnetic field).