Topological Heat Transport and Symmetry-Protected Boson Currents

The study of non-equilibrium properties in topological systems is of practical and fundamental importance. Here, we analyze the stationary properties of a two-dimensional bosonic Hofstadter lattice coupled to two thermal baths in the quantum open-system formalism. Novel phenomena appear like chiral edge heat currents that are the out-of-equilibrium counterparts of the zero-temperature edge currents. They support a new concept of dissipative symmetry-protection, where a set of discrete symmetries protects topological heat currents, differing from the symmetry-protection devised in closed systems and zero-temperature. Remarkably, one of these currents flows opposite to the decreasing external temperature gradient. As the starting point, we consider the case of a single external reservoir already showing prominent results like thermal erasure effects and topological thermal currents. Our results are experimentally accessible with platforms like photonics systems and optical lattices.

with V = − J x,y a † x+1,y a x,y e iθ X x,y + a † x,y+1 a x,y e iθ Y x,y + h.c.
Here a x,y stands for the bosonic operator on the site (x, y) of the N × N lattice and where A denotes a gauge field. The Hamiltonian Eq. (S1) can be written in matrix form as where Ψ = (a 1,1 , . . . , a N,1 , a 1,2 , . . . , a N,2 , . . . , a 1,N , . . . , a N,N ) t and with N × N matrices so that H S is a N 2 × N 2 matrix. By diagonalizing the matrix H S = UDU † we obtain the spectrum of H S which can be visualized as the celebrated Hofstadter butterfly [S1]. The eigenmodes of this Hamiltonian are given by b = U † Ψ, with b = (b 1 , b 2 , . . . , b N 2 ) t . In Fig. S1 we have represented the density of states for an 8 × 8 lattice, the low density regions correspond to energies associated with edge eigenmodes.
In terms of these eigenmodes the system-reservoir Hamiltonian yields where we have taken into account that the mode a x,y is allocated as the component N (y − 1) + x of the vector Ψ and u i,j denotes the components of the unitary matrix U. After some rearrangement, we write for operators of left and right reservoirs, respectively. Now, we obtain the Davies generator of the weak coupling limit [S2] by applying the standard procedure (see, for instance, [S3]). Note that there are not mixed terms between different site reservoirs in the correlation functions, so the second order dissipator becomes whereL j,k (τ ) = y u N (y−1)+1,k [A j,y exp(−iν j,y τ ) + A † j,y exp(iν j,y τ )], and (similarly for)R j,k (τ ), are the operators L j,k and R j,k in the interaction picture (ν j,y is the frequency of the bath mode A j,y ), ω k is the frequency associated to the eigenmode b k , and β h,c = 1/(k B T h,c ) are inverse temperatures of hot and cold baths. Following the usual steps for the derivation [S3], we obtain the following master equation: wheren k (T ) = {exp[ ω k /(k B T )] − 1} −1 denotes the mean number of bosons with frequency ω k and temperature T , and γ k is a constant that depends on the strength of the coupling via the spectral density f (ω) ∼ g 2 j δ(ω j − ω). For the sake of simplicity we have assumed in the main text the same decay rate for each eigenmode γ k = γ, although this is not relevant to our conclusions. Furthermore the constants s k and r k are related to the matrix U via: Since the columns of the matrix U are the coordinates in real space of "one-particle" wavefunctions ψ k (x, y), after reordering, we obtain In the thermal equilibrium situation T h = T c = T , the master equation (S11) becomes withγ k = γ(s k + r k ). This equation describes the dynamics of the system towards thermal equilibrium with the baths, so that the steady state obtained for long times is the Gibbs state at the same temperature as the baths: In the general nonequilibrium case T h > T c , it is usually involved to obtain the steady state at the stationary limit. However, in this case we can find it by noting that the master equation (S11) can be rewritten as with some "effective" temperature T eff k depending on the mode: Therefore, the master equation (S17) is the same as the one describing the dynamics of a collection of N 2 independent modes with different frequencies interacting with N 2 thermal baths with different temperatures T eff k . Hence, we conclude that the steady state is of the form:

CURRENT OPERATORS
External and internal currents are derived from continuity equations. For the external case, using the master equation (S11) we have where J h (J c ) is the current operator that describes the heat flux between system and hot (cold) bath, and L denotes the Liouvillian at Eq. (S11) in the Heisenberg picture. Note that in the usual convention a positive current (meaning an increment of H S ) is associated to energy flowing from outside to the system, whereas a negative current describes energy flowing from system to outside. On the other hand, since and we identify the currents as This identification is standard in the theory of open quantum systems, and it can be proven [S4] that the time-evolution described by the master equation (S11) fulfills the entropy production inequality: where S = −k B Tr(ρ log ρ) is the thermodynamical entropy. In the stationary regime, t → ∞, the system approaches some steady state lim t→∞ ρ(t) = ρ ss , so that d H S dt = 0 and dSss dt = 0, and Eqs. (S20) and (S26) yield Since T h > T c , the above relations impose that J h ss = − J c ss > 0. This is in agreement with the second law of thermodynamics, in particular in the Clausius formulation as if no work is performed on system, the current goes from hot to cold bath [S5]. Note that if all baths are at the same temperature, T h = T c = T , once the stationary limit has been reached, the net external current between system and baths becomes obliviously zero. Coming back to the system of our interest, we can split the current operator corresponding to hot J h and cold J c baths as a sum of currents operators for individual baths. Specifically, using Eq. (S14), is the current operator accounting for the energy flow between the hot bath at position y and the system. Similarly, . To obtain J y h ss and J y c ss , we solve the dynamical equation for b † k b k : Therefore, the external currents at the stationary state are In order to derive internal currents we make use of the exact form of the continuity equation and the Davies' theorem [S2]. Specifically, the exact equation for the population at the site (x, y) is given by Then, the (internal) current operators are identified as: where the subindex (→ x) is a short notation for (x − 1 → x). So that J (→x),y denotes the operator for the current leaving the site (x − 1, y) and entering in (x, y). Similarly for (→ y). The term i [H SB , a † x,y a x,y ] in (S35) is not zero only for x = 1 and x = N , and defines the exact external currents. Of course we cannot compute the exact time derivative d a † x,y a x,y /dt; our approximation to it is given by the master equation (S11). However, if one tries to define internal currents directly from Eq. (S11), one finds a problem because the approximation introduces fictitious dissipative couplings among all oscillators. Nevertheless, the Davies theorem [S2] asserts that the dissipative part of (S11) is actually a weak-coupling approximation of the term Tr B (−i[H SB , ρ]) [S6]. This suggests that, in a weak coupling regime, it is consistent to take the above exact internal currents operators as internal current operators also in the master equation approximation. In this manner, the terms i [H SB , a † 1,y a 1,y ] are intended to be described by J y h as defined above. In a Landau-type gauge taken throughout the main document, we write A = (−|B|y, 0, 0) and internal currents take the form J (→x),y := iJ a † x,y a x−1,y e −2παiy − a † x−1,y a x,y e 2παiy , where α stands for the flux of B per plaquette. These equations, in terms of normal modes, become At the stationary limit the steady state is given by Eq. (S19) and we obtain Note that these internal currents describe the flux of carriers or quanta per time, we do not aim to associate any specific energy with the current from some site to its adjacent independently of what normal mode is excited in the lattice.
In addition, note also that if s k = r k , T eff k is invariant under the exchange T h ↔ T c , (because s k = r k ). This, in particular, implies that the internal currents do not change, neither their absolute value nor their sign, under the exchange T h ↔ T c . This may be surprising at first sight, but it is due to the fact that a simple exchange of T h ↔ T c can be seen as lattice rotation R π , because in that case left and right temperature are exchanged and the magnetic field and lattice properties remain invariant. Under such a rotation, it seems natural that currents do not change their value as they are chiral. However, under a reflection along the temperature gradient direction x ↔ −x, the temperatures are exchanged T h ↔ T c but also the magnetic field changes B ↔ −B. In this situation the currents indeed change their sign. This is showed in Fig. S2.
Finally, in Fig. S3 we depict the edge/bulk current ratio as a function of the temperature for the numerical parameters taken throughout the manuscript. We can distinguish the three phases in Fig. 2 of the main text. For low and high temperatures the bulk current is similar to (or higher than) the edge current. However, for a large intermediate range of temperatures the system remains in a phase with high edge current concentration.

MASTER EQUATION IN THE "LOCAL" APPROACH
In the analysis of internal currents it is also common to consider the sometimes refereed to as the "local" approach to the master equation. In this approach the master equation for our system reads This master equation corresponds to the introduction of dissipation in a kind of "adiabatic" way; by assuming that the dissipative process is not affected when changing J from 0 to some small parameter (see discussion in [S3, S7]). So, this master equation is expected to provide a good description of the dynamics as long as J is small in comparison with the typical time scale of the system when J = 0. In the stationary regime at the limit t → ∞, this equation is not expected to be a good description and it presents some problems from a thermodynamical point of view. Namely, 1. For T h = T c = T , the state at the long time limit provided by this master equation is not a Gibbs state. Of course, if J is small enough both states are very close [S7], however this does not guaranty a correct thermodynamic description [S8].
2. For T h > T c , external currents may violate the second law of thermodynamics [S9], accounting for an unphysical heat flow from cold to hot bath.
Yet, for the sake of comparison we have solved this equation (S44): no chiral currents are found, and the pattern is independent of the values of α. This is in agreement with the results reported on [S10] for a lattice with two rows (a ladder) where no chiral current was obtained.