Quantum Transport Through Open Systems by Jinshan Wu B.Sc., The Beijing Normal University, 1999 M.Sc., The Beijing Normal University, 2002 M.Sc., The University Of British Columbia, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Physics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April, 2011 c
Jinshan Wu 2011 Abstract To study transport properties, one needs to investigate the system of inter- est when coupled to biased external baths. This requires solving a master equation for this open quantum system. Obtaining this solution is very chal- lenging, especially for large systems. This limits applications of the theories of open quantum systems, especially insofar as studies of transport in large quantum systems, of interest in condensed matter, is concerned. In this thesis, I propose three ecient methods to solve the Redeld equation | an example of such a master equation. The rst is an open- system Kubo formula, valid in the limit of weak bias. The second is a solu- tion in terms of Green's functions, based on a BBGKY (Bogoliubov{Born{ Green{Kirkwood{Yvon)-like hierarchy. In the third, the Redeld equation is mapped to a generalized Fokker-Planck equation using the coherent-state representation. All three methods, but especially the latter two, have much better eciency than direct methods such as numerical integration of the Redeld equation via the Runge-Kutta method. For a central system with a d-dimensional Hilbert space, the direct methods have complexity of d3, while that of the latter two methods is on the order of polynomials of log d. The rst method, besides converting the task of solving the Redeld equation to solving the much easier Schrodinger's equation, also provides an even more important conceptual lesson: the standard Kubo formula is not applicable to open systems. Besides these general methodologies, I also investigate transport proper- ties of spin systems using the framework of the Redeld equation and with direct methods. Normal energy and spin transport is found for integrable but interacting systems. This con
icts with the well-known conjecture linking anomalous conductivity to integrability, and it also contradicts the relation- ship, suggested by some, between gapped systems (Jz > Jxy) and normal spin conductivity. I propose a new conjecture, linking anomalous transport to the existence of a mapping of the problem to one for non-interacting particles. ii Preface The work presented in Chapters 2 and 6 are a continuation of the work done for my M.Sc. thesis. This work is reported in the manuscript: \Heat transport in quantum spin chains: the relevance of integrability", Jinshan Wu and Mona Berciu, submitted to Physical Review B and posted as arXiv:1003.1559. The work presented in Chapter 3 is reported in Jinshan Wu and Mona Berciu 2010 Europhysics Letters 92 30003. The rst part of the work presented in Chapter 4 is reported in Jinshan Wu 2010 New J. Phys. 12 083042. The second part, which is on the second- order BBGKY-like method, is currently under preparation for publication. Original work reported in other chapters will be prepared for publication in the near future. All the original work in this thesis was carried out by JW who also developed the conception, scope and methods of this research, with various degrees of consultation and editing support from MB. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Two major theoretical approaches . . . . . . . . . . . . . . . 2 1.2 A big challenge for the system-baths scenario . . . . . . . . . 5 1.3 Our methods and their applications . . . . . . . . . . . . . . 7 1.3.1 Linear response theory . . . . . . . . . . . . . . . . . 7 1.3.2 The BBGKY-like hierarchy . . . . . . . . . . . . . . . 8 1.3.3 The GFPE in the coherent-state representation . . . 9 1.3.4 Thermal conductance of spin chains . . . . . . . . . . 9 1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Chapter 2 The projector technique and the Redeld equation 12 2.1 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 The eective equation of motion for open systems . . . . . . 12 2.2.1 The projector operator technique . . . . . . . . . . . 14 2.2.2 An alternative derivation . . . . . . . . . . . . . . . . 19 2.2.3 QME in terms of creation and annihilation operators 20 2.3 Relaxation towards thermal equilibrium . . . . . . . . . . . . 22 2.3.1 One-site open system . . . . . . . . . . . . . . . . . . 22 2.3.2 Two-site open system . . . . . . . . . . . . . . . . . . 24 2.4 Various other extensions . . . . . . . . . . . . . . . . . . . . 28 iv Table of Contents 2.5 Evolution towards a NESS . . . . . . . . . . . . . . . . . . . 30 2.5.1 NESSs from the RE . . . . . . . . . . . . . . . . . . . 31 2.5.2 NESSs from the LOLE and the MQME . . . . . . . . 34 2.5.3 The general Lindblad form . . . . . . . . . . . . . . . 37 2.6 Major challenge for large systems . . . . . . . . . . . . . . . 38 Chapter 3 Linear Response Theory: Kubo formula for open systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 General Linear Response Theory . . . . . . . . . . . . . . . . 41 3.3 Limitations of the standard Kubo formula . . . . . . . . . . 43 3.4 OKF: the Kubo formula for open nite-size systems . . . . . 47 3.5 Comparison of the two OKFs . . . . . . . . . . . . . . . . . . 50 3.6 Summary and discussion . . . . . . . . . . . . . . . . . . . . 53 Chapter 4 BBGKY-like hierarchy for the Redeld equation 55 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Derivation of the BBGKY-like hierarchy . . . . . . . . . . . 56 4.3 Non-equilibrium Wick's theorem . . . . . . . . . . . . . . . . 59 4.4 Truncation of the hierarchy . . . . . . . . . . . . . . . . . . . 62 4.4.1 Method 1: equilibrium substitution . . . . . . . . . . 62 4.4.2 Method 2: rst-level cluster expansion . . . . . . . . 67 4.4.3 Method 2: second-level cluster expansion . . . . . . . 72 4.5 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . 78 Chapter 5 Solving the Redeld equation using coherent-state representations . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2 The model and its eective equation of motion . . . . . . . . 81 5.3 Exact solution for non-interacting systems . . . . . . . . . . 85 5.4 Interacting systems: BBGKY-like hierarchy derived from the GFPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.5 Approximate numerical solution for interacting systems . . . 90 5.6 Solving the LOLE via the GFPE method . . . . . . . . . . . 94 5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Chapter 6 Application: heat transport in spin chains . . . . 97 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2 The model and its Redeld equation . . . . . . . . . . . . . . 99 6.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . 100 v Table of Contents 6.4 Denitions of the thermal current and local temperatures . . 101 6.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Chapter 7 Conclusions and discussions . . . . . . . . . . . . . 112 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Appendix A Perturbational decomposition of the operators m̂ 128 A.1 Expanded in eigenbasis of HS . . . . . . . . . . . . . . . . . 128 A.2 Perturbational expansion starting from non-interacting sys- tems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Appendix B Second-order cluster expansion: G1, G2 . . . . . 134 B.1 Second-order cluster expansion . . . . . . . . . . . . . . . . . 134 B.2 Second-order equations of G1 and G2 . . . . . . . . . . . . . 135 Appendix C Estimation of convergence of the BBGKY-like methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 C.1 E;(1) 1 from method 1 . . . . . . . . . . . . . . . . . . . . . . 141 C.2 C;(1) 1 from method 2 . . . . . . . . . . . . . . . . . . . . . . 145 Appendix D Coherent states and coherent-state representa- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 D.1 Coherent states for systems with a single mode . . . . . . . . 148 D.2 Coherent states for multi-mode systems . . . . . . . . . . . . 151 Appendix E FPE and Langevin equation . . . . . . . . . . . . 153 E.1 Langevin equation for standard FPEs . . . . . . . . . . . . . 153 E.2 Analytical solution for the Ornstein-Uhlenbeck process . . . 155 E.3 Stochastic dierence equation for truncated GFPEs . . . . . 155 Appendix F Solving the RE and LOLE of bosonic systems with BBGKY-like hierarchy . . . . . . . . . . . . . . . . . . . 157 F.1 Exact solution on non-interacting systems . . . . . . . . . . . 158 F.2 Approximate solution on interacting systems . . . . . . . . . 159 vi List of Figures 1.1 Sketch of a typical setup of heat conduction. . . . . . . . . . 10 2.1 Sketch of a one-site system coupled to a heat bath. . . . . . . 23 2.2 Sketch of a two-site system coupled to a heat bath. . . . . . . 24 2.3 Sketch of a two-site system coupled to two heat baths. . . . . 27 2.4 Dierence between non-equilibrium and equilibrium states . . 33 2.5 Stationary states from the local-operator Lindblad equation and the multi-mode quantummaster equation compared against those from the Redeld equation . . . . . . . . . . . . . . . . 36 3.1 Dierence between (1) (1) and (2) (1) for dierent values of . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2 Dierence between (1) (1) and (2) (1) for dierent values of . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3 Dierence between (1) (1) and (2) (1) for dierent values of V0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1 g E;(1) 1 compared with g Ex 1 for interacting systems at non- equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 g C;(1) 1 compared with g Ex 1 for non-equilibrium interacting sys- tems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3 g C;(2) 1 compared with g Ex 1 for non-equilibrium interacting sys- tems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4 g C;(2) 2 is compared with g Ex 2 for non-equilibrium interacting systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.5 Particle currents calculated from the second-, rst- and the zeroth-order forms are compared . . . . . . . . . . . . . . . . 76 4.6 Heat currents calculated from the second-, rst- and the zeroth- order forms are compared . . . . . . . . . . . . . . . . . . . . 77 5.1 Analytical solution compared against the BBGKY numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 vii List of Figures 5.2 The numerical simulation of the generalized Fokker-Planck equation for interacting systems . . . . . . . . . . . . . . . . . 93 5.3 A similar application to the local-operator Lindblad equation 96 6.1 Typical temperature proles for N = 10 spin chains . . . . . 104 6.2 Current v.s. system size . . . . . . . . . . . . . . . . . . . . . 105 6.3 Normalized proles . . . . . . . . . . . . . . . . . . . . . . . . 106 6.4 Typical temperature prole for the cases of small Jz and small Bz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.5 Temperature proles on Ising spin chains and fermionic chains 109 7.1 A sketch of an exactly solvable decoherence model and its variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 viii Acknowledgements It is a pleasure to thank those who made this thesis possible. I owe my deepest gratitude to my supervisor, Dr. Mona Berciu. This thesis would not have been possible without her great eorts to educate and to help me in so many ways. She taught me not only physics, but also more importantly how to appreciate physics, and how to write up and explain things clearly and simply. She also encouraged and protected my interest in fundamental physical questions. Throughout my years of graduate study at UBC, she has always been a good advisor and a good friend to me. It is an honor for me to thank the many people with whom I discussed the related work and have been beneted greatly from such discussion: Ian Aeck, Keiji Saito, Philip Stamp, Tomaz Prosen, Peter Drummond, Hong Guo and Shouyong Pei. I would like to thank the many great teachers I ever had especially: Yue'e Yi, Jianjiang Zhang, Anshen Qi, Jiaruan Wang, Yongcheng Wang, Zhanru Yang, Shouyong Pei, Canbing Liang, Jingdong Bao, Zengru Di, Qiang Yuan, Andrew DeBenedictis, Leslie Ballentine, Ariel Zhitnitsky, Mark Van Raamsdonk and Ian Aeck. I am indebted to my many colleagues for the stimulating discussion and sometimes even excited and heated argument: Yong Sun, Mandy Wang, Kelly Cheung, Bayo Lau, Arthur James Charbonneau, Jennell Van Dongen, Kory Stevens, Si Chen, Junliang Song, Menghui Li, Liang Gao, Ying Fan, Dahui Wang, Limei Xu, Yougui Wang, Zhangang Han, Honggang Li. Lastly but most importantly, I wish to thank my wife, Qian Feng, for her unconditional love and support; and my daughters, Lindsay Wu and Elysia Wu for the happiness that they have already brought to this family. To them I dedicate this thesis. ix To Qian Feng, Lindsay Wu and Elysia Wu x Chapter 1 Introduction Theories of open quantum systems [1{3] usually start from an eective equa- tion of motion of the central systems of interest, which is usually derived from the exact equation of motion for the composite system | the cen- tral system plus the bath(s). Although various approximations are involved in deriving such eective equations of motion, which are sometimes called quantum master equations, they are believed to be reasonably reliable since they are built on relatively rm grounds, i.e. an equation of motion derived from rst-principles. Theories of transport [4, 5] in solid state physics deal with a specic type of open system, where the central system is coupled to multiple baths held at dierent temperatures, chemical potentials, etc. As a result, study of transport should be a branch of the theories of open quantum systems. In practice, dierent theories based on other approaches have been developed for transport studies, such as the Kubo formula [6, 7], the Landauer-Buttiker formula [8{11] and the non-equilibrium Green's function method [4, 12{16]. One important reason why theories of open quantum systems are not widely used when studying transport is the computational diculty in applying them to large systems. In studies of relaxation times in nuclear magnetic resonance [17], where the theories of open quantum systems provide the standard tools, a typical central system has usually only a few degrees of freedom. In contrast, in discussions of transport, one studies quite often mesoscopic systems whose Hilbert space's dimension easily exceeds 220. It is the primary goal of this thesis to nd new ways of making the theo- ries of open quantum systems amenable for the study of transport in large(r) systems, by nding ecient calculation schemes. For non-interacting sys- tems there are straightforward methods to calculate stationary states and therefore physical quantities, even if the system is large. Thus, the real challenge is for systems with interactions, just as is the case in the study of equilibrium properties and dynamical processes [18]. We will propose various ecient and accurate approximations, and test them. This chapter serves as an introduction to the subject and also an overview of this thesis. Sections x1.1 and x1.2 brie
y review the theories of open sys- 1 1.1. Two major theoretical approaches tems while focusing on their applications to the study of transport. Section x1.3 lists our contributions, including the development of new methodolo- gies and discovery of new physics from application of these methodologies to certain specic systems. The rest of this thesis follows the same line of logic as this overview, but in more details. Chapter 2 reviews the projector technique used in theo- ries of open quantum systems. Examples are given to illustrate the basic procedure. The resulting eective equations of motion are solved via direct methods, which have very poor eciency. Then we present an approximate solution based on linear response theory in Chapter 3. This approxima- tion has better eciency than the direct methods but not good enough to be comparable to the non-equilibrium Green's function method. Chapters 4 and 5 introduce two very ecient and accurate approximations. One is based on a Bogoliubov{Born{Green{Kirkwood{Yvon (BBGKY)-like equa- tion hierarchy and the other is based on the coherent-state representation of density matrices. In Chapter 6 we apply the direct methods to discuss thermal transport of spin systems and the validity of Fourier's Law. At last, in Chapter 7 we summarize these results and list some of the potentially valuable questions to be investigated in future, with our newly developed methods. We also comment on the value and implications of the work presented in this thesis. 1.1 Two major theoretical approaches A rough picture of a typical experimental setup for measuring transport properties is that of a quasi 1-d system coupled to multiple biased baths. We want to calculate currents
owing through the system given the parameters of the baths such as their temperatures, chemical potentials, etc. Such calculations are relevant for comparison against experimental measurements of charge, energy and spin transport, for example through quasi 1-d spin systems [19], quantum wires including nanotubes [20], molecular devices [21{ 24] and quantum dots [25, 26]. In particular, in recent years there have been many experiments looking at heat transport of spin systems. The heat conductance = JhT , where Jh is the heat current and T the applied temperature bias, is measured for a crystal with spin-chain structures. Typical crystals sizes are of roughly a few mm. In order to measure the heat conductance contributed by the spin chains, one has to subtract the conductance due to phonons and charge carriers. Usually this is achieved by using anisotropic spin-chain compounds 2 1.1. Two major theoretical approaches such as (Ca,Sr)14Cu24O41 [27], which contain S = 1 2 two-leg ladders with an intra-chain coupling much stronger than the inter-chain coupling. By measuring the dierence between the conductance parallel to the chain and that perpendicular to the chain, one nds the pure spin conductance. The idea is that the parallel conductance is due to both the spin-chain and other sources, while the perpendicular one only has the contributions from other sources. Experimental results on spin chains in strong magnetic elds have also been reported [28, 29]. However, at present, for spin chains, there is no available data regarding temperature proles and nite size-scaling for the heat currents. This makes comparison with our theoretical results, discussed in Chapter 6, impossible at this time. In a recent measurement of thermal conductance of nanotubes [30], Wang et al. reported the size dependence of heat currents and concluded that it violates Fourier's law. For small, and especially also quantum, systems, it is not reliable to use phenomenological laws such as Fourier's law and Ohm's law to characterize their transport. To understand and model such systems, a theoretical frame- work based on rst principle is necessary. Developing such a framework is not just a matter of theoretical interest, but also has obvious practical value. Various theoretical models have been used so far. Some studies are based on the phenomenological balance equation [31] or rate equation [32, 33], but the most common and systematic approaches are the Kubo for- mula [6, 34, 35], the non-equilibrium Green's function method [13{15], and the open quantum system approach [1{3] with the system-bath scenario [36{ 38]. In the approach based on the standard Kubo formula, the system itself is assumed to be innite and homogeneous. There is no explicit use of the baths but instead the system is treated as being isolated. The current is calculated as a response to a presumed applied electrical potential [6] or temperature gradient drop [39], or uneven mass gradient distribution [40]. In the approach of non-equilibrium Green's function and Landauer formula, the system is modeled by three pieces: the left lead, the central system and the right lead. The leads are taken as half-innite non-interacting systems and the central system is usually nite. There are no explicit baths in- cluded in this picture. The parameters of the baths, such as temperatures and chemical potentials, are attached to the corresponding leads. Finally, the theories of open quantum systems deal with a nite-size central system coupled explicitly to baths, whose degrees of freedom are then integrated out. There are no innite leads in this picture. These three dierent physical pictures are supposed to describe the same problem. Fisher and Lee [41] showed that for non-interacting systems the 3 1.1. Two major theoretical approaches rst two are somewhat equivalent if the same three-piece system is used when applying the Kubo formula. We should note, however, that the ma- jority of studies of conductance via the Kubo formula are using an innite homogeneous system [35, 42], or extrapolate from results on nite-size sys- tems [43, 44]. With that being said, let us ignore the question of whether these methods are equivalent, and focus only on the latter two approaches. Comparisons between them, or even settling how to make them directly comparable, has not been extensively studied. Non-equilibrium Green's functions is a perturbational theory starting from a well-dened piece-wise equilibrium state of the disconnected three- piece system. The left/right lead is assumed to be in its own equilibrium state correspondingly dened by the left/right temperature and chemical potential, while the central system is in an arbitrary equilibrium state [14]. From this, the Green's function of the disconnected system is calculated and later used to construct the full Green's function of the connected sys- tem, with the connection between the central system and the leads treated as perturbation. The standard diagrammatic perturbation theory does not apply here since the initial and nal states are qualitatively dierent. It is an essential assumption of the standard diagrammatic perturbation theory that the states at both 1 in time are the same equilibrium state (or the ground state if T = 0) of the unperturbed system. At this point a non- equilibrium Keldysh formalism [12, 18] is used instead. Conceptually, the condition that the two leads are half-innite is essential. With any nite- size leads, an equilibrium state { instead of the desired non-equilibrium stationary state { will be reached ultimately since the unbalanced energies or chemical potentials will be eventually dissipated over the whole system. However, even for innite leads, theoretically there is still no guarantee that the proper stationary state can be calculated perturbationally [13]. In prac- tice, for interacting systems, the non-equilibrium Green's function method is often used together with the density functional theory [13, 16]. Addi- tional approximations, thus possibly additional errors [45, 46], are involved in dealing with the interactions. Theories of open quantum systems, on the other hand, are based on relatively rmer grounds. They start from the rst-principle equation of motion of the composite system, and arrive at an eective equation of mo- tion of the central system. Various equations can be derived depending on the level of approximations. Using the in
uence functional theory [47, 48], exact time-dependent master equations have been derived for simple cen- tral systems [49]. Examples include a two-level system, one harmonic os- cillator and two coupled harmonic oscillators, which is the current limit of 4 1.2. A big challenge for the system-baths scenario such exact approaches [50]. When the Markovian approximation is used, time-independent master equations have also been derived from the exact ones [51]. For more general central systems, some direct brute-force gener- alization of the master equation of a simple two-site system have also been used in the literature. We will call one kind of such equations the local- operator Lindblad equation for reasons that will be clear later. It is in fact very hard to derive exact master equations from in
uence functional theory for general central systems. Instead, approximate equations can be derived from the projector technique [51] or other similar techniques. Such a master equation is sometimes called the Redeld equation [52{54]. Applications of the Redeld equation to transport studies are in an early but fast develop- ing stage [36{38, 55{57]. Others prefer to use the local-operator Lindblad equation [58, 59], which turns out to be relatively easier to solve numerically. It is usually easier to deal with large systems using non-equilibrium Green's functions, than it is with theories of open quantum systems. In practice, quite often a set of few-particle Green's functions are enough for calculating the physical quantities of interest. The number of single-particle Green's functions, for instance, is much smaller than the number of ma- trix elements of the density matrix, which appears in the Redeld equation. Consider for example an N -site 1-d spinless fermionic system. The set of all single-particle Green's functions has N2 unknowns while a density matrix has 22N unknowns. Therefore, unless an ecient calculation scheme whose complexity scales as N2 can be found, numerically the two approaches are not comparable at all. This explains why the non-equilibrium Green's func- tion method is much more widely used than other methods to study trans- port in solid state systems. To many, it does not even seem necessary and urgent to compare results from the two methods and determine which is more accurate or even which one is correct. This situation may change if ecient methods to solve the Redeld equa- tion are found, as is the goal of this thesis. 1.2 A big challenge for the system-baths scenario Both the general Redeld equation and the general local-operator Lindblad equation can be cast into the form of a generalized Liouville equation, d dt (t) = L (t) ; (1.1) where (t) is the density matrix of the central system and L is a linear superoperator acting on the density matrix. We are only interested in the 5 1.2. A big challenge for the system-baths scenario long-term stationary state and the steady current associated with it. The stationary state is dened by L (1) = 0; (1.2) where (1) is the long-term non-equilibrium stationary state. If the Hilbert space of the central system is d-dimensional, then the density matrix has d2 elements. If we represent the density matrix as a vector, for example, P = [11; 12; ; dd]T ; (1.3) then the linear superoperator L is a d2d2 matrix. The task of nding (1) thus implies nding the zeroth mode of L | the eigenvector corresponding to the zero eigenvalue. One can imagine that this becomes nontrivial when d > 210, as in the example of a fermionic chain with N > 10 mentioned earlier. The central task of this thesis is to nd ecient methods to solve this equation. Before we present our work, we brie
y review the currently avail- able methods. Instead of solving directly for the zero mode of L, we may instead let the system evolve from an arbitrary initial state long enough to arrive at (1). Thus, the rst class of solutions are based on propagators { this includes the direct Runge-Kutta method [37, 60], the Newtonian polynomial propagator [61], the short-iterative-Arnoldi propagator [62] and the short- time Chebyshev polynomial propagator [63]. All these methods have been compared in Ref. [64]. As for their eciency, they are all roughly on the same scale { their complexities scale roughly as d3. A more ecient approach is the stochastic wave-function method [2, 65, 66]. The basic idea is to replace the eective equation of motion by a stochastic Schrodinger equation for wave functions. The method has a complexity scaling as d2. On the local-operator Lindblad equation it has been shown that the stochastic wave-function method leads to accurate re- sults [67]. Breuer et al. showed that it can also be applied to the Redeld equation [65], but its eciency is still not comparable with that of the non- equilibrium Green's functions method. For the local-operator Lindblad equation, in fact, there is another quite ecient method developed recently by Prosen et al. [59]. It is based on the time-dependent density matrix renormalization technique. This method makes it possible to analyze systems roughly of the same size as those that can be studied with the non-equilibrium Green's functions method. It is still 6 1.3. Our methods and their applications an open problem to extend this approach to the Redeld equation. How- ever, we should note that the same central system might behave dierently under evolution described by the local-operator Lindblad equation versus the Redeld equation. In fact, this issue has not been settled because of the lack of comparison between the results of the two approaches. The lack of ecient methods for the Redeld equation has certainly contributed to this state of aairs. With our newly developed methods, such a comparison can nally be undertaken. 1.3 Our methods and their applications In this section, we roughly summarize the ideas, performances and applica- tions of our newly developed methods. Their details will be discussed in the rest of this thesis. We have developed three dierent methods, with dierent strenghts and weaknesses. 1.3.1 Linear response theory The rst method is a linear response theory based on the Redeld equation. The basic idea is very simple. The equilibrium state of the central system (when unbiased) is known to be the Boltzmann distribution. The stationary state of the Redeld equation under appropriate conditions, i.e. all baths at the same temperature and chemical potential, is indeed the Boltzmann distribution. This makes it possible to separate the operator L into a large part L0 and a small part L, where the former part is the operator describ- ing the equilibrium conditions. Then, the non-equilibrium distribution can be treated as a perturbational response to the small part L. This method has eciency comparable to the stochastic wave-function method. Besides the gain in eciency, which is not that remarkable, an impor- tant conceptual conclusion is drawn from this work. The standard Kubo formula [6], as we already knew, does produce a proper rst-order correc- tion to equilibrium states when an additional potential is included into the Hamiltonian. However, our work shows that it fails to describe the non- equilibrium stationary states. If innite-size systems are used as the cen- tral system, there are ways to get around this conceptual problem. For non-equilibrium stationary states in nite-size systems, the standard Kubo formula must be replaced by a Kubo-like formula based on the Redeld equation { the standard Kubo formula is no longer valid. The key dierence is that the standard Kubo formula uses Schrodinger's equation and the sys- tem is treated as isolated while in the new Kubo-like formula, the Redeld 7 1.3. Our methods and their applications equation of an open system is used. Given that there are still many studies of transport based on various forms of the standard Kubo formula, this is a very important conclusion. 1.3.2 The BBGKY-like hierarchy The second method mirrors the ideas of the well-known BBGKY hierarchy. The original BBGKY hierarchy [68] is a set of coupled equations for equilib- rium correlation functions (Green's functions). For non-interacting systems, the hierarchy is decoupled, meaning that the equation of a single-particle Green's function only depends on other single-particle Green's functions, leading to a closed system of equations. Using Wick's theorem, which holds for correlation functions of non-interacting systems, higher order n-particle Green's functions can then be constructed from the single-particle Green's functions. However, when there are interactions in the system, equations for all orders of Green's functions are coupled together. Fortunately, other ap- proximations can be used to truncate the hierarchy and then one solves the truncated hierarchy. These approximations can be performed in such a way that they are equivalent to the diagrammatic perturbation theory [18, 69]. Here, we extend the same approach to non-equilibrium stationary states. A BBGKY-like hierarchy is derived from the Redeld equation. We nd again that for non-interacting systems, the equations decouple and a similar Wick's theorem can be proved. All equations are coupled together for inter- acting systems. Two further approximations are introduced to truncate the hierarchy. The rst one is based on the known equilibrium distribution of the interacting system. The second one is based on the cluster expansion [69]: the observation that higher-order Green's functions can be constructed from the lower-order ones through Wick's Theorem and an additional correlation term. Ignoring these correlation terms at a certain order leads to a trunca- tion of the hierarchy at the corresponding order. Our tests of the two approximations for small systems, where exact solu- tions needed for comparisons are also possible, indicate that they have very good accuracy. Moreover, they can be systematically improved by going to higher orders in the truncation schemes. The complexity of this approach scales only like polynomials of N , not d, therefore this is a huge improvement of the eciency. 8 1.3. Our methods and their applications 1.3.3 The general Fokker-Planck equation in the coherent-state representation The third method we propose is to map and solve the Redeld equation in the coherent-state representation of quantum states [70, 71]. Consider a simple harmonic oscillator for example. The set of all coherent states ji, which are eigenvectors of the annihilation operator a ji = ji, provides an overcomplete basis for quantum states. Expanding a density matrix in such a basis results in a \distribution" function P () over the complex plane C. Here we use quotation marks because of the possibility that P () could be a negative or even complex number. In this way, creation and annihilation operators become dierential operators, and the Redeld equation becomes a stochastic dierential equation of the distribution function over CN for a general coupled interacting system of N harmonic oscillators. Instead of dealing with the huge size of the linear system, like with the direct meth- ods, now we need to solve an N -dimension stochastic dierential equation, which generally can be simulated classically. For non-interacting systems, in fact, the equation can be solved analytically. We have obtained analytical expressions of non-equilibrium stationary states and compared them to the exact solutions. We nd perfect agreement, as expected. For interacting systems, once the classical simulations are performed, the same comparisons can be done. The number of variables increases linearly with the system size so this method is capable of looking at very large systems. A similar procedure [72{74] has been proposed and implemented for the study of equilibrium and pure dynamical systems. There the method is applied successfully for bosonic systems with roughly N = 106 modes [75]. In a way, our method can be regarded as an extension of that method to the study of non-equilibrium stationary states. Although this extension involves additional technical diculties, its eciency remains comparable to that of its counterpart for the study of equilibrium states and pure dynamical processes. 1.3.4 Thermal conductance of spin chains The physical question we discuss using this open-system approach is the validity of Fourier's Law of heat conduction, jh = rT; (1.4) where jh is the heat current density, T is the local temperature and is the heat conductivity. This phenomenological law is found to be obeyed by 9 1.3. Our methods and their applications jh TL TR Figure 1.1: Sketch of a typical setup of heat conduction: left and right baths have dierent temperatures and the conductor is a quasi 1D system. jh is the current density on the conductor. It is usually assumed that there is a linearly decreasing temperature prole on the conductor. macroscopic systems, when measured experimentally in set-ups such as the one sketched in Fig. 1.1. The derivation of the law from the microscopic equation has been an open question for many decades [34, 76]. The question has been investi- gated in both classical [76] and quantum systems [34]. For classical systems, a plausible relation between chaotic behavior and the validity of Fourier's Law (also referred to as normal transport), has been observed. But there is still no consensus on the criteria required for a system to exhibit normal transport. The situation is even less clear-cut for quantum systems. Most experimental [19] and theoretical work [34] have focused on spin and energy transport of spin chains or spin ladders. It has been conjectured, based on investigations using the standard Kubo formula [34, 35], that integrability is the criterion for anomalous conductivity. Other studies based on dier- ent approaches, however, reach dierent conclusions. For example, it has been suggested that the existence of gap leads to anomalous transport [77] and also that spin transport could be qualitatively dierent from energy transport [59]. Experimental results on integrable systems range too from anomalous conductance [78{80] to normal conductance [81]. Most of the previous studies use the standard Kubo formula [34, 44, 82]. Some are based on the local-operator Lindblad equation [56, 59] and some are based on the Redeld equation [36, 37, 83]. However, the spin systems studied via the Redeld equation are either non-interacting (and thus triv- ially integrable) or non-integrable. We denote, for example, the XY chains as non-interacting since via Jordan-Wigner transformation [84] they can be mapped to systems of free fermions. In the same sense we regard the integrable Heisenberg chains to be interacting since by the same transfor- mation the resulting fermionic systems are interacting. Anomalous trans- port was found for (trivially) integrable systems and normal transport for non-integrable systems. Disordered spin systems, which are non-integrable, 10 1.4. Summary have also been investigated and found to have normal transport [37]. Prior to our work, there were no integrable interacting systems investigated via the Redeld equation. We apply the Redeld equation to investigate the energy transport of spin 1=2 chains. We also nd anomalous transport for XY chains, but normal transport for XXZ chains, which are integrable but interacting. Therefore, we conclude that integrability does not necessarily imply anomalous transport. We then investigate a wider class of systems and conjecture another criterion: interacting systems have normal transport. The reliability of this conjecture needs to be checked further. 1.4 Summary In summary, this chapter has presented a brief synopsis of the challenges that we are trying to address in this thesis, and the directions we have used to achieve these goals. The rest of the thesis will now discuss each of these issues in detail. 11 Chapter 2 The projector technique and the Redeld equation 2.1 Outline In this chapter we review the projector technique for central systems coupled to single or multiple baths, and illustrate its application to small systems with only a few degrees of freedom, where the equations can be solved di- rectly. The challenge it faces when applied to large systems, and our solution to this challenge, will be the topics of other chapters. This chapter is organized as follows. Section x2.2 contains a brief re- view of an open system's dynamics based on the system-bath scenario, as a preparation for our method. In Section x2.3, we discuss how systems cou- pled to a single bath, or to multiple baths held at the same temperature and the same chemical potential, arrive at thermal equilibrium. As examples, the Redeld equations of both a single-mode central system and a two-site central system are solved to illustrate the general procedure. Besides the Redeld equation, there are other various forms of open-system quantum master equations. Usually they can be regarded as generalizations of the Redeld equation of a single-mode central system. Two such equations, which are widely used in literature, are presented in Section x2.4. Then, in Section x2.5 we generalize the Redeld equation to systems coupled to two or more baths of dierent temperatures but with the same chemical potential, and discuss their heat conductance. The method is also valid for electrical and spin transport, where the chemical potentials of baths could also be dierent. 2.2 The eective equation of motion for open systems In this section we derive a general eective equation of motion for a central system coupled to a single bath, starting from the Schrodinger's equation for 12 2.2. The eective equation of motion for open systems the total composite system. A similar derivation can be found in Refs. [1, 2, 51, 85]. Unless otherwise stated, in this thesis we use ~ = 1; kB = 1; c = 1. Consider a system S interacting with a single bath/reservoir B. In the Schrodinger picture, the Liouville equation for the evolution of the density matrix of the total closed quantum system S +B is: i @̂T (t) @t = h Ĥ; ̂T (t) i L̂̂T (t) : (2.1) Here we denote the density matrix of the total composite system as ̂T and the reduced density matrix of the central system (the bath) as ̂ (̂B), ̂ (t) = trB (̂T (t)) ; (2.2) ̂B (t) = tr S (̂T (t)) ; (2.3) where trB (trS) stands for the operation of partial trace over the bath (the central system) degrees of freedom. The total Hamiltonian has the generic form: Ĥ = Ĥ0 (S;B) + V̂ (S;B) ; (2.4) where Ĥ0 (S;B) = ĤS (S) Î (B) + Î (S) ĤB (B) (2.5) describes the system, respectively the bath, and V̂ (S;B) = X i L̂i (S) F̂i (B) (2.6) describes their interactions, with L̂i (S) and F̂i (B) being system, respec- tively bath operators. Î (B) and Î (S) are the identity operators in the Hilbert space of the bath and the central system, respectively. If we are interested only in the properties of S, we call it an open quantum system. Typically, we need to calculate the ensemble average of physical quantities depending on operators acting on S only, hAi = tr [ÂS Î (B)]̂T : (2.7) This can be trivially rewritten as: hAi = trS ÂS ̂ : (2.8) Therefore, it is convenient to try to nd ̂ (t) directly, and avoid solving for the total density matrix ̂T (t), which is a much more complicated problem 13 2.2. The eective equation of motion for open systems due to the signicant increase in the number of degrees of freedom. Instead we are looking for the evolution equation of the reduced matrix density in a form similar to Eq. (2.1): @̂ (t) @t = L̂ ̂ (t) ; (2.9) but now L̂ is obviously neither h Ĥ; i of the whole system nor h ĤS ; i of the central system only. The task of this section is to derive the proper Liouville operator starting from knowledge of the total Hamiltonian Ĥ. There are two main types of approaches to discuss the dynamics of an open system: based on dynamic variables or based on density matrices. The rst approach uses equations of motion for dynamic variables such as X̂; P̂ etc. The dynamic variables of the bath are integrated out, resulting in noise-sources terms 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 central system is known. From a quantum mechanics point of view, this approach is similar to the Heisenberg picture. The second approach focuses on distribution functions in the phase space of the system, such as (x; p) for classical systems and the density matrix ̂ 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 the Schrodinger picture. We will use the latter approach. In this section we summarize the general formalism of this method. 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 rst one in reasonable detail and then shortly review the second one. 2.2.1 The projector operator technique and the reduced Liouville equation One way to derive the unknown L̂ in Eq. (2.9) is through the technique of projection. We dene the projector operator on the Hilbert space of the total system S +B, P : H ! H, such that P (̂T ) , trB (̂T ) ̂eqB = ̂ ̂eqB (2.10) where ̂eqB is the density matrix of the isolated bath in equilibrium, although this requirement could be relaxed. It is easy to verify that P is a projector, 14 2.2. The eective equation of motion for open systems i.e. that PP = P. We denote Q = IP, and rewrite the equation of motion for the total system, Eq. (2.1), as: @ @tP ̂T (t) = PL̂P ̂T (t) + PL̂Q̂T (t) @ @tQ̂T (t) = QL̂P ̂T (t) +QL̂Q̂T (t) (2.11) The second equation is integrated formally and its solution is substituted into the rst equation, to obtain: @ @t P ̂T (t) = PL̂P ̂T (t) + PL̂ Z t 0 de(t)QL̂QL̂P ̂T () + PL̂etQL̂Q̂T (0) : (2.12) If we start from the reasonable initial condition ̂T (0) = ̂ (0) ̂eqB ; (2.13) then Q̂T (0) = 0 (2.14) and the last term vanishes identically. Let LS ;LB;LV be the Liouville op- erators corresponding respectively to ĤS (S) ; ĤB (B) and V̂ (S;B). Then, L̂ = L̂S + L̂B + L̂V , L̂0 + L̂V : (2.15) Since L̂B ̂eqB = 0, we can rewrite the dynamic equation as @ @t ̂ (t) ̂eqB = PL̂S ̂ (t) ̂eqB + PL̂1 ̂ (t) ̂eqB (2.16) +PL̂ Z t 0 de(t)QL̂QL̂ ̂ (t) ̂eqB : (2.17) 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.6)] vanishes: hbj F̂i (B) jbi = 0: (2.18) In this case, we have: i hbj L̂1 ̂ (t) ̂eqB jbi =X i hbj h L̂i (S) F̂i (B) ; ̂ (t) ̂eqB i jbi = X i L̂i (S) ̂ (t) hbj F̂i (B) ̂eqB jbi ̂ (t) L̂i (S) hbj ̂eqB F̂i (B) jbi = 0 15 2.2. The eective equation of motion for open systems showing that indeed, PL̂V P = 0. This condition, which is satised for most types of couplings, is not however necessary. If it is not obeyed, there will be a term of the form P i;b eEb Z D b F̂i bE hL̂i; i, which can be added to L̂S . It follows that this assumption is not essential, although all the coupling Hamiltonians we study satisfy it. With these assumptions and after tracing out the bath degrees of freedom, we nd that: @ @t ̂ (t) = L̂S ̂ (t) + trB PL̂ Z t 0 de(t)QL̂QL̂P ̂T (t) : (2.19) Subject to the two assumptions mentioned above, this equation is exact. Clearly, all the eects of the bath on the system's dynamics are contained in the second term on the right-hand side. This term can be further simplied. One can verify that QL̂P = QL̂0P +QL̂V P = QL̂V P = PL̂V P + L̂V P = L̂V P (2.20) where the second assumption has been invoked again. Now we use the rst essential assumption, replacing: e(t)QL̂ ! e(t)QL̂0 : (2.21) 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 sim- plify: e(t)QL̂0 = e(t)(IP)L̂0 = e(t)(L̂0PL̂0) = e(t)L̂0 ; (2.22) where we make use of PL̂0 = PL̂0P + PL̂0Q = 0. Therefore we obtain: @ @t ̂ (t) = L̂S ̂ (t) + trB L̂V Z t 0 de(t)L̂0L̂V ̂ () ̂eqB ; (2.23) where indeed the second term is of second order in the interaction V̂ . This can be further simplied in the interaction picture , where we dene ̂I (t) = eiĤSt̂ (t) eiĤSt; (2.24a) L̂Ij (t) = e iĤStL̂je iĤSt; (2.24b) F̂ Ij (t) = e iĤBtF̂je iĤBt: (2.24c) 16 2.2. The eective equation of motion for open systems In this case, the rst term cancels out, and the dynamic equation becomes: @̂I (t) @t = X jl Z t 0 d nh L̂Ij (t) ; L̂ I l (t ) ̂I (t ) i Gjl () h L̂Ij (t) ; ̂ I (t ) L̂Il (t ) i Glj () o (2.25) where Gjl () = tr B F̂ Ij () F̂ I l (0) ̂ eq B : (2.26) are Green's functions which depend only on the bath characteristics and can be calculated directly. Depending on the specic form of the Hamiltonian, Eq. (2.25) may or may not have a unique long-term stationary solution. For example, it may oscillate between several metastable states. The exact criteria guaranteeing the existence of such a steady-state are not known. Finding such criteria would be an interesting topic for further research. In the following, we take the pragmatic approach and assume that if a stationary steady-state solution has been found, then it is unique. Now comes our second essential assumption, the Markov approximation: ̂I (t ) ̂I (t) for t c; (2.27) where c is the characteristic relaxation time of Gjl () (i.e., of the bath). Since in this thesis we investigate the steady-state solution when t!1, the Markov approximation is likely justied. Using the in
uence functional tech- nique, similar but more complicated equations can be derived if the Markov approximation is relaxed [47{50, 86]. These more complicated equations are necessary in the study of transient process. We believe that for the purpose of study of stationary states, the Markov approximation is accurate enough. With this nal approximation, we nd the desired equation of motion for the reduced density matrix: @̂I (t) @t = X jl Z t 0 d nh L̂Ij (t) ; L̂ I l (t ) ̂I (t) i Gjl () h L̂Ij (t) ; ̂ I (t) L̂Il (t ) i Glj () o : (2.28) It can be transformed back to the Schrodinger picture, @̂ (t) @t = i [HS ; ̂ (t)] X jl Z t 0 d nh L̂j ; L̂ I l () ̂ (t) i Gjl () h L̂j ; ̂ (t) L̂ I l () i Glj () o : (2.29) 17 2.2. The eective equation of motion for open systems This equation can also be projected in the eigenbasis of HS , denoted as fjsig such that HS jsi = Es jsi. In this case, the eective equation of motion becomes d dt Is0s (t) = X m;n Rs0smn (t) I mn (t) ; (2.30) with Rs0smn = eiEs0smnt "X k sn + s0kkm +nss0m nss0m + X k s0m nkks # ; (2.31) where Es0smn = E S s0 ESs + ESn ESm (2.32) and +mkln (t) = X ij hmj L̂i jki hlj L̂j jni Z t 0 ei(E S l ESn)Gij () d (2.33) and mkln (t) = +nlkm (t) . The above eective equation of motion is sometimes called a quantum master equation. If we focus only on the diag- onal terms for Iss (t), then we obtain the secular quantum master equation, which requires the so-called secular approximation [1] in order to be consis- tent with Fermi's Golden Rule. The secular approximation is to ask that ESs0 ESs + ESn ESm = 0; (2.34) so as to eliminate the explicit time dependence of the right hand side of Eq. (2.31). Then, nally, d dt Iss (t) = X m Wsm I mm (t) X m Wms I ss (t) ; (2.35) where Wsm = + mssm + mssm. This is now consistent with Fermi's Golden Rule. While there is a considerable body of work based on such secular quantum master equations, we choose to work directly with the full eective equation of motion Eq. (2.28). 18 2.2. The eective equation of motion for open systems 2.2.2 An alternative derivation: the non-correlation approximation This is a dierent way to derive the above equation of motion for the reduced density matrix, using the so-called non-correlation approximation. Suppose that both ĤS and ĤB are solvable, and their eigenstates are respectively fjsig for ESs and fjbig for EBb . Then, a complete basis for the total Hilbert space can be chosen as fjsi jbig. In the interaction picture I (t) = ei(ĤS+ĤB)t j (t)i, the time evolu- tion of the operators describing the coupling is L̂Ii (t) = e iĤStL̂i (t) e iĤSt (2.36) and F̂ Ii (t) = e iĤBtF̂i (t) e iĤBt: (2.37) The dynamical equation of the whole system is: i @̂IT (t) @t = h V̂ I (t) ; ̂IT (t) i : (2.38) Integrating this formally and substituting the resulting expression for ̂I in the right-hand side of Eq. (2.38), we obtain: @̂IT (t) @t = i h V̂ I (t) ; ̂IT (0) i Z t 0 d h V̂ I (t) ; h V̂ I () ; ̂IT () ii : (2.39) Integrating out the bath, this becomes: @̂I (t) @t = itrB h V̂ I (t) ; ̂IT (0) i Z t 0 dtrB h V̂ I (t) ; h V̂ I () ; ̂IT () ii ; (2.40) which is still exact. Now we begin to introduce several approximations. First, we again assume hbj F̂i (B) jbi = 0; (2.41) so that the rst term in Eq. (2.40) vanishes. As before, this is not an essential assumption. Second , we also assume, ̂IT () = ̂ I () ̂IB () ; (2.42) 19 2.2. The eective equation of motion for open systems such that the second term in Eq. (2.40) can be separated and simplied. This is the essential approximation. In fact this is stronger than the corre- sponding approximation, Eq. Eq. (2.13), used in the previous section, where this factorization is assumed to be true only for the initial state. Third, we take, ̂IB () = ̂ eq B = 1 ZB eĤB : (2.43) With these assumptions, we get a simplied dynamic equation for the re- duced density matrix identical to Eq. (2.25). If we further apply the Markov condition, we arrive at the same eective equation of motion as in Eq. (2.28). 2.2.3 The quantum master equation in terms of creation and annihilation operators From now on, we will work in the Schrodinger picture and directly deal with the reduced density matrix of the central system. In this thesis, we will mostly use the second quantization language. In terms of creation and annihilation operators, Eq. (2.28) can be simplied further. Let us now derive such a simplied version for a general central system, whose th site is coupled to the th bath with bath parameters T ; and coupling constant Vk; and , HS =H aj ; a y j ; (2.44a) HB = X k ! (k) byk;bk; ; (2.44b) HSB = X k; Vk;a y bk; + h:c: : (2.44c) Here H aj ; a y j can be any function of aj and a y j . refers to the strength of the coupling between the central system and the baths. It could be absorbed into the coupling coecients Vk; but we separate it and use it as a relatively small dimensionless number for bookkeeping purpose. Here we assume that every site connected to a bath, is connected to its own bath. If two sites happen to share a bath, then the formula we derive below needs to be adjusted. First we identify L (1) k; = Vk;a ; F (1) k; = b y k; ;L (2) k; = V k;a y ; F (2) k; = bk; (2.45) 20 2.2. The eective equation of motion for open systems and note that the only non-vanishing bath Green's functions are: G1;2k;;k;0 () = trB F (1) k; ()F (2) k;0 (0) B = hnk;i ei!(k);0 ; (2.46a) G2;1k;;k;0 () = trB F (2) k; ()F (1) k;0 (0) B = h1 nk;i ei!(k);0 : (2.46b) The ;0 comes from the assumption of independent baths for each site, which means that Green's functions among operators belonging to baths coupled to dierent sites are zero. By denition, hnk;i hn (! (k) ; T ; )i = e ! (k) T 1 1 (2.47) is the average particle number in state k in bath . For simplicity of notation, hnk;i will be denoted as nk; whenever it is safe to do so. After plugging these Green's functions into Eq. (2.28), we arrive at a quantum master equation @(t) @t = i[HS ; (t)] 2 X nh ay ; m̂(t) i + a ; ̂m(t) + h:c: o ; (2.48) with operators m̂ and ̂m dened as follows: m̂ = X k jVkj2 Z 1 0 da () ei!(k) h1 nk;i ; (2.49a) ̂m = X k jVkj2 Z 1 0 day () ei!(k) hnk;i : (2.49b) Here a () = e iHStae iHSt is expressed in the Heisenberg picture. Note that we have also changed the integral limit from t to 1. Again this is reasonable only if the stationary solution itself is of interest, not the whole evolution process. This quantum master equation in creation and annihila- tion operators is the central equation of this chapter and will be used later throughout the whole thesis. Eq. (2.48) is sometimes called the Redeld equation [52, 53, 87]. In some cases, it is convenient to use the eigenstates of HS as a basis, fjEnig. In that case, using element-wise products in this basis, the expres- sions of the operators m̂ and ̂m can be rewritten as m̂ = a ; ̂m = ay ; (2.50) 21 2.3. Relaxation towards thermal equilibrium where = X m;n jmi hnj [1 n (En Em; T ; )]J (En Em) ; (2.51a) = X m;n jmi hnjn (Em En; T ; ) J (Em En) ; (2.51b) and J () = X k jVkj2 ( !(k)) (2.52) is usually called the spectrum density function [36, 48] of the bath . 2.3 Relaxation towards thermal equilibrium We rst apply this general procedure to the simplest example | a single-site fermionic system coupled to a bath | and show that the system evolves to the correct thermal equilibrium as t ! 1. We then brie
y investigate a corresponding two-site system. 2.3.1 One-site open system Consider a one-state fermionic system in a strong magnetic eld, describing for instance a single quantum dot, dened by a Hamiltonian: HS = a ya (2.53) where a; ay are fermionic operators. The system is coupled to a single bath ( = 1 only), which is regarded as an ensemble of fermions with xed T and , so that both particles and energy are exchanged. Using the general quantum master equation Eq. (2.48), we have @ @t = i h aya; i 2 nh ay; m̂ i + a; ̂m + h:c: o ; (2.54) where m̂ = J () [1 n()] a; (2.55a) ̂m = J ()n()ay; (2.55b) and n() is shorthand for n(; T; ) = e T 1 1 , the average number of particles with energy in the bath (since there is a single bath, there is no 22 2.3. Relaxation towards thermal equilibrium B S Figure 2.1: Sketch of a one-site system coupled to a heat bath. need to explicitly write down T and everywhere). J () can be dened accordingly for this specic system using Eq. (2.52). Now we can solve the dynamical equation and check the stationary solu- tion. The density matrix of the above system can be written in the general form: (t) = 1 + p (t) 2 j0i h0j+ 1 p (t) 2 j1i h1j+ q (t) j0i h1j+ q (t) j1i h0j ; (2.56) so that the trace is preserved. Then we get evolution equations for these parameters from the Redeld equation: d dt p (t) = 2 (1 2n()) 2p (t) ; (2.57a) d dt q (t) = q (t) : (2.57b) This has a stationary solution, independent of the initial state: p (1) = 1 2n() = e () 1 e() + 1 ; (2.58a) q (1) = 0; (2.58b) which leads indeed to the expected equilibrium grand-canonical distribution (1) = 1 e() + 1 j0i h0j+ e () e() + 1 j1i h1j (2.59) for the two-level system. Note that in this example, the value of itself is irrelevant: this equilibrium state will be reached irrespective of the strength of the coupling to the bath. As we show later, this is not true for non- equilibrium stationary states, which depend on the details of the coupling between the central system and the baths. 23 2.3. Relaxation towards thermal equilibrium 1 2 Bath Figure 2.2: Sketch of a two-site system coupled to a heat bath. 2.3.2 Two-site open system Let us now consider a two-site fermionic system coupled to one heat bath, see Fig. 2.2. The system's Hamiltonian is taken to be a simple hopping Hamiltonian: HS = ay1a2 + a y 2a1 ; (2.60) and we assume that only site 1 is coupled to the bath, so that: V = X k Vka1b y k + V k a y 1bk : (2.61) Using the general eective equation of motion Eq. (2.48), we arrive at d dt = i [HS ; ] 2 nh ay1; m̂1 i + a1; ̂m1 + h:c: o ; (2.62) where m̂1 = J () [1 n ()] 2 (a1 + a2) + J () [1 n ()] 2 (a1 a2) ; (2.63a) ̂m1 = J ()n () 2 ay1 + a y 2 + J ()n () 2 ay1 ay2 ; (2.63b) and J () can be calculated similarly as in Eq. (2.52) but now for the eigenenergies . Again, for simplicity we assume J() = 1. Note that when = 0, the operators m̂1 and ̂m1 become m̂1 = J (0) [1 + n (0)] a1; (2.64a) ̂m1 = J (0)n (0) a y 1; (2.64b) 24 2.3. Relaxation towards thermal equilibrium in agreement with results in the previous section. This shows that only if is zero, then the operators m̂1 and ̂m1 depend only on a1 and a y 1. Otherwise, they also involve the creation and annihilation operators for other sites in the system (here, site 2). In fact, unless ay happens to be the creation operator of an eigenmode of HS , which is generally not the case, ̂m involves other operators than a ; a y . We will revisit this issue in Section x2.4. The eigenmodes of this Hamiltonian are, HS = ay+a+ + aya; (2.65) where a = a1 a2p 2 : (2.66) In terms of occupation states of these eigenmodes jn+ni, eigenstates are j00i, j10i, j01i, j11i with eigenvalues respectively 1 = 0, 2 = , 3 = , 4 = 0. A density matrix in this basis = 2664 p11 p12 p13 p14 p21 p22 p23 p24 p31 p32 p33 p34 p41 p42 p43 p44 3775 (2.67) can be cast into a 16-dimensional vector, P = [p11; p12; ; p44]T : (2.68) 25 2.3. Relaxation towards thermal equilibrium In this basis, operators are represented by: a1 = 1p 2 2664 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 3775 ; (2.69) a2 = 1p 2 2664 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 3775 ; (2.70) m̂1 = p 2 2664 0 1 n() 1 n() 0 0 0 0 1 n() 0 0 0 1 + n() 0 0 0 0 3775 ; (2.71) ̂m1 = p 2 2664 0 0 0 0 n() 0 0 0 n() 0 0 0 0 n() n() 0 3775 : (2.72) Plugging the above operators and density matrix into Eq. (2.62), we can derive the explicit form of matrix and rewrite Eq. (2.62) in the following form, d d P = P; (2.73) where the matrix has dimension 1616, or generally 4N4N for an N -site spinless fermionic system. Fortunately, since we are only interested in equilibrium states in this section, we know that the o-diagonal elements of the equilibrium density matrices will vanish. Therefore, we work with P dia = [p11; p22; p44; p44] T in this section and derive and solve only its evolution equation. In this case, the matrix can also be reduced to a matrix dia. Using the above operators a1, a y 1, m̂1, and ̂m1, we nd dia = 2 2664 n() n() 1 n() 1 n() 0 n() n() n() 1 0 1 n() n() 0 n() n() 1 1 n() 0 n() n() n() + n() 2 3775 ; (2.74) 26 2.3. Relaxation towards thermal equilibrium 21 Left bath Right bath Figure 2.3: Sketch of a two-site system coupled to two heat baths. and the reduced equation becomes d d P dia = diaP dia: (2.75) A general time-dependent solution of Eq. (2.75) with initial condition P dia0 is P dia () = e diaP dia0 : (2.76) It converges to a long term stationary state only when all eigenvalues of dia are non-positive. Its stationary state is the eigenvector of dia with eigenvalue equal to zero, diaP dia (1) = 0: (2.77) We veried that dia indeed only has one zero eigenvalue, and that its cor- responding eigenvector P dia (1) leads to the proper thermal equilibrium state, P dia (1) = 1 Z h 1; e(+); e(); e2 iT ; (2.78) where Z = 1 + e(+) + e() + e2. Up to this point, we have illustrated the general procedure for a system coupled to one bath and conrmed that its stationary state is the expected thermal equilibrium. A straightforward extension is to consider a system coupled to multiple baths, for example dierent sites in the system are locally coupled to dierent baths. In fact, our general formula Eq. (2.48) has already been derived to be applicable for such situations. The only change to be made is to allow every bath to have its own operators m̂ and ̂m . 27 2.4. Various other extensions Specically, if as sketched in Fig.2.3, we couple site 2 to a bath 2 with the same bath parameters, then our eective equation of motion reads: d dt = i [HS ; ] 2 nh ay1; m̂1 i + a1; ̂m1 + h:c: o 2 nh ay2; m̂2 i + a2; ̂m2 + h:c: o ; (2.79) where the additional operators are dened as m̂2 = [1 n2 ()] 2 (a1 + a2) [1 n2 ()] 2 (a1 a2) ; (2.80a) ̂m2 = n2 () 2 ay1 + a y 2 n2 () 2 ay1 ay2 : (2.80b) Repeating the above calculation, we nd the same stationary state. It is now natural to ask what kind of stationary state will emerge if the two baths have dierent parameters, for example dierent temperatures and/or chemical potentials. Ideally this should be the non-equilibrium sta- tionary state which can be used to calculate transport properties.1 We will use examples of simple systems to investigate this question in x2.5. Before that, however, we have a subtle issue to settle, in the next section. 2.4 Various other extensions of the single-site quantum master equation There are several forms of the eective equation of motion used in the lit- erature. Besides the multi-bath Redeld equation derived above, a local- operator Lindblad equation and a multi-mode quantum master equation are also used quite often to study transport or relaxation processes. In this section, we will introduce the latter two and discuss their limitations. Another way to write Eq. (2.54) for a single-site system is: @ @t = i [HS ; ] J () n (1 n()) aya+ aya 2aay +n() aay+ aay 2aya o : (2.81) 1What we have argued above shows only that it is natural to link this state to the non-equilibrium stationary state, since the proper equilibrium states are reached when the system is coupled to unbiased baths. However, we note that this is not a real rst- principle derivation: thermal equilibrium of the baths has been explicitly assumed. The question of establishing thermal equilibrium from only rst principles is not the scope of this thesis although many mathematicians and physicists are very much interested in this issue. 28 2.4. Various other extensions It is very tempting to generalize the above form directly to multi-site systems coupled to multiple local baths as follow, @ @t = i [HS ; ] X J () n (1 n()) aya+ a y a 2aay +n() aa y + aa y 2aya o ; (2.82) i.e. as if each new bath adds the contribution it would if the site it is coupled to was isolated from the rest of the central system. Here, is the on-site energy of site and n() is the occupation number for a mode of energy given the parameters of bath . We call this simplied equation the local-operator Lindblad equation. Because of its simplicity, it has been used often for the study of transport properties [58, 59, 88, 89]. We have seen, however, that Eq. (2.48) or more specically Eq. (2.79) can not be cast into this form because the operator m̂ usually involves not only a , but also the operators for all other sites a0 through formulas that depend on the particularHs of the system. Therefore, the two equations are dierent. The question, then, is whether the stationary states of the two equations are qualitatively dierent, and if yes, whether the solution of the local-operator Lindblad equation still captures most features of the non-equilibrium stationary states properly? These questions have not been answered yet. In this section, we will provide some preliminary results. Before that, we also introduce the multi-mode quantum master equation that arises from another extension of the Redeld equation for a single- site system [70]. For multi-site systems, let us assume we can rewrite the Hamiltonian as follows: HS = X q (q) ayqaq; (2.83) We further assume that each eigenmode of the system is coupled to its own individual bath with temperature Tq: V̂ = X k;q h Vqkaq byq;k + h:c: i (2.84) where creation and annihilation operators of bath q satisfy the fermionic commutation relation, fbq;k; byq0;k0g = k;k0q;q0 . Note that this is dierent from our setup of Redeld equation, where some of the sites, not eigenmodes, are coupled each to its own individual bath. 29 2.5. Evolution towards a NESS For this new setup, every eigenmode is independent so that the operators m̂q and ̂mq involve only aq and a y q. Therefore, from Eq. (2.48) the eective equation of motion reads d dt = i X q (q) h ayqaq; i X q Jq ( (q)) n (1 nq( (q))) ayqaq+ a y qaq 2aqayq +nq( (q)) aqa y q+ aqa y q 2ayqaq o ; (2.85) i.e. it is simply a sum over all eigenmodes, consisting of contributions similar to those in Eq. (2.54). We call this equation the multi-mode quantummaster equation. It has been used in quantum optics to study relaxation processes. We should emphasize that the two equations Eq. (2.82) and Eq. (2.85) are dierent from each other, except for a single-site system in which case all three equations are the same. For example, for our two-site system, a in Eq. (2.82) refers to a1; a2 while aq in Eq. (2.85) refers to a = a1a2p2 . In the local-operator Lindblad equation a1; a y 1 does not mix with a2; a y 2, while in the Redeld equation there are terms mixing them. Next, we investigate the dierence between these three equations: the Redeld equation, the local-operator Lindblad equation and the multi-mode quantum master equation. We will argue rstly that while the local-operator Lindblad equation gives both wrong equilibrium states and wrong non- equilibrium stationary states, it may lead to qualitatively correct results for some physical quantities such as certain currents; and secondly that while the proper equilibrium states can be found from the multi-mode quantum master equation, it does not describe any proper non-equilibrium stationary states. In order to do this, we begin by nding the \proper" non-equilibrium stationary state from the Redeld equation of a two-site system coupled to two baths. This solution is discussed in the next section. 2.5 Evolution towards a non-equilibrium stationary state In this section, we rst solve the Redeld equation for a two-site system locally coupled to two baths, and then compare it with solutions from the local-operator Lindblad equation and the multi-mode quantum master equa- tion. The general procedure presented in this section is the general method 30 2.5. Evolution towards a NESS underlying all the later work in this thesis. In a sense, all the later work focuses either on the technical side of this procedure or on applications of this procedure. Of course, these are non-trivial because solving the Redeld equation directly for large systems is a dicult task. 2.5.1 Non-equilibrium stationary states from the Redeld equation We begin from Eq. (2.79). Note that when the bath parameters are dierent: T1 6= T2 and/or 1 6= 2, unlike in the case of thermal equilibrium states, the o-diagonal elements of the density matrix are no longer necessarily zero. As a result, we need to nd all the elements of the density matrix. Correspondingly the matrix has d4 = 28 elements. For such a large linear system we have to solve the problem numerically. We are interested in its stationary solution, i.e. P (1) = 0: (2.86) This can be solved by nding the eigenvector of L corresponding to the zero eigenvalue. Or, it can also be solved as the solution of the following linear system of equations, P (1) = V ; (2.87) where V = [1; 0; ]T and is dened by replacing the rst row of by the normalization condition tr () = 1, i.e. P j Pjd+j = 1 with our notation. The advantage of this is that is a singular matrix but is not. Therefore, while Eq. (2.86) has to be treated as a computationally costly eigenvalue problem, Eq. (2.87) is a non-singular system of linear equations that can be solved more eciently. For simplicity of notation in later chapters, the above two equations dening the non-equilibrium stationary states are also written as L (1) = 0: (2.88) and L (1) = : (2.89) Here instead of superoperator ; and supervector P;V , we use directly the Liouvillian L; L and matrix ; v, which are dened from mapping the supervector P;V in the above supervector representation back into the usual 31 2.5. Evolution towards a NESS density matrix representation. These two sets of equations are equivalent, but the later has simpler notation since it is still formulated in terms of the Liouvillian and density matrices. Let us summarize the general algorithm for this direct method: 1. Find all eigenvalues and eigenstates of HS , and express all the original operators a and a y in this basis of eigenstates; 2. Evaluate the correlation functions for every bath to calculate the ma- trices and dened in Eq. (2.51); 3. From these, derive the expressions of the operators m̂ and ̂m as in Eq. (2.50); 4. Construct or using Eq. (2.48) and the previous results; 5. Solve for the zero eigenvectors of or the linear system dened by in Eq. (2.89). Using this procedure we solve Eq. (2.79). Typical results are shown in the Fig. 2.4 for the specic parameters = 1; = 0:1; 1 = 0 = 2; T1 = T 1 + 2 ; T2 = T 1 2 , and T = 2:0. We have calculated the dierence between the general non-equilibrium stationary state (1) and the thermal equilibrium Eq (T ) at temperature T . The results are presented separately for o-diagonal elements and diagonal elements, ddiag = vuuutPj jj Eqjj 2P ij (ij) 2 ; d off = vuuutPi6=j ij Eqij 2P ij (ij) 2 : (2.90) The charge current operator is dened as J = ie ay1a2 ay2a1 ; (2.91) and we also separate its expectation value into contributions from the diag- onal and o-diagonal parts, Jdiag = X j Jjjjj ; J off = X i6=j Jjiij : (2.92) From Fig. 2.4(a) we see that compared to the diagonal part, the o-diagonal elements of the non-equilibrium stationary states are relatively small but 32 2.5. Evolution towards a NESS 0 0.5 1 1.5 2 δ 0 0.05 0.1 0.15 0.2 d ddiag doff 0 0.5 1 1.5 2 δ 0 0.005 J Jdiag Joff Figure 2.4: Density matrices for non-equilibrium stationary states are com- pared to those of equilibrium states. In (a) ddiag and doff are plotted vs. the temperature bias . Diagonal elements of the density matrices for non- equilibrium stationary states are found to be quite dierent from those of equilibrium states. The dierence between o-diagonal elements of non- equilibrium stationary states and equilibrium density matrices are much smaller. Note that the o-diagonal elements of equilibrium density matrices are zero always. (b) Charge currents, Jdiag and Joff , for the non-equilibrium stationary states, are plotted vs. . We see that only the o-diagonal terms of the density matrix contribute to the current. 33 2.5. Evolution towards a NESS non-zero. However, Fig. 2.4(b) shows that the diagonal part does not con- tribute at all to the current, Jdiag = 0; the current comes entirely from the o-diagonal part. So the seemingly small change in the o-diagonal part makes a big qualitative dierence: o-diagonal elements are essential to de- scribe transport properties. This is not surprising since we know that if the time-reversal symmetry is not broken, the eigenstates of HS do not support currents, Jjj = 0. Next, we will perform a comparison between stationary states from the local-operator Lindblad equation/multi-mode quantum master equation and these non-equilibrium stationary states obtained from the Redeld equation, to see whether they agree or not. 2.5.2 Non-equilibrium stationary states from the local-operator Lindblad equation and the multi-mode quantum master equation We have claimed before that the local-operator Lindblad equation leads to wrong equilibrium states but possibly qualitatively correct non-equilibrium stationary states, while the multi-mode quantum master equation gives proper equilibrium states but not non-equilibrium stationary states. In this section, we will conrm these statements by presenting results for the simple example of a two-site fermionic system coupled to two baths hold at dierent temperatures but the same chemical potential. The Redeld equation for this system has been written down in Eq. (2.79). Following the general form in Eq. (2.85) and Eq. (2.82), the local- operator Lindblad equation and the multi-mode quantum master equation for this specic system are, respectively, @L @t = i[HS ; L] 2 n h1 + n (0; T1; 1)i h ay1; a1 L i + hn (0; T1; 1)i h a1; a y 1 L i + h1 + n (0; T2; 2)i h ay2; a2 L i + hn (0; T2; 2)i h a2; a y 2 L i + h:c: o : (2.93) 34 2.5. Evolution towards a NESS and @M @t = i[HS ; M ] 2 n h1 + n (; T; )i h ay+; a+ M i + hn (; T; )i h a+; a y + M i + h1 + n (; T; )i h ay; a M i + hn (; T; )i h a; a y M i + h:c: o ; (2.94) where n (; T; ) is the average particle number in state with energy in a bath with temperature T and chemical potential . For the multi-mode quantum master equation, in which every eigenmode | not every site | is treated independently, it is impossible to incorporate the local coupling (the fact that only certain local sites couple to their own individual local baths), so we have just set all Tq = T . This makes the multi-mode quantum master equation not applicable to calculation of non-equilibrium stationary states. We will nd so in the following numerical study. Following a similar procedure like in the last section, we solve for the stationary states, L (1) and M (1). We compare them against the non- equilibrium stationary states of the Redeld equation, R (1) using the following measure for distances: dAB = vuuuut P i;j Aij Bij 2 P ij Bjj 2 : (2.95) Typical results are shown in Fig. 2.5. We see that the dierence between the multi-mode quantum master equation and the Redeld equation is rela- tively small compared to the dierence between the local-operator Lindblad equation and the Redeld equation. At equilibrium ( = 0), the multi-mode quantum master equation and the Redeld equation lead to the same sta- tionary state, which is of course the thermal equilibrium. The local-operator Lindblad equation is quite dierent from the Redeld equation for all values of , and in particular it fails to predict the proper thermal equilibrium for unabiased baths. It can be shown that this is always the case, even for more general systems. However, we cannot discard the local-operator Lindblad equation on this basis, when analyzing transport, since we just argued that transport is re- lated only to o-diagonal matrix elements. Fig. 2.5(b) shows that while 35 2.5. Evolution towards a NESS 0 0.5 1 1.5 2 δ 0 0.1 0.2 0.3 0.4 0.5 d dMR dLR 0 0.5 1 1.5 2 δ 0 0.01 0.02 J JR JM JL Figure 2.5: L (1) and M (1) are compared with R (1). In (a) we plot dLR and dMR vs. the temperature bias . We see that at equilibrium ( = 0), the dierence between M (1) and R (1) is zero but the dierence between L (1) and R (1) is not. For all other cases, the dierences are nite. However, as shown in (b), the steady-state charge current JL and JR, re- spectively calculated from L (1) and R (1), show qualitatively similar behavior. On the other hand, JM = 0 shows that the multi-mode quantum master equation does not capture non-equilibrium stationary states at all. 36 2.5. Evolution towards a NESS the charge current from the multi-mode quantum master equation is al- ways zero, the currents from the local-operator Lindblad equation and the Redeld equation are qualitatively similar. The rst fact conrms that the multi-mode quantum master equation captures properly only the equilib- rium distribution, but cannot describe transport through the system, as should be obvious from its structure. The second fact suggests that the local-operator Lindblad equation may be used to study, qualitatively, trans- port properties. However, since the ratio between the two currents changes for dierent parameters, it is impossible to link JR quantitatively to JL. This qualitative similarity, if conrmed for more complicated problems, could still be valuable for qualitative studies, for example in the question of normal v.s. anomalous conductance. More investigation are needed to check how universal this similarity is, but they are hampered by the lack of ecient solutions for the Redeld equation. If conrmed, this similarity would be useful since it is much easier to solve the local-operator Lindblad equation than the Redeld equation [59]. While not a proof, the results shown here for a 2-site system coupled to 2 baths clearly illustrate why the local-operator Lindblad equation and the multi-mode quantum master equation approaches are questionable when studying transport through open systems. This is why we focus on the Redeld equation from now on. To summarize, up to this point we have derived the general eective equation of motion and the more specic Redeld equation, and we have also illustrated the general procedure of solving the Redeld equation by brute force direct methods. Before we develop more ecient approaches in the next chapters, we would also like to discuss its relation to the most general evolution equation, the semigroup Lindblad form. This is the topic of the next section. 2.5.3 The general Lindblad form The evolution operators of a closed quantum system H, U (t) = eiHt, form a group. In general, however, the evolution operators of an open system form a semigroup so as to keep the density matrix always positive and normalized, but not necessary reversible, (t) = U (t) (0) (2.96) such that U (t1 + t2) = U (t1)U (t2) : (2.97) For a closed system, we have the additional condition: U (t) = U1 (t) (2.98) 37 2.6. Major challenge for large systems which insures reversibility. Lindblad discussed the generators of such trace-preserving positive semi- groups [90]. As shown in a theorem by Gorini, Kossakowski and Sudarshan [91], a completely positive semigroup evolution of a quantum system in a n-dimensional Hilbert space can be characterized by a Lindblad operator that governs its evolution d dt (t) = L (t) ; (2.99) where generally: L (t) = i h Ĥ; i + 1 2 n21X ij Aij h Li; L y j i + h Li; L y j i ; (2.100) where the operators fLig form an orthonormal basis of the operators on the system's Hilbert space, and the constants coecient matrix A must be positive. Note that the local-operator Lindblad equation is of the Lindblad form, where Li = a and Aij = Aiij , and = L;R stands for the left and the right baths respectively. This is why, especially in practical usage, some [59, 75] refer to the local-operator Lindblad equation as the Lindblad equation. More importantly, the Redeld equation of Eq. (2.48) is not of the Lind- blad form. In principle, this means that it may lead to non-physical results such as non-positive density matrices, however we have never encountered such problems in any of our simulations. 2.6 Major challenge for large systems As it is clear from the examples presented above, the major challenge in solving the Redeld equation has to do with the dimension of the vector P . This is 22N for a N -site system of spinless fermions, (2S + 1)2N for an N -site chain of spins-S, and innity for bosonic systems. With such exponential growth, it is impossible to deal with large systems by direct methods, in- cluding the Runge-Kutta method, the eigenvector corresponding to the zero eigenvalue, the solution of the linear system and other propagator-based methods [64]. Using these methods, current computational power allows us to deal with systems with up to N = 10 [37, 38]. There are interesting physi- cal systems well within this size which already show rich behavior. However, such limitations on the system's size does not give us a lot of freedom to 38 2.6. Major challenge for large systems separate reliably bulk behavior from boundary eects, for larger systems. In contrast, quite often the non-equilibrium Green's function methods [4, 13] are capable of dealing with system sizes of up to N 200. The various methods mentioned in the introduction are either not e- cient enough, or involve additional approximations. As we mentioned earlier, ecient stochastic wave-function methods [65] have been developed for both the local-operator Lindblad equation and the Redeld equation. This has better eciency then the direct methods and it is capable of dealing with roughly N 20. For the local-operator Lindblad equation only, an ecient method has been developed based on the time-dependent density matrix renormalization [59], capable of dealing with N 100. However, we have argued that the local-operator Lindblad equation's predictions are quantita- tively dierent from those of the Redeld equation. In the rest of this thesis, we will search for more ecient methods to calculate non-equilibrium sta- tionary states of large systems for the Redeld equation, and discuss some preliminary applications. This chapter served as a foundation for our theory. None of it, except for the comparison between the Redeld equation, the local-operator Lind- blad equation and the multi-mode quantum master equation, is our original contribution. These start in the next chapter. 39 Chapter 3 Linear Response Theory: Kubo formula for open systems 3.1 Introduction Consider a nite size one-dimensional system whose ends are coupled to two baths held at dierent temperatures and/or dierent chemical potentials, leading to
ow of energy or charge through the system. The question of interest to us is to nd the stationary state and thus physical quantities, in particular the steady-state charge and heat currents. Since it is hard to solve exactly for the non-equilibrium stationary states of large systems, in this chapter we introduce an approximate solution for the Redeld equation valid for small biases. We call this approach the open Kubo formula as it is somewhat similar to the standard Kubo formula for isolated systems. While the size limitation of direct methods is around N = 10 (measured in qubit), in certain circumstances detailed below the open Kubo formula is capable of dealing with systems of up to N = 20 sites, the same size limitation as with exact diagonalization of Hamiltonians. The open Kubo formula gives the rst-order corrections towards non- equilibrium stationary states for nite-size systems, starting from the cor- responding equilibrium states. Two forms of such Kubo-like formulae are proposed. The rst, while computationally less ecient, gives proper general density matrices. The more ecient second form produces correct average values only for certain special physical quantities. Conditions for these spe- cial cases will also be discussed. Besides eciency, this work is also motived by considerations of the ap- plicability of the standard Kubo formula to nite-size open systems. The standard Kubo formula has been widely used for both innite [42] and nite- size systems [82], usually with periodic boundary conditions. Open bound- ary conditions have also been considered, leading to qualitatively dierent 40 3.2. General Linear Response Theory results [44]. For nite-size systems, the standard Kubo formula leads to a summation of many Dirac -function peaks at dierent energies. Var- ious ways to get around or to smooth out such peaks have been pro- posed [92, 93]. However, there is no guarantee that correct physical quanti- ties can be extracted via such procedures. More importantly, we show that the standard Kubo formula gives the rst-order corrections towards equilibrium states of the perturbed systems, instead of the desired non-equilibrium stationary states. We also prove that the standard formulas for Drude weight and dc conductance only work for innite-size systems, not for nite-size ones. One way to solve these issues is to explicitly take into consideration the coupling to baths. This is why we are interested in linear response theory based on the Redeld equation. This chapter is organized as follows. In Section x3.2 we brie
y review the general linear response theory, and in Section x3.3 we discuss the ap- plicability of the standard Kubo formula. In Section x3.4 we present the open Kubo formula based on the Redeld equation, and nally, in Section x3.5, for small systems, we compare the approximate solution from this open Kubo formula with the direct numerical solution of the Redeld equation, in order to gauge its validity. 3.2 General Linear Response Theory Consider a general linear equation of motion, @ (t) @t = (L0 +L) ; (3.1) where L0 is the \large" term and L is the \perturbation". Explicit forms for L0 and L will be given later. We are particularly interested in station- ary states of the above equation, (1), given by: (L0 +L) (1) = 0: (3.2) Let us assume that 0 is known and it satises L00 = 0: (3.3) Then, to rst order, = 0 is the solution of @ (t) @t = L0+L0; (3.4) 41 3.2. General Linear Response Theory where we have neglected the second order term L. A general solution of the above equation is (t) = Z t t0 deL0(t)L0 + (t0) ; (3.5) Assuming that L was turned on at t0 = 1 when the initial state was 0 so that (t0) = 0, we have the steady-state solution at t = 0 to be, (t = 0) = Z 0 1 dteL0tL0: (3.6) Changing the variable t to t, we then arrive at (1) = Z 1 0 dteL0ttL0: (3.7) Here we have inserted an innitesimal positive number ( ! 0+) to make sure that the expression converges. This is necessary if the real part of any eigenvalues of L0 are zero. The stationary state is (1) = (1)+ 0. This is the basic formulation of linear response theory. It is required that every eigenvalue of L0 must have a negative or zero real part, so that the unperturbed system described by L0 eventually reaches a stationary state. If every eigenvalue of L0 happens to have a negative real part, as will be the case in one of the methods described below, then L0 is invertible and (1) = (L0 )1L0: (3.8) This can also be derived from Eq. (3.2) since L0 is non-singular, 0 = (L0 +L) (0 + (1))) (1) = L10 L0: (3.9) This is a straightforward formula if L0 is indeed invertible but it becomes tricky if L0 has zero eigenvalues. In this latter case, we can still use Eq. (3.7), but not Eq. (3.8). We may transform L0 to a non-singular L0 using the fact that is normalized, tr () = 1 as we have done in rewriting Eq. (2.88) as Eq. (2.89). Following that we nd a corresponding expression similar to Eq. (3.8), more specically (1) = L10 L0: (3.10) 42 3.3. Limitations of the standard Kubo formula 3.3 Limitations of the standard Kubo formula The Liouville-von Neumann equation of motion for the density matrix (t) of a closed system is: @ (t) @t = LH (t) = i[H; ]; (3.11) where H is its Hamiltonian. Long-term stationary solutions of this equation are not unique. A special class of such states are the Boltzmann equilibrium distributions: eq (H) = 1 Z eH ; (3.12) where Z = tr eH , = 1=T . If H = HS + Vext, where HS is the isolated system and Vext is a static weak coupling to an external eld, one can formally use the above general linear response theory to nd a stationary solution (1) = 0+ (1) near a state 0 of the unperturbed system, LHS0 = 0: (1) = Z 1 0 dteL0ttLVext0: (3.13) If 0 = eq(HS) describes the unperturbed system in equilibrium, then this leads to the standard Kubo formula [6, 40]: (1) = i Z 1 0 dtet Vext (t) ; e HS Z ; (3.14) where Vext(t) = e iHStVexte iHSt. We can use the identity [Vext(t); eHS ] = ieHS R 0 d _Vext (t i), where _Vext (t i) = ddtVext (t) jt=ti , to rewrite: (1) = Z 1 0 dtet Z 0 d0 _Vext (t i) : (3.15) In terms of the eigenvectors of HS , HS jni = njni, hmj _Vext (t) jni = i (m n) hm jVextjni ei(mn)t; (3.16) leading to: (1) = X m;n m 6=n em en Z hm jVextjni m n i jmihnj : (3.17) 43 3.3. Limitations of the standard Kubo formula Incidentally, note that there is no contribution from states with n = m, for which hmj _Vext (t) jni = 0. This also follows directly from Eq. (3.14); if we write Vext = V 0 ext+V ? ext, where V 0 ext = P m=n hm jVextjni jm ihnj commutes with HS , then [Vext(t); 0] = [V ?ext(t); 0]. The \diagonal" part V 0ext of Vext does not contribute to (1), and consequently has no in
uence on the static response functions. The lack of diagonal contribution is, however, very puzzling. The well- known Drude weight, which is derived from the standard Kubo formula and is used quite often when discussing charge or thermal transport [42, 44, 82] reads D = L X m;n m=n em Z jhmjĴ jnij2: (3.18) It has contributions only from states with m = n. To understand the reason for this dierence, consider the derivation of Eq. (3.18) from Eq. (3.15), e.g. for spinless fermions in a one-dimensional chain (lattice constant a = 1), described by HS = t X l cyl cl+1 + h:c: + V0 X l nlnl+1; (3.19) where nl = c y l cl, plus a static electric potential Vext = X l Vlnl (3.20) induced by a homogeneous electric eld E = rV . From the continuity equation: _Vext(t) = X l Vl d dt nl(t) = X l Vl[Jl+1(t) Jl(t)]; (3.21) where Jl = it c+l+1cl c+l cl+1 is the local current operator. The sum can be changed to X l [Vl1Jl(t) VlJl(t)] = EJ(t); (3.22) where J(t) = P l Jl(t) is the total current operator. Using _Vext(t) = EJ(t) (or, on a lattice, E = Vl1 Vl) in Eq. (3.15) gives (1) = E Z 1 0 dtet Z 0 d0J (t i) : (3.23) 44 3.3. Limitations of the standard Kubo formula The dc conductivity is then = Z 1 0 dtet Z 0 d hJ (t i) Ji ; (3.24) where hOi = Tr[eq(HS)O]. This expression can be further simplied to arrive at Eq. (3.18). The only questionable step in this derivation, and the one responsible for going from a result with no contributions from states with n = m, to one with contributions only from these states, is the changeX l VlJl+1(t)! X l Vl1Jl(t): (3.25) This is only justied for an innite system (where boundary terms are pre- sumed to be negligible), or a system with periodic boundary conditions and an external eld with the same periodicity. This latter condition can only be achieved for a charge current in a nite system with periodic bound- ary conditions driven by a varying magnetic
ux through the area enclosed by the system. For a nite-size system (even one with periodic boundary conditions) in a static applied electric eld this approach is not valid. The same is true for thermal transport, which cannot experimentally be induced in a system with periodic boundary conditions. In both cases, the physical relevance of the results of Eq. (3.18) are hard to fathom. There is an even more serious conceptual problem with the standard Kubo formula: the resulting distribution ~ = eq(HS) + (1) is the rst order perturbational expansion of eq (H), not the expected non-equilibrium stationary state. This statement is proved below, where for convenience, we assume that the diagonal part V 0ext = 0. If it is not, we simply remove the \diagonal" part V 0ext from Vext and add it to HS . Consider then the eigenstates of the full Hamiltonian, Hj~ni = ~nj~ni, to rst order perturbation in Vext. Since hm jVextjni = 0 for all m = n, we can apply the rst order perturbation theory for non-degenerate states to all the states, whether degenerate or not, to nd ~n = n +O V 2ext (3.26) and j~ni = jni+ X m;m 6=n hmjVextjni n m jmi+O V 2ext : (3.27) 45 3.3. Limitations of the standard Kubo formula This immediately leads to eq (H) = X n 1 ~Z e~n j~nih~nj = eq(HS) + (1) +O(V 2ext); (3.28) where (1) is indeed given by Eq. (3.17). Intuitively this is easy to un- derstand. In a sense, we found the original thermal equilibrium distribution eq(HS) from the Liouville equation with H = HS . Then, it is not so sur- prising that a perturbation based on the same Liouville equation for the full H = HS + Vext leads to eq (H). This veries our previous statement that the standard Kubo formula implies (1)! eq(H). This is rather problematic because normally { for example if invariance to time reversal symmetry is not broken { no currents are generated in a thermal equilibrium state and therefore no steady-state transport through the closed nite system can be described by this approach. However, because we only keep the rst-order perturbational correction, the situation is less clear-cut. In principle, it is not impossible for of Eq. (3.17) to also be a rst order approximation to the true non-equilibrium stationary state, or at least to capture a sizable part from it. This needs to be investigated in more detail, and we do so below. Before that, let us also note that the fact that the use of Eq. (3.18) and its equivalents for nite size systems is problematic can also be seen from the following technical considerations. The spectrum of a nite-size system is always discrete. As a result, any pair of degenerate eigenstates m = n gives a -function contribution to the response functions. Such a singular response is unphysical for nite-size systems (for an innite system, the integration over the continuous spectrum removes these singularities). Various techniques have been proposed to smooth out these singular con- tributions in order to extract some nite values, such as use of imaginary frequencies [93] or averaging (!) over a small range of frequencies ! and then taking ! ! 0 [92]. These dierent approaches may lead to dierent results. Moreover, the order in which the various limits are approached, e.g., taking ! 0 before L!1 or vice versa, also make a dierence [92]. All of these subtleties of the standard Kubo formula are related to the potential divergence whenever m = n in Eq. (3.18). To summarize, using the standard Kubo formula blindly for nite-size systems is fraught with both conceptual and technical problems. Its solution gives the rst order correction of the equilibrium state of the full Hamilto- nian, instead of the expected non-equilibrium stationary state. The derived Drude weight formula is applicable only for innite-size systems; applying 46 3.4. OKF: the Kubo formula for open nite-size systems it directly to nite-size systems involves dealing with unphysical -function peaks. In order to x these diculties, we now discuss how to calculate the non-equilibrium stationary states from the Redeld equation, which takes the connection between system and baths explicitly into consideration. 3.4 Open Kubo formula: linear response theory for open nite-size systems For concreteness, let us assume that the central system is coupled to thermal baths kept at temperatures TL=R = T T2 and investigate the thermal transport in the resulting steady state. If T T , this will lead to a Kubo- like formula which replaces Eq. (3.15). This approach can be generalized straightforwardly to derive a Kubo-like formula for charge transport. The general Redeld equation can be written as: @(t) @t = [LH + LL(TL) + LR(TR)] (t); (3.29) where LH = i[H; ], just like for an isolated system, while LL=R are additional terms that describe the eects of the left/right thermal baths (assumed to be in equilibrium at their corresponding temperatures TL=R) on the evolution of the system. We have seen from Eq. (2.48) how the expressions for LL=R depend on the Hamiltonian H of the system and on its coupling to the baths. Generally speaking, H could also contain an external potential due to coupling to the external eld, Vext, such thatH = HS+Vext. Such a term is expected to appear for charge transport, but not for thermal transport, therefore in the following we set Vext = 0. As we have seen in Eq. (2.48), the temperature and chemical potential of the left (right) baths enter the Redeld equation only via average particle numbers n (; TL; L) (n (; TR; R)). Therefore, if T T , we can Taylor expand the particle number into a large part and a small part. Thus we also expand the LL=R and re-arrange the Redeld equation to read: @(t) @t = [LHS + LB(T ) + LP (T )] (t) = L(t); (3.30) where LB(T ) = LR(T ) + LL(T ) is the contribution from the thermal baths if both are kept at the same temperature, while LP (T ) collects the terms proportional to T . 47 3.4. OKF: the Kubo formula for open nite-size systems We are interested in the t ! 1, stationary state solution (1) of the above equation, satisfying both Eq. (2.88) and Eq. (2.89), which are rewrit- ten in the following form for convenience, L (1) = 0; (3.31) and L (1) = : (3.32) As we have shown in examples in the previous chapter, L has one non- degenerate zero eigenvalue and all its other (transient) eigenvalues have a negative real part, so (1) is unique for any value of T . We have also shown that for small systems, Eq. (3.31) and Eq. (3.32) can be solved numerically. We call this solution ex and we use it to validate the solutions of various approximation schemes. A Kubo-like formula, which is potentially more ecient, can be obtained using linear response theory. The rst step is to separate the Liouvillian L of Eq. (3.30) into a \large" plus a \small" part. There are two possible choices: ( L (1) 0 = LHS + LB(T ) L(1) = LP (T ) (3.33) or ( L (2) 0 = LHS L(2) = LB(T ) + LP (T ) (3.34) We begin with the rst choice. Assume that L (1) 0 has eigenvalues n L (1) 0; o and left/right eigenvectors fjL)g, fjR)g. As discussed, the unique (zero order in perturbation theory) steady-state solution of L (1) 0 0 = 0 is 0 = eq(HS). The deviation (1) K due to the perturbation L (1) is obtained like in Eq. (3.13): (1) (1) = X Z 1 0 dteL (1) 0;tt jR) (LjL(1)0 = X jR) (Lj L (1) 0; L(1)0 = X >0 jR) (Lj L (1) 0; L(1)0: (3.35) Note that the only divergent term, due to L (1) 0;0 = 0, disappears because (L0jL(1)0 = (0jL(1)0 = 0. To see why, we start from Eq. (3.31), 48 3.4. OKF: the Kubo formula for open nite-size systems L(0 + ) = 0, project it on (0j and keep terms only to the rst order, to nd 0 = (0j L (1) 0 +L (1) (0 + ) = (0jL(1)0 since L(1)0 0 = 0. As a result, Eq. (3.35) has only regular contributions. A similar approach has been suggested in Ref. [88], but for the local-operator Lindblad equa- tion [90] instead of the Redeld equation. Eq. (3.35) is dicult to use in practice: nding all eigenstates of L (1) 0 is a hard task unless the system has an extremely small Hilbert space. A computationally simpler solution is obtained from Eq. (3.32), where, in matrix terms, L is dened by replacing the rst row of the equation L (1) = 0 by Tr ( (1)) = 1, so that is a vector whose rst element is 1, all remaining ones being 0. As a result det L 6= 0 while det (L) = 0. We can also solve it to obtain a Kubo-like formula by dividing L = L (1) 0 + L(1). Again, the overbar shows that in matrix terms, L (1) 0 is obtained from L (1) 0 by replacing its rst row with Tr ( (1)) = 1, while L(1) is obtained from L(1) by replacing its rst row with zeros. We nd (1) (1) = [ L(1)0 ]1L(1)0: (3.36) This is much more convenient because inverting the non-singular matrix L (1) 0 is a much simpler task than nding all the eigenvalues and eigen- vectors of L (1) 0 . We have veried that both schemes produce identical re- sults on systems on which both of them can be performed. We denote 0 + (1) (1) = (1) (1). The second option is to take L (2) 0 = LHS and L (2) = LB(T )+LP (T ). In this case, we can still choose the stationary solution associated with L (2) 0 to be the thermal equilibrium state at T , 0 = eq(HS). However, this solution is no longer unique, since any matrix 0 that is diagonal in the eigenbasis of HS satises L (2) 0 0 = 0. Expanding the corresponding analogue of Eq. (3.14) in the eigenbasis of HS , we now nd: (2) (1) = i X n;m hmjL(2)0jni m n i jmihnj: (3.37) We call 0 + (2) (1) = (2) (1). Note that unlike (1) (1) of Eq. (3.36), this solution has divergent contributions from states with n = m. As such, it is analogous to the standard Kubo formula for innite systems. This is not an accident. As discussed, the standard Kubo formula for innite systems always ignores the coupling to the leads. It also takes L0 = LHS and assumes that 0 = eq(HS). Moreover, the driving force for transport 49 3.5. Comparison of the two OKFs is only the potential Vext added to HS , so that L ! LVext . With these assumptions, Eq. (3.37) maps straightforwardly into the standard Kubo formula of Eq. (3.17). 3.5 Comparison of the two open Kubo formulae To see which of these two solutions { the regular solution (1) (1) or the singular solution (2) (1) { gives a proper density matrix, we compare them against the exact numerical solution ex of Eq. (3.31) in the limit T T . We do this for an N -site chain of spinless fermions HS = t N1X l=1 cyl cl+1 + c y l+1cl + V0 N1X l=1 cyl+1cl+1c y l cl: (3.38) coupled to two heat baths, modeled as collections of fermions: HB = X k;=L;R !k;b y k;bk;; (3.39) where indexes the left and right-side baths and we set ~ = 1; kB = 1, the lattice constant a = 1, and hopping t = 1. The system-baths coupling is chosen as: VSB = X k; V k cybk; + cb y k; ; (3.40) where the left (right) bath is coupled to the rst (last) site: cL = c1 and cR = cN . Bath parameters, i.e. the temperature and chemical potential, are chosen to be (TL; ) and (TR; ) with TL=R = T T2 . The corresponding Redeld equation reads from Eq. (2.48), @(t) @t = i[HS ; (t)] 2 X =L;R nh cy; m̂(t) i + c; ̂m(t) + h:c: o ; (3.41) where operators m̂ and ̂m are dened by Eq. (2.49), which is rewritten here for convenience: m̂ = X m;n jmihnjhmjcjni (1 n ( nm)) ; (3.42a) ̂m = X m;n jmihnjhmjcyjnin ( mn) : (3.42b) 50 3.5. Comparison of the two OKFs 0 0.2 0.4 0 0.02 0.04 0.06 0.08 D δT (λ=0.1, V0=0.2) 0 40 80 120 160 D2 D1 (a) 0 0.2 0.4 δT (λ=0.1, V0=0.2) 0 0.002 0.004 0.006 0.008 J Jex J(1) J(2) 0 0.2 0.4 2.5 3 3.5 -900 -450 0 N(2) Nex N(2) (b) Figure 3.1: (a) D1 (squares) and D2 (circles) of Eq. (3.43) vs. = T=2T , for N = 8; t = 1:0; V0 = 0:2; T = 2:0; = 1:0; = 0:1; = 0:00001. (b) the steady-state electric current calculated with ex (triangles) and (1;2) (1) (squares, circles and stars). The inset shows the corresponding particle numbers. Eq. (3.41) is an example of the general Eq. (3.29). Since TL=R enters only in the Fermi-Dirac distributions, it is easy to expand when TL=R = T T2 ;T T to identify LB(T ) and LP (T ). We characterize the distance between the exact numerical solution ex and the two possible Kubo solutions (i) (1), i = 1; 2 using the norm: Di = sX n;m hnjex (i) (1) jmi2: (3.43) For the proper solution, this dierence should be small but still nite, be- cause of higher-order perturbation terms. Results typical of those found in all the cases we investigated are shown for N = 8; V0 = 0:2; = 0:1; = 10 5 in Fig. 3.1(a), where we plot D1;2 vs. = T=2T . We see that D2 (circles, axis on the right) is very large. In fact, because of the singular contributions from n = m states, D2 is divergent, with a magnitude controlled by the cuto . In contrast, D1 (squares, left axis) is small and independent of . Fig. 3.1(b) shows the electric current (J) and the total number of particles (N) calculated with ex, and (1)(1;2)(triangles, squares, respectively circles). Both N (1) and J (1) are very close to the exact values N ex, Jex. However, N (2) is very dierent from N ex while J (2) is close to Jex. These results conrm that (1) (1) of Eq. (3.36) is the proper Kubo solution. They also show that (2) (1) can also be used, but only for quantities A for which hm jAjni = 0 51 3.5. Comparison of the two OKFs 0 0.05 0.1 0.15 0.2 λ(δT=0.1, V0=0.2) 0 1 2 3 4 5 6 D D1 D2 (a) 0 0.05 0.1 0.15 0.2 λ(δT=0.1, V0=0.2) 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 J Jex J(1) J(2) 0 0.05 0.1 0.15 0.2 2.5 3 3.5 -50 -40 -30 -20 -10 0 N(2) Nex N(1) (b) Figure 3.2: (a)D1 (squares) and D2 (circles) of Eq. (3.43) vs. , for N = 8; t = 1:0; V0 = 0:2; T = 2:0; = 1:0; = 0:1; = 0:00001. (b) the steady-state electric current calculated with ex (triangles) and (1;2) (1) (squares and circles). The inset shows the corresponding particle numbers. whenever m = n, so that the divergences in Eq. (3.37) disappear. This explains the previous success of this formula to be somewhat of an accident. Fig. 3.2 shows the same quantities as Fig. 3.1, at a xed bias T as a function of the strength of the system-bath coupling . It conrms again that (1) (1) is the proper approximation of ex, but that for the charge current (2) (1) works well too. This also shows that the non-equilibrium stationary states depend on the coupling strength . This is not surprising for a nite-size system: the intrinsic conductance of the system is added to comparable \contact" contributions from the interfaces between the system and the baths, and experiments measure the total conductance. It follows that quantitative modeling of transport in nite systems will require a careful consideration of the entire experimental set-up. Fig. 3.3 shows again the same quantities as Fig. 3.1, at a xed bias T as a function of the strength of the interaction V0. It conrms again that (1) (1) is the proper approximation of ex. We also nd that values of currents Jex are quite dierent between V0 = 0 and V0 = 1 for example. The charge currents calculated from both (1;2) (1) indeed capture the major part of such dierence. One may note a discontinuity in the data around V0 = 0:6. We think this is due to numerical reasons. The calculation depends very strongly on the values of energy dierences among all energy levels. For dierent values of V0, such energy dierences may either shift continuously or change qualitatively. 52 3.6. Summary and discussion 0 0.2 0.4 0.6 0.8 1 V0(λ=0.1, δT=0.1) 0 0.5 1 1.5 2 2.5 3 D D1 D2 (a) 0 0.2 0.4 0.6 0.8 1 V0(λ=0.1, δT=0.1) 0.00065 0.0007 0.00075 0.0008 J Jex J(1) J(2) 0 0.2 0.4 0.6 0.8 1 -8 -4 0 4 Nex N1 N2 (b) Figure 3.3: D1 (squares) and D2 (circles) of Eq. (3.43) vs. V0, for N = 8; t = 1:0; T = 2:0; = 1:0; = 0:1; = 0:1; = 0:00001. (b) the steady- state electric current calculated with ex (triangles) and (1;2) (1) (squares and circles). The inset shows the corresponding particle numbers. 3.6 Summary and discussion To conclude, we rst showed that the standard Kubo formula fails to provide an approximation for non-equilibrium stationary states, instead approximat- ing the equilibrium state of the whole Hamiltonian. Secondly, applying the standard Kubo formula to nite systems generically leads to divergences, which are unphysical. We then showed that taking into consideration the coupling to baths explicitly solves both problems and leads to a well-behaved Kubo-like formula. Finally, we showed that the improper solution similar to those used in literature can give correct average values but only for certain physical quantities. Although we only considered small systems so as to be able to calculate ex, Eqs. (3.36) and (3.37) can be used for larger systems, with the latter being more ecient but not always valid. In fact, the above results give us enough grounds to question the validity of the standard Kubo formula, where the eect of the baths are taken care of only via an additional potential. We believe that this is not enough. There are two dierent types of \driving forces" responsible for a steady-state current
ow. The rst force is the applied electric eld, which gives rise to an additional interaction term Vext to be included in the total Hamiltonian H = HS + Vext of the system. The second one is the imbalance in the chemical potentials and/or temperatures of the two baths, which results in a non-equilibrium distribution function (1). As we have shown above, the standard Kubo formula calculates the current from eq(H), the thermal equilibrium state of the full Hamiltonian, so it takes into account only the 53 3.6. Summary and discussion eect of Vext, but not the full non-equilibrium distribution function. We believe that this approach is problematic when calculating charge currents. This approach is even less reliable for thermal currents, where there is only one driving force: the non-equilibrium distribution. There is no thermal equivalent for the electric potential. Some additional assumptions such as local equilibrium [39] or use of gravitational potentials [40] have been used to make the situation similar to that of charge transport. Either way, an articial thermal potential term Vext is constructed and added to HS and then one uses eq(H) instead of the proper non-equilibrium stationary state, in direct analogy with the usual approach for charge currents. The work presented in this chapter points to a solution for these concep- tual problems, if the baths are explicitly taken into consideration. However, the resulting open Kubo formula can still only be used for system of up to N 20, which may be too small for transport in solid state systems. In the next two chapters, we present two ecient ways to solve the Redeld equation for larger systems. Furthermore, unlike the open Kubo formula based on linear response theory, the new methods are applicable even when the baths' bias is not small. 54 Chapter 4 BBGKY-like hierarchy for the Redeld equation 4.1 Introduction In this chapter we introduce a much more ecient way to solve the Redeld equation, based on a BBGKY-like hierarchy. Its eciency depends on the level at which the hierarchy is truncated. At rst order, where only single- particle Green's functions are involved, it leads to a linear system with N2 unknowns for a N -site spinless fermionic system. This makes it possible to analyze systems larger than those that can be treated with the non- equilibrium Green's function method. Like the diagrammatic perturbation theory [18], the BBGKY equation hierarchy for systems in equilibrium also provides an equivalent systematic approach to calculate many-particle correlation functions [68, 69]. For non- interacting systems, the hierarchy is decoupled, meaning that the equation of a single-particle Green's function is only related to other single-particle Green's functions, forming a closed system of equations. The same holds for general n-particle Green's functions denoted as Gn. However, when there are interactions in the system, equations of all orders of Green's functions are coupled together. Generically, equations for n-particle Green's functions, Gn, also involve Gn+1, so that the system of BBGKY equations becomes innite. In this case, approximations such as the cluster expansions can be used to truncate and then solve the truncated hierarchy [69, 94]. Here we use exactly the same idea but now for the Redeld equation. Not surprisingly, for non-interacting central systems we nd that the hierarchy is decoupled, just like for systems in equilibrium. In fact, in Ref. [36], the Redeld equation with N 100 has been solved in terms of single-particle Green's functions. However, they did not regard this as a special case of the more general BBGKY-like equation hierarchy. From this point of view, our idea is basically to extend this Green's function based solution of the Redeld equation from Ref. [36] to interacting systems, where equations for 55 4.2. Derivation of the BBGKY-like hierarchy Green's functions of various orders become coupled. Two systematic approximations are proposed to truncate and solve the hierarchy. As a test of reliability of the proposed methods, we apply them to small systems which can also be solved with exact direct methods. Consis- tent results are found. Just as the usual Wick's theorem of non-interacting systems provides a basis for the usual diagrammatic perturbation theory of interacting systems, one of the approximations requires a non-equilibrium Wick's theorem. Such a theorem is proved starting from the Redeld equa- tion of non-interacting systems. This chapter is organized as follow. In Section x4.2, using a specic exam- ple, a BBGKY-like equation hierarchy is derived from the Redeld equation. In Section x4.3, non-interacting systems are discussed and a non-equilibrium Wick's theorem is proved. In Section x4.4, we introduce the two methods to truncate and solve the coupled hierarchy for interacting systems, and com- pare them against exact results. We nd that both methods are signicantly more ecient than the direct methods. The rst is capable of dealing with relatively small systems but with strong interactions, while the second can deal with much larger systems but with relatively weaker interactions. For the latter case, we demonstrate that a systematic increase in the level at which the truncation is performed leads to a systematic improvement of the results, as well as an increase in the range of interactions for which good accuracy is obtained. 4.2 Derivation of the BBGKY-like hierarchy For concreteness in presenting our general formulation, let us start from the same Redeld equation describing an N -site chain of spinless fermions coupled with two fermionic baths. Our system of interest is described by: HS = t N1X l=1 cyl cl+1 + c y l+1cl + V0 N1X l=1 cyl+1cl+1c y l cl = H0 + VS ; (4.1) and its Redeld equation, given by Eq. (3.41), is rewritten here for conve- nience: @(t) @t = i[HS ; (t)] 2 X =L;R nh cy; m̂(t) i + c; ̂m(t) + h:c: o ; (4.2) 56 4.2. Derivation of the BBGKY-like hierarchy where the operators m̂ and ̂m are dened as follows: m̂ = X k jV k j2 Z 1 0 dc () ei!k; h1 n (!k;)i ; (4.3a) ̂m = X k jV k j2 Z 1 0 dcy () ei!k; hn (!k;)i : (4.3b) Here, = L;R refer to the two baths, c correspond to the two sites at the boundaries coupled each to its corresponding bath, !k; are the energies of the modes in bath , and hn (!k;)i are their average occupation numbers. If U (t) = eiHSt is known, then so are c (t) = U y (t) cU (t) and there- fore the operators m̂. Finding explicit forms for m̂ and ̂m thus requires a full diagonalization of HS . Using its eigenvectors, one can perform the above integrals to get these operators. More details are included in Appendix A. Again we are only interested in the long-time steady-state solution (1) and we want to calculate the values of various correlation functions. For a physical quantity of the central system, denoted by the operator A, from Eq. (4.2) we nd: 0 = i h[A;H0]i+ i h[A; VS ]i+ 2 X nDh A; cy i m̂ E + [A; c] ̂m Dm̂y [A; c]E D ̂my hA; cyiEo ; (4.4) where hAi = tr (A (1)) and m̂y( ̂my) is the hermitian conjugate of m̂( ̂m). Here we have used the cyclic property of the trace operator. All equations of Green's functions in the rest of this paper will be derived from this equa- tion. For example, the rst and the second equation of the hierarchy for respectively single-particle Greens functions: G1 my; n = D cymcn E (4.5) and two-particle Greens functions: G2 my; ny;m 0 ; n 0 = D cymc y ncm0 cn0 E (4.6) can be derived by using A = cymcn, respectively A = c y mc y ncm0 cn0 in Eq(4.4): 57 4.2. Derivation of the BBGKY-like hierarchy 0 = it D cym1cn E + it D cym+1cn E it D cymcn+1 E it D cymcn1 E (4.7a) iV0 D cymc y n1cncn1 E + iV0 D cym+1c y mcm+1cn E iV0 D cyn+1c y mcn+1cn E + iV0 D cymc y m1cncm1 E (4.7b) 2 X D mcn ̂m + n ̂m y c y m ncymm̂ mm̂ycn E ; (4.7c) and 0 = it D cym1c y ncm0 cn0 E + it D cym+1c y ncm0 cn0 E + it D cymc y n1cm0 cn0 E +it D cymc y n+1cm0 cn0 E it D cymc y ncm0 cn0+1 E it D cymc y ncm0 cn01 E it D cymc y ncm0+1cn0 E it D cymc y ncm01cn0 E +iV0 D cymc y ncm0 cn0 E m0+1;n0 + m01;n0 m+1;n m1;n (4.8a) iV0 X l=m1;n1 D cyl c y mc y nclcm0 cn0 E +iV0 X l=m01;n01 D cyl c y mc y nclcm0 cn0 E (4.8b) 2 X n m0 D cymc y ncn0 m̂ E n0 D cymc y ncm0 m̂ E +m D cyncm0 cn0 ̂m E n D cymcm0 cn0 ̂m E + n0 D ̂myc y mc y ncm0 E m0 D ̂myc y mc y ncn0 E + n D m̂yc y mcm0 cn0 E m D m̂yc y ncm0 cn0 Eo : (4.8c) Note that since the set of all polynomials of n cl; c y l o form a complete basis of the operator space, operators m̂ are certainly functions of polynomials of n cl; c y l o . Therefore, we expected G1 to depend on G2 from Eq(4.7b), and possibly also G3 or higher Green's functions from Eq(4.7c); while G2 is coupled to G3 from Eq(4.8b), and possibly also G4 or higher Green's functions from Eq(4.8c). Solving this full hierarchy is no easier than directly solving the Redeld equation, unless VS = 0 so that the set of equations for G1 is closed, i.e. not coupled to G2. We may, however, solve these equations by truncating the hierarchy at a certain level using further approximations, such as the molecular-chaos 58 4.3. Non-equilibrium Wick's theorem assumption in the classical Boltzmann equation [51], or replacing higher order Green's functions by a cluster expansion of lower order ones [69, 94, 95]. We use the following two approximate methods: (1) substitution of certain higher-order Green's functions by their values at equilibrium; (2) expressing higher-order Green's functions as combinations of lower-order ones plus a correlation part via cluster expansion, and then ignoring the correlation part at a certain order. A rough estimate of the accuracy of these substitutions at dierent orders can be found in Appendix C. Here we focus on the potential of this BBGKY-like formulation and discuss brie
y the performance of the two approximations. 4.3 Non-interacting systems: proof for a non-equilibrium Wick's theorem In this section, we will prove that when V0 = 0, G2 ky1; k y 2; k3; k4 = G1 ky1; k4 G1 ky2; k3 G1 ky1; k3 G1 ky2; k4 ; (4.9) where G1 and G2 are dened as G1 ky1; k2 = D cyk1ck2 E ; (4.10) respectively G2 ky1; k y 2; k3; k4 = D cyk1c y k2 ck3ck4 E : (4.11) In a general non-equilibrium state, one may have hck3ck4i 6= 0. In that case, the above Wick's theorem should have a more general form. In all of our examples, our choice of the specic system-baths coupling makes such Green's functions vanish. If there are terms involving two creation or two annihilation operators in the coupling between the central system and the baths, the situation is dierent. Here it is more convenient to work in the momentum representation than the position representation, therefore we dene: ck = 1p N NX l=1 sin kl N + 1 cl: (4.12) 59 4.3. Non-equilibrium Wick's theorem For a tight-binding chain with N sites and free ends (c0 = 0 = cN+1), the usual plane-waves of innite-size systems eikl are replaced by the eigenstates sin klN+1 , hence this transformation. In this representation, H0 = NX k=1 kc y kck; (4.13) where k = 2t cos k N + 1 : (4.14) Starting from Eq(4.4) with the above H0 in momentum space and using A = cyk1ck2 and A = c y k1 cyk2ck3ck4 , we nd the equations for G1 ky1; k2 and respectively G2 ky1; k y 2; k3; k4 as follows: 0 = i (k2 k1)G1 ky1; k2 2 2 N + 1 X sin k1l N + 1 sin k2l N + 1 (n (k1) + n (k2)) +2 2 N + 1 X ;k sin k2l N + 1 sin kl N + 1 G1 ky1; k +sin k1l N + 1 sin kl N + 1 G1 ky; k2 (4.15) 60 4.3. Non-equilibrium Wick's theorem and 0 = i (k4 + k3 k2 k1)G2 ky1; k y 2; k3; k4 +2 2 N + 1 X ;k sin k1l N + 1 sin kl N + 1 G2 ky; ky2; k3; k4 +2 2 N + 1 X ;k sin k2l N + 1 sin kl N + 1 G2 ky1; k y; k3; k4 +2 2 N + 1 X ;k sin k3l N + 1 sin kl N + 1 G2 ky1; k y 2; k; k4 +2 2 N + 1 X ;k sin k4l N + 1 sin kl N + 1 G2 ky1; k y 2; k3; k +2 2 N + 1 X sin k2l N + 1 sin k4l N + 1 G1 ky1; k3 (n (k2) + n (k4)) 2 2 N + 1 X sin k2l N + 1 sin k3l N + 1 G1 ky1; k4 (n (k2) + n (k3)) 2 2 N + 1 X sin k1l N + 1 sin k4l N + 1 G1 ky2; k3 (n (k1) + n (k4)) +2 2 N + 1 X sin k1l N + 1 sin k3l N + 1 G1 ky2; k4 (n (k1) + n (k3)) (4.16) Each of these forms a closed linear system of equations and has a unique solution. Therefore, we only need to nd one solution and that must be the unique solution. We rst apply Eq(4.9) to Eq(4.16) to expand G2 into products of G1. It is then easy to prove that the resulting equation is equiv- alent with Eq(4.15), meaning that a solution of Eq(4.15) is also a solution of Eq(4.16). For example, if we collect terms with G1 ky2; k4 together, we will have G1 ky2; k4 n i ( (k3) (k1))G1 ky1; k3 2 2 N + 1 X sin k1l N + 1 sin k3l N + 1 (n (k1) + n (k3)) +2 2 N + 1 X ;k sin k3l N + 1 sin kl N + 1 G1 ky1; k +sin k1l N + 1 sin kl N + 1 G1 ky; k3 ; (4.17) 61 4.4. Truncation of the hierarchy where the term in curly brackets is zero according to Eq(4.15). Therefore, the solutions of Eq(4.15) satisfy Eq(4.16) as well, as long as Eq(4.9), the non-equilibrium Wick Theorem, holds. Since the combined equation has a unique solution, the solution with Eq(4.9) valid is the only possible solution. Therefore, Eq(4.9) is proved. From this, we also know that for non-interacting systems, it is enough to solve only for the single-particle Green's functions, which form a linear system of equations with N2 unknowns. All higher order Green's functions can be factorized in products of various G1. 4.4 Truncation of the hierarchy As already shown, for interacting systems the equations for all Green's func- tions are coupled to one another. In order to solve Eq(4.7) explicitly, we rst have to nd explicit forms of the operators m̂ in terms of the operatorsn cl; c y l o . In Appendix A, we present an exact numerical calculation and a perturbational calculation for these operators. Correspondingly, based on these two methods of nding the operators m̂, we discuss in this section two ways of turning Eq(4.7) into a closed equation by approximating the G2 terms. This is the lowest meaningful level of truncation. We will rst discuss a more accurate method, which works for larger interactions V0, but which is computationally costly: perturbation based on two-particle Green's functions at equilibrium. Then we present a relatively less accurate but computationally much more ecient method, the non- equilibrium cluster expansion. The latter works only for relatively smaller V0 but can be applied to much larger systems. The values G1 my; n obtained from both methods will be compared against GEx1 my; n , the exact solution of Eq(4.2) via the linear system Eq. (2.89). Finally, we generalize the latter method by truncating it at the next level. This improves its accuracy and range of applicability, without sacricing its eciency. 4.4.1 Method 1: using equilibrium Green's functions As explicitly worked out in Eq(A.2) of Appendix A, the operators m̂ can be written in terms of eigenmodes of HS , which can be found from an exact diagonalization of HS . This is a 2 N -dimensional eigenvalue problem. Then, in the language of a super-operator space [88], where operators are treated like vectors { so called super-vectors { m̂ can be expanded under the basis 62 4.4. Truncation of the hierarchy of all polynomials of n cl; c y l o as: m̂ = X l d;lcl + V0D (4.18a) ̂m = X l d;lc y l + V0 D; (4.18b) where using the denition of inner product between super-vectors hhAjBii = tr AyB , we have d;l = 1 2(N1) tr cyl m̂ ; (4.19a) d;l = 1 2(N1) tr cl ̂m : (4.19b) Here 2N1 is a normalization constant to make d;l = 1 when m̂ = cl. Also, d;l and d;l are the expansion coecients at the linear order in cl and cyl ; operators V0D and V0 D are the remaining terms in the expansion of the operators m̂ and ̂m, respectively. V0 is explicitly factorized out because of the fact that when V0 = 0, this remaining part vanishes. With these expressions for m̂, Eq(4.7) becomes, 0 = it D cym1cn E + it D cym+1cn E it D cymcn+1 E it D cymcn1 E +2 X l; D n d;l + d ;l cymcl + m d;l + d ;l cyl cn E (4.20a) 2 X m d;n + n d ;m (4.20b) iV0 D cymc y n1cncn1 E + iV0 D cym+1c y mcm+1cn E iV0 D cyn+1c y mcn+1cn E + iV0 D cymc y m1cncm1 E (4.20c) 2V0 X D mcn D + n D y c y m ncymD mDycn E : (4.20d) Note that every c0; c y 0; cN+1 and c y N+1 that appears in the equation should be set to 0. First let us replace all G2s in Eq(4.20c) by their values at equi- librium, denoted here as G E;(0) 2 where the superscripts indicate the use of thermal equilibrium (denoted by superscript E) as the zeroth order approx- imation (denoted by (0)) of the non-equilibrium G2. Using the rst term as an example, G E;(0) 2 my; n = tr cymc y n1cncn1eq (HS) ; (4.21) 63 4.4. Truncation of the hierarchy where eq (HS) = 1 Z e HS kBT . This requires knowing the eigenstates of HS . Similarly one can dene G E;(0) D from Eq(4.20d), using the rst term as an example, G E;(0) D my; n = mtr cn Deq (HS) : (4.22) Next let us calculate G E;(1) 1 from Eq(4.20), where the superscript (1) means that the approximate calculation takes care of the rst equation of the hi- erarchy, Eq(4.20). We organize all G E;(1) 1 my; n as a vector, g E;(1) 1 = h G1 1y; 1 ; G1 1y; 2 ; ; G1 N y; N iT ; (4.23) and then Eq(4.20) for given values of m;n is the equation occupying the (mN + n)th row and in total there areN2 such equations. After substituting G E;(0) 2 and G E;(0) D for the exact but unknown G2 and GD, the whole set of Eq(4.20) for all m;n becomes a linear system for g E;(1) 1 with dimension N 2, (1)g E;(1) 1 = iV0g E;(0) 2 + 2 + 2V0g E;(0) D ; (4.24) where the vector comes from ordering Eq(4.20b) in the same way as g E;(1) 1 . The same holds for g E;(0) 2 and g E;(0) D correspondingly from ordering Eq(4.20c) and Eq(4.20d). The matrix (1) is extracted from Eq(4.20a). For example, assuming m and n are not at the boundaries, one may read from Eq(4.20), mN+n = X m d;n + n d ;m ; (4.25a) (1) mN+n;(m1)N+n = it: (4.25b) We can calculate single-particle equilibrium Green's functions, G E;(0) 1 , and organize them in the same way into a vector denoted as g E;(0) 1 . We dene the distance between two vectors A and B as, dAB = qP i jAi Bij2qP i jBij2 : (4.26) In order to gauge the accuracy, we compare dE;(0), the dierence between the zeroth order g E;(0) 1 and the exact solution g Ex 1 , and d E;(1), the dierence between the rst order solution above, g E;(1) 1 and the exact solution g Ex 1 . Here GEx my; n = tr cymcnex , where ex is the exact solution from Eq(4.2). 64 4.4. Truncation of the hierarchy Results First, we keep V0 = 0:2 as a constant, and check the accuracy of g E;(1) 1 for dierent values of T . From Fig.4.1(a) we can see that the worst is d(1) 1%. Secondly, we set T = 0:4T as a constant, and check the accuracy of g E;(1) 1 for dierent values of V0. The worst case is d (1) 0:3% as shown in Fig.4.1(b). Overall, the distance dE;(1) is always much smaller than dE;(0). We also compared the particle current calculated from this Green's functions, dened as: J = ie N 1 N1X l=1 (G (l + 1; l)G (l; l + 1)) : (4.27) JE;(0), JE;(1), JEx are calculated respectively from g E;(0) 1 , g E;(1) 1 and g Ex 1 . From Fig.4.1(c) and (d) we see that in both cases, JE;(1) is very close to the exact value, JEx, while JE;(0), the current in the equilibrium state, is always zero. Very high accuracy is found especially for small T . This indicates that the approximation captures the essential part of the non- equilibrium stationary states. It is also worth mentioning that this method generates reasonable results for very large V0. Furthermore, it is likely that the approximation could be further improved, by expansions to higher order polynomials of cl; c y l and substitution of their values at equilibrium for the higher order unknown Green's functions in the higher-order equations of the hierarchy. Stopping the expansion of operators m̂ at linear order of V0 is compatible with solving only the rst equation of the hierarchy. If further equations of the hierarchy are used then one should also expand operators m̂ to higher orders of V0. In order to estimate the accuracy of this level of approximation and also get a rough estimate of the accuracy of higher levels of truncation, we study the leading order of residues in terms of 2 and TT , both of which are assumed to be small in the following. Thus 2V0 V0, therefore we know that gD is smaller than the g2 term so we drop it. Similarly, since 2T T , we drop the 2T term in 2 in Eq(4.24), 2 = 20 (T ) + 2T;T ; (4.28) and keep only the large term, 20 (T ), which is independent of T . Here ;T denotes formally a derivative of T on | ddT . The general idea is then to write down equations for gEx1 and g E;(1) 1 (as shown in Eq. (C.2) of Appendix C), and then compare them to get an estimate for E;(1) 1 = g E;(1) 1 gEx1 . In 65 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 1 ∆T/T(V0=0.2) 0 0.02 0.04 0.06 0.08 d (a) (b) 0 0.5 1 1.5 2 V0(∆T/T=0.4) 0 0.005 0.01 0.015 0.02 dE,(0) dE,(1) 0 0.2 0.4 0.6 0.8 ∆T/T(V0=0.2) 0 0.001 0.002 0.003 0.004 J JE,(0) JEx JE,(1) 0 0.5 1 1.5 2 V0(∆T/T=0.4) 0 0.0005 0.001 0.0015 0.002(c) (d) Figure 4.1: g E;(1) 1 is compared with g Ex 1 for interacting systems at non- equilibrium. dE;(1) (dE;(0)) is the dierence between g E;(1) 1 (g E;(0) 1 ) and g Ex 1 dened using Eq. (4.26). Both dE;(1) and dE;(0) are plotted vs. respectively T in (a) and V0 in (b). In both cases, d E;(1) is much smaller than dE;(0). In (c) and (d) JE;(0) and JE;(1), plotted vs. respectively T and V0, are compared against JEx. We see that JE;(0) is zero while JE;(1) is close to JEx even for relatively large V0. In all these sample calculations, t = 1:0, = 0:1, = 1:0. 66 4.4. Truncation of the hierarchy order to understand how an approximation to the next order improves the accuracy, we also want to compare E;(1) 1 to E;(0) 1 = g E;(0) 1 gEx1 , which is estimated in the same way from the dierence between the equations respectively for gEx1 and g E;(0) 1 , see Appendix C for further details. Here we summarize the results, namely E;(0) 1 = T (1) 0 1 (1) ;T g Ex 1 + iV0 (1) 0 1 E;(0) 2 (4.29) and E;(1) 1 = (1) 0 1 V 20 (2)0 1E;(0)3 +iV0T (2) 0 1 (2) ;T g Ex 2 + iV0 2 (2) 0 1 E;(0) 1 : (4.30) Here E;(0) n = g E;(0) n gExn and E;(1)n = gE;(1)n gExn for general n-particle Green's functions gn. We refer readers to Appendix C for denitions of all matrices. Most importantly here we see that (from the last term in Eq. (4.30)) E;(0) 1 is multiplied by a small number 2V0 and then is included into E;(1) 1 . Furthermore, this relation holds generally for higher-order levels of this approximation. Judging from this it follows that, as long as 2V0 t the method is very reasonable. As for the other two terms, they can be regarded as V 20 g Ex 3 + V0g Ex 2 T . Therefore, the maximum value of V0 where this method is still accurate is determined by gEx2 1 or gEx3 12 . Such limit could be much larger than 1 since roughly gExn = gExn | smaller for larger n. This explains why this method is applicable even for V0 larger than t, as we see from Fig.4.1. Here, we note that since this approximation requires a full diagonaliza- tion of HS , it is necessarily not ecient for large systems. As a result, we now discuss another approach, that can be used for much larger systems. 4.4.2 Method 2: the non-equilibrium cluster expansion | rst level A dierent way to turn Eq(4.7) into a closed equation is to use the cluster expansion. Taking G2 as an example, this leads to: G2 my; ny;m 0 ; n 0 = G1 my;m 0 G1 ny; n 0 +G1 my; n 0 G1 ny;m 0 +G2 my; ny;m 0 ; n 0 : (4.31) 67 4.4. Truncation of the hierarchy If we set the part that describes correlations, G2 ! 0 (4.32) we have an approximative expression of G2 in terms of G1, which allows us to obtain a close set of equations for G1 (for more details, see below). This approach can be applied to higher-order Green's functions, for example using a similar expansion for G3 in terms of G1 and G2 and setting G3 = 0. In fact the non-equilibrium Wick's Theorem proved in Section x4.3 shows that indeed G2 = 0 when V0 = 0. This makes such expansions plausible for non-equilibrium Green's functions. Setting G2 = 0 with V0 6= 0 is similar to using the Hartree-Fock approximation. Depending on the system and the physical problem under investigation, one may need to go to the next level of approximation, i.e. keeping G2 but ignoring G3 to truncate the equation hierarchy at the level instead of the rst one. Here, we will study both levels, beginning rst with the lowest level approximation, i.e. setting G2 ! 0. However, the cluster expansion cannot be applied to the operators D dened in the previous section. Instead, we have to expand the operators m̂ in higher-order polynomials of n cl; c y l o . This can be done as follows. In order to avoid the costly exact diagonalization, the operators m̂ can also be found analytically using perturbation theory (see Appendix A for full details). The basic idea is to start by assuming cl (t) = c (0) l (t) + V0c (1) l (t) +O V 20 ; (4.33) and then derive and solve the equations of motion of c (0) l ; c (1) l from Heisen- berg's equation. In this way, one avoids the direct diagonalization of HS . This simplies the calculation but its accuracy depends on the order of V0 at which the expansion stops. Stopping at the linear order of V0 is compat- ible with the cluster expansion for G2 (rst level of approximation). If the cluster expansion in a higher-order Green's functions is applied, then the operators m̂ should also be expanded to higher orders of V0. Keeping only the rst order, the operators m̂ become m̂ = X m D;mcm + V0 X m1m2m3 D;m1m2m3cm1c y m2cm3 +O V 20 (4.34a) ̂m = X m D;mc y m V0 X m1m2m3 D;m1m2m3c y m3cm2c y m1 +O V 20 ; (4.34b) where the denitions of D;m and D;m1m2m3 are given in Appendix A. 68 4.4. Truncation of the hierarchy With the rst-order cluster expansion and the above expansion of oper- ators m̂ plugged into Eq(4.7), we have 0 = itG1 (m 1)y; n + itG1 (m+ 1)y; n itG1 my; n+ 1 itG1 my; n 1 +2 X l; h n D;l + D ;l G1 my; l +m D;l +D ;l G1 ly; n i (4.35a) +2V0 X ;m1;m2 (D;nm2m1 D;m1m2n)G1 my1;m2 m +2V0 X ;m1;m2 (D;mm2m1 D;m1m2m)G1 my2;m1 n (4.35b) 2 X m D;n + n D ;m (4.35c) +2V0 X ;m1 (D;m1m1nm +D;m1m1mn) (4.35d) iV0G my; (n 1)y; n; n 1 +iV0G (m+ 1)y;my; (m+ 1); n iV0G (n+ 1)y;my; n+ 1; n +iV0G my; (m 1)y; n;m 1 : (4.35e) Here G2 should be interpreted according to Eq. (4.31) with G2 ! 0. This is the closed system of equations for the unknowns G1. Next we dene a vector g C;(1) 1 , where as before superscript C means cluster expansion and (1) symbolizes keeping only the rst equation in the hierarchy. For simplicity we order Eq(4.35e) in the same way and denote it as g C;(1) 2 = g C;(1) 1 , where refers to the nonlinear function | summation of products | of g C;(1) 1 in Eq(4.35e). Then the above equation can be written in matrix form as: (1) 0 + 2V0 (1) D g C;(1) 1 = 20 + 2V01 + iV0g C;(1) 2 ; (4.36) where the ve terms are respectively the ve sub equations in Eq(4.35), for 69 4.4. Truncation of the hierarchy example, (0)mN+n = X m D;n + n D ;m : (4.37) This equation can be solved iteratively g [n+1] 1 = (1) 0 + 2V0 (1) D 1 20 + 2V01 + iV0 g [n] 1 ; (4.38) starting from the initial value g [0] 1 = (1) 0 1 20: (4.39) Here we start from g [0] 1 = g C;(0) 1 , which is the exact solution of Eq(4.36) when V0 = 0. Through the iteration dened above we get the solution g C;(1) 1 = limn!1 g [n] 1 . In practice we stop at large enough n such that g [n] 1 g[n1]1 is small enough. Results Using again Eq. (4.26), we dene dC;(0) as the relative distance between g C;(0) 1 and g Ex 1 , and d C;(1) as the relative distance between g C;(1) 1 and g Ex 1 . First, we keep V0 = 0:2 as a constant, and check the accuracy of g C;(1) 1 for dierent values of T . From Fig.4.2(a) we see that the worst case is d(1) 1%. Secondly, we set T = 0:4T as a constant, and check the accuracy of g C;(1) 1 for dierent values of V0. The worst case is d (1) 2% as shown in Fig.4.2(b). Overall, the dierence dC;(1) is always much smaller than dC;(0). From Fig.4.2(c) we see that for a small V0, J C;(0) already provides most of the total current. However, Fig.4.2(d) shows that as V0 increases, the dierence between JC;(1) and JC;(0) grows. We should also note that for larger V0, JC;(1) starts to deviate from JEx. This indicates that the approximation captures the essential part of the interaction but it is qunatitatively accurate only for small V0. Of course, this can be improved by going to the next level in the trucation, as shown in section 4.4.3. In order to estimate the accuracy of this approximation, let us assume that 2 and V0 are small. We dene C;(0) n = g C;(0) n gExn and C;(1)n = g C;(1) n gExn . Again we start from the equations of the three: gC;(0)1 , gC;(1)1 and gEx1 (as we have done in Eq. (C.14) of Appendix C), and then compare them while ignoring certain higher-order, smaller terms such as those which 70 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 ∆T/T(V0=0.2) 0 0.01 0.02 0.03 0.04 d (a) (b) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0 0.03 0.06 0.09 0.12 0.15 dC,(0) dC,(1) 0 0.2 0.4 0.6 0.8 ∆T/T(V0=0.2) 0 0.001 0.002 0.003 0.004 J JC,(0) JEx JC,(1) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0.001 0.0012 0.0014 0.0016 0.0018 (c) (d) Figure 4.2: g C;(1) 1 is compared with g Ex 1 for non-equilibrium interacting systems. Both dC;(1) and dC;(0) are plotted vs. respectively T in (a) and V0 in (b). In both cases, d C;(1) is always much smaller than dC;(0). In (c) and (d) JC;(0) and JC;(1), plotted vs. respectively T and V0, are compared against JEx. From (c), where V0 = 0:2 and it is relatively small, we see that for a given value of V0, both J C;(0) and JC;(1) are consistent with JEx. From (d) we nd that for relatively larger V0, J C;(1) agrees much better with JEx than JC;(0). 71 4.4. Truncation of the hierarchy are proportional to 2V0, see Appendix C for more details. We arrive at: C;(0) 1 = iV0 (1) 0 1 gEx2 V0 gEx1 2 ; (4.40) and C;(1) 1 = (1) 0 1 iV 20 (2) 0 1 gEx3 +2V0 (2) 0 1 C;(0) 1 V 20 gEx1 3 + 2V 20 gEx1 2 : (4.41) We refer the readers to Appendix C for denitions of all matrices. This estimate agrees with the numerical results that C;(0) 1 is propor- tional to V0 while C;(1) 1 is proportional to V 2 0 . Most importantly, we see again that C;(0) 1 is multiplied by a small number 2V0 and then be- comes a part of C;(1) 1 . Since roughly gExn = gExn, the other term, V 20 g Ex 3 V 20 gEx1 3, is also much smaller than C;(0)1 V0 gEx1 2. However, for large enough V0 the other approximation used in this method, the per- turbational expansion of the operators m̂, becomes invalid. Therefore, as long as V0 t this level of approximation is expected to be very reasonable. It should be noted that this method is capable of dealing with large systems since it does not require a direct diagonalization of a 2N -dimension matrix, like the previous method. Instead, it deals with vectors of dimension N2. 4.4.3 Method 2: second-level cluster expansion In the previous section, we showed that even truncation at rst level leads to quite accurate results, and a very ecient method. Here, we discuss the truncation at the second level of the BBGKY hierarchy. This requires us to keep G2 as unknowns in the equations which now include the equations for G2, while letting G3 = 0 in order to truncate the resulting larger sys- tem. While the rst-level approximation is equivalent to the Hartree-Fock approximation, this second-level form goes beyond that. The cluster expan- sion of G3 [94] expresses it in terms of G1, G2 and G3 (see Appendix B for the explicit expression): G my1;m y 2;m y 3;m4;m5;m6 = X P (1)PG1G1G1 + X P (1)PG1G2 +G3 my1;m y 2;m y 3;m4;m5;m6 : (4.42) 72 4.4. Truncation of the hierarchy Here P P (1)PG1G1G1 is a short-hand notation for the various ways of pairing fmy1;my2;my3;m4;m5;m6g into three groups using G1 and taking care of the anti-commutation relations, for example G1 my1;m6 G1 my2;m5 G1 my3;m4 : (4.43) Similarly P P (1)PG1G2 denotes all dierent ways of pairing these into two groups using G1 and G2, for example G1 my1;m6 G my2;m y 3;m4;m5 : (4.44) Setting G3 ! 0, we aim to derive a close system of equations for G1 and G2 from Eq. (4.7) and Eq. (4.8). In order to be consistent, i.e. so that all terms lower than G3 should be included in both equations, bath operators m̂ should be truncated to second order, i.e. at terms proportional to V 20 . However, including the V 2 0 terms greatly complicates both equations (see Appendix B for details). Therefore, here we use only the expression of operators m̂ truncated at the rst order of V0. In this case, Eq. (4.7) becomes an equation very similar to Eq. (4.35) but with a dierent meaning for G2: G2 should be interpreted according to Eq. (4.31) with explicit G2 included. We denote this, using compressed notations similar to the ones in Eq. (4.24) and Eq. (4.36), as (1)g1 = 20 + 2V01 + iV0 (g1g1) + iV0g2: (4.45) Note that now we have an additional term iV0g2, so this is not a closed equation by itself. Here (g1g1) is the short-hand notation for the combi- nations of products between g1 and g1 the cluster expansions in Eq. (4.31). Meanwhile, Eq. (4.8) becomes, (2)g2 = (2) (g1g1) + 2g1 + 2V0g1 +iV0 (g1g1g1) + iV0 (g1g2) ; (4.46) such that we get a closed system of equations when we combine the two sets. The explicit form of this latter equation can be found in Appendix B. Here (g1g2) is the short-hand notation for the combinations of products between g1 and g2 from Eq. (4.42). Same holds for (g1g1g1). The matrices (1) and (2) come from the following terms (1) = (1) H0 + 2 (1) B0 + 2V0 (1) B1 ; (4.47a) (2) = (2) H0 + iV0 (2) V + 2 (2) B0 + 2V0 (2) B1 : (4.47b) 73 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 1 ∆T/T(V0=0.8) 0 0.05 0.1 d (a) (b) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0 0.05 0.1 0.15 0.2 d1 (0) d1 (1) d1 (2) Figure 4.3: G1 calculated with dierent levels of the approximations are compared. d (2) 1 is plotted in addition to the previous shown d (1) 1 and d (0) 1 . We see that the distance d (2) 1 is smaller than, but rather close to d (1) 1 , while both are much smaller than d (0) 1 . The various s are obtained by considering only the term indicated by their subscripts. For example, (1) H0 is the one derived from H0 only, without using the potential V and bath operators m̂, while 2 (1) B0 (2 (1) B1 ) is related to the component of the bath operators that is of the zeroth (rst) order in V0. The combination of Eq. (4.45) and Eq. (4.46) is a closed system for g1 and g2. Solving it gives us results of the second-level cluster expansion, and provides nite values for G2. In other words, explicit correlations are built in, and thus this method goes well beyond the Hartree-Fock approximation. In the following we are investigate how large are these G2 and also what is the gain in accuracy when compared against exact numerical solutions, and against the rst-level truncation. First, we dene d (2) 1 to be the dierence between the G1 calculated from this second-order cluster expansion and the exact G1, similarly as d C;(1), which is now denoted as d (1) 1 . In Fig.4.3 we nd that d (2) 1 is smaller than but very close to d (1) 1 . This means that the second-order cluster expan- sion has better accuracy, but quantitatively does not lead to a signicant improvement of G1. 74 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 1 ∆T/T(V0=0.8) 0 0.05 0.1 0.15 0.2 0.25 0.3 d (a) (b) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0 0.1 0.2 0.3 0.4 d2 (0) d2 (1) d2 (2) Figure 4.4: G2 calculated with dierent levels of the approximations are compared. d (2) 2 is plotted v.s. T and V0. We see that d (2) 2 is signicantly smaller than d (1) 2 . Second, in the same way, we dene d (2) 2 and compare the dierences be- tween G2s calculated from the various truncations to the exact G2. From Fig.4.4 we nd that the second order truncation leads to a sizable improve- ment in G2. To nd how signicant is this dierence for physical quantities, we com- pare the currents calculated from the approximations and exact calculation. In order to enhance these dierences, we increase the coupling constant to = 0:5 in the following gures. This is because the absolute value of the particle current depends very strongly on the value of (J / 2). For the same purpose we have also set V = 0:8 instead of V = 0:4 when checking the T dependence. Figure 4.5 shows that the agreement between J (2) and JEx is much better than the one between J (1) and JEx, especially for larger V0. This proves that the correlations kept at the second-level truncation are important for an accurate description of the transport properties. Since the second-level approximation allows us to calculate G2 | the correlated part of G2 { we also consider the heat current, a physical quantity 75 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 1 ∆T/T(V0=0.8) 0 0.02 0.04 0.06 0.08 J JEx J(0) J(1) J(2) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0.02 0.025 0.03 (a) (b) Figure 4.5: Particle currents, J (2); J (1); J (0), calculated respectively from the second-, rst- and the zeroth-order truncations are compared against JEx, which is found from the exact diagonalization method. For all the range of V0 shown in this gure, J (2) is very close to JEx, while J (1) is accurate only for small V0. 76 4.4. Truncation of the hierarchy 0 0.2 0.4 0.6 0.8 1 ∆T/T(V0=0.8) 0 0.005 0.01 0.015 0.02 0.025 J Jh Ex Jh (0) Jh (1) Jh (2) 0 0.2 0.4 0.6 0.8 1 V0(∆T/T=0.4) 0.004 0.005 0.006 0.007 0.008 0.009 0.01 (a) (b) Figure 4.6: Heat currents, J (2) h ; J (1) h ; J (0) h are compared against J Ex h . For all range of V0 shown in this gure, J (2) h is close to J Ex h . Jh depends strongly on G2 while particle current J in the previous gure relies on G1. making explicit use of G2: hJhi = 2 N 2 N2X l=2 Im t2G1 (l 1; l + 1) +V tG2 (l 1; l + 1; l; l + 1) + V tG2 (l 1; l + 1; l; l 1)] ; (4.48) which is derived from Jh = N2X l=2 i hBl1; h B l (4.49) and hBl = tcyl cl+1 tcyl+1cl + V cyl clcyl+1cl+1: (4.50) Again, we nd from Fig.4.6 that the agreement between J (2) h and J Ex h is much better than the one between J (1) h and J Ex h , especially for larger V0. This indicates that the nite G2 are essential for getting proper values of heat currents. 77 4.5. Conclusion and discussion Keeping orders of V 20 in the bath operators m̂ is consistent with trun- cation at G3 = 0 and may improve accuracy for the second-level truncation even further. However, including these terms greatly complicates the linear equation for G2, so that it has about 200 dierent terms (see Appendix A for more details). An investigation of their eect on the overall results is left as future work. 4.5 Conclusion and discussion To conclude, a BBGKY-like equation hierarchy is derived from the Redeld equation and two systematic approximations are suggested to truncate and solve the hierarchy. Using the rst-level of the rst method and both the rst- and the second-levels of the second method, non-equilibrium station- ary states of interacting systems are calculated. It is found these are in good agreement with exact results, for small systems where the latter are available. We also estimated the accuracy of the two approximations. We nd that the rst method works for strong interaction V0. However, because it requires an explicit diagonalization of the Hamiltonian, it is less computationally ecient. We imagine this method being used to study, for example, systems of a few strongly coupled quantum dots which have large charging energies. The second method is much more ecient and therefore it can be ap- plied to much larger systems, but its accuracy decreases with increasing interaction strength V0. However, truncation at a higher level is shown to signicantly improve both its accuracy, and its range of validity. Even at the second level, good accuracy for physical quantities is seen up to rather high values of V0 and T , and the results have the correct trends even when they become quantitatively less accurate. It is likely that going to even higher levels of truncation further improves accuracy, although, of course, this also leads to signicantly more involved computations. In future work, one can investigate the generalization to higher level trun- cations. Even more interesting, however, is the application of this method to physical problems. 78 Chapter 5 Solving the Redeld equation using coherent-state representations 5.1 Introduction In quantum optics, a common technique to solve quantum master equations for density matrices is to use the coherent-state representation [70]. Quite of- ten a quantum master equation for a single photon (boson) mode is solved in the coherent-state representation, where the operator-form quantum mas- ter equation becomes a dierential equation of c-numbers and a density operator becomes a \distribution" function. The problem of the Hilbert space's innite dimension is bypassed. For example, a distribution function over N complex numbers is sucient to describe a N -site Bose-Hubbard model [96, 97]. The equations for these distribution functions can be solved analytically or through numerical simulations of certain equivalent stochas- tic dierential equations such as Langevin equations (see Refs. [70, 98] and also Appendix E for the relation between generalized Fokker-Planck equa- tions and Langevin equations). A common coherent-state representation is the P -representation (see Ref [70] and also the Appendix D for details on various forms of coherent- state representations). In the P -representation coherent eigenvalues and are treated as conjugate complex numbers, thus one solves only the equa- tions for . Sometimes it is also useful to work with generalizations of the P -representation [70]. In fact, the generalized P -representation, where co- herent eigenvalues and are treated independently, has been used in sim- ulations of both thermal equilibrium and dynamical evolution of quantum systems (see for example [72] and references therein). It has been shown that this approach is capable of dealing with the pure dynamical evolution of an interacting quantum system with 150000 atoms and 106 momentum modes [75]. The work presented in this chapter can be regarded as an ex- 79 5.1. Introduction tension of this method from simulations of equilibrium and pure dynamical states to simulations of non-equilibrium stationary states. Thus, in this chapter we introduce another ecient method to solve the Redeld equation: using the coherent-state representation, the Redeld equation is mapped onto a generalized Fokker-Planck equation. For non- interacting systems, analytical results are obtained for the non-equilibrium stationary states. For interacting systems, we nd that the resulting equa- tions can be solved eciently via classical simulations. The existence of a general distribution function for systems in thermal equilibrium, which depends only on the Hamiltonian and the temperature, greatly simplies the calculation of equilibrium properties. However, there is no known general distribution for non-equilibrium stationary states. For each problem, therefore, one has to calculate the distribution from a quan- tum master equation, such as the Redeld equation, which is capable of taking into account the role of the baths, besides the central system. However, we have seen in the previous chapters that it is very hard to solve the Redeld equation eciently. The fact that there is no analytical expression for non-equilibrium stationary states may simply be due to the technical diculty in solving the quantum master equation. A more e- cient way of solving this equation might lead to dierent results, or maybe to a proof that there is no general expression describing non-equilibrium stationary states. In this chapter, using the coherent-state representation, we derive an analytical expression for non-equilibrium stationary states of non-interacting systems. The question of nding such non-equilibrium sta- tionary states for non-interacting systems has been numerically investigated via the time-dependent density matrix renormalization [99, 100], however for the local-operator Lindblad equation instead of the Redeld equation. The focus of this chapter, however, is to study the eciency of this approach when applied to interacting systems. The chapter is organized as follows. In section 5.2, we introduce the model system, dene its Redeld equation and then convert the operator- form Redeld equation into a dierential equation in the coherent-state rep- resentation. We nd the mapped equation to be in the form of a generalized Fokker-Planck equation. We then proceed to nd analytical solutions of the generalized Fokker-Planck equation for non-interacting systems in Section 5.3. In Section 5.4, we show that the same BBGKY-like hierarchy can be derived from the generalized Fokker-Planck equation as from the original Redeld equation. This provides further support for the work presented in the previous chapter. In section 5.5, we solve the generalized Fokker- Planck equation via numerical simulations for interacting bosonic systems. 80 5.2. The model and its eective equation of motion The method is also applicable to other quantum master equations such as the local-operator Lindblad equation, as we demonstrate in section 5.6. In section 5.7, we conclude that this new method based on the coherent-state representation leads to analytical solutions for non-interacting systems and comment on its possible further development and applications. 5.2 The model and its eective equation of motion For concreteness, the system of interest in this chapter is a 1-D Bose- Hubbard chain [96, 97] coupled to a bath at each of its ends. The left (right) bath is maintained at temperature TL = T + T 2 (TR = T T2 ). We are interested in the stationary state of the central system. The Hamiltonian of the central system (HS), the baths (HB) and the couplings (HSB) are, respectively: HS = H0 + VS = N1X l=1 ayl al+1 + a y l+1al + U 2 NX l=1 nl (nl 1) ; (5.1a) HB = X k; !k; byk;bk; + 1 ; (5.1b) HSB = X k;=L;R V k a y bk; + h:c: : (5.1c) Note that only the ends are coupled to baths ( = L;R): aL = a1; aR = aN etc. The system of interest, the Bose-Hubbard chain, has been realized experimentally [97]. Neglecting the shift of the spectrum of the central system due to the coupling to the baths, the Redeld equation reads, from Eq. (2.48): @(t) @t = i[HS ; (t)] 2 X =L;R nh ay; m̂(t) i + a; ̂m(t) + h:c: o ; (5.2) 81 5.2. The model and its eective equation of motion where m̂L(m̂R) is related to a1(aN ) while ̂mL( ̂mR) is related to a y 1(a y N ), m̂ = X k jV k j2 Z 1 0 da () ei!k; h1 + n (!k;)i; (5.3a) ̂m = X k jV k j2 Z 1 0 day () ei!k; hn (!k;)i: (5.3b) Here hn (!k;)i = e !k; T + 1 1 is the Bose-Einstein distribution for the bath modes. Note that later we will switch the summation over momentum k into an integral over energy , which involves a density of states D (). Generally both D () and jV k j2 may depend on . However, for simplicity we set the density of state and also jV k j2 to constants and absorb them into the coupling constant . Next we map the Redeld equation into a dierential equation in the coherent-state representation. A similar approach has been used in earlier works mentioned in Ref. [72]. However, they are based on dierent equa- tions of motion, such as the pure dynamical equation or the local-operator Lindblad equation [90] with bath temperature set at zero, while we use the Redeld equation for nite bath temperatures. From a technical point of view, in the coherent-state representation it is easier to deal with the local- operator Lindblad equation than the Redeld equation. However, we have shown in section 2.5 of chapter 2 that in order to get proper non-equilibrium stationary states, it is necessary to use the more complicated Redeld equa- tion. As we demonstrated in Appendix A, for fermions the operators m̂ and ̂m can be derived perturbationally to avoid direct diagonalization of HS . Here we apply the same procedure for bosons. Assuming that the Hubbard U is small, we may solve the Heisenberg equation order by order to get al (t) = X n Una (n) l (t) ; a y l (t) = X n Una y;(n) l (t) ; (5.4) and then carry out the integrals to get the operators m̂ and ̂m. An important observation is that when U = 0 the problem is exactly solvable and al (t) (a y l (t)) is a linear combination of all am (a y m), al (t) = a (0) l (t) = 2 N + 1 X km sin kl N + 1 sin km N + 1 eiktam; (5.5) 82 5.2. The model and its eective equation of motion where k = 2 cos k N+1 . This is because for nite-size systems with free ends, the Fourier transform is discrete and eikl is replaced by the sin klN+1 wave- function, see for example Eq. (A.4). For interacting systems, in principle it is possible to perform this pertur- bational calculation up to a large n, however here we will stop at the rst order, a (1) l (t), a (1) l (t) = i Z t 0 d a y;(0) l () a (0) l () a (0) l () : (5.6) After some tedious but straightforward algebra( see Ref. [101] and also Ap- pendix A for details), we nd, to the linear order in U , that: m̂ = X m D;mam + U X m1m2m3 D;m1m2m3a y m1am2am3 +O U2 (5.7a) ̂m = X m D;ma y m + U X m1m2m3 D;m1m2m3a y m3a y m2am1 +O U2 ; (5.7b) where D;m = 2 N + 1 X k sin kl N + 1 sin km N + 1 [1 + n (k; T; )] ; (5.8a) D;m = 2 N + 1 X k sin kl N + 1 sin km N + 1 n (k; T; ) ; (5.8b) D;m1m2m3 = 2 N + 1 4 X k;m;k1;k2;k2 n ( (k2) + (k3) (k1) ; T; ) n ( (k) ; T; ) (k2) + (k3) (k1) (k) sin kl N + 1 sin km N + 1 sin k1m1 N + 1 sin k2m2 N + 1 sin k3m3 N + 1 sin k1m N + 1 sin k2m N + 1 sin k3m N + 1 : (5.8c) Here n (; T; ) is the average Bose-Einstein occupation number of a state with energy , temperature T and chemical potential . Next, as summarized in Appendix D, the creation and annihilation oper- ators ayl and al become dierential operators in the coherent representation, for example: ayl $ l @ @l P ~ : (5.9) 83 5.2. The model and its eective equation of motion Using the above approximation for the operators m̂ and the transformation from creation/annihilation operators to dierential operators, Eq.(5.2) be- comes, in the coherent P -representation, a generalized Fokker-Planck equa- tion, @P @t = LH0 + ULV + 2L (0) B + 2UL (1) B P: (5.10) where LH0 = N1X l=1 i @ @l+1 l + i @ @l l+1 + c:c: ; (5.11a) LV = NX l=1 i @ 2 @2l 2l 2 + i @ @l l 2 l + c:c: ; (5.11b) L (0) B = X ;m D;m D;m @ @ m + D;m @2 @@m + c:c: ; (5.11c) L (1) B = X ;mi D;m1m2m3 @2 @@m2 m1m3 + @2 @@m3 m1m2 @2 @@m1 m2m3 + c:c: (5.11d) @ 3 @@m2@m3 m1 + c:c: : (5.11e) Here c:c: refers to the complex conjugate and it is performed formally, e.g. i @@l+1 l is the c:c: of i @ @l+1 l. The last term, Eq. (5.11e), of L (1) B is separated from the other terms intentionally for reasons that will become apparent soon. If only equilibrium states are of interest, then the rst and second term can be discarded since equilibrium states commute with HS . However, for non-equilibrium stationary states those terms are important. Setting @@tP = 0, we get the stationary state P ~;1 , which in the following is denoted simply as P ~ . Next, we will calculate the stationary of Eq. (5.10). We rst discuss non-interacting systems and then turn to interacting ones. 84 5.3. Exact solution for non-interacting systems 5.3 Exact solution for non-interacting systems For U = 0 only the terms LH0 and L (0) B contribute, and the stochastic equation is a Fokker-Plank equation [98], in terms of complex variables: @ @t P = X i;j " @ @i ijj + @2 @i@j Dij + c:c: # P: (5.12) For our specic case of a non-interacting system we have: lm = i (l;m1 + l;m+1) X ;m 2lm; (5.13) and Dlm = 1 2 2 X D;lm + D;ml : (5.14) Note that we have veried that Dy = D . This is in fact a standard Fokker-Planck equation (FPE) for the Ornstein- Uhlenbeck process. Such a FPE can be solved analytically [70, 98]. Its stationary solution is found to be a Gaussian function: P ~ = 1 Z e 1 2 y1; (5.15) where Z is a normalization constant and + y +D = 0; (5.16) which has the solution, = X l;m 1 l + m hl jD j mi jlihmj : (5.17) Here ; ji ; hj are respectively the eigenvalues, right eigenvectors and left eigenvectors of and hj ; ji are the complex conjugate of the two cor- responding vectors. This expression is correct for both equilibrium and non-equilibrium states. The only dierence is in the matrix D, in whether TL is the same as, or is dierent from TR. The only needed numerical task is to nd an eigenvector decomposition of , an N -dimensional matrix. Once is known, all correlation functions can be calculated, in particular the single-particle Green's function is G ly;m = hl mi = 2lm. 85 5.3. Exact solution for non-interacting systems Next we work out explicitly the distribution function for a two-site (N = 2) sample system to rstly conrm that it leads to the proper ther- mal equilibrium, and secondly to show the general procedure to nd such a distribution function. It should be emphasized that the dimension of the space grows linearly with system size, so we chose to work with N = 2 just because it is easier to get some intuitive picture, not because of algo- rithm complexity issues. In principle, for two-site systems we need to work with 4-dimensional matrices, but for non-interacting systems they reduce to 2-dimensional matrices. We have for this specic case, = 2 i i 2 ;D = D11 D12 D12 D22 ; (5.18) where D11 = 2 n(1; TL; L) + n(1; TL; L) 2 ; D12 = 2 n(1; TL; L) n(1; TL; L) + n(1; TR; R) n(1; TR; R) 4 ; D22 = 2 n(1; TR; R) + n(1; TR; R) 2 : (5.19) >From this we get the stationary distribution, = 1 42 D11 +D22 2D12 2D12 D11 +D22 + D11 D22 4 (42 + 1) 2 i i 2 : (5.20) For equilibrium states TL = TR = T , L = R = and therefore D11 = D22, so that P ~ / e 1n(1;T;) 1+2p 2 1+2p 2 1 n(1;T;) 12p 2 12p 2 ; (5.21) which is exactly the P -representation form of the Boltzmann distribution in terms of the eigenmodes a1a2p 2 with eigenenergies 1. Note that this is independent of the coupling constant since the second term in vanishes and the 2 in the rst term cancel out. This conrms that the equilibrium state of the central system does not depend on the details of the coupling between the system and the baths. However, this is unique for equilibrium states. A general NESS does depend on . In fact, the non-equilibrium part D1;1D2;2 / 2. In the following calculation, a relatively large is used in order to enhance this non-equilibrium part. 86 5.4. Interacting systems: BBGKY-like hierarchy derived from the GFPE Next we compare correlations from the non-equilibrium distribution with correlations calculated with the BBGKY-like method discussed in the pre- vious chapter [101]. Numerical exact solution for non-interacting systems have been reported in Refs. for the Redeld equation. For non-interacting systems, the BBGKY method can be solved exactly. Here we use it as a reference to check the analytical solution for non-interacting systems. We denote the single-particle correlation functions calculated from the above analytical solution of generalized Fokker-Planck equation, respectively from the numerical solution of the BBGKY-like method as Gre;An1 and G re;B 1 . Similar notations are used for particle currents Jre;An and Jre;B. We mea- sure the dierence between the two correlations, for example the analytical solution and the BBGKY solution, through, dre;AnB = rP l;m Gre;An1 (ly;m)Gre;B1 (ly;m)2rP l;m Gre;B1 (ly;m)2 : (5.22) In Fig.5.1, we compare this analytical solution and the BBGKY-like nu- merical results for the Redeld equation. Since dre 0, the two methods are indeed equivalent. The inset shows that currents calculated from both meth- ods also agree. This validates the fact that the generalized Fokker-Planck equation approach produces accurate results at least for non-interacting sys- tems. In the same gure, we also plot dre;nuB and J re;nu, which are obtained from a numerical solution that will be discussed in detail in section 5.4. There is a small but nite distance between the two correlation functions, however this has no contribution to the current, which is found to be in good agreement with the exact solution. 5.4 Interacting systems: BBGKY-like hierarchy derived from the generalized Fokker-Planck equation For central systems with interactions, Eq.(5.10) is a generalized Fokker- Planck equation which has additional triple derivative terms added to the standard Fokker-Planck equation. It is still possible to solve directly this generalized Fokker-Planck equation via a set of equivalent stochastic dier- ence equations, which are generalizations of the Langevin equation [98, 102]. We will discuss this technique in the next section x5.5. Before doing that, 87 5.4. Interacting systems: BBGKY-like hierarchy derived from the GFPE 0 0.2 0.4 0.6 0.8 1 ∆T/T 0 0.1 0.2 d dre,AnB dre,NuB 0 0.4 0.8 ∆T/T 0 0.3 0.6 J Jre,B Jre,An JRe,Nu Figure 5.1: Analytical solution of the generalized Fokker-Planck equation is compared against the BBGKY numerical solution for non-interacting sys- tems. The red (circles) curve in the main plot shows the comparison between the above analytical solution and the BBGKY-like numerical solution of the Redeld equation. Both Green's functions (dre;AnB , main panel) and particle currents (Jre;An and Jre;B, inset) exactly agree with each other. The blue (square) curve is a comparison between a numerical simulation, to be dis- cussed in x5.4, and the BBGKY solution. Here we nd that the numerical simulation agrees with analytical solution quite well. Other parameters are = 0:5; T = 2:0; = 2:0. 88 5.4. Interacting systems: BBGKY-like hierarchy derived from the GFPE we verify here that the same equation hierarchy as the one presented in the previous chapter can also be derived from this generalized Fokker-Planck equation. The derivation of the BBGKY equations for bosonic systems mir- rors that of the previous chapter, and can be found in Appendix F. We focus again only on the single-particle stationary state Green's func- tions G ly;n , dened as hl ni = R d~l nP ~ . Let us then derive a set of equations for them from the generalized Fokker-Planck equation, Eq. (5.10). For example, consider a rst-order derivative term in Eq. (5.11b). After performing a partial integration, we have,X m Z d~l ni @ @m m 2 mP ~;1 = i X m Z d~P ~;1 m 2 m @ @m (l n) = i Z d~P ~;1 l n 2 n = i l n2n : (5.23) Therefore, setting @@thl ni to zero, we have 0 = ihl n+1i ihl n1i+ ihl+1ni+ ihl1ni (5.24a) +2 X ;m D;m D;m [nhl mi+ lhmni] (5.24b) +2 X D;nl + D;ln (5.24c) +iUhl l lni iUhl nnni (5.24d) +2U X ;mi D;m1m2m3 n hm1m2ilm3 + hm1m3ilm2 +l hm2m1inm3 + hm3m1inm2 : (5.24e) This is exactly the same with Eq. (4.7). This conrms that the BBGKY- like hierarchy method and the generalized Fokker-Planck equation in the coherent-state representation are equivalent, as expected. Next, we discuss how to solve numerically the generalized Fokker-Planck equation. Results from the BBGKY-like method will be used as a reference. 89 5.5. Approximate numerical solution for interacting systems 5.5 Approximate numerical solution for interacting systems The generalized Fokker-Planck equation, Eq.(5.10), has triple derivative terms. The Pawula theorem [98] states that a distribution function with negative values will arise from a generalized Fokker-Planck equation with a nite series of derivatives beyond the second order. This does not, how- ever, rule out the possibility that a generalized Fokker-Planck equation with an innite series of derivatives may lead to a positively valued distribution function. Our Redeld equation is, in principle, such a generalized Fokker- Planck equation with innite derivatives, if all terms in the expansion of the m̂ operators are kept. A truncation at a certain nite order might provide a more accurate result than a truncation at the second order, but the resulting distribution function will have negative values in some small regions. Ref. [98] presented such an example. Of course, new techniques are required to solve such generalized Fokker-Planck equation. Besides the example from Ref. [98], Ref. [102] proposed a set of equivalent stochastic dierence equations for a generalized Fokker-Planck equation with triple derivative terms. The authors showed that a better distribution function than the one obtained after truncating at second derivatives, is given by their technique, which is based on a specic kind of stochastic dierence equation [102]. However, this technique is still under development and discussion. One may totally avoid such triple derivative terms by using generalized Gaussian phase-space method [72] or applying the phase-space method to the whole composite system [103]. Especially with the latter approach, one can simulate the dynamics of the whole system instead of that of the central system only. Of course, in this case, the baths will be treated as large but nite systems. Another benet of this approach is that one could also go beyond Markov approximation. On the other hand, one should also note the possibility that drop- ping such triple derivative terms may be appropriate in certain circum- stances [104]. For example, in Eq. (5.24) | the rst-order BBGKY equa- tion derived above for two-variable (corresponding to single-particle) corre- lations, we see that these third-order terms are not in the equation explicitly. They appear only in equations of three-or-more-variable correlations, for ex- ample in the equation of four-variable correlations, which itself appears in the equation of two-variable correlation. Therefore such third order terms aect two-variable correlations, but indirectly via four-variable correlations. Following this line of logic we have assumed that odd-number-variable cor- 90 5.5. Approximate numerical solution for interacting systems relations vanish. Since the hierarchical cluster expansion method [101] generally relates higher order correlations to lower ones, if one truncates it at a certain order as we do here, one captures certain aspects of the interaction although not its whole eects. This is better than not taking the interaction in consideration at all. In this sense, neglecting the third order terms is equivalent to the rst order of the cluster expansion, where the four-variable correlations are approximated by products of two-variable ones such that eects of the third order derivatives are neglected. To summarize, in the rest of this chapter we will neglect all third order derivatives in the generalized Fokker-Planck equation, Eq.(5.10). The eects of their inclusion will be the topic of further investigations. For our specic two-site interacting system, a classical simulation of the generalized Fokker-Planck equation without the third order derivatives can be casted into the following form, @ @t P = X i;j @ @i (2) ij j @ @i (2) ijj + @2 @i@j D (2) ij + @2 @i@j D (2) ij (5.25) + @2 @i@j D (2) ij + @2 @i@j D (2) ij P; (5.26) where we have replaced i with i and treat them as independent variables. This treatment in fact means that we are now working in the positive-P representation [70, 73]. For non-interacting systems, the diusion matrix D is always positive denite so that we can work with the P representation. For interacting systems, however, D is not positive denite therefore we have to use the more general representation. Here we use the the positive- P representation. We denote the coecients as (2) ij , (2) ij and D (2) ij , D (2) ij , D (2) ij , D (2) ij . In principle, there could also be terms such as @@i (2) ijj , but for this specic problem such terms do not appear. The explicit denitions of those coecients are, (2) ij = ij + iUiiij ; (2) ij = (ij) iUiiij ; (5.27) and D (2) ij = i U 2 2i ij 2U 2 X ;m1;m2 (D;j;m1;m2;i +D;i;m1;m2;j) m1m2 ; (5.28a) 91 5.5. Approximate numerical solution for interacting systems D (2) ij = i U 2 2i ij 2U 2 X ;m1;m2 (D;j;m1;m2;i +D;i;m1;m2;j) m1m2 ; (5.28b) D (2) ij = D (2) ji = Dij + 2U 2 X ;m1;m2 [(D;m1;j;m2 +D;m1;m2;j) ;i +(D;m2;i;m1 +D;m2;m1;i) ;j ] m1m2 : (5.28c) Eq. (5.26) can be solved through the following Langevin equation [98] in positive-P representation [70, 72], di = X j (2) ij jdt+ p 2dwi (t) ; (5.29a) di = X j (2) ijjdt+ p 2d wi (t) ; (5.29b) where dw (t) and d w (t) are stochastic processes satisfying, dwi (t) dwj t0 = D (2) ij tt0dt; (5.30a) dwi (t) d wj t0 = D (2) ij tt0dt; (5.30b) d wi (t) dwj t0 = D (2) ij tt0dt; (5.30c) d wi (t) d wj t0 = D (2) ijtt0dt: (5.30d) We can dene a 2N -dimensional symmetric matrix D, D = 24 D (2) ij NN D (2) ij NN D (2) ij NN D (2) ij NN 35 : (5.31) Then these random variables can be simplied further and organized as a linear transformations of 2N independent Wiener noises d~! [70], d~w = Bd~!; (5.32) where D = BBT : (5.33) This requires a diagonalization of the 2N2N matrix D at every time step, which is the most time consuming part. However, it is still much simpler 92 5.5. Approximate numerical solution for interacting systems 0 0.2 0.4 0.6 0.8 1 ∆T/T 0 0.1 0.2 0.3 0.4 d dre,NuB 0 0.4 0.8 ∆T/T 0 0.3 0.6 J JB Jre,Nu Figure 5.2: The numerical simulation of the generalized Fokker-Planck equa- tion approach is compared against the BBGKY-like numerical solution of the Redeld equation, for interacting systems. The calculated current Jre;Nu is very close to Jre;B while dre;NuB in the main panel is relatively large. Other parameters are set to be = 0:5; U = 0:1; T = 2:0; = 2:0. than solving directly the d2-dimension linear system for the central system with a d-dimensional Hilbert space, since d increases exponentially with N . We have seen in Fig.5.1 that the numerical simulation leads to accu- rate results for non-interacting systems, where analytical solutions are pos- sible. Next, we compare this numerical simulation against the BBGKY-like method for interacting systems, in Fig.5.2. The currents Jre;Nu and J re;B are in agreement with each other up to high temperature biases. However, dre;NuB | the dierence between the two Green's functions { is relatively large. This indicates that a reasonable result can be obtained using this numerical simulation, but there is space its improvement. In the numerical simulation of the generalized Fokker-Planck equations we used the Euler method for the Ito form of the stochastic dierential equation. The Euler method is not always stable, especially for large U . In fact, more stable and ecient methods such as the semi-implicit method for the Stratonovich form have been developed to solve Fokker-Planck equa- 93 5.6. Solving the LOLE via the GFPE method tions [105]. Using such methods might improve stability and accuracy of the results. However, currently there is an additional technical diculty in implementing such methods for the Redeld equation: matrices B are found numerically so it is hard to derive the corresponding Stratonovich form. The relatively large dre;NuB could also be dure to the truncation of the true generalized Fokker-Planck equation to the standard Fokker-Planck equation. This calls for new methods to solve the generalized Fokker-Planck equation. 5.6 Solving the local-operator Lindblad equation via the generalized Fokker-Planck equation method Instead of the Redeld equation, some authors prefer to use the local- operator Lindblad equation in the study of transport [58, 59, 99]. In this section, we demonstrate how to nd NESSs from the local-operator Lindblad equation using the coherent-state representation. The Redeld equation nds that the operators m̂ for the sites coupled to baths typically involve more than just the operator at the coupling sites. The local-operator Lind- blad equation, on the other hand, chooses m̂ / h1 + n (!k;)ia, which results in the much simpler form of the coupling to baths: @(t) @t = i[HS ; (t)] 2 n h1 + n (1; TL; L)i h ay1; a1(t) i + hn (1; TL; L)i h a1; a y 1(t) i + h:c: +h1 + n (N ; TR; R)i h ayN ; aN(t) i + hn (N ; TR; R)i h aN ; a y N(t) i + h:c: o : (5.34) In the case of this specic two-site systems, we have chosen that N = 2 and 1;2 = 0. Note that we get Eq. (5.34) from the Redeld equation Eq. (5.2) by letting D;m = m [1 + n (; T; )] ; D;m = mn (; T; ) ; D;m1m2m3 = 0: (5.35) These relations greatly simplify the resulting generalized Fokker-Planck equation. Following the same procedure as for the Redeld equation, we 94 5.6. Solving the LOLE via the GFPE method derive a similar Langevin equation for Eq. (5.34), dn = " i (n1 + n+1) iU2nn 2 X (n) # dt+B (n) 1 dw (n) ; (5.36a) dn = " i (n1 + n+1) + iUn2n 2 X (n) # dt+B (n) 2 dw (n) ; (5.36b) where the matrix B(n) = B (n) 22 satises, B(n) B(n) T = iU2n 22P n (; T; ) n 22 P n (; T; ) n iU 2 n ; (5.37) and n dw (n) o is a set of 2N independent Wiener noises. Note that similar equations have been used to describe pure isolated dynamics [75] or T = 0 systems driven by the local-operator Lindblad equation [73], for interact- ing bosonic systems. The 2N -dimensional matrices B are decoupled into N 2-dimensional matrices B(n), which can be diagonalized analytically. There- fore, it is technically easier to solve the local-operator Lindblad equation than the Redeld equation and the Stratonovich form is readily supplied into the xmds software [105]. Results from this local-operator Lindblad equation have been plotted in Fig.5.3. For non-interacting systems, ana- lytical expressions can be derived, and one nds that the NESSs from the local-operator Lindblad equation follows the general Gaussian distribution in Eq. (5.15) and that the correlation matrix lole becomes, lole = n (0; TL; L) + n (0; TR; R) 4 1 0 0 1 + 2 [n (0; TL; L) n (0; TR; R)] 4 (42 + 1) 2 i i 2 : (5.38) In Fig.5.3, we have compared these analytical and numerical results against the BBGKY results. Results for both currents and Green's func- tions are found to be consistent for the two approaches. This conrms that our generalized Fokker-Planck equation method works for both the Redeld equation and the local-operator Lindblad equation. 95 5.7. Conclusions 0 0.2 0.4 0.6 0.8 1 ∆T/T 0 0.1 0.2 d dlole,AnB dlole,NuB 0 0.4 0.8 ∆T/T 0 0.2 0.4 J Jlole,B Jlole,An Jlole,Nu 0 0.2 0.4 0.6 0.8 1 ∆T/T 0 0.1 0.2 0.3 0.4 d dlole,NuB 0 0.4 0.8 ∆T/T 0 0.2 0.4 J Jlole,B Jlole,Nu Figure 5.3: Results from the generalized Fokker-Planck equation approach are compared against the BBGKY solution of the local-operator Lindblad equation for non-interacting(a) and interacting systems(b). Results for both calculated currents and the Green's functions are consistent. 5.7 Conclusions To summarize, in this section we map the Redeld equation into a gener- alized Fokker-Planck equation, using the coherent-state representation. For non-interacting systems, the resulting generalized Fokker-Planck equation can be solved analytically and an expression of non-equilibrium stationary states is obtained. This avoids working with the quantum master equation in a space with an exponentially growing (fermions and spins) or innite (bosons) dimension. For interacting systems, rstly this approach is found to be consistent with the BBGKY-like method. Secondly, a preliminary classical simulation via the Langevin equation with 2N continuous complex numbers is discussed. It is found to work well for weakly interacting sys- tems, although more work is needed. The advantage of this method, if it can be made to be accurate at stronger interactions, is that the total number of variables (2N) grows linearly with the system size. The same method is applied to solve the local-operator Lindblad equa- tion. Here the method is even more ecient, since at every time step one needs to solve N 2-dimensional eigenvalue problems instead of a 2N - dimensional one, like for the Redeld equation. Results are found to be consistent with those from the BBGKY-like method. Further development of the techniques for solving the generalized Fokker- Planck equation is required to make full use of the potential of this approach based on the coherent-state representation. 96 Chapter 6 Application: heat transport in spin chains 6.1 Introduction The physical question addressed in this chapter is the validity of the phe- nomenological Fourier's law of thermal transport. A system which obeys (disobeys) Fourier's law is said to have normal (anomalous) transport, so this question is equivalent with identifying the conditions for normal transport. Even for classical systems this issue is not completely settled (see Ref. [76] for a review), although it was found that normal transport is closely related to the onset of chaos. For quantum systems a similar relation has also been proposed [106], where the quantum chaos is characterized by features of level spacing statistics [37, 107{113]. A widely believed conjecture is that anomalous heat transport occurs in systems described by integrable Hamiltonians (see Ref. [34] for a review). In this context, integrability is usually taken to mean the existence of a Bethe ansatz solution [114]. The validity of this conjecture is still not clear for heat transport in quantum spin chains, despite considerable eort [19, 35, 36, 42{ 44, 56, 59, 60, 82, 83, 106, 115]. Most of this work [35, 42{44, 60, 82, 115] studied innite and/or periodic chains, and used the standard Kubo for- mula [6, 40] where nite/zero Drude weight signals anomalous/normal trans- port. For integrable systems, the standard Kubo formula always predicts anomalous heat transport. In fact, full integrability is not even necessary, all that is needed is commutation of the heat current operator with the total Hamiltonian [115]. However, in Chapter 3, we have shown that the validity of the standard Kubo formula is questionable in the study of ther- mal transport, and that one needs to consider explicitly a system connected to baths. We note that the \integrability" of the chain connected to baths may be lost even if the isolated chain is integrable, since the terms describing the coupling to the baths lead to a non-vanishing commutator between the heat current operator and the total Hamiltonian [116{118]. Both these facts 97 6.1. Introduction might invalidate the main argument for anomalous transport of integrable systems, based on the standard Kubo formula. For simplicity of terminol- ogy, here we call a system integrable if it is so when isolated, since the term is used in this sense in the above conjecture. Studies which explicitly consider the eects of the baths, while fewer, also give contradictory results. Using the local-operator Lindblad equation, Prosen et al. showed that an integrable gaped (Jz > Jxy) XXZ chain has normal spin conductivity while its energy transport is anomalous [59, 119], while using the same approach, Michel et al. claimed that when Jz > 1:6Jxy, spin chains have normal energy transport [56]. Particularly, a recent work in Ref. [119] shows that results from the local-operator Lindblad equation and from the standard Kubo formula are consistent. Such con
icting results exist not only in theoretical studies, but also in experiments. Anomalous heat transport observed experimentally in systems described by integrable models, such as (Sr,Ca)14Cu24O41, Sr2CuO3 and CuGeO3, [79, 120, 121] seems to validate this conjecture, but Ref. [81] nds normal transport in Sr2CuO3 at high temperatures. Even if there was no con
ict between results based on the local-operator Lindblad equation [56, 59], as discussed in Chapter 2 this approach is less reliable than the Redeld equation. For example, only the latter leads to the proper Boltzmann distribution if the baths are not biased. This is why here we investigate heat transport in nite spin chains coupled to ther- mal baths, using the Redeld equation. There are already several other such studies [36, 83]. However, these are for spin systems that are either non-interacting (thus trivially integrable), or are non-integrable. (We call a spin chain \non-interacting" if it can be mapped, for example through a Jordan-Wigner transformation, to a Hamiltonian for non-interacting spinless fermions). Anomalous transport is found for the former and normal trans- port is found for the latter. Disordered spin systems have also been found to have normal transport [37]. There are no examples of clean, interacting but integrable systems investigated via the Redeld equation. We study such systems here. We nd no direct relation between inte- grability and anomalous transport. Instead, based on our systematic inves- tigations, we propose a new conjecture: anomalous conductivity is observed in systems that can be mapped to non-interacting fermions. Those that map to interacting fermion models, whether integrable or not, show normal transport. This is in qualitative agreement with the demonstration by Sirker et al. that diusion is universally present in interacting 1D systems [42]. It is important to note that, chronologically, this was the rst project undertaken once the formalism of the Redeld equation was adopted. This 98 6.2. The model and its Redeld equation explains why the results shown below are derived using (inecient) direct methods to solve the Redeld equation, and therefore are restricted to rather small systems. This is also why our conjecture needs to be further examined, for larger systems. It was this computational ineciency that drove the subsequent research on more ecient methods, presented in the previous chapters. In future, we plan to use these ecient methods to extend this study to larger systems. The chapter is organized as follows. Section 6.2 denes the system and its Redeld equation. In Sections 6.3 we discuss the numerical methods and in Section 6.4 we dene the thermal current and local temperatures. We then proceed to present numerical results in Section 6.5, and then conclude. 6.2 The model and its Redeld equation We consider an N -site chain of spins-12 described by the Hamiltonian: HS = N1X i=1 Jxs x i s x i+1 + Jys y i s y i+1 + Jzs z i s z i+1 ~B NX i=1 ~si (6.1) while the heat baths are collections of bosonic modes: HB = X k; !k;b y k;bk; (6.2) where = R=L indexes the right/left-side baths. The system-baths coupling is taken as: V = X k; V () k s y byk; + bk; (6.3) where syL = s y 1 and s y R = s y N , i.e. the left (right) thermal bath is only coupled to the rst (last) spin and can induce its spin-
ipping. This is because we choose ~B ~ey = 0 while j ~Bj is nite, meaning that spins primarily lie in the x0z plane so that sy acts as a spin-
ip operator. From Eq. (2.29), the resulting Redeld equation for (t) is: @(t) @t = i[HS ; (t)] 2 X =L;R ([sy; m̂(t)] + h:c:) (6.4) where m̂ = s y . Here, () refers to the element-wise product of two matrices, hnja bjmi = hnjajmihnjbjmi. The bath matrices L;R are dened 99 6.3. Numerical methods in terms of the eigenstates of the system's Hamiltonian HS jni = Enjni as: = X m;n jmihnj h ( mn)n ( mn)D ( mn) jV ()kmn j2 +( nm) (1 + n ( nm))D ( nm) jV ()knm j2 i where mn = EmEn = nm and kmn is dened by !kmn; = mn, i.e. a bath mode resonant with this transition. Furthermore, (x) is the Heaviside function, n( ) = e 11 is the Bose-Einstein equilibrium distribu- tion for the bosonic modes of energy at the bath temperature T = 1=, and D( ) is the bath's density of states. The product D ( mn) jV ()kmn j2 is the bath's spectral density function. For simplicity, we take it to be a constant independent of , m and n. 6.3 Numerical methods It is straightforward to use the Runge-Kutta method to integrate the Red- eld equation, Eq. (6.4), starting from some initial states. The memory cost is proportional to 22N for an N -spin chain, but it takes a very long time for the integration to converge to the stationary state. This is not surprising since, in principle, the stationary state is reached only as t!1. As already discussed, another approach [58] is to solve directly for the stationary state from: L(1) = 0; (6.5) i.e. to nd the eigenstate for the zero eigenvalue of the 2N -dimensional matrix L. However, in this case, the memory cost is proportional to 42N , which is much worse than for the Runge-Kutta method. The method which has better memory eciency than this eigenvalue problem and also better time eciency than the Runge-Kutta method, is to convert the equation into a linear system and solve it via matrix-free methods like the Krylov space methods, which requires only matrix-vector multiplication but not explicitly the matrix. The eigenvalue problem can be rewritten as a linear system of equations after explicitly using the normal- ization condition tr () = 1 such that L(1) = ; (6.6) where = [1; 0; ]T and L is found from L by replacing the rst row by tr(). Then we use for example the generalized minimal residual method 100 6.4. Denitions of the thermal current and local temperatures (GMRES) [122], which requires only the matrix-vector multiplication rule. This new way of solving the Redeld equation has memory cost of 22N and time eciency of a linear system with dimension 22N . It is still a di- rect method so its eciency is not comparable with methods such as the BBGKY-like approach and the coherent-state representation method that we discussed previously, however, unlike them it gives an exact result. This is important until one can understand whether the further approximations made in the more ecient methods can, for example, destroy the integra- bility of the system. Also, because the commutation relations between spin operators are quite dierent from those of fermions and bosons, the more ecient methods are not readily generalizable. Developing and applying them to this problem will be the topic of future work. 6.4 Denitions of the thermal current and local temperatures We rewrite HS = PN1 i=1 hi;i+1 + PN i=1 hi, where hi;i+1 is the exchange be- tween nearest-neighbor spins and hi is the on-site coupling to the magnetic eld. We can then dene a local site Hamiltonian h (S) i = 1 2 hi1;i + hi + 1 2 hi;i+1 (6.7) with h0;1 = hN;N+1 = 0, and a local bond Hamiltonian h (B) i = 1 2 hi + hi;i+1 + 1 2 hi+1 (6.8) such that HS = NX i=1 h (S) i = N1X i=1 h (B) i : (6.9) The local bond Hamiltonians can be used to derive the heat current operator from the continuity equation ĵi!i+1 ĵi1!i = rĵ = @h (B) i @t = i h HS ; h (B) i i : (6.10) This results in ĵi!i+1 = i h h (B) i ; h (B) i+1 i (6.11) 101 6.4. Denitions of the thermal current and local temperatures for i = 1; ; N 2. As expected, in the steady state we nd hĵi!i+1i = tr ĵi!i+1 (1) = J to be independent of i. Similarly, we dene the spin polarization hszi i and the spin current for XXZ model, Js = Jxy syi s x i+1 sxi syi+1 : (6.12) Knowledge of the steady state heat current J , as such, is not enough to decide whether the transport is normal or not. Consider an analogy with charge transport in a metal connected to two biased leads. What shows whether the transport is anomalous or not is the prole of the electric potential, not the value of the electric current. In anomalous transport (for clean, non-interacting metals) all the voltage drop occurs at the ends of the sample, near the contacts. Away from these contact regions, electrons move ballistically and the electric potential is constant, implying zero intrinsic resistance. For a dirty metal, scattering takes place everywhere inside the sample and the electric potential decreases monotonically in between the contact regions, i.e. the sample has nite intrinsic resistivity. In principle, the scaling of the current with the sample size, for a xed eective bias, also reveals the type of transport: for anomalous transport, the current is independent of the sample size once its length exceeds the sum of the two contact regions, while for normal transport it decreases like inverse length. The problem with this approach is that one needs to keep constant the eective bias, i.e. the dierence between the applied bias and that in the contact regions. Furthermore, since we can only study relatively short chains, the results of such scaling may be questionable. It is therefore desirable to use the equivalent of the electric potential for heat transport and to calculate its prole along the chain, in order to determine the type of transport. This, of course, is the \local temperature", which is a dicult quantity to dene properly. One consistency condition for any denition is that if TL = TR = T , i.e. the system is in thermal equilibrium at T , then all local temperatures should equal T . We dene local site temperatures Ti which fulll this condition in the following way. Since we know all eigenstates of HS , it is straightforward to calculate its equilibrium density matrix at a given T , eq (T ) = 1 Z P n e En jnihnj, where Z = P n e En . Let then hh(S)i ieq (T ) = tr[eq (T )h(S)i ]. We dene Ti to be the solution of the equation: hh(S)i ieq (Ti) = tr[ (1)h(S)i ]: (6.13) In other words, the steady-state value of the energy at that site equals the energy the site would have if the whole system was in equilibrium at Ti. 102 6.5. Results Of course, we can also use other \local" operators such as h (B) i to calculate a local bond temperature Ti+ 1 2 . We nd that when these denitions are meaningful, the results are in very good agreement no matter what \local" operator is used. This type of denition of Ti is meaningful only if a large magnetic eld B is applied. For small B, the expectation values hh(S)i ieq (T ) are very weakly T -dependent, so that tiny numerical errors in the steady-state value can lead to huge variations in Ti. Addition of a large B is needed to obtain hh(S)i ieq (T ) which varies fast enough with T for values of interest, so that a meaningful Ti can be extracted. Since we could not nd a meaningful denition for Ti when ~B ! 0, we cannot investigate such cases. Note, however, that most integrable models remain integrable under addition of an external eld ~B = Bêz. 6.5 Results In all our calculations, we take Bz = 1 and the exchange J 0:1. Temper- atures TL=R = T (1 =2) should not be so large that the steady state is insensitive to the model, or so small that only the ground-state is activated. Reasonable choices lie between min(Jx; Jy; Jz) and NB, which are roughly the smallest, respectively the largest energy scales for the N -site spin chain. In Fig. 6.1 we show typical results for (a) local temperature proles Ti; Ti+ 1 2 and (b) local spin polarization proles Szi . We apply a large bias = (TL TR)=T = 0:4 for clarity, but we nd similar results for smaller . For these values, it seems that the \contact regions" include only the end spins. The prole of the central part of the chain is consistent with anomalous transport (
at prole) for the XX chain (Jx = Jy, Jz = 0) and shows normal transport (roughly linear prole) in all the other cases. XY chains with Jx 6= Jy behave similarly with the XX chain. We nd similar results for ferromagnetic couplings. We also nd that the ratio between the eective temperature dierence, T2 TN1, and the applied temperature dierence TL TR, is not a constant for dierent system sizes N . Therefore, in our numerical calculation, it is not possible to keep the eective temperature dierence as a constant by applying the same temperature dierence on the ends while changing the system size N . This implies that the dependence of J v.s. N can not be used as an indicator of normal or anomalous transport. At most, it can provide a very rough qualitative picture. We have plotted several of the J vs. N in Fig.6.2. We see that J is independent of N for XX 103 6.5. Results 1 2 3 4 5 6 7 8 9 10 i 1.8 2 2.2 Ti XX XXZ(0.05) XXX XXZ(0.2) XYZ 1 2 3 4 5 6 7 8 9 10 i 0.1 0.12 0.14 Szi Figure 6.1: Plot of (a) local temperature prole and (b) local spin polariza- tion prole for chains with N = 10. In all cases TL = 2:4; TR = 1:6; = 0:1; Jx = 0:1; Bz = 1:0. Other parameters are Jy = 0:1; Jz = 0:0(XX); Jy = 0:1; Jz = 0:05 (XXZ0.05); Jy = Jz = 0:1 (XXX); Jy = 0:1; Jz = 0:2 (XXZ0.2) and Jy = 0:2; Jz = 0:3 (XYZ). Only the XX chain shows
at Ti and Szi proles. All other models have almost linearly decreasing proles of both temperature and Sz. 104 6.5. Results 4 5 6 7 8 9 10 N 0 0.0005 0.001 0.0015 0.002 J XX XYZ J=0.001*(N-1.9)-0.89 XXX J=0.0017*(N-0.23)-0.37 Figure 6.2: Heat currents are plotted as a function of system size for three models: XX, XY Z and XXX. J is independent of N for XX chains, and it decreases with N for the other two cases. Due to the uncertainty about the eective applied temperature dierence, the tting curves should only be regarded as a guide for eyes. The same parameters as those in Fig.6.1 are used. chains while it decreases with increasing N for XY Z and XXX chains. Another way to examine size-dependence is to compare the temperature and Sz proles for all available sizes. We normalize all the proles for various models with dierent sizes by harvesting only the data in the central regions and converting them into dimensionless values in [0; 1], for example, from (i; Ti) to ( i2 N3 ; TiTN1 T2TN1 ); i = 2; 3; : : : ; N 1. In Fig.6.3, we plot such normalized temperature and Sz proles. We see that curves for all values of N collapse onto one single straight line and furthermore, there is no dierence between XY Z and XXX chains. Currently we can only study small systems, but from these normalized proles, no qualitative dierence has been found for dierent sizes. All these are integrable models. The XX model is special because it can be mapped to non-interacting spinless fermions with the Jordan-Wigner transformation [84]. A nite Jz leads to nearest-neighbor interactions be- tween fermions. Eigenmodes for models with Jz 6= 0 can be found using Bethe's ansatz, but they cannot be mapped to non-interacting fermions. In order to investigate where this transition between anomalous and 105 6.5. Results 0 0.2 0.4 0.6 0.8 1 (i-1)/(N-2) 0 0.5 1 (T i- T N -1 )/( T 2 - T N -1 ) N=10,XXX9 8 7 6 5 4 10, XYZ 0 0.2 0.4 0.6 0.8 1 (i-1)/(N-2) 0 0.5 1 (S z i - Sz 2)/ (S z N -1 - Sz 2) Figure 6.3: Normalized temperature and Sz proles are plotted. Curves for dierent values of N collapse onto a single straight line and we see no dierence between XY Z and XXX chains. normal transport occurs, we plot temperature proles for systems with small Jz in Fig.6.4(a). We see that even for very small Jz the above observation still holds. For Jz = 0:01 there is a slight qualitative dierence, namely that the contact re- gions seem larger than just the end spins, as can be seen from the normalized proles. We think this may be due merely to limitations of the numerical accuracy. In Fig.6.4(b), we investigate what is the minimum value of Bz that we can use with condence. As we pointed out before, our denition of local temperature and also the idea of local a spin polarization Szi are only meaningful for suciently large Bz. We nd that roughly this implies Bz > 0:3. This is reasonable, since for this limiting value the energy scale due to coupling to the local Bz eld (0:3=2) is comparable with the energy related to exchange (3 0:1=4). We found this generic behavior for a wide range of parameters. When 2 [0:03; 0:2], T 2 [0:3; 30:0] and 0:01, the spin chain has normal conductivity when Jz 2 [0:02; 0:5] and anomalous conductivity when Jz = 0. These results lead us to conjecture that it is the presence or absence of interactions, rather than integrability, that determines whether or not the heat transport is normal. We nd no dierence between thermal transport 106 6.5. Results 2 4 6 8 10 i 1.8 1.9 2 2.1 2.2 B z =0.1 B z =0.3 0 0.4 0.8 0 0.5 1 2 4 6 8 10 i 1.8 1.9 2 2.1 2.2 T i J z =0.00 J z =0.01 J z =0.02 J z =0.03 J z =0.04 (a) (b) Figure 6.4: (a) Temperature proles for small Jz values are plotted. Here Jx = Jy = 0:1; Bz = 1:0. The eective temperature drop becomes smaller for smaller Jz, but the slope is still nite if Jz 0:02. The inset plots normalized proles, and shows no obvious dierences for models with Jz 0:02. The case of Jz = 0:01 shows a nite slope but the normalized prole indicates a slight dierence: the contact regions seem to be more extended than in the other cases. (b) Temperature proles of XXX chains with J = 0:1 and small values of Bz. A linear temperature drop is observed for Bz = 0:3, but a
at prole for Bz = 0:1. Other parameters are the same as in Fig.6.1. 107 6.5. Results and spin transport, unlike the results of Prosen et al. [59], based on the local-operator Lindblad equation. The conjecture may be tested further on various systems. One candidate is the Ising model in a transverse eldBx. It maps to non-interacting spinless fermions [123] and if we add a Bz eld, it becomes interacting. Another closely related model is a system of spinless fermions on a tight-binding chain, with the nearest-neighbor interaction: HS = NX l=1 cyl cl t N1X l=1 cyl cl+1 + c y l+1cl + V0 N1X l=1 cyl+1cl+1c y l cl; (6.14) HB = X k;=L;R !k;b y k;bk;; (6.15) VSB = X k; cybk; + cb y k; : (6.16) The XXZ chains map exactly into this HS model, after using the Jordan- Wigner transformation. However, the coupling to the baths in this fermionic system is dierent from that in the spin systems. For a spin system, the Jordan-Wigner transformation maps the y operator into an operator much more complicated than the cy or c used in this fermionic model. Therefore, although the two models are closely related, they are not identical. Results on this model can be interpreted as another test of our conjecture or at least a check of whether the observations reported above are in
uenced by the specic model for the system-bath coupling. In Fig.6.5, we plot local temperature proles along the Ising spin chains in panel (a) and for fermionic chains in panel (b). The results support our conjecture: interactions lead to normal transport. Note that is set to be much larger than t, to mirror the condition Bz J . When is comparable to t, both non-interacting and interacting systems show almost
at temperature proles. In summary, the rst conclusion we draw from these results is that in- tegrability is not sucient to guarantee anomalous transport: several inte- grable models show normal heat transport, in agreement with other studies [42, 58, 59]. The second conclusion is that only models that map onto Hamil- tonians of non-interacting fermions exhibit anomalous heat transport. This is a reasonable sucient condition, since once inside the sample (past the contact regions) such fermions propagate ballistically. However, we cannot, at this stage, demonstrate that this is a necessary condition as well. We therefore can only conjecture that this is the criterion determining whether 108 6.5. Results 1 2 3 4 5 6 7 8 9 10 i 1.8 2 2.2 V0=0.0 V0=0.2 1 2 3 4 5 6 7 8 9 10 i 1.6 1.8 2 2.2 T i Ising x Ising xz (a) (b) Figure 6.5: Temperature proles of (a) Ising spin chains and (b) the V -t model of fermionic chains. The Ising spin chains have anomalous transport when Bz = 0 (Isingx, circles) and normal transport when Bz = 1:0 (Isingxz, squares). The fermionic chains show a
at temperature prole when V0 = 0 (circles) but a linear temperature drop when V0 = 0:2 (squares). Other parameters are Jz = 0:1; Bx = 1:0; = 1:0; t = 0:1; TL = 2:4; TR = 1:6; = 1:0; = 0:1; N = 10. 109 6.5. Results the heat transport is anomalous. If we believe that our results are not ar- tifacts due to mainly the boundary eects of the limited system sizes, then this conjecture is the only consistent qualitative conclusion that we may draw from all above results. In this context, it is important to emphasize again the essential role played by the connection to the baths. In its absence, an isolated integrable model is described by Bethe ansatz type wavefunctions. Diusion is impos- sible since the conservation of momentum and energy guarantees that, upon scattering, pairs of fermions either keep or interchange their momenta. For a system connected to baths, however, fermions are continuously exchanged with the baths, and the survival of a Bethe ansatz type of wavefunction becomes impossible. In fact, even the total momentum is no longer a good quantum number. We believe that this explains why normal transport in systems mapping to interacting fermions is plausible. Normal transport is also possible for non-interacting fermions, if they are subject to elastic scattering on disorder. This can be realized, for example, by adding to the XY model a random eld Bz at various sites. In that case, as expected, normal conductance has been found in disordered systems [37]. On the other hand, anomalous transport can also occur in models which map to homogeneous interacting fermions if the bath temperatures are very low. Specically, consider the XXZ models. Because of the large Bz we use, the ground-state of the isolated chain is ferromagnetic with all spins up. The rst manifold of low-energy eigenstates have one spin
ipped (single magnon states), followed by states with two spins
ipped (two magnon states), etc. The separation between these manifolds is roughly Bz, although because of the exchange terms each manifold has a fairly considerable spread in energies and usually overlaps partially with other manifolds. If both TL; TR Bz, only single-magnon states participate in the trans- port. We can then study numerically very long chains by assuming that the steady-state matrix elements nm vanish for all other eigenstates (Sz;tot is a good quantum number for these models). In this case we nd anomalous transport for all models, whether integrable or not. This is reasonable, since the lone magnon (fermion) injected on the chain has nothing else to interact with, so it must propagate ballistically. We can repeat this restricted calculation by including the two-magnon, three-magnon, etc. manifolds in the computation. As expected, the results agree at low TL; TR, but dierences appear for higher TL; TR, when these higher-energy manifolds become thermally activated. In such cases, the transport becomes normal for the models mapping to interacting fermions as soon as the probability to be in the two (or more) magnon sector be- 110 6.6. Conclusions comes nite. In other words, as soon as multiple excitations (fermions) are simultaneously on the chain, and inelastic scattering between them becomes possible. These results may explain the heat transport observed experi- mentally in compounds such as Sr2CuO3 [81], where at low temperature anomalous transport was found while at high temperature normal transport was reported. 6.6 Conclusions Based on an extensive study of quantum spin chains using the Redeld equa- tion, we propose a new conjecture for what determines the appearance of anomalous heat transport at all temperatures in spin chains. Unlike pre- vious suggestions linking it to the integrability of the Hamiltonian or the existence of energy gaps, we propose that, for clean systems, the criterion is the mapping of the Hamiltonian onto a model of non-interacting fermions. Our results showing normal energy and spin transport in systems with - nite Jz contradict the major conclusions of Ref. [56] and Ref. [59], which, however, are based on the less trustworthy local-operator Lindblad equation. This conjecture should be checked for larger systems, where more reliable information on the relation between J and N can be extracted. The more ecient methods presented in the previous two chapters might allow us to undertake such a study. For discussion of the thermal conductance, which is related to higher order Green's functions | G2 and G3, not only G1 | extension of the methods to higher orders may be needed. 111 Chapter 7 Conclusions and discussions We began this project with a study of the thermal transport of spin chains in the framework of theories of open quantum systems. At that time, simi- lar questions had been investigated using either the standard Kubo formula, the non-equilibrium Green's functions method or the local-operator Lind- blad equation. However, the idea of using the Redeld equation for these questions had just been proposed and it had been applied only to non- interacting systems or very small interacting but non-integrable systems. We rst developed and applied the direct methods to the Redeld equation for spin systems. From the preliminary results for short chains we found a new conjecture, linking anomalous transport to existence of mapping onto non-interacting fermions. This is the main result of Chapter 6. It then became necessary to study larger interacting systems to either conrm or reject this observation. However, at the time there were no avail- able methods capable of solving the Redeld equation for large interacting systems. Thus, we turned to investigate such methods. We rst tried to make direct use of the standard Kubo formula. However, we found that the standard Kubo formula is not applicable to open systems explicitly cou- pled to baths and instead we had to derive a similar Kubo formula from the Redeld equation using linear response theory. We call the result the open Kubo formula. This is the content of Chapter 3. Several forms of such OKFs were proposed and tested. The proper approach was shown to not suer from singularities, like the the standard Kubo formula, however it is not very ecient. A more ecient OKF was also identied, however only for certain types of quantities. Neither approach is ecient enough to allow the study of systems of the same size as permitted by the non-equilibrium Green's functions method. Next, inspired by the non-equilibrium Green's function method[13{15] and also by Saito's work on non-interacting systems[36], we focused on Green's functions instead of the reduced density matrices and developed the BBGKY-like hierarchical method presented in Chapter 4. Further ap- proximations are needed to truncate this hierarchy, if there are interactions. We identied two possible approaches. The rst replaces, at a certain level of 112 Chapter 7. Conclusions and discussions the hierarchy, all higher-order Green's functions by their equilibrium values. This method is computationally less ecient but accurate up to large values of the interaction. The second method uses a cluster expansion, basically factorizing the higher-order Green's functions in terms of lower-order ones. This approach is applicable to much larger systems. While the simplest version (truncation at rst level of the hierarchy) is accurate only for rather small interactions, we showed that systematic increase of the truncation level leads to both a systematic improvement of its accuracy, and a larger range of interactions where this holds. Also, truncation at second level gen- erates correlated two-particle Green's functions, beyond the Hartree-Fock approximation. In Chapter 5 we proposed another method: solving the Redeld equation in the coherent-state representation. The resulting dierential equation has the form of the generalized Fokker-Planck equation, so we call this method the generalized Fokker-Planck equation method. Similar methods have been applied to pure dynamical evolution or equilibrium states, where studies of large interacting systems with roughly 105000000-dimensional Hilbert spaces have been performed[72{74]. We extend the method to apply it to the Red- eld equation and study also non-equilibrium stationary states. We have conrmed that for non-interacting systems exact results can be found ana- lytically with this new approach, even for non-equilibrium stationary states. Other advantages are that, from a numerical point of view, in this approach interacting systems are handled similarly like non-interacting systems. Also, Green's functions of all orders can be calculated all at once. The challenge, still not fully settled, is to nd ecient and stable ways to solve these dif- ferential equations. We showed some preliminary results, based on classical simulations, that are promising yet require more work. The next step will naturally be to apply these methods to the studies of transport properties of various systems, this time for large systems. The question of validity of Fourier's Law is not yet answered denitely. The second-order form of the BBGKY-like method seems to have a good poten- tial there. However, the study of spin systems requires dierent kinds of BBGKY-like equations, and also to understand whether the cluster approx- imation interferes with the integrability of a model. Therefore, although the methods we proposed are in principle applicable to this problem, spe- cic forms of hierarchical equations need to be re-derived. Future possible projects, thus, are to develop the BBGKY-like method for spin systems, or to use the current BBGKY-like method to discuss thermal transport of bosonic or fermionic systems, instead of spin systems. Alternatively, we can also use the generalized Fokker-Planck equation method. Currently the simulation 113 Chapter 7. Conclusions and discussions (a) giσ z i σ z c LBL LBR Lc λL λR (b) Figure 7.1: (a) A sketch of exactly solvable decoherence model: a single spin (~c) is coupled to a bath of N spins f~ig. (b) A variation of (a) in that the central system is coupled to the baths only via the ends. is not stable enough and it works only for bosonic systems. Improving this, and also developing similar methods for fermionic systems could be another future project. There are also many small systems showing rich physics, such as sys- tems of a few coupled quantum dots or molecular devices coupled to reser- voirs. Application of the proposed methods to such systems should be very straightforward. The methods can also be applied to more general physical systems, not only to study transport but also other physical questions, such as deco- herence control and classical simulations of quantum computations. In all models studied in this thesis, the baths are assumed to be collections of eigenmodes j! (k)i, connected to the central systems via V k (see for ex- ample Eq. (2.44)). We then assume V k to have certain properties. A more realistic and specic setup is to connect the baths locally to the central system, as sketched in Fig7.1(b). In that case, V k is implicitly dened. The system sketched in Fig.7.1(a) is in fact a famous exactly solvable model of quantum decoherence [85, 124]: a single spin-12 is connected to a 114 Chapter 7. Conclusions and discussions bath of N spins-12 via Ising coupling with strength gi to bath i [125]. The dynamics of the whole composite system is exactly solvable and Zurek et al. showed that the o-diagonal elements (in z basis), thus the coherent part, die out very quickly. A related but rather ambitious question is to study a counterpart with a thermal bath and to ask the following questions: whether or not the reduced density matrix of the central system converges to the expected thermal equilibrium state; whether or not the o-diagonal elements disappear and if so, which \diagonal" and \o-diagonal" parts are relevant, and how that depends on the specic forms of interaction? In fact, such a natural selection of the basis is called, in theories of quantum measurement, eiselection of the pointer states [85, 124]. Furthermore, how do the answers to the above questions depend on various parameters of the composite systems, such as size of the bath, coupling strength between the central system and the bath, and temperature of the bath and so on. We can ask all the questions for a similar setup but now with two baths, biased or not, connected to the systems locally, as sketched in Fig.7.1(b). These questions can be investigated using dierent approaches: starting from the pure dynamics of the whole composite system, which usually how- ever is only applicable to non-interacting systems (both the central system and the baths); through the Redeld equation, which assumes the Markovian approximation, weak coupling and idealized innite-size baths; and nally by the non-equilibrium Green's function method, which is typically applied for congurations similar to the variation in Fig.7.1(b). In fact, there is another approach | the exact master equation or the in
uence functional method, which relaxes the assumption of the Markovian approximation and the weak coupling limit, but still takes the baths to be innitely large. How- ever, it currently works only for non-interacting systems since it is hard to solve such an exact master equation for interacting systems. So let us focus on the above three methods. The generalized Fokker-Planck equation method makes it possible to im- plement the rst approach even for large interacting systems. The second approach can be realized through either the BBGKY-like method or the generalized Fokker-Planck equation method. Implementations of the third method are already well developed. Note, however, that if we set the tem- peratures of the baths to be equal then the non-equilibrium Green's function arrives at the corresponding equilibrium state of the whole composite sys- tem, which does not guarantee that the reduced density matrix of the central system follows the proper equilibrium distribution. On the other hand, the reduced density matrix generally converges to an equilibrium distribution in the theories of open quantum systems approach. Therefore, it will not 115 Chapter 7. Conclusions and discussions be much of a surprise if the non-equilibrium Green's function leads to qual- itatively dierent results. If we want to compare them, a proper model, on which the non-equilibrium Green's function method and the Redeld equa- tion ideally should produce the same result, is required. Comparison be- tween the rst two might shed some light on dierences between the Marko- vian and the non-Markovian approximation, the validity of the innite-size bath assumption, etc. Comparison between the latter two might help us understand them better or even improve both. Lastly, our ecient methods for the Redeld equation can be applied as well to solving the local-operator Lindblad equation. This should make the local-operator Lindblad equation a more powerful tool in those cases where it is applicable. As we mentioned earlier, there are still many groups that use the local-operator Lindblad equation instead of the Redeld equation. It should be very straightforward to implement the methods for the local- operator Lindblad equation. 116 Bibliography [1] G. Mahler and V.A. Weberru. Quantum Networks: Dynamics of Open Nanastructures. Springer-Verlag, 1998. [2] H.-P. Breuer and F. Petruccione. The theory of open quantum systems. Oxford, 2002. [3] U. Weiss. Quantum Dissipative Systems. World Scientic, 1999. [4] S. Datta. Quantum Transport: Atom to Transistor. Cambridge Uni- versity Press, 2005. [5] D.K. Ferry and S. M. Goodnick. Transport in Nanostructures. Cam- bridge Univ. Press, 1997. [6] R. Kubo. Statistical-mechanical theory of irreversible processes 1. general theory and simple applications to magnetic and conduction problems. Journal of the Physical Society of Japan, 12(6):570{586, 1957. [7] R. Kubo, M. Yokota, and S. Nakajima. Statistical-mechanical theory of irreversible processes 2. response to thermal disturbance. Journal of the Physical Society of Japan, 12(11):1203{1211, 1957. [8] R. Landauer. Electrical resistance of disordered one-dimensional lat- tices. Philosophical Magazine, 21(172):863{867, 1970. [9] R. Landauer. Spatial variation of currents and elds due to local- ized scatterers in metallic conduction. IBM Journal of Research and Development, 1(3):223{231, 1957. [10] M. Buttiker. 4-terminal phase-coherent conductance. Physical Review Letters, 57(14):1761{1764, 1986. [11] M. Buttiker, Y. Imry, R. Landauer, and S. Pinhas. Generalized many- channel conductance formula with application to small rings. Physical Review B, 31(10):6207{6215, 1985. 117 Bibliography [12] H. Haug and A. P. Jauho. Quantum Kinetics in Transport and Optics of Semiconductors. Springer-Verlag, 1996. [13] M. Koentopp, C. Chang, K. Burke, and R. Car. Density functional calculations of nanoscale conductance. Journal of Physics-Condensed Matter, 20(8):083203, 2008. [14] C. Caroli, R. Combesco, P. Norzieres, and D. Saintjam. Direct calcu- lation of tunneling current. Journal of Physics C-Solid State Physics, 4(8):916{929, 1971. [15] R. Combesco. Direct calculation of tunneling current 3: Eect of localized impurity states in barrier. Journal of Physics C-Solid State Physics, 4(16):2611{2622, 1971. [16] J. Taylor, H. Guo, and J. Wang. Ab initio modeling of quantum transport properties of molecular electronic devices. Physical Review B, 63(24):245407, 2001. [17] A. G. Redeld. Relaxation theory: density matrix formulation. In D. M. Grant and R. K. Harris, editors, Encyclopedia Nuclear Magnetic Resonance, pages 4085{4092. Wiley, 1996. [18] G. D. Mahan. Many-Particle Physics. Plenum Press, 1990. [19] A. V. Sologubenko and H. R. Ott. Energy transport in one-dimensional spin systems. In D. Baeriswyl and L. Degiorgi, editors, Strong Interac- tions in Low Dimensions, chapter 12, pages 383{417. Springer, 2005. [20] A. Yacoby. Transport in quatnum wire. In D. Baeriswyl and L. De- giorgi, editors, Strong Interactions in Low Dimensions, chapter 10, pages 321{346. Springer, 2005. [21] A. Nitzan and M. A. Ratner. Electron transport in molecular wire junctions. SCIENCE, 300(5624):1384{1389, 2003. [22] A Salomon, D Cahen, S Lindsay, J Tomfohr, VB Engelkes, and CD Frisbie. Comparison of electronic transport measurements on or- ganic molecules. Advanced Materials, 15(22):1881{1890, 2003. [23] N. J. Tao. Electron transport in molecular junctions. Nature Nan- otechnology, 1(3):173{181, 2006. [24] W. Goepel and C. Ziegler. Nanostructures Based on Molecular Mate- rials. Wiley-VCH, 1992. 118 Bibliography [25] L. P. Kouwenhoven, C. M. Marcus, P. L. Mceuen, S. Tarucha, R. M. Westervelt, and N. S. Wingreen. Electron transport in quantum dots. In L. L. Sohn, L. P. Kouwenhoven, and G. Schon, editors, Mesoscopic Electron Transport, volume 345 of NATO ADVANCED SCIENCE IN- STITUTES SERIES, SERIES E, APPLIED SCIENCES, pages 105{ 214, 1997. [26] D. Natelson. Handbook of Organic Electronics and Photonics. Amer- ican Scientic Publishers, 2006. [27] A. V. Sologubenko, T. Lorenz, H. R. Ott, and A. Freimuth. Thermal conductivity via magnetic excitations in spin-chain materials. Journal of Low Temperature Physics, 147(3-4):387{403, 2007. [28] A. V. Sologubenko, K. Berggold, T. Lorenz, A. Rosch, E. Shimshoni, M. D. Phillips, and M. M. Turnbull. Magnetothermal transport in the spin-12 chains of copper pyrazine dinitrate. Phys. Rev. Lett., 98(10):107201, Mar 2007. [29] A. V. Sologubenko, T. Lorenz, J. A. Mydosh, A. Rosch, K. C. Short- sleeves, and M. M. Turnbull. Field-dependent thermal transport in the haldane chain compound nenp. Phys. Rev. Lett., 100(13):137202, Apr 2008. [30] C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, and A. Zettl. Breakdown of fourier's law in nanotube thermal conductors. Phys. Rev. Lett., 101(7):075903, Aug 2008. [31] B. Soree, W. Magnus, and W. Schoenmaker. Energy and momen- tum balance equations: An approach to quantum transport in closed circuits. Physical Review B, 66(3):035318, 2002. [32] T. H. Stoof and Y. V. Nazarov. Time-dependent resonant tunneling via two discrete states. Physical Review B, 53(3):1050{1053, 1996. [33] O. N. Jouravlev and Y. V. Nazarov. Electron transport in a double quantum dot governed by a nuclear magnetic eld. Physical Review Letters, 96(17):176804, 2006. [34] X. Zotos and P. Prelovsek. Transport in one dimensional quanutm systems. In D. Baeriswyl and L. Degiorgi, editors, Strong Interactions in Low Dimensions, chapter 11, pages 347{382. Springer, 2005. 119 Bibliography [35] X. Zotos, F. Naef, and P. Prelovsek. Transport and conservation laws. Physical Review B, 55(17):11029{11032, 1997. [36] K. Saito, S. Takesue, and S. Miyashita. Energy transport in the in- tegrable system in contact with various types of phonon reservoirs. Physical Review E, 61(3):2397{2409, 2000. [37] Y. Yan, C.-Q. Wu, G. Casati, T. Prosen, and B. Li. Nonballistic heat conduction in an integrable random-exchange ising chain studied with quantum master equations. Physical Review B, 77(17):172411, 2008. [38] J. Wu and M. Berciu. Heat transport in quantum spin chains: the relevance of integrability. arXiv:1003.1559, submitted to PRB, 2010. [39] H. Mori. Statistical-mechanical theory of transport in
uids. Physical Review, 112(6):1829{1842, 1958. [40] J. M. Luttinger. Theory of thermal transport coecients. Physical Review, 135(6):A1505{A1514, 1964. [41] D. S. Fisher and P. A. Lee. Relation between conductivity and trans- mission matrix. Physical Review B, 23(12):6851{6854, 1981. [42] J. Sirker, R. G. Pereira, and I. Aeck. Diusion and ballistic trans- port in one-dimensional quantum systems. Physical Review Letters, 103(21):216602, 2009. [43] F. Heidrich-Meisner, A. Honecker, and W. Brenig. Thermal transport of the xxz chain in a magnetic eld. Physical Review B, 71(18):184415, 2005. [44] M. Rigol and B. S. Shastry. Drude weight in systems with open bound- ary conditions. Physical Review B, 77(16):161101, 2008. [45] C. Toher, A. Filippetti, S. Sanvito, and K. Burke. Self-interaction er- rors in density-functional calculations of electronic transport. Physical Review letters, 95(14), 2005. [46] F Evers, F Weigend, and M Koentopp. Conductance of molecular wires and transport calculations based on density-functional theory. Physical Review B, 69(23), 2004. [47] R. P. Feynman and F. L. Vernon. the theory of a general quantum system interacting with a linear dissipative system. Annals of Physics, 24(1):118{173, 1963. 120 Bibliography [48] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger. Dynamics of the dissipative 2-state system. Reviews of Modern Physics, 59(1):1{85, 1987. [49] B. L. Hu, J. P. Paz, and Y. H. Zhang. Quantum brownian-motion in a general environment - exact master equation with nonlocal dissipation and colored noise. Physical Review D, 45(8):2843{2861, 1992. [50] C.-H. Chou, T. Yu, and B. L. Hu. Exact master equation and quantum decoherence of two coupled harmonic oscillators in a general environ- ment. Physical Review E, 77(1):011112, 2008. [51] R. Kubo, M. Toda, and N. Hashitsume. Statistical physics II: nonequi- librium statistical mechanics. Springer, 2nd edition, 1991. [52] A. G. Redeld. On the theory of relaxation processes. IBM Journal of Research and Development, 1(1):19{31, 1957. [53] A. G. Redeld. The theory of relaxation processes. Adv. Magn. Reson., 1:1{32, 1965. [54] W. T. Pollard, A. K. Felts, and R. A. Friesner. The redeld equation in condensed-phase quantum dynamics. Advances in Chemical Physics, 93:77{134, 1996. [55] J. M. Jean, R. A. Friesner, and G. R. Fleming. Application of a multi- level redeld theory to electron-transfer in condensed phases. Journal of Chemical Physics, 96(8):5827{5842, 1992. [56] M. Michel, O. Hess, H. Wichterich, and J. Gemmer. Transport in open spin chains: A monte carlo wave-function approach. Physical Review B, 77(10):104303, 2008. [57] S. Andergassen, V. Meden, J. Splettstoesser, H. Schoeller, and M. R. Wegewijs. Charge transport through single molecules, quantum dots and quantum wires. Nanotechnology, 21(27), 2010. [58] M. Michel, M. Hartmann, J. Gemmer, and G. Mahler. Fourier's law conrmed for a class of small quantum systems. European Physical Journal B, 34(3):325{330, 2003. [59] T. Prosen and M. Znidaric. Matrix product simulations of non- equilibrium steady states of quantum spin chains. Journal of Statistical Mechanics-Theory and Experiment, page P02035, 2009. 121 Bibliography [60] K. Saito. Strong evidence of normal heat conduction in a one- dimensional quantum system. Europhysics Letters, 61(1):34{40, 2003. [61] W. Huisinga, L. Pesce, R. Koslo, and P. Saalfrank. Faber and newton polynomial integrators for open-system density matrix propagation. Journal of Chemical Physics, 110(12):5538{5547, 1999. [62] W. T. Pollard and R. A. Friesner. Solution of the redeld equation for the dissipative quantum dynamics of multilevel systems. Journal of Chemical Physics, 100(7):5054{5065, 1994. [63] H. Guo and R. Q. Chen. Short-time chebyshev propagator for the Liouville-Von Neumann equation. Journal of Chemical Physics, 110(14):6626{6634, 1999. [64] I. Kondov, U. Kleinekathofer, and M. Schreiber. Eciency of dierent numerical methods for solving redeld equations. Journal of Chemical Physics, 114(4):1497{1504, 2001. [65] H. Wichterich, M. J. Henrich, H.-P. Breuer, J. Gemmer, and M. Michel. Modeling heat transport through completely positive maps. Physical Review E, 76(3):031115, 2007. [66] H.-P. Breuer and F. Petruccione. Stochastic dynamics of quantum jumps. Physical Review E, 52(1):428{441, 1995. [67] P Saalfrank. Stochastic wave packet vs direct density matrix solu- tion of liouville-von neumann equations for photodesorption problems. Chemical Physics, 211(1-3):265{276, 1996. [68] D. Y. Petrina. Mathematical foundations of quantum statistical me- chanics. Kluwer Acad. Publ., 1995. [69] L. P. Kadano and P. C. Martin. Theory of many-particle systems II: superconductivity. Physical Review, 124(3):670{697, 1961. [70] C. W. Gardiner and P. Zoller. Quantum Noise. Springer, 2000. [71] K. E. Cahill and R. J. Glauber. Density operators for fermions. Phys- ical Review A, 59(2):1538{1555, 1999. [72] P. D. Drummond, P. Deuar, and J. F. Corney. Quantum many-body simulations using gaussian phase-space representations. Optics and Spectroscopy, 103(1):7{16, 2007. 122 Bibliography [73] P. Deuar and P. D. Drummond. First-principles quantum dynamics in interacting bose gases: I. the positive p representation. Journal of Physics A-Mathematical and General, 39(5):1163{1181, 2006. [74] P. Deuar and P. D. Drummond. First-principles quantum dynamics in interacting bose gases ii: stochastic gauges. Journal of Physics A-Mathematical and General, 39(11):2723{2755, 2006. [75] P. Deuar and P. D. Drummond. Correlations in a bec collision: First- principles quantum dynamics with 150 000 atoms. Physical Review Letters, 98(12):120402, 2007. [76] S. Lepri, R. Livi, and A. Politi. Thermal conduction in classical low- dimensional lattices. Physics Reports-Review Section of Physics Let- ters, 377(1):1{80, 2003. [77] K. Damle and S. Sachdev. Spin dynamics and transport in gapped one- dimensional heisenberg antiferromagnets at nonzero temperatures. Physical Review B, 57(14):8307{8339, 1998. [78] A. V. Sologubenko, S. M. Kazakov, H. R. Ott, T. Asano, and Y. Ajiro. Diusive energy transport in the s=1 haldane-gap spin chain compound agvp2s6. Journal of Magnetism and Magnetic Materials, 272(Suppl. 1):E689{E690, 2004. [79] A. V. Sologubenko, K. Gianno, H. R. Ott, A. Vietkine, and A. Revcolevschi. Heat transport by lattice and spin excitations in the spin-chain compounds srcuo2 and sr2cuo3. Physical Review B, 64(5):054412, 2001. [80] N. Hlubek, P. Ribeiro, R. Saint-Martin, A. Revcolevschi, G. Roth, G. Behr, B. Buechner, and C. Hess. Ballistic heat transport of quantum spin excitations as seen in srcuo2. Physical Review B, 81(2):020405, 2010. [81] K. R. Thurber, A. W. Hunt, T. Imai, and F. C. Chou. O-17 nmr study of q=0 spin excitations in a nearly ideal s=1/2 1d heisenberg antiferro- magnet, sr2cuo3, up to 800 k. Physical Review Letters, 87(24):247202, 2001. [82] F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brenig. Zero- frequency transport properties of one-dimensional spin-1/2 systems. Physical Review B, 68(13):134436, 2003. 123 Bibliography [83] K. Saito, S. Takesue, and S. Miyashita. Thermal conduction in a quantum system. Physical Review E, 54(3):2404{2408, 1996. [84] P. Jordan and E. Wigner. ber das paulische quivalenzverbot. Zeitschrift fr Physik A Hadrons and Nuclei, 47(9-10):631, 1928. [85] M. A. Schlosshauer. Decoherence and the Quantum-To-Classical Tran- sition. Springer, 2007. [86] J.-H. An and W.-M. Zhang. Non-markovian entanglement dynamics of noisy continuous-variable quantum channels. Physical Review A, 76(4):042127, 2007. [87] K. Blum. Density Matrix Theory and Applications: Physics of Atoms and Molecules. Springer, 2nd edition, 1996. [88] M. Michel, J. Gemmer, and G. Mahler. Heat conductivity in small quantum systems: Kubo formula in liouville space. European Physical Journal B, 42(4):555{559, 2004. [89] M. Michel, G. Mahler, and J. Gemmer. Fourier's law from schrodinger dynamics. Physical Review Letters, 95(18):180602, 2005. [90] G. Lindblad. Generators of quantum dynamical semigroups. commu- nications in mathematical physics, 48(2):119{130, 1976. [91] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan. Completely posi- tive dynamical semigroups of N-level systems. Journal of Mathematical Physics, 17(5):821{825, 1976. [92] U. Brandt and M. Moraweck. Rigorous results on the dc kubo- greenwood conductivity of disordered-systems. Journal of Physics C- Solid State Physics, 15(25):5255{5268, 1982. [93] D. J. Thouless and S. Kirkpatrick. Conductivity of the disordered linear-chain. Journal of Physics C-Solid State Physics, 14(3):235{245, 1981. [94] M. Kira and S. W. Koch. Many-body correlations and excitonic ef- fects in semiconductor spectroscopy. Progress in quantum electronics, 30(5):155{296, 2006. [95] P. C. Martin and J. Schwinger. Theory of many-particle systems i. Physical Review, 115(6):1342{1373, 1959. 124 Bibliography [96] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller. Cold bosonic atoms in optical lattices. Physical Review Letters, 81(15):3108{3111, 1998. [97] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch. Quantum phase transition from a super
uid to a mott insulator in a gas of ultracold atoms. NATURE, 415(6867):39{44, 2002. [98] H. Risken. The Fokker-Planck Equation: Methods of Solutions and Applications. Springer, 1989. [99] T. Prosen and B. Zunkovic. Exact solution of markovian master equations for quadratic fermi systems: thermal baths, open xy spin chains and non-equilibrium phase transition. New Journal of Physics, 12:025016, 2010. [100] S. R. Clark, J. Prior, M. J. Hartmann, D. Jaksch, and M. B. Plenio. Exact matrix product solutions in the heisenberg picture of an open quantum spin chain. New Journal of Physics, 12:025005, 2010. [101] J. Wu. Non-equilibrium stationary states from the equation of motion of open systems. New Journal of Physics, 12:083042, 2010. [102] L. I. Plimak, M. K. Olsen, M. Fleischhauer, and M. J. Collett. Beyond the fokker-planck equation: Stochastic simulation of complete wigner representation for the optical parametric oscillator. Europhysics Let- ters, 56(3):372{378, 2001. [103] Benoit Palmieri, Yuki Nagata, and Shaul Mukamel. Phase-space algo- rithm for simulating quantum nonlinear response functions of bosons using stochastic classical trajectories. Physical Review E, 82:046706, 2010. [104] S.E. Homann, J.F. Corney, and P. D. Drummond. Hybrid phase- space simulation method for interacting bose elds. Physical Review A, 78:013622, 2008. [105] J. et al. Hope. Xmds: extensible multi-dimensional simulator, http://xmds.sourceforge.net. [106] C. Mejia-Monasterio, T. Prosen, and G. Casati. Fourier's law in a quantum spin chain and the onset of quantum chaos. Europhysics Letters, 72(4):520{526, 2005. 125 Bibliography [107] M. Distasio and X. Zotos. Connection between low-energy eective- hamiltonians and energy-level statistics. Physical Review Letters, 74(11):2050{2053, 1995. [108] Y. Avishai, J. Richert, and R. Berkovits. Level statistics in a heisenberg chain with random magnetic eld. Physical Review B, 66(5):052416, 2002. [109] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong. Random-matrix physics - spectrum and strength
uc- tuations. Reviews of Modern Physics, 53(3):385{479, 1981. [110] T. C. Hsu and J. C. A. Dauriac. Level repulsion in integrable and almost-integrable quantum spin models. Physical Review B, 47(21):14291{14296, 1993. [111] K. Kudo and T. Deguchi. Unexpected non-wigner behavior in level- spacing distributions of next-nearest-neighbor coupled xxz spin chains. Physical Review B, 68(5):052510, 2003. [112] K. Kudo and T. Deguchi. Level statistics of xxz spin chains with discrete symmetries: Analysis through nite-size eects. Journal of the Physical Society of Japan, 74(7):1992{2000, 2005. [113] K. Kudo and T. Deguchi. Level statistics of xxz spin chains with a random magnetic eld. Physical Review B, 69(13):132404, 2004. [114] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin. Quantum Inverse Scattering Method and Correlation Functions. Cambridge University Press, 1993. [115] A. Klumper and K. Sakai. The thermal conductivity of the spin- 1/2 xxz chain at arbitrary temperature. Journal of Physics A- Mathematical and General, 35(9):2173{2182, 2002. [116] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin. Quantum Inverse Scattering Method and Correlation Functions. Cambridge University Press, 1993. [117] M. Shiroishi and M. Wadati. Integrable boundary conditions for the one-dimensional hubbard model. Journal of the Physical Society of Japan, 66(8):2288{2301, 1997. 126 [118] E. K. Sklyanin. Boundary-conditions for integrable quantum-systems. Journal of Physics A-Mathematical and General, 21(10):2375{2389, 1988. [119] T. Prosen. Open xxz spin chain: Nonequilibrium steady state and strict bound on ballistic transport. arXiv:1103.1350v2, 2011. [120] A. V. Sologubenko, E. Felder, K. Gianno, H. R. Ott, A. Vietkine, and A. Revcolevschi. Thermal conductivity and specic heat of the linear chain cuprate sr2cuo3: Evidence for thermal transport via spinons. Physical Review B, 62(10):R6108{R6111, 2000. [121] Y. Ando, J. Takeya, D. L. Sisson, S. G. Doettinger, I. Tanaka, R. S. Feigelson, and A. Kapitulnik. Thermal conductivity of the spin-peierls compound cugeo3. Physical Review B, 58(6):R2913{R2916, 1998. [122] Y. Saad and M.H. Schultz. Gmres - a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientic and Statistical Computing, 7(3):856{869, 1986. [123] R. V. Jensen and R. Shankar. Statistical behavior in deterministic quantum-systems with few degrees of freedom. Physical Review Let- ters, 54(17):1879{1882, 1985. [124] M. Schlosshauer. Decoherence, the measurement problem, and in- terpretations of quantum mechanics. Review of Modern Physics, 76(4):1267{1305, 2004. [125] F. M. Cucchietti, J. P. Paz, and W. H. Zurek. Decoherence from spin environments. Physical Review A, 72(5):052113, 2005. [126] G. Hu. Non-equilibrium Statistical Physics. unpublished, 2002. 127 Appendix A Perturbational decomposition of the operators m̂ In this appendix we will derive expressions of operators m̂ and ̂m in two ways: using the eigenvectors of HS = H0+ VS and expanding the operators perturbationally in orders of V0 based on expressions of the operators for the corresponding non-interaction H0. A.1 Expanded in eigenbasis of HS From its denition in Eq(4.3), in the representation of eigenvalues fEmg and eigenstates fjmig of HS , the operator m̂ can be written as hm jm̂jni = (m̂)mn = (c)mn X k jV k j2 Z 1 0 dei(EnEm!k;) h1 n (!k;)i = (c)mn Z d!D (!) jV (!) j2h1 n (!)i ( mn !) = (c)mn D (En Em) jV (En Em) j2h1 n (En Em)i; (A.1) where we have used R1 0 de i! = (!)+iP 1 ! and neglected the principal value part. We have also assumed that it is possible to perform a change of variable on V k such that it becomes V (knm), where kmn is dened by !kmn; = mn = Em En = nm, i.e. a bath mode resonant with this transition. This limits the possible forms of V k and !k;. For example, for each given energy mn, there should be a unique value of V knm . In all the work in this thesis, we take V k as a constant so this condition is satised. 128 A.2. Perturbational expansion starting from non-interacting systems D(!) is the bath's density of states. We arrive at m̂ = X m;n jmihnjhmjcjni (1 n ( nm))D ( nm) jV knm j2; (A.2a) ̂m = X m;n jmihnjhmjcyjnin ( mn)D ( mn) jV kmn j2: (A.2b) We furthermore set V knmD(!) as a constant and absorb it into 2. Sim- ilar expressions can be derived for bosonic systems with the Fermi-Diract distribution replaced by the Bose-Einstein distribution and (1 n ( nm)) replaced by (1 + n ( nm)). This approach involves a direct diagonalization of the isolated system HS . One can avoid the diagonalization by nding such operators m̂ perturbationally. A.2 Perturbational expansion starting from non-interacting systems To present the general idea concretely, let us work on a specic system, the Hamiltonian HS = H0 + VS dened in Eq. (4.1), which we rewrite here for convenience, HS = t N1X l=1 cyl cl+1 + c y l+1cl + V0 N1X l=1 cyl+1cl+1c y l cl = H0 + VS : (A.3) Next assuming V0 is small and treating it perturbationally, we want to ex- press the operator m̂ in terms of polynomials of fcmg operators and in various orders of V0. When V0 = 0 the system is a tight-binding open chain. The following basis transformation | discrete Fourier transform for open chains ck = 1p N NX l=1 sin kl N + 1 cl; (A.4) diagonalizes H0, H0 = NX k=1 kc y kck; (A.5) where k = 2t cos k N + 1 : (A.6) 129 A.2. Perturbational expansion starting from non-interacting systems Therefore, c (t) is a linear function of all cm, cl (t) c(0)l (t) = 2 N + 1 X km sin kl N + 1 sin km N + 1 eiktcm: (A.7) Hence, m̂ is also a linear combination of all cm operators. One can imagine that for small V0, m̂L should not be too far from such a linear combination. Let us call cl (t) when V0 = 0 as c (0) l (t). Treating this as the zeroth order form for the full dynamical cl (t), and expanding cl (t) = X n V n0 c (n) l (t) ; (A.8) we derive a perturbation equation of c (n) l (t): _c (n) l = it c (n) l1 + c (n) l+1 id(n1)l ; (A.9) where we dene, for non-negative integers n; n1; n2; n3, d (n) l = X n1;n2;n3P i ni=n n c (n1) l c y;(n2) l1 c (n3) l1 + c (n1) l c y;(n2) l+1 c (n3) l+1 o : (A.10) The solution of the above equation can be written as c (n) l (t) = i Z t 0 d 2 N + 1 X km sin kl N + 1 sin km N + 1 eik(t)d(n1)m () : (A.11) Here the initial condition c(n) (0) = 0(8n 1) is used. As the expression of c (0) l is given in Eq. (A.7), explicitly the expression of c (1) l is c (1) l (t) = 2 N + 1 4 X km kimi ei((k2)(k1)(k3))t ei(k)t (k1) + (k3) (k2) (k) sin kl N + 1 sin k1l N + 1 sin km N + 1 sin k1m1 N + 1 sin k2m2 N + 1 sin k3m3 N + 1 sin k2 (l 1) N + 1 sin k3 (l 1) N + 1 + sin k2 (l + 1) N + 1 sin k3 (l + 1) N + 1 cm1cym2cm3 : (A.12) 130 A.2. Perturbational expansion starting from non-interacting systems Similarly we can derive c (2) l , which will be shown explicitly later when we discuss terms proportional to V 20 in bath operators m̂. Plugging this gen- eral solution into Eq(4.3), and after straightforward but tedious algebra, we arrive at the decomposition of m̂ and ̂m in Eq(4.34) with expansion coecients dened as follows, D;m = 2 N + 1 X k sin kl N + 1 sin km N + 1 [1 n (k; T)] ; (A.13a) D;m = 2 N + 1 X k sin kl N + 1 sin km N + 1 n (k; T) ; (A.13b) for the zeroth order expansion, D;m1m2m3 = 2 N + 1 4 X k;m;k1;k2;k2 n ( (k) ; T) n ( (k1) + (k3) (k2) ; T) (k1) + (k3) (k2) (k) sin kl N + 1 sin k1m1 N + 1 sin k2m2 N + 1 sin k3m3 N + 1 sin km N + 1 sin k1m N + 1 sin k2 (m+ 1) N + 1 sin k3 (m+ 1) N + 1 + sin k2 (m 1) N + 1 sin k3 (m 1) N + 1 ; (A.14) 131 A.2. Perturbational expansion starting from non-interacting systems for the rst order expansion. For the second order expansion we make use of c (2) l and we have D;m1m2m3m4m5 = 2 N + 1 7 X k;m;k0;m0;k1;k2;k3;k4;k5 sin km N + 1 5i=0 sin kimi N + 1 sin kl N + 1 sin k1l N + 1 1 (k1) + (k3) (k2) (k0) n ( (k)) n ( (k5) (k0) (k4)) (k0) + (k4) (k5) + (k) n ( (k)) n ( (k1) + (k3) + (k5) (k2) (k4)) (k2) + (k4) (k1) (k3) (k5) + (k) sin k0l N + 1 X =1 sin k2 (l + ) N + 1 sin k3 (l + ) N + 1 X =1 sin k4 (l + ) N + 1 sin k5 (l + ) N + 1 + 1 (k2) + (k4) (k3) (k0) n ( (k)) n ( (k0) + (k1) + (k5)) (k) (k0) (k1) (k5) n ( (k)) n ( (k1) + (k3) + (k5) (k2) (k4)) (k2) + (k4) (k1) (k3) (k5) + (k) X =1 sin k0 (l + ) N + 1 sin k4 (l + ) N + 1 k5 (l + ) N + 1 sin k2l N + 1 sin k3l N + 1 + sin k2 (l + 2) N + 1 sin k3 (l + 2) N + 1 + 1 (k3) + (k5) (k4) (k0) n ( (k)) n ( (k1) (k0) (k2)) (k) + (k0) (k1) + (k2) n ( (k)) n ( (k1) + (k3) + (k5) (k2) (k4)) (k2) + (k4) (k1) (k3) (k5) + (k) X =1 sin k0 (l + ) N + 1 sin k2 (l + ) N + 1 sin k3 (l + ) N + 1 sin k4l N + 1 sin k5l N + 1 + sin k4 (l + 2) N + 1 sin k5 (l + 2) N + 1 : (A.15) Here we introduce a short-hand notation = 1 in (l + ) to stand for the neighbor sites of l. In the last expression we have denoted n ( (k) ; T) as 132 A.2. Perturbational expansion starting from non-interacting systems n ( (k)) to make it t into the page width. Therefore up the the order of V 20 , the bath operators can be written down as m̂ = X m D;mcm + V0 X m1m2m3 D;m1m2m3cm1c y m2cm3 +V 20 X m1m2m3m4m5 D;m1m2m3m4m5cm1c y m2cm3c y m4cm5 +O V 30 (A.16a) ̂m = X m D;mc y m V0 X m1m2m3 D;m1m2m3c y m3cm2c y m1 V 20 X m1m2m3m4m5 D;m1m2m3m4m5c y m5cm4c y m3cm2c y m1 +O V 30 : (A.16b) The same procedure is applicable for bosonic systems and similar ex- pressions can be derived. This has been used in Chapter 5 to calculate the operators m̂ and ̂m, and nd non-equilibrium stationary states for bosonic systems using the BBGKY-like method. 133 Appendix B Second-order cluster expansion and second-order equations of G1, G2 In this chapter we present explicitly forms of the second-order cluster ex- pansion and the second-order equations for G1, G2. B.1 Second-order cluster expansion First-order cluster expansion as in Eq. (4.31) expresses G2 in terms of G1 and neglecting G2. For convenience, we rewrite it here, G2 my; ny;m 0 ; n 0 = G1 my;m 0 G1 ny; n 0 +G1 my; n 0 G1 ny;m 0 +G2 my; ny;m 0 ; n 0 : (B.1) 134 B.2. Second-order equations of G1 and G2 A second-order cluster expansion similarly writes G3 in terms of G1; G2 and neglecting G3[94], G my1;m y 2;m y 3;m4;m5;m6 = G my1;m6 G my2;m5 G my3;m4 +G my1;m5 G my2;m4 G my3;m6 +G my1;m4 G my2;m6 G my3;m5 G my1;m6 G my2;m4 G my3;m5 G my1;m5 G my2;m6 G my3;m4 G my1;m4 G my2;m5 G my3;m6 (B.2a) +G my1;m6 G my2;m y 3;m4;m5 G my1;m5 G my2;m y 3;m4;m6 +G my1;m4 G my2;m y 3;m5;m6 +G my2;m5 G my1;m y 3;m4;m6 G my2;m4 G my1;m y 3;m5;m6 G my2;m6 G my1;m y 3;m4;m5 +G my3;m4 G my1;m y 2;m5;m6 G my3;m5 G my1;m y 2;m4;m6 +G my3;m6 G my1;m y 2;m4;m5 (B.2b) +G my1;m y 2;m y 3;m4;m5;m6 : (B.2c) The same equation has been written down in short-hand notation in Eq. (4.42). B.2 Second-order equations of G1 and G2 Using the perturbational expansion of bath operators Eq. (A.13), Eq. (A.14) and Eq. (A.15) together with cluster expansions Eq. (B.1) and Eq. (B.2), 135 B.2. Second-order equations of G1 and G2 Eq. (4.7) of G1 | the rst equation of the equation hierarchy becomes 0 = itG (m 1)y; n + itG (m+ 1)y; n itG my; n+ 1 itG my; n 1 +2 X l; h n D;l + D ;l G my; l +m D;l +D ;l G ly; n i +2V0 X ;m1;m2 (D;nm2m1 D;m1m2n)G my1;m2 m +2V0 X ;m1;m2 (D;mm1m2 D;m2m1m)G my1;m2 n +2V 20 X ;m1;m2;m3 [D;m3m2m1m3n +D;m3m3nm2m1 +D;nm2m3m3m1 D;m1m2m3m3n D;m3m3m1m2n D;m3m2nm3m1 ]G my1;m2 m +2V 20 X ;m1;m2;m3 [D;mm1m3m3m2 +D;m3m3mm1m2 +D;m3m1m2m3m D;m3m1mm3m2 D;m3m3m2m1m D;m2m1m3m3m]G my1;m2 n (B.3a) 2 X m D;n + n D ;m +2V0 X ;m1 (D;m1m1nm +D;m1m1mn) +2V 20 X ;m1;m2 (D;m1m1m2m2nm +D;m1m1m2m2mn) (B.3b) iV0G my; (n 1)y; n; n 1 + iV0G (m+ 1)y;my; (m+ 1); n iV0G (n+ 1)y;my; n+ 1; n + iV0G my; (m 1)y; n;m 1 +2V 20 X ;m1;m2;m3;m4 n G my1;m y 2;m3;m4 [m (D;m1m3nm4m2 +D;nm3m2m4m1 D;m2m4m1m3n) + n (D;m3m1mm2m4 D;mm1m3m2m4 D;m4m2m3m1m)]g (B.3c) 136 B.2. Second-order equations of G1 and G2 Using the condensed matrices notation introduced in Chapter 4, we can denote this equation as (1) H0 + 2 (1) B0 + 2V0 (1) B1 + 2V0 (1) B2 g1 = 20 + 2V01 + 2V 20 2 + iV0 + 2V 20 (g1g1) + iV0 + 2V 20 g2: (B.4) Similarly we can derive an equation for g2 truncated at the order of V 2 0 using these bath operators. It is a 10-pages long tedious equation so we here just write down its formal form in condensed matrices notation, (2) H0 + 2 (2) B0 + 2V0 (2) B1 + 2V0 (2) B2 (g2 + (g1g1)) = 2g1 + 2V0g1 + 2V 20 g1 + iV0 + 2V 20 (g1g1g1) + iV0 + 2V 20 (g1g2) : (B.5) In the calculation illustrated in the main text, we truncated the bath oper- ators at the order of V0 instead of V 2 0 . In this case, the second equation of 137 B.2. Second-order equations of G1 and G2 the hierarchy becomes, 0 = itG my; (n 1)y ;m0 ; n0 + itG my; (n+ 1)y ;m 0 ; n 0 +itG (m 1)y ; ny;m0 ; n0 + itG (m+ 1)y ; ny;m 0 ; n 0 itG my; ny;m 0 ; n 0 + 1 itG my; ny;m 0 ; n 0 1 itG my; ny; m 0 + 1 ; n 0 itGmy; ny;m0 1 ; n0 +iV0G my; ny;m 0 ; n 0 m0+1;n0 + m01;n0 iV0G my; ny;m 0 ; n 0 (m+1;n + m1;n) (B.6a) iV0 X l=m1;n1 G ly;my; ny; l;m 0 ; n 0 +iV0 X l=m01;n01 G ly;my; ny; l;m 0 ; n 0 (B.6b) 2 X n m h G ny;m 0 D;n0 Gny; n0 D;m0i n h G my;m 0 D;n0 Gmy; n0 D;m0i +m0 h G my; n 0 D;n Gny; n0 D;mi n0 h G my;m 0 D;n Gny;m0 D;mio (B.6c) +2 X ;l D;l + D;l h mG ly; ny;m 0 ; n 0 + nG my; ly;m 0 ; n 0 +m0G my; ny; l; n 0 + n0G my; ny;m 0 ; l i (B.6d) to be continued 138 B.2. Second-order equations of G1 and G2 continued +V0 2 X ;m1 n m0 h G my; n 0 D;m1m1n G ny; n 0 D;m1m1m +G my1; n 0 (D;nm1m D;mm1n) i n0 h G my;m 0 D;m1m1n G ny;m 0 D;m1m1m +G my1;m 0 (D;nm1m D;mm1n) i +m h G ny;m 0 D;m1m1n0 G ny; n 0 D;m1m1m0 +G ny;m1 D;n0m1m0 D;m0m1n0 i n h G my;m 0 D;m1m1n0 G my; n 0 D;m1m1m0 +G my;m1 D;n0m1m0 D;m0m1n0 io (B.6e) +V0 2 X ;m1;m2 n m0 h D;m1n0m2G my; ny;m1;m2 (D;m1m2m D;mm2m1)G my2; n y;m1; n 0 +(D;m1m2n D;nm2m1)G my2;m y;m1; n 0i n0 h D;m1m0m2G my; ny;m1;m2 (D;m1m2m D;mm2m1)G my2; n y;m1;m 0 +(D;m1m2n D;nm2m1)G my2;m y;m1;m 0i + m h D;m1nm2G my1;m y 2;m 0 ; n 0 + D;m1m2n0 D;n0m2m1 G ny;my1;m 0 ;m2 D;m1m2m0 D;m0m2m1 G ny;my1; n 0 ;m2 i n h D;m1mm2G my1;m y 2;m 0 ; n 0 + D;m1m2n0 D;n0m2m1 G my;my1;m 0 ;m2 D;m1m2m0 D;m0m2m1 G my;my1; n 0 ;m2 io (B.6f) Here G3 should be interpreted according to Eq. (4.42) but neglecting G3. 139 B.2. Second-order equations of G1 and G2 However, for certain systems with relatively stronger interaction, it might be necessary to use these new equations at the order of V 20 . 140 Appendix C Estimation of convergence of the two methods of truncating the BBGKY-like hierarchy In this chaprter we present our estimation of the leading order of G1 such as E;(1) 1 and C;(1) 1 , dened in Chapter 4. We will see that E;(1) 1 in fact involves E;(0) 2 , which in turn needs Eq. (4.8) | the equation of G2. A similar equation is needed for estimation of C;(1) 1 . C.1 E;(1) 1 from method 1 In order to estimate the accuracy of the rst-order form of this approxima- tion and also to get an estimate of accuracy for higher orders, we study the leading order of residues in terms of 2 and TT , both of which are assumed to be small in the following. Thus 2V0 V0, therefore we know that in Eq. (4.24) the gD term is smaller than the other g2 term so we drop it. Similarly since 2T T , we drop the 2T term in 2 in Eq(4.24), 2 = 20 (T ) + 2T;T ; (C.1) and keep only the large term, 20 (T ), which is independent of T . Here ;T denotes formally a derivative of T on | d dT . The general idea is then to write down respectively equations for gEx1 and g E;(1) 1 , and then compare the two to estimate E;(1) 1 = g E;(1) 1 gEx1 . In order to understand how an approximation to the next order improves the accuracy, we also want to compare E;(1) 1 to E;(0) 1 = g E;(0) 1 gEx1 , which is estimated in the same way from the dierence between the equations respectively for gEx1 and g E;(0) 1 . More generally we dene E;(0) n = g E;(0) n gExn and E;(1)n = gE;(1)n gExn for n-particle Green's functions gn. 141 C.1. E;(1) 1 from method 1 After dropping the gD term and the term which is proportional to 2T , and keeping only up to the linear order of T , gEx1 ; g E;(0) 1 and g E;(1) 1 respec- tively satisfy (1) 0 + (1) ;T T gEx1 = iV0g Ex 2 + 20; (C.2a) (1) 0 g E;(0) 1 = iV0g E;(0) 2 + 20; (C.2b) (1) 0 + (1) ;T T g E;(1) 1 = iV0g E;(0) 2 + 20; (C.2c) where (1) 0 + (1) ;T T is the zeroth and rst order in T from (1) of Eq(4.24). (1) ;T denotes formally a derivative of T on (1) | ddT (1). To consider E;(0) 1 , one may use Eq(C.2a) and Eq(C.2b), E;(0) 1 = T (1) 0 1 (1) ;T g Ex 1 + iV0 (1) 0 1 E;(0) 2 : (C.3) E;(1) 1 can be estimated from Eq(C.2a) and Eq(C.2c), E;(1) 1 = iV0 (1) 0 + (1) ;T T 1 E;(0) 2 ; (C.4) where E;(0) 2 is required. We nd that roughly speaking E;(1) 1 takes the second term of E;(0) 1 but drops the rst term. Therefore, next we only need to show that the second term, iV0 (1) 0 1 E;(0) 2 , is much smaller than the rst term, or equivalently, smaller than the whole E;(0) 1 . Estimation of E;(0) 2 involves the second equation of the hierarchy, i.e. equation ofG2, which can be derived from substituting Eq(4.18) into Eq(4.8): 142 C.1. E;(1) 1 from method 1 0 = it D cymc y ncm0 cn0+1 E + it D cymc y ncm0 cn01 E +it D cymc y ncm0+1cn0 E + it D cymc y ncm01cn0 E it D cymc y n1cm0 cn0 E it D cymc y n+1cm0 cn0 E it D cym1c y ncm0 cn0 E it D cym+1c y ncm0 cn0 E +2 X l; n n0d;l D cymc y ncm0 cl E + m0d;l D cymc y nclcn0 E +n d;l D cymc y l cm0 cn0 E + m d;l D cyl c y ncm0 cn0 E +nd ;l D cymc y l cm0 cn0 E + md ;l D cyl c y ncm0 cn0 E +n0 d;l D cymc y ncm0 cl E + m0 d;l D cymc y nclcn0 Eo +iV0 D cymc y ncm0 cn0 E m0+1;n0 + m01;n0 m+1;n m1;n (C.5a) iV0 X l=m1;n1 D cyl c y mc y nclcm0 cn0 E +iV0 X l=m 01;n01 D cyl c y mc y nclcm0 cn0 E (C.5b) 2 X h n d;m0 D cymcn0 E + m d;n0 D cyncm0 E +m0 d;n D cymcn0 E + n0 d;m D cyncm0 E n d;n0 D cymcm0 E m d;m0 D cyncn0 E n0 d;n D cymcm0 E m0 d;m D cyncn0 Ei (C.5c) 2V0 X n m0 D cymc y ncn0D E n0 D cymc y ncm0D E +m D cyncm0 cn0 D E n D cymcm0 cn0 D E + n0 D Dyc y mc y ncm0 E m0 D Dyc y mc y ncn0 E + n D Dyc y mcm0 cn0 E m D Dyc y ncm0 cn0 Eo : (C.5d) This can be written in a compact matrix form as (2)gEx2 = iV0g Ex 3 + 2gEx1 + 2V0g Ex D3; (C.6) 143 C.1. E;(1) 1 from method 1 where vector gEx2 is dened as follows, gEx2 = h G2 1y; 1y; 1; 1 ; G2 1y; 1y; 1; 2 ; ; G2 N y; N y; N;N iT : (C.7) Note that some of the elements of gEx2 are naturally zero but we still keep them in this vector. gEx3 , g Ex 1 and g Ex D3 come from ordering respectively Eq(C.5b), Eq(C.5c) and Eq(C.5d) in the same way as gEx2 . The matrix (2) can be read o from Eq(C.5a), for example assuming m;n;m 0 ; n 0 are all dierent and not equal to 1 or N , we have (2) mN3+nN2+m 0 N+n 0 ;mN3+nN2+m 0 N+n 0 +1 = it: (C.8) Since 2V0 V0 we drop the gExD3 term in the following estimation of the accuracy. Therefore, the exact solution and the equilibrium solution satisfy respectively, (2) 0 + (2) ;T T gEx2 = iV0g Ex 3 + 2gEx1 ; (C.9a) (2) 0 g E;(0) 2 = iV0g E;(0) 3 + 2g E;(0) 1 ; (C.9b) where (2) 0 and (2) ;T T stand for the zeroth and rst order in T in (2). Now E;(0) 2 can be analyzed: E;(0) 2 = iV0 (2) 0 1 E;(0) 3 + (2) 0 1 (2) ;T Tg Ex 2 + (2) 0 1 2 E;(0) 1 : (C.10) Here, in fact (1) and (2) have dierent dimensions. However, in this esti- mation of order of magnitudes, we ignore this dierence and furthermore the matrices are regarded as constants with order 1. For the moment, let us focus on the last term of E;(0) 2 . Recall that we want to compare iV0 (1) 0 1 E;(0) 2 against E;(0) 1 . Focusing only on the last term, we have iV0 E;(0) 2 iV02E;(0)1 ; (C.11) which is much smaller than E;(0) 1 as long as V0 2 1. All together we have arrived at: E;(0) 1 = T (1) 0 1 (1) ;T g Ex 1 + iV0 (1) 0 1 E;(0) 2 (C.12) 144 C.2. C;(1) 1 from method 2 while E;(1) 1 = (1) 0 1 V 20 (2)0 1E;(0)3 + iV0T (2)0 1 (2);T gEx2 +iV0 2 (2) 0 1 E;(0) 1 : (C.13) Most importantly here we see that E;(0) 1 is multiplied by a small number 2V0 and then becomes a part of E;(1) 1 . Judging from this it follows that, as long as 2V0 is a small number compared with t the method is very reasonable. As for the other two additional terms, they can be regarded as V 20 g Ex 3 + V0g Ex 2 T . Therefore, the limit of maximum value of V0, where this method is still applicable, is determined by gEx2 1 or gEx3 12 . Such limit could be much larger than 1 since roughly gExn = gExn | smaller for larger n. This explains why this method is applicable even for V0 larger than t, as conrmed from Fig.4.1. C.2 C;(1) 1 from method 2 In order to estimate the accuracy of this approximation, let us assume that 2 and V0 are small. We dene C;(0) n = g C;(0) n gExn and C;(1)n = gC;(1)n gExn . Again we start from the equations of the three: g C;(0) 1 , g C;(1) 1 and g Ex 1 , and then compare the three equations while ignoring certain higher-order terms such as terms which are proportional to 2V0. In this case, after dropping the (1) D term and terms which are propor- tional to 2V0, we nd that g Ex 1 , g C;0 1 and respectively g C;(1) 1 satisfy the following equations: (1) 0 g Ex 1 = iV0g Ex 2 + 20; (C.14a) (1) 0 g C;(0) 1 = 20; (C.14b) (1) 0 g C;(1) 1 = iV0g C;(1) 2 + 20: (C.14c) From Eq(C.14b) and Eq(C.14a) we nd C;(0) 1 = iV0 (1) 0 1 gEx2 : (C.15) 145 C.2. C;(1) 1 from method 2 Comparing Eq(C.14c) and Eq(C.14a), we get C;(1) 1 = iV0 (1) 0 1 g C;(1) 2 gEx2 : (C.16) Note that the magnitude of C;(1) 2 = g C;(1) 2 gEx2 is in fact smaller than the magnitude of C;(0) 2 = g C;(0) 2 gEx2 . Both involve the second equation of the hierarchy, i.e. the equation for G2. So we may analyze the latter to get an upper bound of the former. In this case, we substitute Eq(4.34) into Eq(4.8). The resulting equation will have the same structure as Eq(C.5) but every d;l and d;l are replaced respectively by D;m and D;m, and a similar substitution for D and D. Ignoring terms which are proportional to 2V0, it follows that g Ex 2 and g C;(0) 2 are respectively the solutions of (2) 0 g Ex 2 = iV0g Ex 3 + 2gEx1 ; (C.17a) (2) 0 g C;(0) 2 = 2g C;(0) 1 : (C.17b) Comparing these two equations, we nd that C;(0) 2 = iV0 (2) 0 1 gEx3 2 (2) 0 1 C;(0) 1 : (C.18) Focusing only on the last term, we have iV0 C;(0) 2 iV02C;(0)1 ; (C.19) which is much smaller than C;(0) 1 as long as V0 2 1. We arrived at: C;(0) 1 = iV0 (1) 0 1 gEx2 V0 gEx1 2 ; (C.20) and C;(1) 1 = (1) 0 1 iV 20 (2) 0 1 gEx3 + 2V0 (2) 0 1 C;(0) 1 V 20 gEx1 3 + 2V 20 gEx1 2 : (C.21) This agrees with the numerical results that C;(0) 1 is proportional to V0 while C;(1) 1 is proportional to V 2 0 . Most importantly, we see again that C;(0) 1 is multiplied by a small number 2V0 and then becomes a part of C;(1) 1 . Since 146 C.2. C;(1) 1 from method 2 roughly gExn = gExn, the other term, V 20 gEx3 V 20 gEx1 3, is also much smaller than C;(0) 1 V0 gEx1 2. However, for large enough V0 the other approximation used in this method, the perturbational expansion of the operators m̂, becomes invalid. Therefore, as long as V0 is a small number compared to t the method is very reasonable. It should be noted that this method is capable of dealing with large systems since it does not require a direct diagonalization of a 2N -dimension matrix. 147 Appendix D Coherent states and coherent-state representation To provide a convenient reference, in this appendix, we brie
y describe the general idea of coherent-state representation. More details can be found in for example Ref. [70, 71]. Here we will only discuss bosonic coherent states. Fermionic ones can be dened similarly with the usual complex numbers replaced by Grassman numbers[71]. In the xD.1, coherent states for single-mode bosonic systems is intro- duced and the simplest P -representation is discussed. In the xD.2, we discuss correspondingly the coherent-state representations for multi-mode bosonic systems. We also brie
y review the more general positive P -representation. D.1 Coherent states for systems with a single mode A single-mode bosonic system is dened by a Hamiltonian H a; ay , where a; ay are respectively annihilation and creation operators of a particle with the particular mode. Its Hilbert space, expanded in occupation-number representation, is fjnig, where n = 0; 1; 2; . A coherent state ji is dened to be an eigenstate of annihilation oper- ator, a ji = ji (D.1) with eigenvalue . Operator a is not a hermitian operator so is not a real but a complex number (c-number). Let us dene a displacement operator D () = ea ya; (D.2) which equals to D () = e jj2 2 ea y e a; (D.3) 148 D.1. Coherent states for systems with a single mode utilizing the Baker-Campbell-Hausdor formula[70]. From there one can verify that ji = D () j0i ; (D.4) where j0i is the vacuum state. This can be seen in the following. Let us express D () j0i in terms occupation-number states jni and note jni = (a y) n p n! j0i, D () j0i = e jj 2 2 ea y j0i = e jj2 2 1X n=0 n n! ay n j0i = e jj2 2 1X n=0 np n! jni : (D.5) Therefore aD () j0i = e jj 2 2 1X n=0 np (n 1)! jn 1i = e jj2 2 1X m=0 mp (m)! jmi = D () j0i : (D.6) Inner product between two coherent states ji and ji is, h j i = e jj 2 2 jj2 2 ; (D.7) so jh j ij2 = ejj2 ; (D.8) which is close to zero when is very dierent from . This implies that although no pair of the coherent states are orthogonal to each other but the overlap tends to be small if their corresponding eigenvalues are very dierent. This suggests that maybe the set of coherent states can be used as a basis to expand operators and vectors. Due to the non-orthogonal relation, such forms of basis are not unique. One particular convenient and important choice is I = 1 Z d2 ji hj : (D.9) 149 D.1. Coherent states for systems with a single mode This can be veried by considering arbitrary jni such that jni = 1 Z d2 ji h jni ; (D.10) which becomes jni = X m 1 Z d2 jmi hm ji h jni = X m jmi 1 Z d2ejj 2 mp (m)! ()np (n)! = X m jmi mn = jni : (D.11) Here we have used 1 Z d2ejj 2 mp (m)! ()np (n)! = mn; (D.12) which is 1 only when m = n and zero otherwise. With this basis, a density matrix (or a more general physical observable operator A) can always be expanded as[126] = 1 2 ZZ d2d2e jj2 2 jj2 2 R (; ) ji hj : (D.13) where R (; ) = e jj2 2 + jj2 2 hj ji : (D.14) This distribution function R (; ) holds generally for all and it is called the R representation. However, since the overcompleteness and the non- uniqueness of the coherent-state representation, there are other expansions, which hold not for all but a large class of density matrices[70]. For exam- ple, the P representation and the Positive P representation are respectively dened as = Z d2P (; ) ji hj ; (D.15) and = ZZ d2d2P (; ) e jj2 2 + jj2 2 ji hj : (D.16) 150 D.2. Coherent states for multi-mode systems Next we convert operator a and ay to be operators acting on such distribution function P , in for example, the P representation. For a, a = Z d2P (; ) a ji hj = Z d2P (; ) ji hj ; (D.17) so a! . For ay, noticing ay ji = dd ji ay = Z d2P (; ) d d ji hj = Z d2 d d P (; ) ji hj ; (D.18) so ay ! dd . Similarly we need to consider also a and ay. Summarize all those operators, we have a$ P (; ) ; ay $ P (; ) ; ay$ @ @ P (; ) ; a$ @ @ P (; ) : (D.19) D.2 Coherent states for multi-mode systems For a multi-mode (countably many, denoted as N) system, a coherent basis can be dened similarly as~E = D ~ j0i = e~~ay~~a j0i = ej~j22 e~~ay j0i ; (D.20) where ~ = (1; 2; ; N )T and ~ay = ay1; a y 2; ayN T . ayl is the creation operator at site l, l is a c number and j0i is the vacuum state. Coherent states are over complete. One convenient choice of representation is the P -representation[70], which decomposes the identity as I = 1 Z d2~ ~ED~ ; (D.21) 151 D.2. Coherent states for multi-mode systems and it maps a density matrix (a physical quantityA) to a function P ~; ~ (A ~; ~ ) such that = Z d2~P ~; ~ ~ED~ ; (D.22) hAi = Z d2~P ~; ~ A ~; ~ : (D.23) If A is in normal order of al; a y l . In this representation, operators a; ay become dierential operators on P ~; ~ , al$ lP ~; ~ ; ayl $ l P ~; ~ ; ayl $ l @ @l P ~; ~ ; al $ l @ @l P ~; ~ : (D.24) In this way, the operator-form QME can be mapped to a dierential equa- tion. 152 Appendix E Fokker-Planck equation and Langevin equation A standard Fokker-Planck equation with time-independent coecients reads, @ @t P ~; t = 24X j @ @j Aj ~ + 1 2 X i;j @2 @i@j Dij ~ 35P ~; t : (E.1) When there are higher-order derivative terms besides the above rst two terms, it is called a generalized Fokker-Planck equation. The Pawula theorem[98] states that in order to keep the distribution function (P ~; t ) non-negative, there has to be either only 2 terms or innity terms. In the later case, unless the equation is exactly solvable usually the generalized Fokker-Planck equa- tion is truncated rst and then solved. However, while truncation at the second order keeps the positiveness, a further truncation beyond the second order which destroys such positiveness in small regions, could result a more accurate distribution function[98]. In this appendix we discuss solutions of the Fokker-Planck equation and the generalized Fokker-Planck equation. Note here ~ is a vector of real random variables. It is not directly applicable for the generalized Fokker- Planck equation in c-numbers we derived in the previous chapters. E.1 Langevin equation for standard Fokker-Planck equations There are several approaches to solve the standard Fokker-Planck equation besides brute-force numerical solution in the whole space of n ~ o . Convert- ing it to a set of equivalent Langevin equation[98] with white noise and then numerically simulate the Langevin equation is a common and general method. A standard Langevin equation states, _i = hi ~ + gij ~ wj (t) ; (E.2) 153 E.1. Langevin equation for standard FPEs where we only consider the time-independent coecientsAj ~ andDij ~ , and ~w (t) is a white noise, hwi (t)i = 0; hwi (t)wj t0 i = ij t t0 : (E.3) Using Itô calculus[98], one can dene a correspondence between the standard Fokker-Planck equation and the Langevin equation. Assuming Dij ~ is positive semidenite such that there is a matrix G satisfying D ~ = G ~ Gy ~ ; (E.4) then gij ~ = Gij ~ ; (E.5) hi ~ = Ai ~ : (E.6) In numerical simulation, one works with directly di = hi ~ dt+ gij ~ dWj (t) ; (E.7) where hdWi (t)i = 0; hdWi (t) dWj t0 i = ijdtt;t0 : (E.8) Here t;t0 is the Kronecker not the Dirac . Another way to solve the Fokker-Planck equation is to convert it into equations of correlations by considering for example, hxli and hxlxmi etc.. Let us work out the rst derivative term for the above one-variable average, d dt hli = hl X j @ @j Aii = h X j @ @j lAii+ hAi X j @ @j li = hAli: (E.9) Here we have used the following property: for arbitrary F ~ , as long as limj!1 F ~ does not diverge, we have h @ @j F ~ i = Z d2~ @ @j h F ~ P ~ i = 0: (E.10) In this way, one can derive equations for correlations of any number of variables. However, then one has to solve such possibly coupled equations if many-variable correlations are of interest. 154 E.2. Analytical solution for the Ornstein-Uhlenbeck process E.2 Analytical solution for the Ornstein-Uhlenbeck process For a special case { the Ornstein-Uhlenbeck process, there is an analytical solution. An Ornstein-Uhlenbeck process has a constant diusion matrix D and a linear drift term ~A = ~. We further require is positive denite in or- der to have a non-divergent distribution function. Although time-dependent solution is also possible, here we focus on stationary solution, P = (2) N 2 (Det) 1 2 e 1 2 ~T 1~; (E.11) where = X l;m 1 l + m D vl jD j vm E julihumj ; (E.12) and = X l l ulihvl : (E.13) The last expression denes vl (juli) as the left (right) eigenvector of with eigenvalue l. Out of all above approaches, only the simulation via Langevin equation is generally applicable, even for generalized Fokker-Planck equa- tions. E.3 Stochastic dierence equation for truncated generalized Fokker-Planck equations A generalized Fokker-Planck equation with time-independent coecients truncated at the third-order derivatives reads, @ @t P ~ = 24X j @ @j Aj ~ + 1 2 X i;j @2 @i@j Dij ~ + 1 6 X i;j;k @3 @i@j@k Dijk ~ 35P ~ : (E.14) Now let us solutions of this generalized Fokker-Planck equation. For such generalized Fokker-Planck equation an equivalent stochastic dierence equation[102] for numerical molecular dynamical simulations can be found. In that case, usually the positive P -representation[72] or the gauge P -representation is 155 E.3. Stochastic dierence equation for truncated GFPEs used. The basic idea is to use proper random variables to mimic the eect of the third order term just like what was done for the second derivative terms. For example, @ @i Ai ! hii = AiT; (E.15) @2 @i@j Dij ! hiji = DijT; (E.16) @3 @i@j@k Dijk ! hijki = DijkT: (E.17) Therefore, there will be terms like ijk (T ) 1 3 appearing in the stochastic dierence equation to make sure the above average holds, i = hit+ gi (t) 1 2 + ki (T ) 1 3 : (E.18) just like gi happens to be gijwj such that Eq. (E.16) is satised, ki has to be dened properly to obey Eq. (E.17). Generally relation between ki and Dijk is still missing, but specic forms for example equations can be found in Ref.[102]. 156 Appendix F Solving the Redeld equation and the local-operator Lindblad equation of bosonic systems with BBGKY-like hierarchy The BBGKY-like method is based on the equation of motion of the Green's functions. The method starts from turning the Redeld equation, Eq. (4.2) into equation of motion of Green's functions, which is dened by single- particle Green's function G ly;n = D ayl an (1) E (G1 for short), two-article Green's function G iy; jy; l; n = D ayia y jalan (1) E (G2 for short) and so on. Here again the density matrix (1) refers to a stationary state and it is dened by L (1) = 0. Then generally for interacting systems the stationary equation becomes an equation hierarchy, in the rst order of which are equations for G1. But such G1 is coupled to G2. The next order in the hierarchy are equations for G2 and such G2 is coupled to G1 and G3, and so on. The hierarchy can be truncated and solved if one takes certain further approximations such as cluster expansion[101]. Since the dierence between the Redeld equation and the local-operator Lindblad equation is only at the bath operators, the following calculation is in fact valid for both equations. Of course, all those coecients D
;m should then be substituted properly for the local-operator Lindblad equation. The operator-form Redeld equation reads @ @t = LH0+ ULV + 2L0B+ 2UL1B; (F.1) 157 F.1. Exact solution on non-interacting systems where LH0 = i N1X l=1 h ayl al+1 + a y l+1al; i (F.2a) LV = i 2 NX l=1 h ayl al ayl al 1 ; i (F.2b) L(0)B = X
=L;R n D
;m h ay
; am i + D
;m h a
; a y m i + h:c: o (F.2c) L(1)B = X
;m1;m2;m3 D
;m1m2m3 nh ay
; a y m1am2am3 i + h a
; a y m3a y m2am1 i + h:c: o : (F.2d) Let us now transform the above Redeld equation in the operator form into equations of Green's functions by multiplying ayl an at the RHS of Eq.(3.30) and then performing a trace, 0 = ihayl an+1i ihayl an1i+ ihayl+1ani+ ihayl1ani (F.3a) +2 X
;m D
;m D
;m h n
hayl ami+ l
haymani i (F.3b) +2 X D
;nl
+ D
;ln (F.3c) +iUhayl ayl alani iUhayl aynanani (F.3d) +2U X
;mi D
;m1m2m3 h n haym1am2ilm3 + haym1am3ilm2 +l haym2am1inm3 + haym3am1inm2 i : (F.3e) F.1 Exact solution on non-interacting systems Note that when U = 0 only Eq. (F.3a), Eq. (F.3b) and Eq. (F.3c) are there so that this equation is closed. In this section we consider only the U = 0 case. This is an N2-dimension linear system for unknown G (0) 1 , where (0) refers to U = 0, H0g (0) 1 + 2B0g (0) 1 = 2(0): (F.4) Here we arrange all unknownG1 into a vector g1 = (G1 (1; 1) ; G1 (1; 2) ; ; )T . Elements of the two matrices and the vector can be read o from Eq. (4.7). 158 F.2. Approximate solution on interacting systems For example, (H0)lN+n;lN+n+1 = i; (F.5) and when l = 1 (or l = N , or n = 1; N), (0) 1N+n = 2 D1;n; (F.6) otherwise, it is zero. Solving this N2-dimension linear system we nd all G1. F.2 Approximate solution on interacting systems Note that when U 6= 0 this equation is not closed. This means that one need to consider the next order equation in terms of G2 together to solve the whole Eq.(3.30). One can derive such equation by multiplying ayia y jalan and then performing a trace. In fact, if one do so, one will nd that equations of G2 is coupled to G3 and so on[101]. Therefore, one has to solve the whole hierarchy to solve Eq.(3.30) exactly, which is as hard as dealing with directly Eq. (3.30). So we have to come up with some further approximation to truncate the hierarchy and then solve it. In the present work, we will truncate the hierarchy at the rst order by using the cluster expansion[101], hayl ayl alani = 2hayl alihayl ani+ o (U) : (F.7) So we will do not need explicitly the second order equations here. After taking the approximation of the cluster expansion, Eq. (F.3d) becomes, iU2hayl alihayl ani iU2hayl anihaynani: (F.8) Therefore, Eq. (F.3) becomes a closed equation, H0g (1) 1 + 2B0g (1) 1 + 2UB1g (1) 1 = 2(0) + iU g (1) 1 ; (F.9) where g (1) 1 lN+n = 2 g (1) 1 lN+n g (1) 1 nN+n 2 g (1) 1 lN+l g (1) 1 lN+n : (F.10) 159 F.2. Approximate solution on interacting systems This equation can be solved iteratively, g [n+1] 1 = H0 + 2 (0) D + 2U (1) D 1 h 2(0) + iU g [n] 1 i ; (F.11) with initial solution g[0] = g (0) 1 for the corresponding non-interacting system. Then the solution is g (1) 1 = limn!1 g [n] 1 ; (F.12) which in practice stops at large enough n such that g [n] 1 g[n1]1 is small enough. 160