Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Non-equilibrium evolution of quantum systems connected to multiple baths Wu, Jinshan 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2006-0719.pdf [ 4.39MB ]
JSON: 831-1.0084847.json
JSON-LD: 831-1.0084847-ld.json
RDF/XML (Pretty): 831-1.0084847-rdf.xml
RDF/JSON: 831-1.0084847-rdf.json
Turtle: 831-1.0084847-turtle.txt
N-Triples: 831-1.0084847-rdf-ntriples.txt
Original Record: 831-1.0084847-source.json
Full Text

Full Text

Non-equilibrium evolution of quantum systems connected to multiple baths by Jinshan Wu B .Sc , The Beijing Normal University, 1999 M . S c , The Beijing Normal University, 2002 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Master of Science in The Faculty of Graduate Studies (Physics) The University Of British Columbia October 2006 © Jinshan Wu, 2006 i i A b s t r a c t In this thesis we develop a systematic formalism that allows the calculation of transport properties for small quantum systems coupled to multiple baths with different temperatures and/or chemical potentials. Our approach is based on a generalization of the projector operator technique, previously used to study the evolution towards equilibrium of systems weakly coupled to a single bath. Applying this technique to a system weakly coupled to multiple baths, we find a dynamical equation for the reduced density matrix of the system (the re-duced density matrix is obtained by tracing out from the total density matrix the baths' degrees of freedom). Integration of this dynamical equation gives the time evolution of the biased system. After a time interval long enough compared to the characteristic time-scales of the problem, the system arrives to a non-equilibrium steady-state. The density matrix in this non-equilibrium steady state can also be calculated directly, avoiding the time-integration of the dynamical equation. In the second part of the thesis we study non-equilibrium steady states and heat transport through short spin chains coupled to two thermal baths with different temperature. Two new definitions for the local temperature in non-equilibrium are proposed, based on local spin and energy at different sites. The definition based on the local energy seems to be more reliable. We find that for small biases, the heat and spin currents increase proportionally to the bias, however non-linear behavior is observed for large biases. The temperature profile of the chain depends strongly on the parameters of the system. i i i C o n t e n t s Abstract i i Contents i i i List of Figures iv Acknowledgements v 1 Introduction 1 2 Dynamics of open systems 4 2.1 General formalism: system coupled to a single bath 5 2.1.1 Projector operator technique and the reduced Liouville equation 6 2.1.2 Alternative derivation: the non-correlation approximation 9 2.1.3 Lindblad form: general equation for the reduced density matrix 10 2.2 Applications: thermal equilibrium 11 2.2.1 One-state open systems 11 2.2.2 Multi-state open systems 15 2.3 The local coupling scheme is wrong 17 3 Spin chains coupled to thermal baths 20 3.1 Single spin coupled to a heat bath 20 3.2 Ising spin chain coupled to a heat bath 22 3.3 X Y chain coupled with one heat bath: thermal equilibrium . . . 25 3.3.1 Sample calculation for a 2-site X Y chain 25 3.4 X Y chain coupled to two baths 28 3.4.1 Sample calculation for a 2-site X Y chain 29 3.5 Definition of the local temperature and heat current 32 3.6 Current and local temperature profiles in non-equilibrium steady-states 38 4 Conclusion 50 Bibliography 52 iv L i s t o f F i g u r e s 1.1 Temperature profile of a classical harmonic chain 2 2.1 A two-site chain coupled to a heat bath 17 3.1 3-spin X Y chain coupled with two heat baths 28 3.2 Time evolution of off-diagonal elements of the density matrix for a biased system 32 3.3 (<7JZ)T as function of the inverse temperature, for (a) two-site chain, and (b) three-site chain 34 3.4 Expectation value of the "local Hamiltonian", (hi)T as a function of inverse temperature, for (a) two-site chain and (b) three-site chain 36 3.5 Temperature profile for a three-site chain in the stationary state, when a bond is being removed 38 3.6 Temperature profile for a two-site chain in the stationary state, when the coupling to the right bath is being removed 39 3.7 Spin current in the non-equilibrium steady-state of the two-site chain, as a function of the temperature of the left bath 40 3.8 Heat and spin currents in the non-equilibrium steady-state of the three-site and four-site chains, as a function of the temperature of the left bath 41 3.9 Steady-state temperature profile for a 2-site spin chain as a func-tion of T i 43 3.10 Steady-state temperature profile for 3 and 4-site spin chains as a function of TL , 44 3.11 Steady-state temperature profile for 2, 3 and 4-site chains when TL = 4,TR = 1 45 3.12 Steady-state temperature profile for 3 and 4-site chains as a func-tion of TL, for JXY = 4 and bz = 2 47 3.13 Steady-state temperature profile for 3 and 4-site chains as a func-tion of TL, for JXY = 2 and bz = 0.1 48 V A c k n o w l e d g e m e n t s With a deep sense of gratitude, I wish to express my sincere thanks to my supervisor Mona Berciu, for her guidance to me into this new field, for her support without any reservation even when I diverged from the topic for a long time and by a long distance, for her constantly grabbing back onto the right direction when sometimes (or quite.often) my too wide interests are consuming me, and for her great help and stimulating discussions on some other projects I did just for fun. It must have been a hard task for her to find such a good balance to make me stick on one problem and let me wonder free. I also want to thank all my teachers, classmates for always being so nice to me, even when sometimes I was so annoying and stubborn. I asked all kind of good and bad questions in classroom and office-rooms, but they have never refused to discuss with me and always provided comments when I asked for them. Instead of listing all of them, I especially want to mention James Charbonneau, Mandy Wong and Kory Stevens, Ian Affleck, Ariel Zhitnitsky, Mark Van Raamsdonk, Moshe Rozali, Leslie Ballentine, Andrew DeBenedictis, Zhanru Yang, Shouyong Pei and Zengru Di . I am very grateful to my wife Qi an Feng, for her love and patience during the extraordinary long Master period, which includes this two-year Master program in U B C , one year of P H D program in S F U and a three-year Master program back in B N U . Without her belief and encouragement I do not think I could have made it. Finally I have to say that I owe everything to my mother Yundi Huang, who did not get the chance to share any of the happiness of this graduation, or even of my bachelor degree and first master degree in B N U . It is her belief in me that drives me forward when sometimes I feel lost. It is her fight with her illness that gives me the power to conquer my own difficulties. It is her hospitality, generosity and caring for other people that made me a popular person so that I did not feel lonely. M y mom has been a great treasure for me, forever. 1 C h a p t e r 1 I n t r o d u c t i o n The long term goal of this work is to develop a formalism suitable for inves-tigating transport properties (of charge, spin, heat, etc.) and more generally, for studying the evolution of interacting quantum systems far from equilibrium. We first focus on systems with only a few degrees of freedom that are coupled to one or more baths (environments), although in the longer term the aim is to study large systems (bulk materials). The study of such "small" systems is con-venient for technical reasons, as development of numerical code is simpler and run times are conveniently short, so that the emphasis is on finding the correct framework. Moreover, transport in quantum systems with just a few degrees of freedom, such as, for example, a few coupled quantum dots, is currently of large experimental and theoretical interest. There are several ways to approach the computation of transport proper-ties, including the Kubo formalism [1], the Landauer-Biittiker formalism [2, 3] and the more general non-equilibrium Keldysh formalism [4, 5]. Compar-isons between these approaches have also been pursued, see for example Ref. [6]. However, the Kubo formula is. not suitable for systems far from equilib-rium/stationarity, while the Landauer-Biittiker approach is not suitable for in-teracting systems. While in principle exact, the Keldysh formalism involves (the usually impossible task of) summation of infinite numbers of Feynman diagrams, thus approximations are in fact required to perform any calculations. In the study of open systems, there is another well-known approach, the so-called projector operator technique. For a system S coupled to an environ-ment/bath B, the general idea is to start from the full equation of motion for the total system T = S + B, and by doing a projection, to derive a reduced equation of motion for the evolution of S alone. Our idea is to make use of such projection techniques for the''situation where the system of interest is coupled to two (or more) environment/baths, which is the typical case for any transport measurement. Once we obtain the reduced equation of motion of the central system, we can, for instance, find its stationary solution and therefore calcu-late transport properties in the stationary limit, and their dependence on the various system, baths, and/or system-baths coupling parameters. We can also investigate short-time evolution towards this stationary state and other non-equilibrium or time-dependent properties. Of course, the results we obtain with this approach should and will be gauged against those of other methods, where appropriate, and against experimental results. In the work described in this Thesis, we focus on heat transport as a first test case for our formalism. One reason for this is that thermal baths can be Chapter 1. Introduction 2 T(i) Figure 1.1: A typical temperature profile for a classical harmonic chain coupled to left and right baths at TL and T+ respectively. We thank the authors of Ref. [7] for their permission to use this figure. simply modeled as collections of many bosonic degrees of freedom equilibrated at the desired temperature. By contrast, to study electrical transport, one needs baths/reservoirs of electrons, and both the temperatures and chemical potentials become relevant parameters. The heat-transport problem is also interesting on its own, since the question of under what kind of circumstances a system has normal thermal conductivity, i.e. it follows Fourier's Law, is still open. There are lot of works on thermal conductivity (for a review, see Ref. [7]), with most studies investigating classical systems, although there are also several studies of quantum spin systems [8, 9]. It is well known that a purely harmonic classical chain does not exhibit a normal thermal conductivity. Instead, it has an almost flat temperature profile so that its conductivity diverges. Fig. 1.1, taken from Ref. [7], shows this typical behavior. It is believed that either higher order (inharmonic) interactions or random scatterers are needed to lead to a normal conductivity. Recently, several groups have used the concept of self-consistent baths, where each site is coupled to its own thermal bath, with a temperature adjusted self-consistently so that no energy flows in or out of the local bath. It was shown that in this case, normal conductivity can be recovered even for pure harmonic systems [10]. This is not so surprising, since these local baths effectively model incoherent scattering at each site. On the other hand, behavior in their presence is obviously very different than the behavior of the system alone, so questions remain as to how reasonable this approach is. Our general approach to this problem is to follow these steps: First, we choose reasonable models for the heat baths and for their (weak) coupling to the system of interest. For simplicity, the systems of interest in this Thesis will be short chains of spins \ , coupled to one another through anisotropic Heisenberg exchanges, and possibly also placed in an external magnetic field. Secondly, we use the projection technique to obtain the reduced equation of motion of the central system. Thirdly, we solve this time-dependent reduced equation of motion to find the reduced density matrix at all times, in particular its long-term (stationary) limit. Finally, we use the definition of the thermal current operator and of the local temperature to compute the current and temperatures Chapter 1. Introduction 3 profiles. This allows us tp compute global and local thermal conductivities and also compare our results to those from other techniques. Generalizations to other types of transport (primarily, of electrical charge) should follow the same framework, with suitable modifications. The thesis is organized as follows. Chapter 2 contains a brief review of open-system dynamics as a preparation for our method. Some relevant and hopefully interesting comments are included. In Chapter 3, we discuss first the thermal equilibrium of systems coupled to a single bath, or to multiple baths at the same temperature. We then generalize to systems coupled to two or more baths of different temperatures, and discuss their heat conductance. At last, in Chapter 4 we summarize our results and thoughts on future application of this approach to calculating electrical conductivity. C h a p t e r 2 D y n a m i c s o f o p e n s y s t e m s Consider a system S interacting with a single bath/environment B. In the Schrodinger picture, the Liouville equation for the evolution of the density ma-trix of the total closed quantum system S + B is: H,p(t) . (2.1) dt The total Hamiltonian has the generic form: H = H0(S,B) + V{S,B), (2.2) where H0 (S, B) = Hs [S) ®I{B)+ I (S) ® HB (B) (2.3) describes the system and the bath, and V(S,B) = J2Li(S)®Fi(B) (2.4) i describes their interactions, with Li (S) and Fi (B) being system, respectively bath operators. When we are interested only in the properties of S, we call such a system an open quantum system. Typically, we need to calculate the ensemble average of physical quantities depending on operators acting on S only, (A)=tr([As®i(B)]p). (2.5) This can be trivially rewritten as: (A) = trs (AsRs) , (2.6) where Rs (t) = trB (p(t)) (2-7) is the reduced density matrix of the system 5. Therefore, it is convenient to try to find Rs {t) directly, and avoid solving for the total density matrix p (t), which is a much more complicated problem due to the significant increase in the number of degrees of freedom. Note that Eq. (2.1) can be written in a unified form suitable for both classical Chapter 2. Dynamics of open systems 5 and quantum systems by introducing the Liouville operator dp{t) dt Cp(t), (2.8) where for quantum systems, C , while for classical systems it is a Poisson bracket, C = - j . We would like to write the evolution of the reduced matrix density in a form similar to Eq. (2.8): M = ^ S ( t ) , . (2.9) although now C is obviously no longer the same as for the whole system. The task of this chapter is to derive such a general Liouville operator - in fact called a Lindblad operator - starting from knowledge of the total Hamiltonian H. In Section 2.1 we briefly review the two techniques most used to deal with this problem when the system is coupled to a single bath, and in Section 2.2 we apply this to a few simple models, and show that the correct equilibrium is reached in the steady-state limit. Finally, in Section 2.3 we comment on the local coupling approach, which has been used by several authors to generalize this scheme to systems coupled to multiple baths. 2.1 General formalism: system coupled to a single bath There are two main types of approaches to discuss the dynamics of an open system: dynamic variable based or density matrix based. The first approach uses equations of motion for dynamic variables such as X, P for a classical system and X, P for quantum system. The dynamic variables of the bath are integrated out, resulting in noise sources in the Langevin-type equations for the dynamic variables of the system. If one can solve these Langevin equations, the full time dependence of the system is known. From a Quantum Mechanics point of view, this approach is similar to Heisenberg's picture. The second approach focuses on distribution functions in the phase space of the system, such as p (x,p) for classical systems and p for quantum systems. Equations of motion are generated for them, the bath is projected out, and then the equation for the reduced density matrix is solved. From a Quantum Mechanics point of view, this approach is more similar to Schrodinger's picture. Since our approach is density matrix based, in this section we summarize this general formalism. There are two well-known techniques to derive the equation of motion of the reduced density matrix, namely the projector operator technique and the non-correlation approximation. We discuss the first one in reasonable detail and then shortly review the second one. Chapter 2. Dynamics of open systems 6 2.1.1 Projector operator technique and the reduced Liouville equation One way to derive the unknown C in Eq. (2.9) is through the technique of projection. We define the projector operator on the Hilbert space of the total system S + B, V : TL —> TC, such that (2.10) where p^ is the density matrix of the isolated bath in equilibrium, although this requirement could be relaxed. It is easy to verify that V is a projector, i.e. that TT = V. We denote Q = I — V, and rewrite the equation of motion for the total system, E q . (2.8), as: §-tVp(t) = VCVp(t)+VCQp(t) df mQp(t) = QCVp(t) + QCQp(t) (2.11) The second equation is integrated formally and its solution is substituted into the first equation, to obtain: d_ di Vp(t) = VCPp{t) + VC [ dTe{t-T)QtQCVp{T)+VCetQtQp{Q). (2.12) Jo If we start from the reasonable initial condition p{0) = Rs (0)®,3f, then N Qp(0) = 0 (2.13) (2.14) and the last term vanishes identically. Let £ 5 , LQ, L\ be the Liouville operators corresponding respectively to Hs (S), HQ (B) and V (S, B). Then, C = Cs + CB + C-i — Co + C\. Since CBP"B = 0, we can rewrite the dynamic equation as (2.15) ^Rs {t)®pe§=VCs [Rs {t)®peg + VCi +VC / dre^-^QC Rs (t) ® p£ Rs (t) ® Pei (2.16) The second term is identically zero if we require that for any eigenstate in the bath's Hilbert space, the expectation value of the bath operators involved in the bath-system coupling [see Eq. (2.4)] vanishes: (6|F i(B)|6)=0. (2.17) Chapter 2. Dynamics of open systems 7 In this case, we have: tfi(6|£i \RS (t)®PB« b) = (b\ \li (S) ® Ft (B), Rs (t) ® /5 e B ' l |6) £ (% (5) A 5 (*) (b\ Pi (B) pe* \b) - Rs (t) U (5) (b\ ffiFi (B) |6>) = 0 showing that in this case, Vt{P = 0. This condition, which is satisfied for most types of couplings, is not however necessary. If it is not obeyed, there wil l be term of the form J2i Pi ^ Li, which can be added to Cs- It follows that this assumption is not essential, although our coupling Hamiltonians wil l satisfy it, in the following. With these assumptions and after tracing out the bath degrees of freedom, we find that: di Rs(t)=£sRs(t) + trB |T?£ j f dre^-^QCVpit)^ . (2.18) Subject to the two (reasonable) assumptions mentioned above, this equation is exact. Clearly, all the effects.of the bath on the system's dynamics are contained in the second term on the right-hand side. This term can be further simplified. One can verify that Q£V = QC0V + QC{P = QdV = -VCXV + C{P = L{P (2.19) where the second assumption has been invoked again. Now we use the first essential assumption, replacing: O(t-T)QC0 (2.20) This is valid when the system-bath interaction is weak, since, as shown in the following, it is equivalent to second order perturbation in the "coupling strength" between system and bath. This follows because we can then simplify: e(t-r)QC0 = e(t-T)(X-P)£0 o(t-r)(C0-VC0) = (t-T)£0 (2.21) where we make use of VC0 = VCQV + VCQQ = 0. Therefore we obtain: ^RS (t) = CSRS (t) + trB JA j f dre^-^C.Rs (r) ® ^  j , (2.22) where indeed the second term is of second order in the interaction V. This can be further simplified in the interaction picture, where we define Rrs (t) = eifl^Rs{t)e-if1^, L]{t) = e ^ ^ e " ^ ' , Pj (t) = e^^P^"^. In this Chapter 2. Dynamics of open systems 8 case, the first term cancels out, and the dynamic equation becomes: dR's (t) dt The two commutators can be re-organized dR's (*) ^ L f ( r ) F / ( - i + r ) , ^ (r) p% . I (2.23) dt where / d r { [ L j ( * ) , £ f ( t - T ) £ £ ( * -ji J o - [ L j ( 0 , ^ ( t - r ) L f ( t - r ) Gjt {r)=trB (FJ{T)Fl(0)pi) Glj(-r) = trB (F/(0)FJ(T)p%«) GI G0x (r) ( - T ) } (2.24) (2.25) are quantities (Green's functions) which depend only on the bath characteristics and can be calculated directly. Now comes our second essential assumption, the Markov approxi-mation: PJS it - r) « PJS (t) for t > r c , (2.26) where r c is the characteristic relaxation time of Gji (r) (i.e., of the bath). Since in this thesis we investigate predominantly the steady-state solution when t —» oo, the Markov approximation is likely justified. It can certainly be relaxed if the equation is integrated numerically, and we plan to investigate its consequences in future work. Wi th this final approximation, we find the desired equation of motion for the reduced density matrix: dR's (t) dt -h ^ J o d T { lLIj { t ) ' l ' { t ~ T ) k ' s {t)]Gjl ( r ) L J ( t ) , ^ ( t ) L f ( t - r ) ] G y ( - r ) } . (2.27) After Rg(t) is solved for, we transform back to the Schrodinger picture Rs{t) = e~~lHstRIs{t)elHst, and then calculate all physical quantities. Alternatively, we can transform this dynamic equation back into an equation for Rs(t), and di-rectly solve that. Note that in some cases, the evolution of the system is such that the reduced density matrix is always diagonal in the eigenbasis of H, Rg, Hs = 0, so that Rs(t) = Rs(t)- For a system coupled to a single bath, the (equilibrium) steady-state limit corresponds to a diagonal reduced density matrix, so certainly we have Rs{t —> oo) = Rs{t —> oo). However, this is generally not true for systems Chapter 2. Dynamics of open systems 9 coupled to mul t ip le baths. 2.1.2 Alternative derivation: the non-correlation approximation Th i s is a different way to derive the above equat ion of mot ion for the reduced density ma t r i x , using the so-called non-correlat ion approx imat ion . Suppose that bo th Hs (S) and HB (B) are solvable, and their eigenstates are respectively {\s)} wi th Ef and {\b)} w i t h E B . Then , a complete basis for the to ta l H i l be r t space can be chosen as {\s) \b)}. In the interact ion picture \ipT (r.)) = elHo{s,B)t |^ ^ Q ^[me e v o l u t i o n of the operators descr ibing the coupl ing is L\ (t) = e l H s ^ t L i (t) e~lHs(S)t a n c j p? (t) = e i H B ^ t F i it) e-iHB{B)t^ T h e d y n a m i c s equat ion of the whole system is: -dp1 (t) ih-dt V1 (t)rP' (t) (2.28) Integrat ing this formal ly and subst i tut ing the result ing expression for p1 i n the right-hand side of E q . (2.28), we obta in : dp1 (t) ^ W , P 7 ( 0 ) ^jyr[v'(t),[v'(T),p'(r)]]. (2.29) dt h Integrat ing out the ba th , this becomes: dJ§r - r B ( [ * " < ' ) ^ ' < » ) ] ) - # I ' d * f l ( [* '»> . [ O ' M . A ' M ] ] ) . (2.30) which is s t i l l exact. Now we begin to introduce several approx imat ions . . F i r s t , we again assume (b\ ^  (B) \b) = 0, (2.31) so that the first t e rm i n E q . (2.30) vanishes. A s before, th is is not an essential assumpt ion. Second, we also assume, p1 (r) = RIs(r)®pIB(r): (2.32) such that the second term i n E q . (2.30) can be separated and s impl i f ied. Th i s is the essential approx imat ion , basical ly equivalent to the second-order per tur -bat ion theory i n the projector operator technique. T h i r d , we take, pIB{r)=peBq = f-Be-pflB. (2.33) W i t h these, we get a s impl i f ied dynamic equat ion for the reduced density ma t r i x ident ica l to E q . (2.24). If we further apply the Markov condition, we arrive at the same equat ion of mot ion as in E q . (2.27). Th i s equat ion can be projected i n the basis. In this case, the reduced Chapter 2. Dynamics of open systems 10 equation of motion becomes /^Rs'smn (t) Rlnn (t) (2.34) -L. r+ _|_ r -kkm 1 nss'm 1 nss'm (2.35) (2.36) where R, = e * i £ - ' " » " * - V c 5 r+ lljs smn ° / y "sn1 s ' where and r+ Hn H ^ <ZI ^ l"> f e " * ^ - ^ ) " ^ (r) dr (2.37) a n d rmfcin (0 = (Tnikm (*))*• I f w e f o c u s o n l Y o n t h e diagonal terms for RTSS (t), we obtain the Quantum Master Equation, which requires a further approxima-tion in order to be consistent with Fermi's Golden Rule. The so-called secular a p p r o x i m a t i o n [11] is to require Ef, - E: E" 0,. (2.38) as to eliminate the explicit time dependence of the right hand side of Eq. (2.35). Then finally -Ri. (t) = ] T WsmRlnm (t) - WmsRL (t), (2.39) where Wsm = T + s s m + T ~ M m . This is now consistent with Fermi's Golden Rule. While there a considerable body of work based on such Quantum Master Equations, we choose to work directly with the dynamic equation (2.27). 2.1.3 Lindblad form: general equation for the reduced density matrix Evolution operators of a closed quantum system form a group, but generally evolution operators of an open system form a semigroup so as to keep the density matrix always positive and normalized to 1. Lindblad discussed the generators of such trace-preserving positive semigroups [12]. As shown in a theorem of Gorini, Kossakowski and Sudarshan [13], a completely positive semigroup evolution of a quantum object in a n-dimensional Hilbert space, p(t)=U(t)p(0) (2.40) Chapter 2. Dynamics of open systems 11 such that W ( t i + t 2 ) = t f ( t i ) t f ( t 2 ) (2.41) can be characterized by a Lindblad operator such that its evolution follows -p(t)=C(t)p(t), where generally, - Tl — 1 £(t)p(t) = - \ [H,P] + - £ A n ([LupL + (2.42) (2.43) The operators S = | l / j | u | - ^ / | must form a complete basis of the operator Tr (U]V) = SuvyU, V e B . (2.44) space, Eq. (2.43) is similar to Eq. (2.27), except that its coefficients (t) may generally depend on t. Thus, the dynamic equation for the reduced density matrix that we derived is one example of the general Lindblad form. This general Lindblad form allows for non-Hermitian evolution. For example, some of the terms in Eq. (2.43) can be re-organized into a form, however, the resulting O is no longer necessarily Hermitian. 2.2 Applications: thermal equilibrium We first apply this general procedure to the simplest examples of a fermionic or bosonic system coupled to a bath, and show that the system evolves to the correct thermal equilibrium as t —» oo. We then briefly investigate a more general case and discuss the so-called "local coupling" scheme. Besides being a useful exercise, this will also demonstrate the subtle point that this local coupling, which is used in literature [8, 9], is not appropriate; 2.2.1 One-state open systems Consider a one-state system, defined by a Hamiltonian: HS = wa+a (2.45) where w > 0 and a, could be either bosonic or fermionic operators. For fermionic operators, this is a two-level type system that can accommodate either zero or 1 particle. For bosonic operators, the system is a harmonic oscillator. The system is coupled to a bath. For the harmonic oscillator, the bath is regarded as an ensemble of phonons (with p = 0) thermalized at a given temperature T, so that only heat is exchanged. For the two-level fermionic Chapter 2. Dynamics of open systems 12 system, the bath is regarded as an an ensemble of fermions with a fixed T and ^ 0, so that both particles and energy are exchanged. We take H = H0 + V, (2.46) H0 = HS®IB+Is® HB, (2.47) HB = YJt(k)blbk. (2.48) k The bk operators are bosonic or fermionic, depending on whether the bath is made of phonons, or of thermalized fermions. For phonons, we always have e(k) > 0. The system-bath coupling is taken to have the general form: V = J2 (A^ a ® bt + A V f c * a t ® 6fc) • ( 2 - 4 9 ) k which describes an exchange of particles between the system and the bath, for the two-level fermionic system, or the excitation of the system upon absorption of phonons from the bath, for the harmonic oscillator. Note that this coupling is not the most general possible, since it excludes terms of the form a ® bk and a* ®bk. For the fermionic system, instead of the particle hopping process described in the former two terms, these later terms correspond to the process where both the central system and bath gain/lose one particle. For the harmonic oscillator, the former two terms describe energy exchange processes between the central system and the bath, however, the later two terms describe processes where both the central system and the bath gain/lose energy. For simplicity, we ignore such processes here; in section 3.1 we demonstrate that indeed they do not contribute to the dynamical equation we obtained within the approximations we described. For concreteness, we assume that initially the system is in its ground-state, and the bath is in equilibrium: p(0) = Rs (0) tg. pe« = |0) <0| ® = , (2.50) where ZB = trB (e-^HB~^N^) = T[k C£nk e-^(fc)-")"*). We use Eq. (2.27) to describe the system's dynamics. First, we identify L[l) = \Vka,Fil) = &t ;4 2 ) = A V > t , i f > = bk. (2.51) The only nonzero bath Green's functions are found to be: Glkl (±r) = trB (F^ (±r) i f > (0) pB) = (nk) e ± ^ G*;jt (±r) = trB [F12) (±r) i f > (0) pB) = (1 ± nk) e^<k)\ ^ where in (l±nk), the — corresponds to a particle bath and the + to a phonon Chapter 2. Dynamics of open systems 13 bath. Then the reduced equation of motion for R!s (t) becomes: [XVkoe-^, WtaU^-^Rtg (t)] G1* (r) 23kW = i T r ' d T J +Wk^^\\Vkae-^-^RIs{t)}Gll{r) dt W^khaTS _ [xvkae-^, R's (t) XV^e^-^] G^ (-r) - [XV^aU^, R!s (t) XVkoe-W-r)] G1^ ( - T ) { <nfc> [a,a^Rrs (t)] e-^Teie^T + (l±nk) [a^aR's (t)] e^e'^ -(l±nk) [a,R's (<)at] e - * " - e + * ( * > --(nk) [a\RJs (t)a] e^e^^ As a last step, for long-time evolution we replace /0* dr by /0°° dr, and use the relation dre-l"T = TT5{LO) - iP (-) . (2.53) / Jo The second term gives an overall energy shift and can be absorbed into u> (or, we can assume that the system was coupled to the bath at t = —oo instead of t — 0). Then, in the interaction picture, we find the equation of evolution for the reduced density matrix of such an open system to be: ^R's (t) = n [ 2 a i ^ (t) at - 0 t o i i / ( t) _ ^ ( t) a t a ] +r 2 [2a +i?^ (t) a - aa f i ?^ (*) - (t) aa f ] , ^ A 2 £ IVfcl2 (1 ± n f c) <5(W - e(k)) = r [1 ± n(w)], (2.55) fc p A 2 |V f c | 2 (n fc) <J(w - e(k)) = rn(w), (2.56) where »"2 = inhere fc is a known quantity depending on the bath-system coupling, and n(u>) = —r-,—~r > £ / 3 ( W - M ) zf l is the average number of phonons/particles with energy w in the bath (for phonons, p. = 0). Now we can solve the dynamical equation and check the stationary solution. Let us first consider the case of the fermionic two-level system. First, note that if the initial reduced density matrix is diagonal (as we assumed it to be, Rs(0) = |0)(0|), then this equation never gives rise to off-diagonal terms. In Chapter 2. Dynamics of open systems 14 this case, Rs(t) = RJs(t), and therefore we can work directly in the Schrodinger picture. We set: Rs (t) = 1-±p$- |o> (0| + ) 1 } ( 1 , t ( 2 . 5 7 ) so that the trace is preserved. The equation for p (t) is then straightforward to extract from Eq. (2.54): jtP (t) = 2 ( n - r 2 ) - 2 ( n + r 2 ) p (t), (2.58) p(0) = 0 . This has the stationary solution *<°°> = r 7 T 3 ' ( 2 ' 5 9 ) which for the fermionic bath is: p(oo) = This leads indeed to the expected equilibrium grandcanonical distribution i e-P(.u-v) Rs(oo) = — s 10> <0| + ^ 7 . 11)(1| for the two-level system. If there are off-diagonal terms when the coupling is turned on, Eq. (2.54) can still be solved by more complex techniques, such as using the coherent state representation [14], with the same end result for the stationary state. In fact, if we assume R's (t) = 1~±^ |0) <0| + |1) (1| + q |0) <1| + q* |1) (0|, (2.60) besides Eq. (2.58), we wil l also get ±q(t) = -q(t), (2.61) which has only q = 0 as its stable stationary solution. As expected, in the stationary state (which in this case is equilibrium with the bath) the reduced density matrix is diagonal in the eigenstate basis. Note that we could also consider this two-level system as being a spin-^ in a magnetic field. In this case, we can couple it to a thermal bath of phonons. This case is discussed in the next chapter. The harmonic oscillator also gives the correct answer. Again, if we start in any eigenstate (or more generally, with a diagonal reduced density matrix Rs(0))t the reduced density matrix remains diagonal at all later times and Chapter 2. Dynamics of open systems 15 there is no distinction between Schrodinger and interaction pictures. We take: oo Rs(t) = J2P"(t)\n)(n\ n=0 with Yl™=oPn(t) ~ 1- Introducing this into Eq. (2.54) leads to: ^Pn(t) = 2 [(n + l)(ripn+i(t) - r2pn{t)) - n(ripn(t) - r2pn-i(t))} so the stationary state solution requires that rxpn = r2pn-i, (V)n. Using the expressions of r\ and r2 for the phonon thermal bath and the normalization condition, this gives e-nf3ui 1 Pn{00) = ^—,Z 1 - e-l3" as expected in thermal equilibrium. Note that in both examples, the value of r itself is irrelevant: if the system is coupled to a bath so that r ^  0, the same equilibrium state wil l be reached sooner or later. 2.2.2 Multi-state open systems For multi-state systems, we generalize the system Hamiltonian to: Hs = Y,uj(q)alaq, (2.62) where we take these to be fermionic operators, for simplicity (a simple specific example is given in the next section). Assume first that each of these system's eigenstates is coupled to its own, independent bath: ^ = A J2 [V"kai ® bU + h-c] (2-63) k,q where {bq,k,b]q, k,} = 5k,k'Sq,q'. In this case, all Green's functions involving operators from different baths vanish, and the extra terms in the dynamic equation of motion jtR's = Cdiag Rs (2.64) are simply a sum over all systems's eigenstates, consisting of contributions sim-ilar to those in Eq. (2.54) coming from each mode: Cdias = J2(ri(q)Cl+r2(q)C2q). (2.65) Chapter 2. Dynamics of open systems 16 /here ri (q) = Ek £A 2 \ V q k \ 2 (1 - nk) 5 (e (fc) - w (q)) (2.66) •^2 (q) = Ek ^ A 2 |V g f c | 2 (nk) 5 (e (fc) - w (q)) and we introduce the short-hand notation: £qRs = [aq, Rsa\] ~ [aj, a 9-Rs] = 2aqRsa\ - a\aqRs - Rsa\aq, (2.67) £2/?s = [aj, Rsaq] - [aq, aqRs] = 2aqRsaq - aqaqRs - Rsaqa]q. (2.68) In this case, the density matrix again remains diagonal, and following the same arguments as before, one finds that each mode equilibrates with its own bath. However, if the system is coupled to a single bath, V = A ]T Wqkaq ®bl + /i.e.] , (2.69) k,q the full dynamic equation also has terms involving different eigenmodes such as [aqi,al2Rs]e-i"-i^^\ (2.70) These come from processes where a particle vacates the qi mode and enters the bath, after which it leaves the bath and occupies the qi mode in the system. In this case, the dynamic equation is determined by the total operator C: C = C d i a g +Coff, (2.71) where C°"RS = [ e - ' A ( 9 l ' 9 2 ) t (n (q2) [aqi,Rsal2] - r 2 (q2) [aqi,aq2Rs]) + eiA^^ ( r 2 (q2) [a^Rsa,,] - n (q2) , aq2Rs])] , (2.72) where A (qi,q2) = w (qi) — UJ (q2). It is easy to see that the qi — q2 terms are included in Cdiag. The off-diagonal terms are explicitly time dependent as long as the energy difference w (qi) — to (q2) ^ 0. If we are interested only in the steady-state (equi-librium) limit, these off-diagonal terms can be ignored. The reason is that the terms multiplied by such explicit time-dependent expressions, like e* 1 ^ 9 , 1 ' 9 2 ) ' must vanish in the steady-state limit. This can be shown to simply be equiv-alent with requesting that the reduced matrix becomes diagonal in this limit. The diagonal values are then seen to be the expected Boltzmann distributions, like in the case where each mode is thermalized by its own bath. However, note that if the system is coupled to more baths which have dif-ferent parameters (different temperatures, chemical potentials, etc), the steady-Chapter 2. Dynamics of open systems 17 1 Bath Figure 2.1: A two-site chain coupled to a heat bath. state limit is not an equilibrium state, and there is no a propri reason why the reduced density should become diagonal in this limit. In fact, as we show in the next section, in the generic situation the density state must have off-diagonal terms in the steady-state limit. In such cases, these off-diagonal terms become very important and must be carefully included in the calculation. To conclude, the general procedure to study the evolution towards equilib-rium of a system coupled to a single bath is: 1. find the eigenstates of the central system, so that its hamiltonian is written in diagonal form like in Eq. (2.62). 2. rewrite the coupling between bath and central system in this new basis, so that it has the general form of Eq- (2.69). 3. find the dynamical equation of the reduced density matrix Eq. (2.64), keeping only the time independent operator Cdlag. The steady-state limit is then found to be the equilibrium state. If the time evolution is of interest, the off-diagonal terms should be included as well, however the resulting steady-state will remain the same. 2.3 The local coupling scheme is wrong Before ending this section, we need to address the so-called "local coupling scheme", which is quite extensively used in literature [8, 9, 11] to write the general form of the operator C directly. To illustrate how this scheme works, consider a two-site electron chain system coupled to a heat bath, see Fig. 2.1. The system's Hamiltonian is taken to be a simple hopping Hamiltonian: (2.73) and we assume that only site 1 is coupled to the bath, so that: Chapter 2. Dynamics of open systems 18 In accordance with our previous discussion, if 7 = 0 so that site 1 is coupled to the bath while site 2 is completely independent, then In this case, the reduced density is diagonal in the basis ci,C2 (which are the eigenstates, if 7 = 0). Site 1 will equilibrate while site 2 will remain forever in its initial state. The "local coupling scheme" assumes that the dynamic equation is described by Eq. (2.74) even if 7 ^ 0. This assumption is also generalized to multi-site systems coupled to multiple independent baths, where for each site coupled to a bath, a term of the general form of Eq. (2.74) is added in the dynamic equation. The "argument" for using this scheme is as follows: in Eq. (2.27) we only need to know the interaction Hamiltonian, i.e. which system operators (gen-erally called Lj, but here representing site creation or annihilation operators) are coupled to which baths, and this fully determines the Lindblad form. Since for the two-site system this interaction Hamiltonian is the same whether 7 = 0 or not, the Lindblad form should be the same in both cases. Generalization to more complicated systems should follow in the same way. The argument is wrong, however, because the dependence on the system operators in Eq. (2.27) is through their form in the interaction representation, and this depends explicitly on the system's Hamiltonian. In particular, for the two-site system, and using these in Eq. (2.27) will also give rise to terms like C]. ,C\ as well as cross-terms, in the dynamic equation. Thus, the Lindblad terms kept in the "local coupling" scheme are a subset of all terms that arise in the general formal-ism. There is no general justification why the neglected terms are unimportant to the dynamics of the system, and therefore the relevance of the local coupling scheme is questionable. The correct dynamic equation can be obtained using the general approach we sketched in the last section. We first diagonalize the Hamiltonian of the system, d_ dt Rs(t)= [nic^Cl +r2{c1)Cli]Rs. (2.74) (2.75) (2.76) where c i + c 2 c i - c 2 (2.77) V2 ;&2 V2 1 Chapter 2. Dynamics of open systems 19 In this basis, the coupling between system and bath becomes: v = Ek(xvkClb{ + \vk*c\bk) (2.78) Here we see that the two eigenmodes 0,1,0,2 are coupled with the same bath. This is exactly the type of case we pointed out before, where we need to also consider the operator involving both 01,02. With the notation introduced in the previous section, and after going through all the steps, we find C = Cdiag + £°ff, where CdiagR3 = r1(a1)C1ai + r2(a1)C2ai + n(a 2)£i a + r2(a2)C (2.79) and £°ffRs = e- i 2 7 t (ri(a 2) [ai,i?s4] - r 2(a 2) [ai,oJi2s]) +e i 2 7 t (r 2(a 2) [a{,i?5o2] -ri(o 2) [al,a 2i? s]) +ei27* (ri(ai) a2,RSa[ - r2(a1) a2,a\Rs ^ • f e - " 7 4 ^r 2(ai) a]2, Rsa^ - n(ai) 4, ai#s])-(2.80) It is now straightforward to check that the correct equilibrium state is reached by the system. Rewriting this Lindblad operators in terms of C\ and c2 operators, one can verify that this dynamic equation is different from that of Eq. (2.74). We con-clude that the "local coupling" scheme is wrong, and that the terms appearing in the Lindblad form depend explicitly on the Hamiltonian of the system. C h a p t e r 3 20 S p i n c h a i n s c o u p l e d t o t h e r m a l b a t h s Although several types of models (electrons on a lattice, coupled harmonic oscil-lators, etc.) could be used as test cases, in the following we will focus on chains of spins ,^ for concreteness. Such systems can also be used to model ensembles of coupled two-level systems (qubits). The Hamiltonian of such systems cannot generally be reduced to a diagonal form in creation and annihilation operators, as in the examples discussed before, and some generalizations are necessary. Also, for such systems, a natural definition for the local temperature and even for the heat current, is not yet available in current literature. We propose such definitions here. We first discuss the process of relaxation towards thermal equilibrium for a single spin, in Section 3.1. In Section 3.2 we consider an Ising chain, and show that Ising coupling is generally not sufficient to allow the system to evolve to its thermal equilibrium state. We then consider chains coupled to multiple thermal baths. In Section 3.3 we show that for baths set at the same temperature, the system reaches equilibrium in the steady state. In Section 3.4 we analyze chains coupled to baths set at different temperatures, and calculate the resulting heat currents and temperature profiles. 3.1 Single spin coupled to a heat bath The systems discussed as examples in the previous chapter have Hamiltonians expressed in terms of fermionic/bosonic operators with simple anti/commutative algebras, allowing for relatively straightforward diagonalization. In contrast, spin operators have more complicated commutation relations [<r_, a+] — — 4<r2 ^ 1 so that generally their Hamiltonians cannot be written as a diagonal form of annihilation and creation operators. One can use various bosonic representa-tions (such as Schwinger bosons, for instance) to do so, but in this case one gen-erally needs to add supplementary constraints to keep the states in the physical sector. If the constraints are built-in, such as for the Holstein-Primakoff repre-sentation, then Hamiltonians become very complicated. Of course, there are a few cases, such as the ID X Y model, where a diagonal form in terms of fermionic operators can be found using the Jordan-Wigner transformation [15]. However, in the more general cases this is not possible, and as a result we wil l work di -rectly with spin operators. We begin with investigating a single spin coupled Chapter 3. Spin chains coupled to thermal baths 21 to a thermal bath, and deriving the equation of motion for its reduced density matrix, following the formalism introduced in the previous section. We take the Hamiltonian of this system to be: Hs = -bz<rz, (3.1) where bz can be positive or negative, depending on the direction of the applied external magnetic field. The heat bath is modelled as an ensemble of bosons HB — y^vkblbk, k where uik > 0. The system is coupled to the bath through V = 5] AV f ccr x ® (b\ + 6 fc) . (3.2) it This coupling has all four possible types of terms, a±bk and <J±bk. This is necessary so that we can treat simultaneously both signs of bz. If bz > 0, the ground state is the | j) state, and processes allowed by conservation of energy are those described by <7_6fc (excitation with absorption of a phonon) and <J+b\ (relaxation with emission of a phonon). If bz < 0, the ground-state is | J.) and the allowed processes are described by cr+bk and cr-b\.. For either case, as we show in the following, the other two corresponding terms give no contribution to the dynamic equation, and it is therefore expedient to include all of them in the interaction Hamiltonian. We can then identify: Lk = Wkax,Fk =b[ + bk. (3.3) We now find the interaction representation for these operators. For the spin operators, these are: . aI±{t)=e*i2b>t<j+, (3.4) so that ^ ( t ) = ^ [ e - i 2 6 - V + + e i 2 6 « t a _ ] . (3.5) Note that we set H = 1 for simplicity. For the phonon operators, we have bfc(0 = e - ^ ' b f c , {blY{t) = eiuJktbl and the only finite Green's function is Gfcfc(r) = trB [(e^bl + e'1^bk){bk + b{)pe^] = e^T(nk) + e ' ^ ^ l + (n fc» (3.6) where (nk) = {el3uk - l ) " 1 . The dynamical equation (2.27) for the reduced density matrix becomes, in this case (we again let t —> oo in the upper integration limit, since we are Chapter 3. Spin chains coupled to thermal baths 22 interested only the in the steady state solution) d&s (t) ^ A 2 | V f c | 2 / dT([ai(t),oi(t-T)Rs(t)}Gkk 1 k J o . (r) dt h? ^ • - J0 • ' " ' ' (3.7) -[^(*),^(*)^(t-r)]G f c f c ( -T)) The time integrals can now be performed directly. We again keep only the real part, which are 5-functions. Given the T-dependence of cr^t — r ) and of Gfcfc(±r), we obtain only terms proportional to either S(2bz — wk) or S(2bz +u>k). If bz > 0, the latter type of terms are identically zero, whereas if bz < 0, the former type of terms are identically zero. These are the terms that correspond to processes that violate conservation of energy, so their disappearance from the dynamical equation is expected. For specificity, let us assume that bz > 0. In this case, after defining h2 and r = ^2J2\Vk\2S(u,k-2bz) k r1=r[l+n(2bz)},r2 = rn(2bz) where • is the average number of phonons of energy cok = 2bz, we finally find: — = r i [2a_Rs (t) a+ - v+*_R's (t) - Rs (t) ( g g ) +r2 [2o-+Rs (t) cr_ - <T-O-+Rs (t) - Rs (t) . It is again straightforward to see that if the system starts in either of its eigenstates, the reduced density matrix remains diagonal at all times so that Rg = Rs- Then, this equation has the stationary solution ' o ( n o ] _ n II) (II +r-2 IT) (Tl _ U) (II + |T) (Tl .... S ( j " n+r2 - 2cosh(/?6 2) • ( 3 ' 9 ) This is indeed the correct equilibrium state for a spin thermalized at this tem-perature. The evolution towards equilibrium can also be investigated. 3.2 Ising spin chain coupled to a heat bath In this section we will demonstrate that pure Ising coupling between the spins in the system is not sufficient to allow the system to evolveu to equilibrium. Chapter 3. Spin chains coupled to thermal baths 23 Consider a chain with N spins coupled through Ising exchange: Hs = J ((Ti,za2,z + h CTJV-I.ZO-JV.Z) • (3.10) Addition of external magnetic field along the z axis wil l not change the remaining arguments. We assume that only spin 1 is coupled to the spin bath, through the same interaction used in the previous section. Since the interaction is the same, the derivation follows the same steps until we arrive at: 9JaT1 = 4 E A W r d r ^ i ^ ^ t - T W ^ G ^ T ) d t h k Jo (3.11) - K x ( 0 , ^ ( * ) ^ . x ( t - T ) ] G f c f c ( - r ) ) Before continuing with this specific problem let us simplify this equation to the final form we wil l use in this and following sections. Consider the first of the 4 operator products appearing in the integrand. Using the relationship between the Schrodinger and interaction pictures, it can be rewritten as: < r U ' K * ( * - T)Rls(t) = eiH8t<TltXa{tX(-T)Rs(t)e-iH^ and the 3 other terms can be manipulated similarly. We now multiply the whole equation by e~lHst from the left, and by elHst from the right, to find in the end: dRs (!) _ i. . „ — d t ~ - - h [ H s ' R s { t ) ] 1 f°° " M E A 2 l V f c | 2 / d T ([^i,x(-T)Rs(t) ~ <d,x(-r)Rs(t)<TllX]Gkk(T) k J° - [<ri,xRs{t)<Tl<x(-T) - Rs{t)*{iX(-T)<TliX]Gkk(-T)) (3.12) Note that this equation is very general - it applies to any system that is con-nected to a single bath through spin 1, with the type of interaction Hamiltonian considered. The details about the specific nature of the system and its effects on the extra Lindblad terms are all embodied in cr[x{— r ) . For this particular system, one can find by direct calculation that: °{,X(-T) = \ [ e - f f l V p + eiQtam] (3.13) where we introduced the shorthand notation: Cp = CT1]X + iaity<J2,z Cm = 0\,x — iV\,yaiyz and fi = 2 J . The r-integrals can now be performed as before. We first take Chapter 3. Spin chains coupled to thermal baths 24 t —» oo in the upper limit, since we are interested only in the steady-state limit. The real part of the integrals are simply 5-functions. If we assume JZ,Q > 0, the final dynamic equation is: dRs (t) _ %.„ p , v, r -^—---[Hs,Rs(t)}--x [<rliX (n( f i )a p + (1 + n(Q))om) Rs(t) - (n(Q)ap + (1 + n(Q))om) Rs{t)irltX - o-hxRs(t) (n{n)am + (1 + n(Q))ap) Rs(t) (n(fi)crm + (1 + n{Q))ap) ahx] (3.14) The first observation is that the extra Lindblad terms are now very different than for the single-spin system, showing again that a "local coupling" type approximation, which has Lindblad terms only containing spin matrices for spin 1, but not 2, is wrong. The second observation is that the operators of spins 3 , 4 , N do not appear in the extra Lindblad terms for this Ising coupling, and as a result those spins remain forever in the state they were initially prepared in. Spin 2 is not changed either, because all the Lindblad terms in the above equation depend only on 02,2, so states where the spin 2 is "up" are never coupled to states where it is down. The only spin that can be changed as a result of the coupling to the bath is spin 1. It is again straightforward to see that in the t -> 00 limit, Rs becomes diagonal, and the coefficients multiplying the various diagonal terms can be found. If we assume that at t = 0 the system has o~2,Z = +1 with probability p, then one finds in the end that Rs (00) = o w o T \ I Tl)(Tl I + 7, • , A T N I I I T2) {T2 2cosh(pJ z ) 2cosh(pJ 2 ) ' + ( 1 - p ) ( 1 T l ) ( T l 1 + w T O 1 1 1 1 i 2 ) ( i 2 Rs,N-2(0) (3.15) where RS,N-2(0) represents the initial state of the other N — 2 spins, and can be obtained from Rs(0) by tracing out spins 1 and 2. This shows that only spin 1 is thermalized to the correct temperature, and that spin 2 simply acts as an effective "magnetic field" for spin 1. Spins 2, N remain essentially in their initial state. This is obviously due to the fact that no term in Hs allows these spins to be flipped, so that they can change their state and evolve to equilibrium. One way to achieve this is to couple each spin to its individual bath. Clearly, each of them wil l equilibrate to the temperature of its own bath, with Boltzmann factors dependent on the "effective" magnetic field created by neighboring spins. No realistic study of heat transport is possible in such a case. In order to study heat transport, we need to have a system Hamiltonian that allows for thermal equilibration to begin with. This can be achieved by allowing for xy coupling between neighboring spins. Chapter 3. Spin chains coupled to thermal baths 25 3.3 X Y chain coupled with one heat bath: thermal equilibrium A spin-^ X Y chain has the following Hamiltonian, Hs = 53 xy\Vi,x&i+l,x + &i,y&i+l,y) , (3.16) i which can be written in terms of spin raising and lowering operators as Hs = 53 §Jxy (o"t ,_0 - i + i i + + Cj,+cr i + l i _), (3.17) i with the usual commutation relations [<Ti,+,^,_] = 4 a i i Z J < \ [ a M , a J - i ± ] = ±2(7 i i±J«. (3.18) Although for a pure XY-chain we could use the Jordan-Wigner transformation to reduce its Hamiltonian to a simple non-interacting fermionic model, we prefer to continue working directly with spin operators here. One reason is that in terms of Jordan-Wigner operators, the coupling of a spin to its bath can become quite complicated. More generally, we want to be able to easily generalize to cases with external magnetic field and/or anisotropic Heisenberg interactions. In this section we study X Y chains and present results for chains with up to N = 4 spins. The calculations shown were performed in Maple (hence the rather small chain size; larger systems can be treated similarly, but one needs to write programs to solve the resulting differential equations numerically). The model studied in this section assumes again that only the first spin is coupled to a heat bath Hg = J2kwkb\.bk through the usual coupling: V = ^2xvk^i,x®(bl + bky (3.19) fc As discussed, this leads to the general Lindblad form of Eq. (3.12). The chal-lenge now is to find the expression of o[ x (r). 3.3.1 Sample calculation for a 2-site X Y chain In this section, we use the 2-site chain as an example to show a typical calcu-lation. Larger systems are treated similarly, but notation becomes too cumber-some to discuss here. First, we perform the direct diagonalization of the central system. We rewrite the Hamiltonian in a matrix representation, using the 2 2 usual base Chapter 3. Spin chains coupled to thermal baths 26 vectors {|TT>, |U>, |1T>, | U » , H = 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 (3.20) Let Q = Jxy/2. After diagonalization we find the eigenvalues and eigenvectors, E = [0,ft,-f2,0] , 1 0 0 0 0 V2 -y/2 0 0 y/2 \ / 2 0 0 0 0 1 (3.21) (3.22) The evolution operator U (t) is, in the eigenbasis representation, U(t) -iHst 1 0 0 0 0 <e-itn 0 0 0 0 eitn 0 0 0 0 1 (3.23) This allows us to compute the interaction picture representations of the various spin operators in the eigen-basis, in particular o(x(t) = U (t) aiiXSU^ (t), where the representations of these operators in the initial basis. For example, we find 0 eiClt 0 0 0 e-iQt 0 0 _ e - i Q t 0 g — iQt _eiQt 0 (3.24) While this example is trivial, direct generalization of this approach to larger systems and more complex Hamiltonians is straightforward. We can now use this expression of cr / x (—r) in the dynamical Eq. (3.12): 9 R S (0 _ * rrr n M ] i r°° ¥ £ A2|^!2 / d T (Wi,*°i,x(-TWt) - olx(-r)Rs(t)ahx}Gkk(r) i. Jo - [oltXRs(t)alx(-T) - Rs(t)<rltX{-T)crliX]Gkk(-T)), Using now the usual Green's functions Gkk(T) = exu>kT{nk) + e~iuJkT(l + {nk)), we see that the r integrals can now be performed. If we assume that Q. > 0, then only the terms proportional to 5(u>k — ft) contribute. The dynamical equation Chapter 3. Spin chains coupled to thermal baths 27 then reduces to: dRs (t) dt = --[Hs,Rs(t)] -r [crilXmxRs(t) - miRs(t)cr±iX - crltXRs(t)m2 + Rs(t)m2critX] (3.25) where m i m 2 71 J _ 72 0 n(fi) 1 +n( f i ) 0 0 1 + n(£2) 0 1 +n(f i ) 0 0 l+n(Q) n(fi) 0 0 1 0 0 - n ( f i ) f n ( Q ) 0 0 n(Q) - l - n ( f i ) 0 n(ft) - 1 - n ( f i ) 0 0 1 + n(il) - n ( f i ) 0 (3.26) (3.27) In other words, the matrices m i and m 2 are similar to <?[ X ( — T ) , with e i r 2 T —> 1 + n(fi), e _ l f 2 r —> n(fi) for m i , and viceversa for m 2 . A l l the other operators must also be represented in the eigen-basis, in this dynamical equation. As before, the rate is: k r = h2 Q) and n(ft) = After performing the matrix multiplications, the equation of motion can be solved by direct integration. If we are only interested in the thermal equilibrium state, we expect the reduced density matrix to be diagonal in this eigen-basis representation, RS(oo) P i 0 0 0 0 Pi 0 0 0 0 P3 0 0 0 0 Pi (3.28) with the normalization condition YltPi = 1- Since now [Hs, Rs (oo)] = 0, it follows that the steady-state solution is obtained by requesting that the term multiplied by r in the above dynamical equation becomes identically zero. After performing the matrix multiplications, the solution is found to be: 1 Pi = Pi = 2 + 2cosh(/3fi)' x '" 2 + 2cosh(/?ft) i.e. the expected thermal equilibrium. ;P3 = 2 + 2 cosh(/3fi) : Chapter 3. Spin chains coupled to thermal baths 28 Right Batli Figure 3.1: A 3-spin X Y chain is coupled to two heat baths with temperature respectively TL and TR. We want to discuss the heat current on such chain. Note that the time-evolution can also be obtained from this dynamical equa-tion, by direct integration. In the long-term limit the off-diagonal matrix el-ements vanish, and the diagonal ones become equal to the Boltzmann proba-bilities. The only significant evolution is seen close to t = 0, after which the system quickly arrives at equilibrium. However,'in that limit, extending the integral JQ dr —> f0°° dr is not justified, so the relevance of these results is some-what unclear. Of course, one could solve the full integro-differential equation (with and without the Markovian approximation) to understand the evolution towards equilibrium/ steady-state. We postpone such studies for future work. 3.4 X Y chain coupled to two baths In order to study transport properties, the system of interest must be coupled to more than just one bath. Here we analyze short xy chains with and without external magnetic fields applied along the z axis. The iV-site chains have the first and last sites coupled to two different thermal baths, as sketched in Fig. (3.1). Let TL and TR be the temperatures of the two baths. If TL — TR we should recover the thermal equilibrium in the steady state limit, however if the baths have different temperatures, the steady state solution wil l be a non-equilibrium state in which heat should flow through the system. Since the two baths are different, the Green's functions involving operators from both of them are zero, so that the Lindblad terms are simply the sum of terms for each individual bath. In other words, the dynamical equation is now [compare with Eq. (3.12)]: 5 ^ - > . . « . ( « „ - ^ £ A 2 I ^ I 2 / ^ ( K X * ( - T ) ^ - [<Ti,xRs(t)<r{,x(-T) ~ Rs(t)alx(-r)aliX}GkLk\-T)) - WN,xRs(tWNiX(-T) - Rs(tyN<x(-T)aN,x]Gif(-T)) (3.29) Chapter 3. Spin chains coupled to thermal baths 29 For simplicity, we assumed that both baths have the same phonon spectrum uik and coupling constants Vk, and use XL and XR to modulate the overall rates ri, and rR for the two baths. The Green's functions are as before, but each is therrrialized at its own temperature: GiLK/R\T)=e^(nk)L/R + e-^{l + (nk)L/R) -where As discussed, the information about the specific system considered is all embedded in cr{x(—r) and cr^ x (—r) . These are calculated in the same way shown for the single bath case. We quickly sketch the main steps for the 2-spin chain as an illustration, and to explain the differences in the calculation of the steady-state solution. 3.4.1 Sample calculation for a 2-site X Y chain The diagonalization of Hs and the resulting expression for cr[x(—r) are the same as previously. After extending the upper integration limit to infinity, the T integrals for the terms coming from the left bath are performed asn before and one obtains terms similar to the ones in Eq. (3.25)", except that r is replaced by r/, = X2Lr/X2, we rename m^j —> a n d n(f2) —> n i ( f i ) = (exp(/?i$7) — l ) - 1 . We deal similarly with the terms coming from the right bath. In this case, 1 71 0 e-int _£int 0 gifit 0 0 e-iQt 0 0 0 e-iSlt 0 (3.30) The r integrals can now be performed, and as before result in replacing the exponentials with either nR(Q,) or l+nR(fl), where nR(Q) = (exp(/?2fi) — l ) - 1 -The final dynamical equation is: c^.xTOi^-RsM - wi J i ; i?s(*)o '2,x - o-2,xRs(t),m)2n> + Rs(t)mK2n'<J2,. AR) (R) AR), (3.31) where J _ 71 0 l+nR(Q) -nR(Q) 0 n f l ( f i ) 0 0 nR(Q) -(l + nR(Q)) 0 0 lnR(Q) 0 l + n f i (Q) n f l (Q) 0 (3.32) Chapter 3. Spin chains coupled to thermal baths 30 and m 2 K ) is obtained from m[H- by exchanging TIR(Q,) <-> 1 + nn(Q). For longer chains, the structure of the dynamical matrix remains as in Eq. (3.31), but the matrices become 2N x 2^ dimensional. They can be generated efficiently, however, and then one can perform the needed matrix multiplications. One essential difference for longer chains, or even for a two site chain with finite Jz or bz terms, is that in these cases there are many different excitation energies, not just one value ft. Phonons can be absorbed/emitted in exciting or relaxing the system by such transitions, so in general there will be multiple effective rates to each bath, of general form r{5E) — TrA 2 /^ 2 ]Cfc \Vk\2S((Jk — SE). In the following we assume that all rates of coupling to a bath are the same, irrespective of the value of 5E, i.e. r(SE) — r. This is true if the density of states of phonons is flat and |Vfc| =const. More realistic models will be considered in future work. The excitation, energy 6E of different transitions still enters the dynamical equation through the phonon occupation numbers n(SE), so it remains relevant. Unless specifically stated otherwise, we also assume that the baths are identical so that = TR. In the steady-state limit, we request that ARQ^ —• 0, and thus we need to solve the system of 22N coupled linear equations described by the right-hand-side of Eq. (3.31). First, let us briefly analyze the case TL = TR. In this case, we know that taking either TL —> 0 or TR —> 0 gives the correct equilibrium diagonal density matrix. It is then clear that having both , rR 0 wil l also admit the correct equilibrium diagonal density matrix as a solution. Going through the actual calculation shows that this is the only possible solution. For TR ^  TL, however, we do not know that the density matrix should be diagonal. In fact, as we argue now, in such a case the density matrix must have of f -d iagonal terms. The reason is simple. Suppose that j is some current operator (of heat, or more generally charge, spin, etc). Since we are working in the system's eigen-basis, it follows that if the reduced matrix is diagonal, the expectation value of the current operator in the steady state is: (j) =trs(jRs) = 5 >nHJ | n ) n where \ri) label the different eigenstates of Hs- However, (n\j\n) = 0 in this case, since there cannot be any finite current in any eigenstate of the isolated Hamiltonian Hs (degenerate states can always be chosen so as to satisfy this requirement). It follows that if the steady-state density matrix were to be di -agonal, transport of any kind would be impossible, irrespective of the system considered or the bias applied to it. This is obviously an unacceptable con-clusion, and therefore we must have a steady-state density matrix with finite off-diagonal matrix elements if the system is biased. Thus, we must solve the full Eq. (3.31). For the time-dependent case (as-suming that the extension of the upper limit of the integral is acceptable) this gives 22N coupled first order differential equations, which can be integrated by usual means. In fact, if we define P(t) as the 22N vector containing all matrix Chapter 3. Spin chains coupled to thermal baths 31 elements of Rs{t), the dynamical equation can be recast in the standard form: ±P(t) = rP(t) The matrix T is non-hermitian. If we denote by Ui its right-hand eigenvectors, Tui = 7»ttt, it follows that if P(0) = YJ» ctiUi, then P(t) = J2i a.ielitUi. Clearly, the steady-state limit corresponds to the eigevector of the zero eigenvalue 7 = 0. A l l other eigenstates must have negative real part. While we have not attempted to understand on general grounds why this matrix should only have one zero eigenvalue while all other ones have negative real parts, our calculations always show this to be the case. As mentioned before, in this thesis we show results for chains with up to N = 4 spins, for which the dimension of T is 256x256. The reason is that the calculations were performed with Maple, which is not optimized for such applications. However, given the nature of this calculation, one can expect to be able to quite easily find the steady-state for systems with up to at least 14 spins, since this implies finding the eigenstate corresponding to the highest eigenvalue (known to be precisely zero) of a sparse matrix of dimension 16, 384 x 16, 384. This is a quite straightforward task if one uses the optimized libraries available for such problems. It is to be hoped that such sizes are large enough to allow for some understanding of bulk-limit properties, if that is what one is interested in. Also, adding couplings to extra baths is straightforward and does not increase the computational effort significantly, so one could just as easily study (shorter) chain ladders, for example. As a final observation, we see that if the system is biased, the steady-state density matrix becomes dependent on the ratio of r^ / r^ . It obviously goes to one or the other equilibrium solutions if this ratio is either zero or infinity, however for finite ratios one can study the dependence of transport properties on the coupling to the baths. As an empirical demonstration that the steady-state for a biased system is not diagonal, we plot in Fig. 3.2 the time evolution of some of the off-diagonal matrix elements of Rs(t), assuming that the system was initially in the ground-state of the isolated chain. Once the steady-state reduced density matrix is known, we can proceed to calculate various quantities of interest. Various current operators can be defined in the usual way, however if we are interested in heat transport and Fourier's law, we should also be able to calculate the "local temperature" of each spin. A generally accepted definition of this quantity seems to be absent in the current literature. In the next section we discuss our approach to this problem. Chapter 3. Spin chains coupled to thermal baths 32 Figure 3.2: Time-dependence of some of the off-diagonal matrix elements of Rs{t), for a 2-spin chain (left) and for a 3-spin chain (right). Parameters are Jxy = 4, bz = —0.1. The chains are coupled to baths with TL = 5,TR — 1 and TL = TR = 1. Clearly, the steady-state density matrix has finite off-diagonal matrix elements in this case. 3.5 Definition of the local temperature and heat current For classical particle systems, the temperature of a particle is related to (p 2 ) /2mj , allowing one to determine it easily from a knowledge of the distribution function. Unfortunately, for a spin chain we do not have such a simple and intuitive defi-nition. One may calculate the temperature Ti of the ith spin in the steady-state limit by finding the solution of the equation trs(aiiZRs(oo)) = tanh(ft&z) (3.33) if a magnetic field is applied along the z axis. This definition is used in literature, for example in Refs. [8, 9], however it seems arbitrary to expect the spin to behave as if it is free, since supplementary bias from neighboring spins is likely to occur in ordered chains. (Also, this quantity is not well defined if bz = 0, although one could calculate it from interpolating bz —> 0). The biggest problem with this definition is that it does not predict Ti = T if one sets both baths at the same temperature TL = TR = T, since at equilibrium, the spin expectation values on an interacting chain are different from those of free spins. We have demonstrated that if TL = TR = T the spins are in equilibrium at temperature T, therefore we have to reject this approach as unsuitable. As mentioned before, another approach used in literature is to couple every spin to a bath, and adjust the temperatures of the extra baths until no heat is exchanged between a bath and its own spin [10]. This approach is appealing in that it mimicks a real temperature measurement. However, one must make sure that the couplings to the "thermometer" baths are so weak as to not significantly Chapter 3. Spin chains coupled to thermal baths 33 modify the behavior of the system. A quantitative definition for this condition is not simple (in fact, in literature this approach seems to be used primarily to introduce significant incoherent scattering for each spin, and thus to modify the properties of the chain considerably. It is therefore unclear whether one learns anything about the properties of the isolated chain coupled to only two leads, in this way). Finding the "self-consistency" condition of no heat flow to and from the "thermometers" is also a non-trivial task, especially for longer chains. Our choice of definition for the local temperature is to follow the idea of Eq. (3.33), however to use the proper expectation values expected for a system in equilibrium. For example, if we add interaction with an external magnetic field Hext = —hz CTi,2 to the xy Hamiltonian, then for the two-spin xy chain equilibrated at temperature T, one expects, by symmetry: (O-\,Z)T = {O2,Z)T = , , o f l , ^ , • , „ m (3.34) cosh(2po2) + cosh(piZ) which is quite different from the tanh(/362) of Eq. (3.33). Of course, these equilibrium expressions differ depending on the number N of spins on the chain, and, for chains with N > 2, also on the position of the spin along the chain. However, since we diagonalize Hs as part of the overall calculation, finding the equilibrium expectation values of each K 2 ) T = ! * r 5 ( < 7 M e - ^ ) (3.35) is a straightforward task. For non-equilibrium steady-state systems, we then use the condition trs{oitZRs (oo)) = (o-iiZ)Ti as one way to calculate the temperature Tj of that spin. This equilibrium expectation values are shown in Fig. 3.3, for 2 and 3 spin chains placed in small, respectively large magnetic fields. For the large magnetic field (right panel) the differences between the equilibrium values of Eq. (3.35) and those given by Eq. (3.33) are not so significant, although they are non-negligible. As a result, one expects that using Eq. (3.35) instead of Eq. (3.33) wil l not change things dramatically. However, for small magnetic fields, the differences are very significant unless temperatures are very large. In particular, we see that (o ri, 2)r is very small at some sites, and is not a monotonic function of /?, so that the solution of for the local temperature is not unique. In such cases, we choose the solution closest to the temperature of the bath near which the spin is located. This approach is at least certain to predict the correct temperature profile in equilibrium. However, if there is no external magnetic field, all these expectation values vanish. Also, there is no obvious reason to favor the spin operators <Tj ) 2 to cal-culate the local temperatures; any other operator could be used just as well. Given the xy symmetry, the other two Pauli matrices c r i x , o i y would give a dif-Chapter 3. Spin chains coupled to thermal baths 34 Figure 3.3: ( C T J I Z } T as function of the inverse temperature, for (a) two-site chain, and (b) three-site chain. The parameters are bz = —2, fl = 1. In both cases, the tanh(/?62) expression used in Eq. (3.33) is also shown (black line). Chapter 3. Spin chains, coupled to thermal baths 35 ferent equation for than Eq. (3.35), however these operators would only have finite expectation values if one applies an in-plane magnetic field bx. The next possible choice involves products of Pauli matrices from different sites, which can be defined in many possible ways. Since we are interested in heat transport, we introduce a second way of calculating the local temperature, using the "local Hamiltonian" as the test operator. Specifically, we divide the Hamiltonian of the total system in "local Hamiltonians": N #s = X > (3-36) i=i where each hi contains all terms of Hs that depend on operators acting at site i only, as well as one-half of the terms involving bonds between this site and other sites it interacts with. Specifically, for an xy chain in external magnetic field bz, we have: hi = —bzGiiZ + 2~^~ {ai'x (°"*+i>* + Ct-i.x) + ai,v (°»+i,!/ + °i-i,3/)] • (3.37) For sites i — 1 and i = N, there is coupling only to spins 2, respectively N — 1. As before, we now calculate the expectation values of these operators if the system is in equilibrium, (hi)T = ^trS(hie-PHs). (3.38) Again, this is a straightforward task. Curves of ( / I J ) T as a function of T are shown in Fig. 3.4 for 2 and 3-site chains. These are always monotonical func-tions, unlike the (cri ) 2)r values. As before, in a non-equilibrium steady-state, we use the condition: trs(hiRs{oo)) = (hi)Ti to find the local temperature Ti, which is now unique. This approach also gives the correct solution in the equilibrium state. Moreover, it works even if bz = 0, which is a definite advantage. The only problem is if the system is very cold, since then it is difficult to find T, accurately. In practice, we calculate the local temperatures by both these approaches. As we show in the next section, if bz is not too small, there is rather good agreement between the two predictions, which suggests that this is a reasonable way to define the non-equilibrium temperature profile. For small bz values, the approach based on (hi)x always gives reasonable results, but the values based on the (oi,z)r a r e sometimes clearly wrong. We plan to follow up on the idea of adding tiny "thermometer" baths to measure local temperatures by yet another approach, in future work. The definition of the local heat current is also not straightforward. If we were to measure a charge current in an electronic system, we could trivially define the density of charge operator at any site, and use an equation of motion Chapter 3. Spin chains coupled to thermal baths 36 Figure 3.4: Expectation value of the "local Hamiltonian", (/IJ)T as a function of inverse temperature, for (a) two-site chain and (b) three-site chain. The parameters are .... Chapter 3. Spin chains coupled to thermal baths 37 and the conservation of charge to derive the expression of the electric currents flowing on each link between sites. The trace of this operator weighted by the reduced density matrix gives the local current. We try to proceed similarly here. We start from the "local" energy at a site, as defined by hi, and use a continuity equation for conservation of energy Yt + V / = 0. (3.39) For a linear chain, we find: - j l ^ i = ~ ^ = i [HsX] • (3-40) where we assume a lattice constant a = 1. This allows us to identify the heat currents for different bonds.. For a pure xy model (bz = 0), the heat current is: •h _  1 I r \2 / Ji_>i-|_l — ^ \Jxy) \0~i,-0~i+\tzO~i+2,+ — 0~ii + CTi+iizO~i+2,-+0~i-l,-0~i,z°~i+l,+ — 0 i - l , + < 7 t , z 0 ' t + l , - ) For i = 1 or N — 1, terms including co,±, &N+2,± should be removed. In these cases there are also terms involving bath operators. We simply set these to zero. If the baths were truly in equilibrium, these terms would vanish since they are proportional to b + tf combinations. However, if the baths are in equilibrium there would be no overall energy flow into/out of the system. This suggests that currents on the outer links may not be exact, but the ones calculated for inner links, where bath operators do not appear, should be correct. In the stationary state we should find (Ji_i+1) = So one could also define an average heat current as {jh) = E^(kli±A. ( 3 . 4 2 ) Here we exclude the heat current on the end bonds between the end spins and their baths. One peculiar aspect of this definition is that for a 2-site chain and bz = 0, we cannot define a heat-current, since in this case h\=h2= 0.5H$, and therefore the commutator [Hs, hi] = 0. This is also apparent from the expression of Ji-^i+i which involves only products of three neighboring spin operators. This is rather surprising, since one expects some heat to flow from the hotter bath towards the colder one, however it is not clear how this can be measured. If bz is finite, there is a second contribution to Ji_i+i that is proportional to bz, which we call the "spin-current". This describes the rate of change of each <7j i2 and is derived similarly to the heat current. (3.41) Chapter 3. Spin chains coupled to thermal baths 38 S If a. <l A it * * * Figure 3.5: Temperature profile for a three-site chain in the stationary state, as a function of the temperature TL of the left bath, for TR = 1 respectively T R = 3. The xy coupling between sites 1 and 2 is J\2 = 0.2, while the xy coupling between sites 2 and 3 is 0, 0.01, 0.1 and 0.2, respectively. The applied magnetic field is bz = 2. 3 . 6 Current and local temperature profiles in non-equilibrium steady-states In this section we present some of the steady-state results we have obtained so far. In all cases shown here JZ = 0; there is only xy coupling between neighbor-ing spins and the applied magnetic field bz is uniform. A full investigation of the entire parameter space is still yet to be completed, and study of longer chains still needs to be implemented numerically. The conclusions we draw from these results may, therefore, change somewhat when more data becomes available. We begin with some simple tests to verify that our programs give reasonable results. One such check is that if we remove the interactions between two neighboring spins, then the left/right side of the chains should thermalize to the temperatures of their respective baths, TL/R. That this is true is shown in Fig. 3.5. Each of the four panels contains two sets of curves, showing the local temperature Tj as a function of the temperature TL, as TL is increased from 1 to 2, respectively from 3 to 4. In all cases, TR is kept at the lowest initial value of T R , i.e. TL = 1 or 3. If the exchange between sites 2 and 3 is removed (first panel) we see that T 3 = TR irrespective of the value of TL , whereas T\ = T2 = TL in all cases. This clearly demonstrates that the two halves are thermalized with their respective baths. If the coupling between sites 2 and 3 is slowly turned on (panels 2, 3 and 4), we see that T2 gradually moves away from T\. For TL = TR UD Chapter 3. Spin chains coupled to thermal baths 39 • Tl _ 1 « T2 / L L Figure 3.6: Temperature profile for a two-site chain in the stationary state, as a function of the temperature TL of the left bath, for TR = 1 respectively TR = 3. Parameters are Jxy — 0.2, bz = 2. The effective couplings to the baths are T\£, = 1 and rR = 1,0.01 and 0.0001, respectively. As the coupling to the right bath goes towards zero, the system thermalizes at TL-the entire chain is in equilibrium, but as TL increases, in the steady state the temperature varies monotonically across the chain, TL ~ T\ > T2 > T 3 w TR. As a second test, we change the ratio of the effective couplings to the two baths, rL/rR. As discussed in a previous section, the steady state reduced density matrix depends only on the ratio of these two quantities. In the limit when one of these couplings is zero, the chain is coupled only to the other bath and its steady-state should be equilibrium at that particular temperature. This is indeed what Fig. 3.6 reveals. The results are for a 2-site chain with Jxy = 0.2, bz = 2. As before, the temperature profile is shown for two fixed values oi TR — 1, respectively 3, as a function of TL- When the effective coupling to the two baths is equal (first panel) the local temperature of each site is equal to that of the bath it is coupled to. When = lOOr/j, the local temperature of site 2 starts to scale with TL, showing that the whole system, not only site 1, heats up. In this case, the right-side bath simply cannot cool site 2 fast enough. Finally, for ri = 10,000r/j, the system is effectively only coupled to the left bath, and both sites are thermalized at TL = T\ ~ T2. In the following, all results shown correspond to ri = rR. Finally, a third check (not shown) is that if we change the sign of the mag-netic field, bz —> — bz, the temperature profiles remain the same, only the sign of the expectation values {ai%z) changes. Chapter 3. Spin chains coupled to thermal baths 40 Figure 3.7: Spin current ( j f_ 2 ) f ° r a two-site chain in the steady state, as a function of the temperature of the left bath TL- The right bath is at TR = 1 (red line) or TR = 3 (black line). The parameters are JXY = 0.2, bz — —2. The heat current is 0'i_ > 2) = - ^ ( i i - , 2 ) -We now present results for the heat and spin currents in the non-equilibrium steady state through chains with N = 2, N = 3 and N = 4, in Figs. 3.7 and 3.8, respectively. For the two-site chain we only show the spin current (jf_, 2), but in this case the heat current is (jh) — —bz(js). In all cases we take Jxy = 0.2, bz = 2. The right-hand bath is kept at a constant temperature TR — 1 or 3, and the average heat current is shown as a function of the temperature of the left reservoir, for TL > TR. Of course, when TL = TR the system is in equilibrium and no current flows through the system. As the temperature difference increases, the currents increase first proportionally to TL — TR, as expected from linear response. However, for larger TL —TR the currents saturate and then begin to decrease. The reason for this is two-fold. First of all, these are finite, small systems, with few eigenstates. Therefore, there must be an upper limit to how much current can be transported through the system. Moreover, if the temperature of one bath is much larger than the temperature of the other bath, the system starts to thermalize with the hotter bath. The reason is that the coupling to a bath depends on products of the type rn(5E), where n(5E) is the number of phonons with energy SE in that bath, and 5E is the splitting between eigenstates connected by the coupling to the bath. Even if TL = TR, the very hot bath will have many more phonons available and will therefore be more effective in heating the entire system, than the colder bath with few phonons is to cool it. Thus, for a large temperature differential we expect the Chapter 3. Spin chains coupled to thermal baths 41 / -ti. / • • / — • — • • / A i i • •  • • A •f . --• m.._ • L '/  t* ' T V / L••• © i A A 4 *'-~Ir-...» • • • • • • Figure 3 . 8 : Heat and spin currents (jh) and {j"), for a three-site (upper panel) and four-site (lower panel) chain in the steady state, as a function of the tem-perature of the left bath Tj_. The right bath is at TR = 1 or TR = 3 . The parameters are Jxy = 0 . 2 , bz = — 2 . Currents on different links are equal to one another, as expected in the steady-state. Chapter 3. Spin chains coupled to thermal baths 42 steady state to be more and more s imi lar to therma l equ i l ib r ium w i t h the hot ba th , and of course the heat current th rough the system becomes smaller and smaller, in this case. In the l im i t i ng case TL > > T R , the entire system w i l l be at temperature TL and one expects heat current to flow on ly between site iV and the r ight ba th . (Note that such a s i tuat ion wou ld not violate the conservation of energy, since the " loca l H a m i l t o n i a n " of the end sites i n fact also includes the coupl ing to the leads, the var ia t ion of wh ich gives rise to therma l currents over and above those inc luded i n j h , and wh ich can affect the end bonds) . We cannot calculate the current on the outside l inks at the current stage, since that impl ies knowing how the coupl ing to the system changes the density matr ices of the heat baths away from their equ i l ib r ium values. For the longer chains, we plot the various currents on several bonds. A s expected i n the steady-state, these are equal to one another, be ing guaranteed to be so by conservation of tota l energy at the sites wh ich are not end sites. Compar i son of the results i n the two panels of F i g . 3.8 shows that these currents are rather s imi lar i n aspect and absolute magni tude to one another. Th i s suggests that f inite size effects are already rather smal l , and that these values are already rather close to the bu lk values that should be establ ished in longer chains. If we define a global "heat conductance" such that (jh) = GH{TL — Tn), these results suggest that G # is a constant for smal l A T , as expected f rom the Four ier law, and then it slowly decreases for larger A T . The second conclusion is that GH is independent of the length of the chain. Th i s is somewhat reminiscent of the electr ical conductance of a perfect meta l (for example, a chain w i t h no impuri t ies) wh ich is given by the contact conductance f rom the coupl ing to the leads, but there is no " in t r ins i c " cont r ibut ion f rom the meta l i tself since any electron that enters the meta l sample at one end is t ransmi t ted w i t h 100% p robab i l i t y to the other end. However, i n this case the conductance does vary strongly w i t h L, the length of the metal l ic cha in , i f the Fe rmi level is swept through the spectrum of the meta l l i c chain. T h e reason is construct ive and destruct ive interference between electrons mov ing back and forth between the ends of the sample, before they can tunne l th rough the contacts out into the leads. To obta in the rough equivalent of this i n our case, we should calculate G # | A T _ O as a funct ion of TL. F r o m the slope of the curves displayed i n F igs . 3.7 and 3.8 when A T —> 0, i t is clear that this quant i ty varies w i t h T L , however i t seems to be the same for a l l three chain lengths, therefore dependence on the chain length seems unl ikely. A second difference w i t h respect to electr ical t ransport , is that the heat conductance i n the steady state is independent of the magni tude of the effective coupl ing to the baths, on ly their rat io TLITR appear ing i n the equat ion that gives the steady-state density mat r ix . O n the other hand, e lectr ical conductance depends very strongly on the tunne l ing ma t r i x between leads and system; the smaller this is, the less l ike ly electrons are to tunne l into the system, and the smaller the electr ical conductance becomes. These results suggest that there are some simi lar i t ies but also some differ-ences between heat and charge t ransport . Th i s is to be expected, especia l ly for the pure xy mode l , wh ich can be mapped into a non-interact ing spinless fermion Chapter 3. Spin chains coupled to thermal baths 43 -• a • a o ma - • • --ii::: n _ • 0 ° n >lo . »0<q>.00 , Figure 3 . 9 : Temperature profile for a 2-site chain in the stationary state, as a function of the temperature TL of the left bath, for T R = 1 respectively T R = 3 . Parameters are Jxy = 0 . 2 , bz = 2 . A l l three measurements of the local temperature T, agree, and find TL ~ T i and T 2 « T R . model. However, the heat transport must be investigated in much more detail before a definite picture can emerge. We now investigate in more detail the temperature profiles in non-equilibrium steady state situations. We continue with a model with Jxy = 0 . 2 , bx = —2, T R either 1 or 3 , and different values of TL- Results for 2 , 3 and 4 spin-chains are shown in Figs. 3 . 9 and 3 . 1 0 . The local temperatures T, were calculated from using the three methods discussed in the previous section, namely from matching (i) trs(ai<zRs) = tanh(/3i&2); (ii) trs(aitZRs) = (&I,z)T, and (iii) trs(hiRs) = (hi)r- In all cases, the results agree for the three methods. The agreement with the tanh fit is somewhat surprising but not totally unexpected, since the magnetic field is rather large and if the spins are significantly po-larized, the differences between the tanh-function and the full {<Ji,z)T become small. The situation is much different for small bz, as we show next. Fig. 3 . 9 shows that each spin on the two-chain site stays close to the tem-perature of its own bath. This is to be expected, since the effective couplings to the baths are equal, and the temperature differential is not too large. The top panel of Fig. 3 . 1 0 shows that for a 3-site chain, the outside sites are still at roughly the bath temperatures, TL ~ T\ and T R ~ T 3 . The temperature of the middle site increases linearly with TL — TR. For the lower overall temperatures when T R = 1, we find T 2 « 1 / 2 ( T R + TL), however for the higher temperatures Chapter 3. Spin chains coupled to thermal baths 44 • m rj • CT * a • CT • • A • . ! * _• _o I^ffWvw I . I . I Figure 3.10: Temperature profile for a 3-site (upper panel) and 4-site (lower panel) chain in the stationary state, as a function of the temperature Ti of the left bath, for T R = 1 respectively T R = 3. Parameters are Jxy = 0.2, bz = 2. A l l three measurements of the local temperature Tj give the same results. Chapter 3. Spin chains coupled to thermal baths 45 • ff _ ; \ . . . » - •— a.. „ - \ . . .» _ 1 • • • — i - *••.., ' 9 • » • • • • • ••-•=: • • \ m-- \ X •- --\ ft - - • <T -— a. • » - 1 1 1 V 1 Figure 3.11: Steady-state temperature profile for 2, 3 and 4-site chains when TL = 4, T R = 1. The two data sets correspond to models with JXY = 4 and bz = 2, respectively Jxy = 0.2 and bz = 2. In the former case, the results of the three methods of estimating tt disagree, although there is convergence for the longer chain between estimates from (/IJ)T and from (oitZ)T- For the latter case, all three estimates agree with one another. when T R = 3, TZ is clearly less than the average bath temperature. Very similar behavior is seen for the four-site chain. In this case the middle sites have the same temperature T 2 = T 3 , and this is very close to the temperature of the middle site of the three-site chain. Based on this very limited data, extrapo-lation to longer chains would suggest a constant temperature along the inner part of the chain, with all temperature differential occuring at the ends of the chains. This is in direct analogy with the expected voltage fall across a sample that is a perfect conductor, in which case all voltage difference occurs across the contacts to the leads, not inside the sample. However, again one needs to do more investigations, in particular for longer chains, before this assertion can be supported. The situation becomes more difficult to analyze if the magnetic field is smaller, or the ory-coupling is larger. Fig. 3.11 shows the temperature pro-files for two, three and four-site chains, for a large temperature differential TL = 4,TR = 1, Jxy = 4 and bz = 2. This is contrasted with the results for the Chapter 3. Spin chains coupled to thermal baths 46 previous case with Jxy =0.2 and bz = 2. As for the previous data shown, for the latter case the three methods of measuring T give the same results. The outside sites are close to the temperatures of the baths, while the inner sites achieve a constant intermediary temperature. The situation for the chains with Jxy = 4 and bz = 2 is very different, however. Since a large Jxy favors a singlet eigenstate and the overall temperatures are not that large, one expects a rather small spin expectation value in these cases, and therefore significant differences between the tanh and the ( c i l Z ) r fits. This is certainly the case. Somewhat more surprising is the significant difference between the predictions from fits to ( ° I , Z ) T a n d to {hijT, observed for the smaller chains. The agreement between the two improves as the chain becomes longer, and for the four-site chain it is quite reasonable. Because of the relatively small bz (compared to the other energy scales), one expects that here the values obtained from (hi)? may be more trustworthy, although for the three-spin chain even these are rather pe-culiar. We ascribe the non-monotonic profile for the 3-spin chain to the large temperature differential. We have checked that if we lower TL to 3.9 instead of 4, the spin temperatures are monotonically decreasing, and they remain so for all smaller temperature differentials (see next figure). More data for three and four-site chains with Jxy = 4 and bz = 2 is shown in Fig. 3.12. The right-bath temperature is kept fixed at TR = 3, while TL is varied between 3 and 4. The results from the tanh-fit (green symbols) are clearly wrong, since even when TL = TR they predict different temperatures for Tj. The values extracted from fits to ( o i ) Z ) r and to (/IJ)T both predict linear increases of the local temperatures of the entire chain with TL- The increase is faster for sites closer to the hotter bath. For the three-site chain, the temperature of the middle site 2 is in good agreement for both methods, however, (hi)T predicts lower/higher temperature for T1/T3 than those obtained from (<7I IZ)T-The same trend is true for the 4-site chain. Both fits to (hi)T and to (CJ ) Z)T predict that T i and T2 are rather close to one another, as are T3 and T4. The whole chain heats as TL increases, but again the temperature of different sites increases at different rates. However, in all cases there is a monotonic decrease of the temperatures across the chain, TL > T\ > T2 > T 3 > T 4 > TR. As alluded to before, this situation may change if there is a very large temperature differential across the chain. • Finally, in 3.13 we show temperature profiles for three and four-site chains with Jxy = 2 and bz = 0.1, as a function of TL when TR is fixed at 3. The magnetic field is very small, and this combined with the tendency of the system towards antiferromagnetic ordering leads to even smaller (o~iiZ)T values. As dis-cussed, in such cases some of the (o~iiZ)T expectation values are non-monotonic functions of T, leading to multiple solutions for Tj. In this case, we plot the value of Tj closest to the temperature of the bath near which the spin is located. Even so, fits from this are very peculiar for the three-site chain, as they suggest that T2 and T3 cool below TR, and that for large TL, the middle site is coldest than T i . This must be wrong, since it is impossible to have the chain colder than the baths it is connected to. For the four site chain, the fits to ( a i , z ) r Chapter 3. Spin chains coupled to thermal baths 47 • • m • o • CT • I A 1 • i i A 1 « * CT T CT • • • * • M * • : " i i * t l • • t • • • A f • ... • • .. m • Figure 3.12: Steady-state temperature profile for 3-site (upper panel) and 4-site (lower panel) chains as a function of TL, for TR = 3, JXY = 4 and bz = 2. Estimates from the tanh fit (green symbols) are completely wrong. Estimates from fits of (/I ,)T (blue symbols) and (cri<z)T (red symbols) are in qualitative but not quantitative agreement. Chapter 3. Spin chains coupled to thermal baths 48 • I t Figure 3.13: Steady-state temperature profile for 3-site (upper panel) and 4-site (lower panel) chains as a function of TL, for TR — 3, Jxy = 2 and bz = 0.1. Estimates from the tanh fit (green symbols) are completely wrong, and estimates from (o"j)Z)x (red symbols) are also rather peculiar, especially for the three-site chain. Estimates from fits of (hi)r (blue symbols) show behavior similar to that observed in the previous case. Chapter 3. Spin chains coupled to thermal baths 49 predict T i > T2 > T 4 > T 3 , but at least they are now all in between Tr, and TR. However, the values obtained from ( / IJ ) T are still reasonable and in agreement with previous data, and also with the results for the four-site chains. Clearly more investigations are needed before a coherent picture of the heat transport in spin chains emerges. First, we still need to understand how reason-able are the definitions of the local temperatures in cases where there is signif-icant difference between their predictions. It is to be hoped that using a small "thermometer" to measure local temperatures wil l help solve this issue. The dependence on chain length and on the various parameters Jz, Jxy, bz, TR, TL must be investigated thoroughly. For Jz = 0 this task may be aided by map-ping to non-interacting fermions and the possibility to use known results from there. Better modeling of the coupling to the leads (so that r is not the same irrespective of the energy of the transition) is also of interest. Finally, different types of baths may be analyzed. We already have results for "fermionic" baths, for the same type of quantities described here. These results are obtained by simply replacing n(E) —> (e / 3 ( £ - ' J ) + l ) - 1 , in the occupation numbers. Such baths could be thought of as collections of independent two-level systems (for instance, modeling coupling to nuclear spins) that are biased by their own local field, fixed by p. The results we have (not shown) reveal some similarities but also clear differences from the results for the phonon baths. These shall also be further investigated in future work. 50 Chapter 4 Conclusion In this thesis we have made some initial steps in the effort of developing a formalism for studying transport properties in systems biased by coupling to several baths. While the approach we propose if quite general, we have chosen short spin-chains as specific models in this work. Our approach is to generalize the projection operator technique to cases with multiple baths. The key approximation made is that the coupling to the baths is weak, so that second-order perturbation is sufficient to describe it. We have also made a Markov approximation to simplify the calculation. This approximation should be reasonable for the steady-state case that we focus on, but its effects on short term evolution remain to be analyzed. This requires numerical integration of some differential equations and it therefore seems feasible, at least for not too large systems. If only the steady-state is of interest, the calculation simplifies even more. In this case there is no need to solve the entire time evolution, instead the eigenstate corresponding to the zero eigenvalue of a certain matrix must be found. Using Maple, we could investigate systems with up to 4 spins, however we estimate that up to 14 or 15 spin chains could be studied quite straightforwardly if we write our own programs. We have shown that in a non-equilibrium steady-state, the reduced density matrix is non-diagonal, as one would expect since otherwise there would never be any currents flowing through the system. This observation raises some doubts about the relevance of studying Quantum Master Equations to learn something about non-equilibrium behavior of quantum systems, for cases where the steady-state is not an equilibrium state. One particular problem in the investigation of heat properties, which are our focus here, is the definition of the local temperature when the system is in a non-equilibrium steady state. For reasons explained in the text, we disagree with the simple definition used in the literature. We have introduced two other methods to calculate the local temperatures. In some cases, the predictions of the two methods agree with one another, and one can be reasonably certain that they are correct. When the two definitions disagree (as seen for small magnetic fields), the results obtained from the fit of the local energy to its equilibrium value are still reasonable, however it would be re-assuring to have a third method to measure the local temperatures that would give similar results. Such work is in progress, based on the idea of using a small "thermometer" coupled to individual sites to measure their temperatures. Of course, all this assumes that the notion of "local temperature" is meaningful. This also remains to be tested, but one expects it to be so at least in cases where the bias across the system is Chapter 4. Conclusion 51 not too large. System-size and available time constraints have severely limited the number of results we could show in this thesis. The results we show suggest that inter-esting physics may be uncovered for spin chains, but clearly significantly more effort is needed before a coherent picture emerges. For xy chains, in particular, one may also use mapping to non-interacting spinless fermions to try to make more progress, although the meaning of the coupling to the baths after such a transformation may be rather peculiar. In general, for all the results we obtain, we plan to also attempt the calculation of similar quantities using other appro-priate techniques, such as linear response for small biases, or non-equilibrium Keldysh formalism for large biases. Such comparisons will also be undertaken, soon. They should help define the limit of validity of this approach. Longer term, we propose to apply this formalism to understand charge trans-port through small systems, like collections of quantum dots, which may or may not also contain spins and/or bosonic degrees of freedom. The class of systems of interest is very large, and clear results and formalisms are currently very few. It is to be hoped that any systematic approach to such problems, as the one we are proposing here, will lead to some clear progress in understanding the non-equilibrium physics in such systems, and allow for a detailed interpretation of experimental results. 52 Bibliography M. Toda R. Kubo and N. Hashitsume. Statistical Physics II, Nonequilib-rium Statistical Mechanics. Springer-Verlag, Berlin Heidelberg New York, 1995. R. Landauer. Electrical resistance of disordered one-dimensional lattices. Phil. Mag., 21, 1970. R. Landauer M . Biittiker, Y . Imry and S. Pinhas. Generalized many-channel conductance formula with application to small rings. Phys. Rev. B. 31, 1985. C. Caroli, R. Combrescot, P. Nozieres, and D. Saint-James. Direct cal-culation of the tunneling current. J. Phys. C: Solid St. Phys, 4:916-929, 1971. G. D. Mahan. Many-Particle Physics. Plenum Press, New York, 1990. D. S. Fisher and P.A. Lee. Relation between conductivity and transmission matrix. Phys. Rev. B, 23, 1981. R. L iv i S. Lepri and A. Politi . Thermal conduction in classical low-dimensional lattices. Physics Report, 377, 2003. M. Michel, M. Hartmann, J . Gemmer, and G. Mahler. Fouriers law con-firmed for a class of small quantum systems. Euro. Phys. J. B, 34:325-330, 2003. M. Michel, J . Gemmer, and G. Mahler. Heat conductivity in small quantum systems: Kubo formula in liouville space. Euro. Phys. J. B, 42:555-559, 2004. J .L . Lebowitz F. Boetto and J . Lukkarinen. Fourier's law for a harmonic crystal with self-consistent stochastic reserviors. J. Stat. Phys., 116, 2003. G.Mahler and V .A . Weberrufi. Quantum Networks: Dynamics of Open Nanastructures. Springer-Verlag, Berlin Heidelberg New York, 1998. G. Lindblad. On the generators of quantum dynamical semigroups. Comm. Math. Phys., 48:119-130, 1976. V. Gorini, A. Kossakowski, and E .C.G. Sudarshan. Completely positive dynamical semigroups of n-level systems. J. Math. Phys., 17:821-825, 1976. Bibliography 53 [14] Hu. Gang. Lecture notes for Non-equilibrium Statistical Physics. Physics Department in Beijing Normal University, unpublished, 2001. [15] P. Jordan and E. Wigner. Uber das paulische Aquivalenzverbot. Z. Physik, 47, 1928. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items