Steady state thermodynamics of two qubits strongly coupled to bosonic environments

When a quantum system is placed in thermal environments, we often assume that the system relaxes to the Gibbs state in which decoherence takes place in the system energy eigenbasis. However, when the coupling between the system and the environments is strong, the stationary state is not necessarily the Gibbs state due to environment-induced decoherence which can be interpreted as continuous measurement by the environments. Based on the einselection proposed by Zurek, we postulate that the Gibbs state is projected onto the pointer basis due to the continuous measurement. We justify the proposition by exact numerical simulation of a pair of coupled qubits interacting with boson gases. Furthermore, we demonstrate that heat conduction in non-equilibrium steady states can be suppressed in the strong coupling limit also by the environment-induced decoherence.


I. INTRODUCTION
The laws of thermodynamics and the principles of statistical mechanics tell us that every system eventually reaches a stationary state known as the Gibbs state, which is the hallmark of thermal equilibrium. The density operator of the Gibbs state is notably a function of only the system Hamiltonian and is thus diagonal in the energy eigenbasis. The coherence between energy eigenstates is completely destroyed. Therefore, thermalization to the Gibbs state must involve decoherence between energy eigenstates, presumably induced by the environments surrounding the system. Such a decoherence process toward the Gibbs state has been investigated under the weak coupling limit [1]. In fact, quantum master equations based on the Born-Markovian approximation are known to converge to the Gibbs state. [2]. However, it has been shown that the non-Markovian dynamics does not necessarily reach the Gibbs state. [3][4][5][6][7][8] For a system strongly coupled to the environments, its equilibrium state cannot be expressed with the system Hamiltonian alone, and an effective Hamiltonian based on the potential of mean force has been developed. [9][10][11][12][13][14][15][16] The resulting stationary state is no longer diagonal in the system energy eigenstates.
Environment-induced decoherence has been intensively investigated in the context of quantum measurement theory and quantum computing. [17] In those theories, the environment does not necessarily induce decoherence in the energy eigenbasis. Zurek [18,19] showed that the decoherence takes place among so-called "pointer states" determined by the coupling Hamiltonian between a system and environments. In general, the system density operator becomes diagonal in the pointer basis under the strong coupling limit. This einselection [18] can be considered as a consequence of continuous measurement of the system by the environment. A similar argument can be used for the thermalization processes, and there have been investigation of thermalization under continuous measurement. [20,21] We investigate thermalization and heat conduction in the strong coupling regime based decoherence in the pointer basis.

II. THERMALIZATION IN THE POITER BASIS
Consider a system in the Gibbs state ρ g s = e −βHs /Z s under the weak coupling, where H s β, and Z s are system Hamiltonian, inverse temperature and a partition function. When the coupling energy becomes significantly larger than the system energy, the Gibbs state is continuously measured by the environments and thus projected to the pointer basis. Our main proposition is that under the strong coupling limit a system tends to relax to a  (1). The Gibbs state on the convex hull Σe is projected onto another convex hull Σp. As the coupling strength increases, the steady state deviates from the Gibbs state (G) along the projection line toward the pointer limit (P ). The maximal entropy state (I) is located on the intersection of the two convex hulls. Noting that P is closer to I than G, the entropy increases as the steady state moves toward the pointer limit. stationary state given by where |p i is the i-th pointer state which we define below. Figure 1 illustrates this proposition. Consider the con- in the Liouville space. The corners of the hull represent the pure states. Any density operator that is diagonal in the energy eigenbasis |e i is in Σ u , including the Gibbs state (G in Fig 1). Similarly, the convex hull contains all possible density operators that are diagonal in the pointer basis |p i . The density operators in the intersection of the two convex hulls are diagonal in both basis sets. A special point I in the figure corresponds to ρ = 1 ds I s where I s is an identity operator and d s is the dimension of the system Hilbert space. The entropy of the system reaches its maximum value ln d s at I. As the coupling gets stronger, the steady state deviates from the Gibbs state (G) toward the pointer limit (P ) along the projection line (GP ). The projection line is "perpendicular" to Σ p , meaning that the diagonal elements in the pointer basis are invariant along the projection line.

III. MODEL AND NUMERICAL SIMULATION
We justify the proposition by numerically investigating the exact dynamics of a simple spin-boson model. Following the standard open quantum system approach [2], we consider an isolated system consisting of a small subsystem H s and environments H b . The unitary evolution of the total system follows the Liouville-von Neumann equation where H b is the Hamiltonian of environment. For simplicity, we assume that the coupling Hamiltonian takes a bilinear form where X and Y are operators in H s and H b , respectively. Furthermore, we assume that [X k , X ] = 0 so that all X share the same eigenkets |p j which we shall call pointer states. If there are degenerate subspaces, we choose a particular basis in the subspace such that the steady state becomes diagonal in the pointer basis. The state of the system is represented by reduced density ρ s = Tr b ρ sb which obeys the equation of motion where we introduced a new operator, Note that the time evolution of the system needs only limited information on the state of the environments through η . In order to demonstrate the proposition, we consider a simple model consisting of a pair of identical qubits S 1 and S 2 whose Hamiltonian is given by where σ z,± , ( = 1, 2) are usual Pauli matrices for theth qubit, and ω 0 and λ s the qubit excitation energy and the internal coupling strength, respectively. We write the energy eigenstates as |e j , (j = 1, · · · , 4) with eigenvalue e j starting from the ground state. Each qubit S is coupled to its own environment B . [22] The environments are assumed to be ideal Bose gases whose Hamiltonians are given by H b = k ω (k) a † (k)a (k), where a † (k) and a (k) are creation and annihilation operators for the k-th mode in B . The interaction Hamiltonian between S and B assumes a simple bilinear form X ⊗ Y where X = σ x and Y = k (k) a † (k) + a (k) . The coupling strength between the system and the k-th mode in B is denoted as (k).
The pointer states in this model are the simultane-ous eigenkets of X 1 and X 2 and denoted as |p 1 = |0 0 , |p 2 = |0 1 , |p 3 = |1 0 , and |p 4 = |1 1 , where |0 and |1 are the eigenkets of σ x . When the coupling is weak, the stationary state is the Gibbs state where ρ e jj = e −βej / i e −βei . Under the strong coupling limit, Proposition (1) claims that the stationary state density is given by where ρ p jj = p j |ρ g s |p j /Z s can be explicitly expressed as Now we show the transition from the Gibbs limit (7) to the pointer limit (8) by numerically solving Eq. (4). Assuming that the total system is initially in a product state we obtain a formally exact expression of the system density in the interaction picture [23] where the super operator K j is defined by (11) with anti(+) and regular(-) commutators S ± = [X , ·] ± .
The dissipation kernel K (d) (t) and noise kernel K (n) (t) are respectively the real and imaginary part of the correlation function C (t) = Y b (t)Y b (t 0 ) 0 where the expectation value is taken with the initial environment state ρ b (t 0 ). The time ordering operator ← − T in Eq. (10) chronologically orders the super operators S ± (t).
Kato and Tanimura [24] showed that Eq (10) can be numerically evaluated if the spectral density of environments is of the Drude-Lorentz type: where γ and λ are the response rate of environment and the overall coupling strength between qubit S and environment B , respectively. Then, the environmental correlation can be expressed with reasonable accuracy as [25] C (t) ≈ λ c e −γ + 2∆ δ(t) where c = 2/β − γ ∆ − iγ and ∆ = γ β /6. Following Kato and Tanimura [24], we introduce a set of auxiliary operators where Index n associated with environment B runs from 0 through infinity. Only the first three lowest order auxiliary operators are needed for ρ s (t) = ζ 0,0 (t), However, the dynamics of auxiliary operators is determined by an infinite set of coupled ODEs or so-called hierarchical equations of motion [24] (16) with the initial condition ζ n1,n2 (t 0 ) = 0 except for ζ 0,0 (t 0 ) = ρ s (t 0 ). The infinite hierarchy is truncated at depth d = 50 such that higher depth auxiliary operators do not significantly contribute to the first two depths.

IV. RESULTS AND DISCUSSION
First, we investigate the equilibrium situation where the initial states of the two environments are identical (λ 1 = λ 2 ≡ λ b , T 1 = T 2 ≡ T , γ 1 = γ 2 ≡ γ). We tried more than ten different initial densities, and all converged to the same stationary state. In Fig. 2, the matrix elements of the stationary state density are plotted as a function of the coupling strength λ b using the energy eigenbasis and the pointer basis. The density matrix in the energy eigenbasis shows that the Gibbs state is realized only at the weak coupling limit. The diagonal elements deviate from the Gibbs state as the coupling increases. The off-diagonal elements indicate that the superposition of eigenstates |e 1 and |e 4 grows rapidly and thus decoherence does not fully take place in the energy eigenbasis. Both the diagonal and off-diagonal elements approach the pointer limit predicted by Eq. (8).
When the matrix elements of the same density operator are evaluated in the pointer basis, all of the offdiagonal elements tend to vanish as the coupling strength increases, suggesting that full decoherence takes place in the pointer basis. The diagonal elements are remarkably insensitive to the coupling strength and in good agreement with Eq. (9) regardless of the coupling strength. The invariance of the diagonal elements confirms that the projection is perpendicular to the convex hull Σ p . (See Fig. 1.) In Fig. 3 the deviation of the steady state from the Gibbs state and its approach to the pointer limit are measured by fidelity F (ρ, ρ ) = tr{ √ ρ ρ √ ρ} 2 . At λ b = 4, the distance between the steady state and the pointer limit nearly vanishes. Through the continuous measurement, the environments gain information of the system and the system loses information. Accordingly, the entropy of the system increases. [18,26] As the coupling gets stronger, more information is expected to be lost and thus the entropy goes up monotonically. Figure 3 confirms the increase of the von Neumann entropy which converges to the pointer limit (8) at the strong coupling limit.
As further evidence of continuous measurement by environments, we also investigated a non-equilibrium steady state. When different temperatures are used, heat flows through the system. Heat from the environment B to the system can be computed as Figure 4 shows the steady state heat current as a function of the coupling strength. In the weak coupling regime, the current increases linearly as expected from the linear response theory. However, the heat current reaches its maximum and dies off rather quickly as the coupling becomes stronger. This suppression of heat is predicted earlier as a consequence of the quantum zeno effect [27] based on a heuristic argument and is observed by Kato-Tanimura [24]. Notably, the heat current vanishes at the strong coupling limit. The lower panel shows the decoherence in the pointer basis for λs = 1.55. The dotted lines are the equilibrium density matrix at the effective temperature T = (T1 + T2)/2 = 1.5. The deviation from the equilibrium density is seen only around λb = 1 where the heat current reaches its maximum.
The present results show that indeed the decoherence due to environments is responsible for the suppression of heat. The off-diagonal elements of the steady state density look almost identical to those in the stationary state at a single effective temperature T = (T 1 + T 2 )/2 However, there is small but significant difference where the heat current is strong. The elements ρ 13 and ρ 24 deviate from ρ 12 and ρ 34 due to the difference in decoherence power between the two environments. In general a higher temperature environment causes stronger decoherence. [28] However, it also depends on the coupling strength as well. When the coupling strength overcomes the asymmetry in temperature, the decoherence power of the two environments becomes nearly identical and eventually the asymmetry in the off-diagonal element responsible for the heat conduction vanishes.
In conclusion, we claim that the "thermal equilibrium" of a small quantum system is not the Gibbs state when the coupling to the environments is strong. Due to continuous measurement by the environment, the stationary state loses the coherency between the pointer states and thus the density is diagonal in the pointer basis rather than in the energy eigenbasis. We further claim that the the stationary state density at the strong coupling limit is the Gibbs state projected onto the pointer basis. The diagonal elements in the pointer basis appear to be insensitive to the coupling strength. We have demonstrated this proposition by exact numerical calculation using the hierarchical equations of motion. This strong coupling limit can be used as a bench mark test for analytic models such as the Hamiltonian of mean force.