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 interest when coupled to biased external baths. This requires solving a master equation for this open quantum system. Obtaining this solution is very challenging, 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 efficient methods to solve the Redfield equation — an example of such a master equation. The first is an opensystem Kubo formula, valid in the limit of weak bias. The second is a solution in terms of Green’s functions, based on a BBGKY (Bogoliubov–Born– Green–Kirkwood–Yvon)-like hierarchy. In the third, the Redfield 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 efficiency than direct methods such as numerical integration of the Redfield 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 first method, besides converting the task of solving the Redfield equation to solving the much easier Schrödinger’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 properties of spin systems using the framework of the Redfield equation and with direct methods. Normal energy and spin transport is found for integrable but interacting systems. This conflicts with the well-known conjecture linking anomalous conductivity to integrability, and it also contradicts the relationship, 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 first 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 secondorder 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Preface List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Dedication x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . 1.1 Two major theoretical approaches . . . . . . . . . . . 1.2 A big challenge for the system-baths scenario . . . . . 1.3 Our methods and their applications . . . . . . . . . . 1.3.1 Linear response theory . . . . . . . . . . . . . 1.3.2 The BBGKY-like hierarchy . . . . . . . . . . . 1.3.3 The GFPE in the coherent-state representation 1.3.4 Thermal conductance of spin chains . . . . . . 1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 5 7 7 8 9 9 11 Chapter 2 The projector technique and the Redfield equation 2.1 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The effective equation of motion for open systems . . . . . . 2.2.1 The projector operator technique . . . . . . . . . . . 2.2.2 An alternative derivation . . . . . . . . . . . . . . . . 2.2.3 QME in terms of creation and annihilation operators 2.3 Relaxation towards thermal equilibrium . . . . . . . . . . . . 2.3.1 One-site open system . . . . . . . . . . . . . . . . . . 2.3.2 Two-site open system . . . . . . . . . . . . . . . . . . 2.4 Various other extensions . . . . . . . . . . . . . . . . . . . . 12 12 12 14 19 20 22 22 24 28 iv Table of Contents 2.5 2.6 Evolution towards a NESS . . . . . . 2.5.1 NESSs from the RE . . . . . . 2.5.2 NESSs from the LOLE and the 2.5.3 The general Lindblad form . . Major challenge for large systems . . . . . . . . . . . . MQME . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 3 Linear Response Theory: Kubo formula for systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 3.2 General Linear Response Theory . . . . . . . . . . . . . 3.3 Limitations of the standard Kubo formula . . . . . . . 3.4 OKF: the Kubo formula for open finite-size systems . . 3.5 Comparison of the two OKFs . . . . . . . . . . . . . . . 3.6 Summary and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . open . . . . . . . . . . . . . . . . . . . . . Chapter 4 BBGKY-like hierarchy for the Redfield equation 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Derivation of the BBGKY-like hierarchy . . . . . . . . . . . 4.3 Non-equilibrium Wick’s theorem . . . . . . . . . . . . . . . . 4.4 Truncation of the hierarchy . . . . . . . . . . . . . . . . . . . 4.4.1 Method 1: equilibrium substitution . . . . . . . . . . 4.4.2 Method 2: first-level cluster expansion . . . . . . . . 4.4.3 Method 2: second-level cluster expansion . . . . . . . 4.5 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . Chapter 5 Solving the Redfield equation using coherent-state representations . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 The model and its effective equation of motion . . . . . . . . 5.3 Exact solution for non-interacting systems . . . . . . . . . . 5.4 Interacting systems: BBGKY-like hierarchy derived from the GFPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Approximate numerical solution for interacting systems . . . 5.6 Solving the LOLE via the GFPE method . . . . . . . . . . . 5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 6 Application: heat transport 6.1 Introduction . . . . . . . . . . . . . 6.2 The model and its Redfield equation 6.3 Numerical methods . . . . . . . . . in spin . . . . . . . . . . . . . . . chains . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 31 34 37 38 40 40 41 43 47 50 53 55 55 56 59 62 62 67 72 78 79 79 81 85 87 90 94 96 . 97 . 97 . 99 . 100 v Table of Contents 6.4 6.5 6.6 Definitions of the thermal current and local temperatures . . 101 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 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 systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 methods . . . . . . . . . . . . . . . . . . . E,(1) C.1 ∆1 from method 1 . . . . . . . . . C,(1) C.2 ∆1 from method 2 . . . . . . . . . of . . . . . . the BBGKY-like . . . . . . . . . . . 141 . . . . . . . . . . . 141 . . . . . . . . . . . 145 Appendix D Coherent states and coherent-state representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 . . . . . . . . . . . . E.1 Langevin equation for standard FPEs . . . . . . . . . . . . . E.2 Analytical solution for the Ornstein-Uhlenbeck process . . . E.3 Stochastic difference equation for truncated GFPEs . . . . . Appendix F Solving the RE and LOLE of bosonic with BBGKY-like hierarchy . . . . . . . . . . . . . . F.1 Exact solution on non-interacting systems . . . . . . F.2 Approximate solution on interacting systems . . . . 153 153 155 155 systems . . . . . 157 . . . . . 158 . . . . . 159 vi List of Figures 1.1 Sketch of a typical setup of heat conduction. 2.1 2.2 2.3 2.4 2.5 Sketch of a one-site system coupled to a heat bath. . . . . . . Sketch of a two-site system coupled to a heat bath. . . . . . . Sketch of a two-site system coupled to two heat baths. . . . . Difference between non-equilibrium and equilibrium states . . Stationary states from the local-operator Lindblad equation and the multi-mode quantum master equation compared against those from the Redfield equation . . . . . . . . . . . . . . . . 3.1 3.2 3.3 4.1 4.2 4.3 4.4 4.5 4.6 5.1 . . . . . . . . . Difference between ρ(1) (∞) and ρ(2) (∞) for different values of δ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Difference between ρ(1) (∞) and ρ(2) (∞) for different values of λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Difference between ρ(1) (∞) and ρ(2) (∞) for different values of V0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 23 24 27 33 36 51 52 53 E,(1) g1 compared with g1Ex for interacting systems at nonequilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . C,(1) g1 compared with g1Ex for non-equilibrium interacting systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C,(2) g1 compared with g1Ex for non-equilibrium interacting systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C,(2) g2 is compared with g2Ex for non-equilibrium interacting systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Particle currents calculated from the second-, first- and the zeroth-order forms are compared . . . . . . . . . . . . . . . . Heat currents calculated from the second-, first- and the zerothorder forms are compared . . . . . . . . . . . . . . . . . . . . Analytical solution compared against the BBGKY numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 71 74 75 76 77 88 vii List of Figures 5.2 5.3 6.1 6.2 6.3 6.4 6.5 7.1 The numerical simulation of the generalized Fokker-Planck equation for interacting systems . . . . . . . . . . . . . . . . . A similar application to the local-operator Lindblad equation 93 96 Typical temperature profiles for N = 10 spin chains . . . . . 104 Current v.s. system size . . . . . . . . . . . . . . . . . . . . . 105 Normalized profiles . . . . . . . . . . . . . . . . . . . . . . . . 106 Typical temperature profile for the cases of small Jz and small Bz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Temperature profiles on Ising spin chains and fermionic chains 109 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 efforts 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 benefited greatly from such discussion: Ian Affleck, 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 Affleck. 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 effective equation of motion of the central systems of interest, which is usually derived from the exact equation of motion for the composite system — the central system plus the bath(s). Although various approximations are involved in deriving such effective equations of motion, which are sometimes called quantum master equations, they are believed to be reasonably reliable since they are built on relatively firm grounds, i.e. an equation of motion derived from first-principles. Theories of transport [4, 5] in solid state physics deal with a specific type of open system, where the central system is coupled to multiple baths held at different temperatures, chemical potentials, etc. As a result, study of transport should be a branch of the theories of open quantum systems. In practice, different 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 difficulty 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 find new ways of making the theories of open quantum systems amenable for the study of transport in large(r) systems, by finding efficient calculation schemes. For non-interacting systems 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 efficient and accurate approximations, and test them. This chapter serves as an introduction to the subject and also an overview of this thesis. Sections §1.1 and §1.2 briefly review the theories of open sys1 1.1. Two major theoretical approaches tems while focusing on their applications to the study of transport. Section §1.3 lists our contributions, including the development of new methodologies and discovery of new physics from application of these methodologies to certain specific 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 theories of open quantum systems. Examples are given to illustrate the basic procedure. The resulting effective equations of motion are solved via direct methods, which have very poor efficiency. Then we present an approximate solution based on linear response theory in Chapter 3. This approximation has better efficiency 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 efficient and accurate approximations. One is based on a Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY)-like equation 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 flowing 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 Jh at heat transport of spin systems. The heat conductance κ = ∆T , 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)14 Cu24 O41 [27], which contain S = 12 two-leg ladders with an intra-chain coupling much stronger than the inter-chain coupling. By measuring the difference between the conductance parallel to the chain and that perpendicular to the chain, one finds 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 fields have also been reported [28, 29]. However, at present, for spin chains, there is no available data regarding temperature profiles and finite 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 framework based on first 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 formula [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 infinite 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-infinite non-interacting systems and the central system is usually finite. There are no explicit baths included 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 finite-size central system coupled explicitly to baths, whose degrees of freedom are then integrated out. There are no infinite leads in this picture. These three different 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 first two are somewhat equivalent if the same three-piece system is used when applying the Kubo formula. We should note, however, that the majority of studies of conductance via the Kubo formula are using an infinite homogeneous system [35, 42], or extrapolate from results on finite-size systems [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-defined piece-wise equilibrium state of the disconnected threepiece system. The left/right lead is assumed to be in its own equilibrium state correspondingly defined 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 system, 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 final states are qualitatively different. It is an essential assumption of the standard diagrammatic perturbation theory that the states at both ±∞ in time are the same equilibrium state (or the ground state if T = 0) of the unperturbed system. At this point a nonequilibrium Keldysh formalism [12, 18] is used instead. Conceptually, the condition that the two leads are half-infinite is essential. With any finitesize 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 infinite leads, theoretically there is still no guarantee that the proper stationary state can be calculated perturbationally [13]. In practice, for interacting systems, the non-equilibrium Green’s function method is often used together with the density functional theory [13, 16]. Additional 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 firmer grounds. They start from the first-principle equation of motion of the composite system, and arrive at an effective equation of motion of the central system. Various equations can be derived depending on the level of approximations. Using the influence functional theory [47, 48], exact time-dependent master equations have been derived for simple central systems [49]. Examples include a two-level system, one harmonic oscillator 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 generalization 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 localoperator Lindblad equation for reasons that will be clear later. It is in fact very hard to derive exact master equations from influence 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 Redfield equation [52–54]. Applications of the Redfield equation to transport studies are in an early but fast developing 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 matrix elements of the density matrix, which appears in the Redfield equation. Consider for example an N -site 1-d spinless fermionic system. The set of all single-particle Green’s functions has N 2 unknowns while a density matrix has 22N unknowns. Therefore, unless an efficient calculation scheme whose complexity scales as N 2 can be found, numerically the two approaches are not comparable at all. This explains why the non-equilibrium Green’s function method is much more widely used than other methods to study transport 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 efficient methods to solve the Redfield equation are found, as is the goal of this thesis. 1.2 A big challenge for the system-baths scenario Both the general Redfield equation and the general local-operator Lindblad equation can be cast into the form of a generalized Liouville equation, d ρ (t) = Lρ (t) , dt (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 defined by Lρ (∞) = 0, (1.2) where ρ (∞) 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 d2 ×d2 matrix. The task of finding ρ (∞) thus implies finding 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 find efficient methods to solve this equation. Before we present our work, we briefly review the currently available 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 ρ (∞). Thus, the first 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 shorttime Chebyshev polynomial propagator [63]. All these methods have been compared in Ref. [64]. As for their efficiency, they are all roughly on the same scale – their complexities scale roughly as d3 . A more efficient approach is the stochastic wave-function method [2, 65, 66]. The basic idea is to replace the effective equation of motion by a stochastic Schrödinger 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 results [67]. Breuer et al. showed that it can also be applied to the Redfield equation [65], but its efficiency is still not comparable with that of the nonequilibrium Green’s functions method. For the local-operator Lindblad equation, in fact, there is another quite efficient 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 Redfield equation. However, we should note that the same central system might behave differently under evolution described by the local-operator Lindblad equation versus the Redfield 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 efficient methods for the Redfield equation has certainly contributed to this state of affairs. With our newly developed methods, such a comparison can finally be undertaken. 1.3 Our methods and their applications In this section, we roughly summarize the ideas, performances and applications of our newly developed methods. Their details will be discussed in the rest of this thesis. We have developed three different methods, with different strenghts and weaknesses. 1.3.1 Linear response theory The first method is a linear response theory based on the Redfield 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 Redfield 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 describing the equilibrium conditions. Then, the non-equilibrium distribution can be treated as a perturbational response to the small part ∆L. This method has efficiency comparable to the stochastic wave-function method. Besides the gain in efficiency, which is not that remarkable, an important conceptual conclusion is drawn from this work. The standard Kubo formula [6], as we already knew, does produce a proper first-order correction to equilibrium states when an additional potential is included into the Hamiltonian. However, our work shows that it fails to describe the nonequilibrium stationary states. If infinite-size systems are used as the central system, there are ways to get around this conceptual problem. For non-equilibrium stationary states in finite-size systems, the standard Kubo formula must be replaced by a Kubo-like formula based on the Redfield equation – the standard Kubo formula is no longer valid. The key difference is that the standard Kubo formula uses Schrödinger’s equation and the system is treated as isolated while in the new Kubo-like formula, the Redfield 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 equilibrium 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 approximations 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 Redfield equation. We find again that for non-interacting systems, the equations decouple and a similar Wick’s theorem can be proved. All equations are coupled together for interacting systems. Two further approximations are introduced to truncate the hierarchy. The first 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 truncation of the hierarchy at the corresponding order. Our tests of the two approximations for small systems, where exact solutions 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 efficiency. 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 Redfield equation in the coherent-state representation of quantum states [70, 71]. Consider a simple harmonic oscillator for example. The set of all coherent states |ξi, which are eigenvectors of the annihilation operator a |ξi = ξ |ξi, 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 differential operators, and the Redfield equation becomes a stochastic differential 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 methods, now we need to solve an N -dimension stochastic differential 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 find 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 difficulties, its efficiency 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 = −κ∇T, (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 different 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 profile 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 investigated 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 different approaches, however, reach different 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 different 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 Redfield equation [36, 37, 83]. However, the spin systems studied via the Redfield equation are either non-interacting (and thus trivially 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 transformation the resulting fermionic systems are interacting. Anomalous transport 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 Redfield equation. We apply the Redfield equation to investigate the energy transport of spin 1/2 chains. We also find 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 Redfield 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 directly. 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 §2.2 contains a brief review of an open system’s dynamics based on the system-bath scenario, as a preparation for our method. In Section §2.3, we discuss how systems coupled 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 Redfield equations of both a single-mode central system and a two-site central system are solved to illustrate the general procedure. Besides the Redfield equation, there are other various forms of open-system quantum master equations. Usually they can be regarded as generalizations of the Redfield equation of a single-mode central system. Two such equations, which are widely used in literature, are presented in Section §2.4. Then, in Section §2.5 we generalize the Redfield equation to systems coupled to two or more baths of different 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 different. 2.2 The effective equation of motion for open systems In this section we derive a general effective equation of motion for a central system coupled to a single bath, starting from the Schrödinger’s equation for 12 2.2. The effective 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 Schrödinger picture, the Liouville equation for the evolution of the density matrix of the total closed quantum system S + B is: i ∂ ρ̂T (t) h i = Ĥ, ρ̂T (t) ≡ L̂ρ̂T (t) . (2.1) ∂t 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) S (2.3) ρ̂B (t) = tr (ρ̂T (t)) , 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) Ĥ0 (S, B) = ĤS (S) ⊗ Iˆ (B) + Iˆ (S) ⊗ ĤB (B) (2.5) where describes the system, respectively the bath, and X V̂ (S, B) = L̂i (S) ⊗ F̂i (B) (2.6) i describes their interactions, with L̂i (S) and F̂i (B) being system, respectively bath operators. Iˆ (B) and Iˆ (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 ⊗ Iˆ (B)]ρ̂T . (2.7) This can be trivially rewritten as: hAi = trS ÂS ρ̂ . (2.8) Therefore, it is convenient to try to find ρ̂ (t) directly, and avoid solving for the total density matrix ρ̂T (t), which is a much more complicated problem 13 2.2. The effective equation of motion for open systems due to the significant 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) = Lˆρ̂ (t) , (2.9) ∂t h i h i but now Lˆ is obviously neither Ĥ, · of the whole system nor ĤS , · 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 first 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 Schrödinger 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 first 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 define the projector operator on the Hilbert space of the total system S + B, P : H → H, such that eq P (ρ̂T ) , trB (ρ̂T ) ⊗ ρ̂eq B = ρ̂ ⊗ ρ̂B (2.10) where ρ̂eq B 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 effective equation of motion for open systems i.e. that PP = P. We denote Q = I −P, and rewrite the equation of motion for the total system, Eq. (2.1), as: ∂ ∂t P ρ̂T (t) = P L̂P ρ̂T (t) + P L̂Qρ̂T (t) (2.11) ∂ ∂t Qρ̂T (t) = QL̂P ρ̂T (t) + QL̂Qρ̂T (t) The second equation is integrated formally and its solution is substituted into the first equation, to obtain: Z t ∂ P ρ̂T (t) = P L̂P ρ̂T (t) + P L̂ dτ e(t−τ )QL̂ QL̂P ρ̂T (τ ) + P L̂etQL̂ Qρ̂T (0) . ∂t 0 (2.12) If we start from the reasonable initial condition ρ̂T (0) = ρ̂ (0) ⊗ ρ̂eq B, (2.13) Qρ̂T (0) = 0 (2.14) then and the last term vanishes identically. Let LS , LB , LV be the Liouville operators 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 ρ̂eq B = 0, we can rewrite the dynamic equation as ∂ eq eq ρ̂ (t) ⊗ ρ̂eq B = P L̂S ρ̂ (t) ⊗ ρ̂B + P L̂1 ρ̂ (t) ⊗ ρ̂B ∂t Z t +P L̂ dτ e(t−τ )QL̂ QL̂ ρ̂ (t) ⊗ ρ̂eq B . (2.16) (2.17) 0 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: hb| F̂i (B) |bi = 0. (2.18) In this case, we have: i X h eq i hb| L̂1 ρ̂ (t) ⊗ ρ̂eq |bi = hb| L̂ (S) ⊗ F̂ (B) , ρ̂ (t) ⊗ ρ̂ i i B B |bi = X i eq L̂i (S) ρ̂ (t) hb| F̂i (B) ρ̂eq |bi − ρ̂ (t) L̂ (S) hb| ρ̂ F̂ (B) |bi i B B i i =0 15 2.2. The effective equation of motion for open systems showing that indeed, P L̂V P = 0. This condition, which is satisfied for most types of couplings, is not however necessary. E h Ifi it is not obeyed, there will be P −βE D a term of the form i,b e Z b b F̂i b L̂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 find that: Z t ∂ B (t−τ )QL̂ ρ̂ (t) = L̂S ρ̂ (t) + tr P L̂ dτ e QL̂P ρ̂T (t) . (2.19) ∂t 0 Subject to the two assumptions mentioned above, this equation is exact. Clearly, all the effects of the bath on the system’s dynamics are contained in the second term on the right-hand side. This term can be further simplified. One can verify that QL̂P = QLˆ0 P + QLˆV P = QLˆV P = −P LˆV P + LˆV P = LˆV P (2.20) where the second assumption has been invoked again. Now we use the first essential assumption, replacing: e(t−τ )QL̂ → e(t−τ )QL0 . ˆ (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 simplify: e(t−τ )QL0 = e(t−τ )(I−P)L0 = e(t−τ )(L0 −P L0 ) = e(t−τ )L0 , ˆ ˆ ˆ ˆ ˆ (2.22) where we make use of P Lˆ0 = P Lˆ0 P + P Lˆ0 Q = 0. Therefore we obtain: Z t ∂ eq B (t−τ )Lˆ0 ˆ ˆ ρ̂ (t) = L̂S ρ̂ (t) + tr LV dτ e LV ρ̂ (τ ) ⊗ ρ̂B , (2.23) ∂t 0 where indeed the second term is of second order in the interaction V̂ . This can be further simplified in the interaction picture , where we define ρ̂I (t) = eiĤS t ρ̂ (t) e−iĤS t , (2.24a) L̂Ij (t) = eiĤS t L̂j e−iĤS t , (2.24b) F̂jI (t) = eiĤB t F̂j e−iĤB t . (2.24c) 16 2.2. The effective equation of motion for open systems In this case, the first term cancels out, and the dynamic equation becomes: i X Z t nh ∂ ρ̂I (t) =− dτ L̂Ij (t) , L̂Il (t − τ ) ρ̂I (t − τ ) Gjl (τ ) ∂t 0 jl h i o − L̂Ij (t) , ρ̂I (t − τ ) L̂Il (t − τ ) Glj (−τ ) (2.25) where Gjl (τ ) = trB F̂jI (τ ) F̂lI (0) ρ̂eq B . (2.26) are Green’s functions which depend only on the bath characteristics and can be calculated directly. Depending on the specific 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 → ∞, the Markov approximation is likely justified. Using the influence functional technique, 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 final approximation, we find the desired equation of motion for the reduced density matrix: i X Z t nh ∂ ρ̂I (t) dτ L̂Ij (t) , L̂Il (t − τ ) ρ̂I (t) Gjl (τ ) =− ∂t 0 jl h i o − L̂Ij (t) , ρ̂I (t) L̂Il (t − τ ) Glj (−τ ) . (2.28) It can be transformed back to the Schrödinger picture, i X Z t nh ∂ ρ̂ (t) = −i [HS , ρ̂ (t)] − dτ L̂j , L̂Il (−τ ) ρ̂ (t) Gjl (τ ) ∂t 0 jl h i o − L̂j , ρ̂ (t) L̂Il (−τ ) Glj (−τ ) . (2.29) 17 2.2. The effective equation of motion for open systems This equation can also be projected in the eigenbasis of HS , denoted as {|si} such that HS |si = Es |si. In this case, the effective equation of motion becomes X d I ρs0 s (t) = Rs0 smn (t) ρImn (t) , dt m,n (2.30) with Rs0 smn = −ei∆Es0 smn t # " X X + − − + δs0 m Γnkks , · δsn Γs0 kkm − Γnss0 m − Γnss0 m + (2.31) k k where S ∆Es0 smn = EsS0 − EsS + EnS − Em (2.32) and Γ+ mkln (t) = X ij Z hm| L̂i |ki hl| L̂j |ni t e−i(El S −E S n )τ G (τ ) dτ ij (2.33) 0 ∗ + and Γ− (t) = Γ (t) . The above effective equation of motion is mkln nlkm sometimes called a quantum master equation. If we focus only on the diagonal terms for ρIss (t), then we obtain the secular quantum master equation, which requires the so-called secular approximation [1] in order to be consistent with Fermi’s Golden Rule. The secular approximation is to ask that S EsS0 − EsS + EnS − Em = 0, (2.34) so as to eliminate the explicit time dependence of the right hand side of Eq. (2.31). Then, finally, X X d I ρss (t) = Wsm ρImm (t) − Wms ρIss (t) , dt m m (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 effective equation of motion Eq. (2.28). 18 2.2. The effective equation of motion for open systems 2.2.2 An alternative derivation: the non-correlation approximation This is a different 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 {|si} for EsS and {|bi} for EbB . Then, a complete basis for the total Hilbert space can be chosen as {|si |bi}. In the interaction picture ψ I (t) = ei(ĤS +ĤB )t |ψ (t)i, the time evolution of the operators describing the coupling is L̂Ii (t) = eiĤS t L̂i (t) e−iĤS t (2.36) F̂iI (t) = eiĤB t F̂i (t) e−iĤB t . (2.37) and The dynamical equation of the whole system is: i i ∂ ρ̂IT (t) h I = V̂ (t) , ρ̂IT (t) . ∂t (2.38) Integrating this formally and substituting the resulting expression for ρ̂I in the right-hand side of Eq. (2.38), we obtain: h i Z t h h ii ∂ ρ̂IT (t) I I = −i V̂ (t) , ρ̂T (0) − dτ V̂ I (t) , V̂ I (τ ) , ρ̂IT (τ ) . ∂t 0 (2.39) Integrating out the bath, this becomes: h i ∂ ρ̂I (t) = −itrB V̂ I (t) , ρ̂IT (0) ∂t Z t h h ii − dτ trB V̂ I (t) , V̂ I (τ ) , ρ̂IT (τ ) , (2.40) 0 which is still exact. Now we begin to introduce several approximations. First, we again assume hb| F̂i (B) |bi = 0, (2.41) so that the first 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 effective equation of motion for open systems such that the second term in Eq. (2.40) can be separated and simplified. This is the essential approximation. In fact this is stronger than the corresponding 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, 1 −β ĤB ρ̂IB (τ ) = ρ̂eq e . (2.43) B = ZB With these assumptions, we get a simplified dynamic equation for the reduced density matrix identical to Eq. (2.25). If we further apply the Markov condition, we arrive at the same effective 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 Schrödinger 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 simplified further. Let us now derive such a simplified 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†j , (2.44a) X (2.44b) HBν = ω ν (k) b†k,ν bk,ν , k HSB X † =λ Vk,ν aν bk,ν + h.c. . (2.44c) k,ν Here H aj , a†j can be any function of aj and a†j . λ refers to the strength of the coupling between the central system and the baths. It could be absorbed into the coupling coefficients 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 ∗ † Lk,ν = λVk,ν aν , Fk,ν = b†k,ν ; Lk,ν = λVk,ν aν , Fk,ν = bk,ν (1) (1) (2) (2) (2.45) 20 2.2. The effective equation of motion for open systems and note that the only non-vanishing bath Green’s functions are: (1) (2) iων (k)τ G1,2 δν,ν 0 , k,ν,k,ν 0 (τ ) = trB Fk,ν (τ ) Fk,ν 0 (0) ρB = hnk,ν i e (2.46a) (2) (1) −iων (k)τ G2,1 δν,ν 0 . k,ν,k,ν 0 (τ ) = trB Fk,ν (τ ) Fk,ν 0 (0) ρB = h1 − nk,ν i e (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 different sites are zero. By definition, −1 ων (k)−µν T ν hnk,ν i ≡ hn (ων (k) , Tν , µν )i = e −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) = −i[HS , ρ(t)] ∂t nh i o X ˆ ν ρ(t) + h.c. , −λ2 a†ν , m̂ν ρ(t) + aν , m̄ (2.48) ν ˆ defined as follows: with operators m̂ and m̄ Z ∞ X 2 |Vk | dτ aν (−τ ) e−iων (k)τ h1 − nk,ν i , m̂ν = k ˆν = m̄ 0 X k Z 2 |Vk | ∞ dτ a†ν (−τ ) eiων (k)τ hnk,ν i . (2.49a) (2.49b) 0 Here aν (τ ) = eiHS t aν e−iHS t is expressed in the Heisenberg picture. Note that we have also changed the integral limit from t to ∞. 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 annihilation operators is the central equation of this chapter and will be used later throughout the whole thesis. Eq. (2.48) is sometimes called the Redfield equation [52, 53, 87]. In some cases, it is convenient to use the eigenstates of HS as a basis, {|En i}. In that case, using element-wise products in this basis, the expresˆ can be rewritten as sions of the operators m̂ and m̄ ˆ ν = a†ν · Σ̄ν , m̂ν = aν · Σν , m̄ (2.50) 21 2.3. Relaxation towards thermal equilibrium where Σν = π X m,n Σ̄ν = π |mi hn| [1 − n (En − Em , Tν , µν )] Jν (En − Em ) , X |mi hn| nν (Em − En , Tν , µν ) Jν (Em − En ) , (2.51a) (2.51b) m,n and Jν () = X |Vk |2 δ( − ων (k)) (2.52) k is usually called the spectrum density function [36, 48] of the bath ν. 2.3 Relaxation towards thermal equilibrium We first 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 → ∞. We then briefly investigate a corresponding two-site system. 2.3.1 One-site open system Consider a one-state fermionic system in a strong magnetic field, describing for instance a single quantum dot, defined by a Hamiltonian: HS = a† a (2.53) where a, a† are fermionic operators. The system is coupled to a single bath (ν = 1 only), which is regarded as an ensemble of fermions with fixed T and µ, so that both particles and energy are exchanged. Using the general quantum master equation Eq. (2.48), we have h i nh i o ∂ ˆ ρ + h.c. , ρ = −i a† a, ρ − λ2 a† , m̂ρ + a, m̄ (2.54) ∂t where m̂ = πJ () [1 − n()] a, † ˆ = πJ () n()a , m̄ (2.55a) (2.55b) −1 −µ , the average number of and n() is shorthand for n(, T, µ) = e T − 1 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 defined accordingly for this specific system using Eq. (2.52). Now we can solve the dynamical equation and check the stationary solution. The density matrix of the above system can be written in the general form: 1 + p (t) 1 − p (t) ρ (t) = |0i h0| + |1i h1| + q (t) |0i h1| + q ∗ (t) |1i h0| , (2.56) 2 2 so that the trace is preserved. Then we get evolution equations for these parameters from the Redfield equation: d p (t) = 2 (1 − 2n()) − 2p (t) , dt d q (t) = −q (t) . dt This has a stationary solution, independent of the initial state: p (∞) = 1 − 2n() = eβ(−µ) − 1 , eβ(−µ) + 1 q (∞) = 0, (2.57a) (2.57b) (2.58a) (2.58b) which leads indeed to the expected equilibrium grand-canonical distribution e−β(−µ) |1i h1| (2.59) e−β(−µ) + 1 e−β(−µ) + 1 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 nonequilibrium stationary states, which depend on the details of the coupling between the central system and the baths. ρ(∞) = 1 |0i h0| + 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 = − a†1 a2 + a†2 a1 , (2.60) and we assume that only site 1 is coupled to the bath, so that: X V = λVk a1 b†k + λVk∗ a†1 bk . (2.61) k Using the general effective equation of motion Eq. (2.48), we arrive at nh i o d ˆ 1 ρ + h.c. , ρ = −i [HS , ρ] − λ2 a†1 , m̂1 ρ + a1 , m̄ dt (2.62) where m̂1 = πJ () [1 − n ()] πJ (−) [1 − n (−)] (a1 + a2 ) + (a1 − a2 ) , (2.63a) 2 2 πJ () n () πJ (−) n (−) † ˆ1 = m̄ a1 + a†2 + a†1 − a†2 , (2.63b) 2 2 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 ˆ 1 become when = 0, the operators m̂1 and m̄ m̂1 = J (0) [1 + n (0)] a1 , (2.64a) (0) n (0) a†1 , (2.64b) ˆ1 = J m̄ 24 2.3. Relaxation towards thermal equilibrium in agreement with results in the previous section. This shows that only if ˆ 1 depend only on a1 and a†1 . Otherwise, is zero, then the operators m̂1 and m̄ they also involve the creation and annihilation operators for other sites in the system (here, site 2). In fact, unless a†ν happens to be the creation operator ˆ ν involves other of an eigenmode of HS , which is generally not the case, m̄ † operators than aν , aν . We will revisit this issue in Section §2.4. The eigenmodes of this Hamiltonian are, HS = −a†+ a+ + a†− a− , (2.65) where a± = a1 ± a2 √ . 2 (2.66) In terms of occupation states of these eigenmodes |n+ n− i, eigenstates are |00i, |10i, |01i, |11i with eigenvalues respectively 1 = 0, 2 = −, 3 = , 4 = 0. A density matrix in this basis p11 p12 p13 p14 p21 p22 p23 p24 ρ= (2.67) p31 p32 p33 p34 p41 p42 p43 p44 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: 0 π 0 m̂1 = √ 2 0 0 ˆ1 m̄ 0 1 1 0 1 0 0 0 1 , a1 = √ 2 0 0 0 −1 0 0 0 0 0 1 −1 0 1 0 0 0 −1 , a2 = √ 2 0 0 0 −1 0 0 0 0 1 − n(−) 1 − n() 0 0 0 1 − n() , 0 0 −1 + n(−) 0 0 0 0 0 0 0 π n(−) 0 0 0 . =√ 0 0 0 2 n() 0 n() −n(−) 0 (2.69) (2.70) (2.71) (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 P = ΓP, dτ (2.73) where the matrix Γ has dimension 16×16, or generally 4N ×4N for an N -site spinless fermionic system. Fortunately, since we are only interested in equilibrium states in this section, we know that the off-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 ˆ 1 , we find a1 , a†1 , m̂1 , and m̄ −n(−) − n() 1 − n(−) 1 − n() 0 n(−) n(−) − n() − 1 0 1 − n() , = λ2 π n() 0 n() − n(−) − 1 1 − n(−) 0 n() n(−) n() + n(−) − 2 (2.74) Γdia 26 2.3. Relaxation towards thermal equilibrium 1 2 Right bath Left bath Figure 2.3: Sketch of a two-site system coupled to two heat baths. and the reduced equation becomes d dia P = Γdia P dia . dτ (2.75) A general time-dependent solution of Eq. (2.75) with initial condition P0dia is dia τ P dia (τ ) = eΓ P0dia . (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, Γdia P dia (∞) = 0. (2.77) We verified that Γdia indeed only has one zero eigenvalue, and that its corresponding eigenvector P dia (∞) leads to the proper thermal equilibrium state, P dia (∞) = 1 h β(µ+) β(µ−) β2µ iT 1, e ,e ,e , Z (2.78) where Z = 1 + eβ(µ+) + eβ(µ−) + eβ2µ . Up to this point, we have illustrated the general procedure for a system coupled to one bath and confirmed that its stationary state is the expected thermal equilibrium. A straightforward extension is to consider a system coupled to multiple baths, for example different sites in the system are locally coupled to different 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 Specifically, if as sketched in Fig.2.3, we couple site 2 to a bath 2 with the same bath parameters, then our effective equation of motion reads: nh i o d ˆ 1 ρ + h.c. ρ = −i [HS , ρ] − λ2 a†1 , m̂1 ρ + a1 , m̄ dt nh i o ˆ 2 ρ + h.c. , −λ2 a†2 , m̂2 ρ + a2 , m̄ (2.79) where the additional operators are defined as π [1 − n2 (−)] π [1 − n2 ()] (a1 + a2 ) − (a1 − a2 ) , (2.80a) 2 2 πn () πn2 (−) † 2 ˆ2 = m̄ a1 + a†2 − a†1 − a†2 . (2.80b) 2 2 Repeating the above calculation, we find the same stationary state. It is now natural to ask what kind of stationary state will emerge if the two baths have different parameters, for example different temperatures and/or chemical potentials. Ideally this should be the non-equilibrium stationary state which can be used to calculate transport properties.1 We will use examples of simple systems to investigate this question in §2.5. Before that, however, we have a subtle issue to settle, in the next section. m̂2 = 2.4 Various other extensions of the single-site quantum master equation There are several forms of the effective equation of motion used in the literature. Besides the multi-bath Redfield equation derived above, a localoperator 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: n ∂ ρ = −i [HS , ρ] − J () (1 − n()) a† aρ + ρa† a − 2aρa† ∂t o +n() aa† ρ + ρaa† − 2a† ρa . (2.81) 1 What 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 firstprinciple derivation: thermal equilibrium of the baths has been explicitly assumed. The question of establishing thermal equilibrium from only first 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, − X ν ∂ ρ = −i [HS , ρ] ∂t n Jν (ν ) (1 − nν (ν )) a†ν aν ρ + ρa†ν aν − 2aν ρa†ν o +nν (ν ) aν a†ν ρ + ρaν a†ν − 2a†ν ρaν , (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 simplified 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 specifically 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 aν 0 through formulas that depend on the particular Hs of the system. Therefore, the two equations are different. The question, then, is whether the stationary states of the two equations are qualitatively different, 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 Redfield equation for a singlesite system [70]. For multi-site systems, let us assume we can rewrite the Hamiltonian as follows: X (q) a†q aq , (2.83) HS = q We further assume that each eigenmode of the system is coupled to its own individual bath with temperature Tq : i Xh V̂ = λ Vqk aq ⊗ b†q,k + h.c. (2.84) k,q where creation and annihilation operators of bath q satisfy the fermionic commutation relation, {bq,k , b†q0 ,k0 } = δk,k0 δq,q0 . Note that this is different from our setup of Redfield 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 ˆ q involve only aq and a†q . Therefore, from Eq. (2.48) the effective m̂q and m̄ equation of motion reads h i X d ρ = −i (q) a†q aq , ρ dt q n X Jq ( (q)) (1 − nq ( (q))) a†q aq ρ + ρa†q aq − 2aq ρa†q − q o +nq ( (q)) aq a†q ρ + ρaq a†q − 2a†q ρaq , (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 quantum master 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 different 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ν 2 in Eq. (2.82) refers to a1 , a2 while aq in Eq. (2.85) refers to a± = a1√±a . In 2 the local-operator Lindblad equation a1 , a†1 does not mix with a2 , a†2 , while in the Redfield equation there are terms mixing them. Next, we investigate the difference between these three equations: the Redfield equation, the local-operator Lindblad equation and the multi-mode quantum master equation. We will argue firstly that while the local-operator Lindblad equation gives both wrong equilibrium states and wrong nonequilibrium 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 finding the “proper” non-equilibrium stationary state from the Redfield 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 first solve the Redfield 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 equation. 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 Redfield equation directly for large systems is a difficult task. 2.5.1 Non-equilibrium stationary states from the Redfield equation We begin from Eq. (2.79). Note that when the bath parameters are different: T1 6= T2 and/or µ1 6= µ2 , unlike in the case of thermal equilibrium states, the off-diagonal elements of the density matrix are no longer necessarily zero. As a result, we need to find 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 (∞) = 0. (2.86) This can be solved by finding 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 (∞) = V , (2.87) where V = [1, 0, · · · ]T and Γ̄ is defined byP replacing the first row of Γ by the normalization condition tr (ρ) = 1, i.e. j Pj∗d+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 efficiently. For simplicity of notation in later chapters, the above two equations defining the non-equilibrium stationary states are also written as Lρ (∞) = 0. (2.88) L̄ρ (∞) = ν. (2.89) and Here instead of superoperator Γ, Γ̄ and supervector P, V , we use directly the Liouvillian L, L̄ and matrix ρ, v, which are defined 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†ν in this basis of eigenstates; 2. Evaluate the correlation functions for every bath to calculate the matrices Σν and Σ̄ν defined in Eq. (2.51); ˆ ν as in 3. From these, derive the expressions of the operators m̂ν and m̄ 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 defined by Γ̄ in Eq. (2.89). Using this procedure we solve Eq. (2.79). Typical results are shown in the Fig. 2.4 for the specific parameters = 1, λ = 0.1, µ1 = 0 = µ2 , T1 = T 1 + 2δ , T2 = T 1 − 2δ , and T = 2.0. We have calculated the difference between the general non-equilibrium stationary state ρ (∞) and the thermal equilibrium ρEq (T ) at temperature T . The results are presented separately for off-diagonal elements and diagonal elements, v v uP 2 2 uP u u Eq Eq u j ρjj − ρjj u i6=j ρij − ρij ddiag = t (2.90) , dof f = t . P P 2 2 (ρ ) ij ij ij (ρij ) The charge current operator is defined as J = −ie a†1 a2 − a†2 a1 , (2.91) and we also separate its expectation value into contributions from the diagonal and off-diagonal parts, X X J diag = Jjj ρjj , J of f = Jji ρij . (2.92) j i6=j From Fig. 2.4(a) we see that compared to the diagonal part, the off-diagonal elements of the non-equilibrium stationary states are relatively small but 32 2.5. Evolution towards a NESS 0.2 diag J off J 0.005 J 0.15 d diag d off d 0.1 0.05 0 0 0.5 1 δ 1.5 2 0 0 0.5 1 δ 1.5 2 Figure 2.4: Density matrices for non-equilibrium stationary states are compared to those of equilibrium states. In (a) ddiag and dof f are plotted vs. the temperature bias δ. Diagonal elements of the density matrices for nonequilibrium stationary states are found to be quite different from those of equilibrium states. The difference between off-diagonal elements of nonequilibrium stationary states and equilibrium density matrices are much smaller. Note that the off-diagonal elements of equilibrium density matrices are zero always. (b) Charge currents, J diag and J of f , for the non-equilibrium stationary states, are plotted vs. δ. We see that only the off-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 contribute at all to the current, J diag = 0; the current comes entirely from the off-diagonal part. So the seemingly small change in the off-diagonal part makes a big qualitative difference: off-diagonal elements are essential to describe 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 Redfield 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 confirm these statements by presenting results for the simple example of a two-site fermionic system coupled to two baths hold at different temperatures but the same chemical potential. The Redfield 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 localoperator Lindblad equation and the multi-mode quantum master equation for this specific system are, respectively, n h i ∂ρL = −i[HS , ρL ] − λ2 π h1 + n (0, T1 , µ1 )i a†1 , a1 ρL ∂t h i + hn (0, T1 , µ1 )i a1 , a†1 ρL h i + h1 + n (0, T2 , µ2 )i a†2 , a2 ρL h i o + hn (0, T2 , µ2 )i a2 , a†2 ρL + h.c. . (2.93) 34 2.5. Evolution towards a NESS and n h i ∂ρM = −i[HS , ρM ] − λ2 π h1 + n (−, T, µ)i a†+ , a+ ρM ∂t h i + hn (−, T, µ)i a+ , a†+ ρM h i + h1 + n (, T, µ)i a†− , a− ρM h i o + hn (, T, µ)i a− , a†− ρM + h.c. , (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 find so in the following numerical study. Following a similar procedure like in the last section, we solve for the stationary states, ρL (∞) and ρM (∞). We compare them against the nonequilibrium stationary states of the Redfield equation, ρR (∞) using the following measure for distances: v 2 uP u A − ρB ρ u i,j ij ij u . (2.95) dA B =t 2 P B ij ρjj Typical results are shown in Fig. 2.5. We see that the difference between the multi-mode quantum master equation and the Redfield equation is relatively small compared to the difference between the local-operator Lindblad equation and the Redfield equation. At equilibrium (δ = 0), the multi-mode quantum master equation and the Redfield equation lead to the same stationary state, which is of course the thermal equilibrium. The local-operator Lindblad equation is quite different from the Redfield 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 related only to off-diagonal matrix elements. Fig. 2.5(b) shows that while 35 2.5. Evolution towards a NESS 0.02 M R d 0.5 J M J L J R L d R J d 0.4 0.01 0.3 0.2 0.1 0 0 0.5 1 δ 1.5 2 0 0 0.5 1 δ 1.5 2 Figure 2.5: ρL (∞) and ρM (∞) are compared with ρR (∞). In (a) we plot dL R and dM R vs. the temperature bias δ. We see that at equilibrium (δ = 0), the difference between ρM (∞) and ρR (∞) is zero but the difference between ρL (∞) and ρR (∞) is not. For all other cases, the differences are finite. However, as shown in (b), the steady-state charge current J L and J R , respectively calculated from ρL (∞) and ρR (∞), show qualitatively similar behavior. On the other hand, J M = 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 always zero, the currents from the local-operator Lindblad equation and the Redfield equation are qualitatively similar. The first fact confirms that the multi-mode quantum master equation captures properly only the equilibrium 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, transport properties. However, since the ratio between the two currents changes for different parameters, it is impossible to link J R quantitatively to J L . This qualitative similarity, if confirmed 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 efficient solutions for the Redfield equation. If confirmed, this similarity would be useful since it is much easier to solve the local-operator Lindblad equation than the Redfield 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 Redfield equation from now on. To summarize, up to this point we have derived the general effective equation of motion and the more specific Redfield equation, and we have also illustrated the general procedure of solving the Redfield equation by brute force direct methods. Before we develop more efficient 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) = e−iHt , 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) U (t1 + t2 ) = U (t1 ) U (t2 ) . (2.97) such that For a closed system, we have the additional condition: U (−t) = U −1 (t) (2.98) 37 2.6. Major challenge for large systems which insures reversibility. Lindblad discussed the generators of such trace-preserving positive semigroups [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 ρ (t) = Lρ (t) , dt (2.99) where generally: n −1 i h i h 1 X Lρ (t) = −i Ĥ, ρ + Aij Li , ρL†j + Li ρ, L†j , 2 h i 2 (2.100) ij where the operators {Li } form an orthonormal basis of the operators on the system’s Hilbert space, and the constants coefficient matrix A must be positive. Note that the local-operator Lindblad equation is of the Lindblad form, where Li = aα and Aij = Ai δij , 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 Redfield equation of Eq. (2.48) is not of the Lindblad 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 Redfield 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 infinity for bosonic systems. With such exponential growth, it is impossible to deal with large systems by direct methods, including 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 physical 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 effects, 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 efficient enough, or involve additional approximations. As we mentioned earlier, efficient stochastic wave-function methods [65] have been developed for both the local-operator Lindblad equation and the Redfield equation. This has better efficiency then the direct methods and it is capable of dealing with roughly N ∼ 20. For the local-operator Lindblad equation only, an efficient 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 quantitatively different from those of the Redfield equation. In the rest of this thesis, we will search for more efficient methods to calculate non-equilibrium stationary states of large systems for the Redfield equation, and discuss some preliminary applications. This chapter served as a foundation for our theory. None of it, except for the comparison between the Redfield equation, the local-operator Lindblad 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 finite size one-dimensional system whose ends are coupled to two baths held at different temperatures and/or different chemical potentials, leading to flow of energy or charge through the system. The question of interest to us is to find 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 Redfield 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 first-order corrections towards nonequilibrium stationary states for finite-size systems, starting from the corresponding equilibrium states. Two forms of such Kubo-like formulae are proposed. The first, while computationally less efficient, gives proper general density matrices. The more efficient second form produces correct average values only for certain special physical quantities. Conditions for these special cases will also be discussed. Besides efficiency, this work is also motived by considerations of the applicability of the standard Kubo formula to finite-size open systems. The standard Kubo formula has been widely used for both infinite [42] and finitesize systems [82], usually with periodic boundary conditions. Open boundary conditions have also been considered, leading to qualitatively different 40 3.2. General Linear Response Theory results [44]. For finite-size systems, the standard Kubo formula leads to a summation of many Dirac δ-function peaks at different energies. Various ways to get around or to smooth out such δ peaks have been proposed [92, 93]. However, there is no guarantee that correct physical quantities can be extracted via such procedures. More importantly, we show that the standard Kubo formula gives the first-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 infinite-size systems, not for finite-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 Redfield equation. This chapter is organized as follows. In Section §3.2 we briefly review the general linear response theory, and in Section §3.3 we discuss the applicability of the standard Kubo formula. In Section §3.4 we present the open Kubo formula based on the Redfield equation, and finally, in Section §3.5, for small systems, we compare the approximate solution from this open Kubo formula with the direct numerical solution of the Redfield equation, in order to gauge its validity. 3.2 General Linear Response Theory Consider a general linear equation of motion, ∂ρ (t) = (L0 + ∆L) ρ, ∂t (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 stationary states of the above equation, ρ (∞), given by: (L0 + ∆L) ρ (∞) = 0. (3.2) Let us assume that ρ0 is known and it satisfies L0 ρ0 = 0. (3.3) Then, to first order, δρ = ρ − ρ0 is the solution of ∂δρ (t) = L0 δρ + ∆Lρ0 , ∂t (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 Z t δρ (t) = dτ eL0 (t−τ ) ∆Lρ0 + δρ (t0 ) , (3.5) t0 Assuming that ∆L was turned on at t0 = −∞ when the initial state was ρ0 so that δρ (t0 ) = 0, we have the steady-state solution at t = 0 to be, Z 0 δρ (t = 0) = −∞ dte−L0 t ∆Lρ0 . Changing the variable t to −t, we then arrive at Z ∞ δρ (∞) = dteL0 t−ηt ∆Lρ0 . (3.6) (3.7) 0 Here we have inserted an infinitesimal 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 ρ (∞) = δρ (∞) + ρ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 δρ (∞) = − (L0 − η)−1 ∆Lρ0 . (3.8) This can also be derived from Eq. (3.2) since L0 is non-singular, 0 = (L0 + ∆L) (ρ0 + δρ (∞)) ⇒ δρ (∞) = −L−1 0 ∆Lρ0 . (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 L̄0 using the fact that ρ is normalized, tr (ρ) = 1 as we have done in rewriting Eq. (2.88) as Eq. (2.89). Following that we find a corresponding expression similar to Eq. (3.8), more specifically δρ (∞) = −L̄−1 0 ∆L̄ρ0 . (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) = LH ρ (t) = −i[H, ρ], ∂t (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 −βH e , Z (3.12) where Z = tr e−βH , β = 1/T . If H = HS + Vext , where HS is the isolated system and Vext is a static weak coupling to an external field, one can formally use the above general linear response theory to find a stationary solution ρ (∞) = ρ0 +δρ (∞) near a state ρ0 of the unperturbed system, LHS ρ0 = 0: Z ∞ dteL0 t−ηt LVext ρ0 . (3.13) δρ (∞) = 0 If ρ0 = ρeq (HS ) describes the unperturbed system in equilibrium, then this leads to the standard Kubo formula [6, 40]: Z ∞ e−βHS −ηt δρ (∞) = −i dte Vext (−t) , , (3.14) Z 0 where Vext (t) = eiHS t Vext e−iHS t . We can use the identity [Vext (−t), e−βHS ] = Rβ d −ie−βHS 0 dτ V̇ext (−t − iτ ), where V̇ext (−t − iτ ) = dt Vext (t) |t=−t−iτ , to rewrite: Z Z ∞ δρ (∞) = − 0 β dte−ηt dτ ρ0 V̇ext (−t − iτ ) . (3.15) 0 In terms of the eigenvectors of HS , HS |ni = n |ni, hm|V̇ext (t) |ni = i (m − n ) hm |Vext | ni ei(m −n )t , (3.16) leading to: δρ (∞) = X e−βm − e−βn hm |Vext | ni |mihn| . Z m − n − iη m,n (3.17) m 6=n 43 3.3. Limitations of the standard Kubo formula Incidentally, note that there is no contribution from states with n = m , for which hm|V̇ext (t) |ni = 0. This also follows directly from Eq. (3.14); if we P 0 + V ⊥ , where V 0 = write Vext = Vext ext ext m =n hm |Vext | ni |m ih n| commutes ⊥ 0 of with HS , then [Vext (−t), ρ0 ] = [Vext (−t), ρ0 ]. The “diagonal” part Vext Vext does not contribute to δρ (∞), and consequently has no influence on the static response functions. The lack of diagonal contribution is, however, very puzzling. The wellknown 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 πβ X e−βm ˆ 2. |hm|J|ni| (3.18) D= L m,n Z m =n It has contributions only from states with m = n . To understand the reason for this difference, 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 X † X nl nl+1 , (3.19) HS = −t cl cl+1 + h.c. + V0 l l where nl = c†l cl , plus a static electric potential X Vext = Vl n l (3.20) l induced by a homogeneous electric field E = −∇V . From the continuity equation: X X d Vl [Jl+1 (t) − Jl (t)], V̇ext (t) = Vl nl (t) = − (3.21) dt l l + where Jl = it c+ l+1 cl − cl cl+1 is the local current operator. The sum can be changed to X (3.22) − [Vl−1 Jl (t) − Vl Jl (t)] = −EJ(t), P l where J(t) = l Jl (t) is the total current operator. Using V̇ext (t) = −EJ(t) (or, on a lattice, E = Vl−1 − Vl ) in Eq. (3.15) gives Z ∞ Z β −ηt δρ (∞) = E dte dτ ρ0 J (−t − iτ ) . (3.23) 0 0 44 3.3. Limitations of the standard Kubo formula The dc conductivity is then Z ∞ dte σ= −ηt Z β dτ hJ (−t − iτ ) Ji , (3.24) 0 0 where hOi = T r[ρeq (HS )O]. This expression can be further simplified 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 change X X Vl Jl+1 (t) → Vl−1 Jl (t). (3.25) l l This is only justified for an infinite system (where boundary terms are presumed to be negligible), or a system with periodic boundary conditions and an external field with the same periodicity. This latter condition can only be achieved for a charge current in a finite system with periodic boundary conditions driven by a varying magnetic flux through the area enclosed by the system. For a finite-size system (even one with periodic boundary conditions) in a static applied electric field 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 ) + δρ (∞) is the first order perturbational expansion of ρeq (H), not the expected non-equilibrium stationary state. This statement is proved below, where for convenience, we 0 = 0. If it is not, we simply remove the assume that the diagonal part Vext 0 “diagonal” part Vext from Vext and add it to HS . Consider then the eigenstates of the full Hamiltonian, H|ñi = ˜n |ñi, to first order perturbation in Vext . Since hm |Vext | ni = 0 for all m = n , we can apply the first order perturbation theory for non-degenerate states to all the states, whether degenerate or not, to find 2 (3.26) ˜n = n + O Vext and |ñi = |ni + X m,m 6=n hm|Vext |ni 2 |mi + O Vext . n − m (3.27) 45 3.3. Limitations of the standard Kubo formula This immediately leads to ρeq (H) = X 1 2 e−β˜n |ñihñ| = ρeq (HS ) + δρ (∞) + O(Vext ), Z̃ n (3.28) where δρ (∞) is indeed given by Eq. (3.17). Intuitively this is easy to understand. In a sense, we found the original thermal equilibrium distribution ρeq (HS ) from the Liouville equation with H = HS . Then, it is not so surprising that a perturbation based on the same Liouville equation for the full H = HS + Vext leads to ρeq (H). This verifies our previous statement that the standard Kubo formula implies ρ (∞) → ρ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 finite system can be described by this approach. However, because we only keep the first-order perturbational correction, the situation is less clear-cut. In principle, it is not impossible for δρ of Eq. (3.17) to also be a first 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 finite size systems is problematic can also be seen from the following technical considerations. The spectrum of a finite-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 finite-size systems (for an infinite system, the integration over the continuous spectrum removes these singularities). Various techniques have been proposed to smooth out these singular contributions in order to extract some finite values, such as use of imaginary frequencies [93] or averaging σ (ω) over a small range of frequencies δω and then taking δω → 0 [92]. These different approaches may lead to different results. Moreover, the order in which the various limits are approached, e.g., taking η → 0 before L → ∞ or vice versa, also make a difference [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 finite-size systems is fraught with both conceptual and technical problems. Its solution gives the first order correction of the equilibrium state of the full Hamiltonian, instead of the expected non-equilibrium stationary state. The derived Drude weight formula is applicable only for infinite-size systems; applying 46 3.4. OKF: the Kubo formula for open finite-size systems it directly to finite-size systems involves dealing with unphysical δ-function peaks. In order to fix these difficulties, we now discuss how to calculate the non-equilibrium stationary states from the Redfield equation, which takes the connection between system and baths explicitly into consideration. 3.4 Open Kubo formula: linear response theory for open finite-size systems For concreteness, let us assume that the central system is coupled to thermal baths kept at temperatures TL/R = T ± ∆T 2 and investigate the thermal transport in the resulting steady state. If ∆T T , this will lead to a Kubolike formula which replaces Eq. (3.15). This approach can be generalized straightforwardly to derive a Kubo-like formula for charge transport. The general Redfield equation can be written as: ∂ρ(t) = [LH + LL (TL ) + LR (TR )] ρ(t), ∂t (3.29) where LH ρ = −i[H, ρ], just like for an isolated system, while LL/R are additional terms that describe the effects 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 field, Vext , such that H = 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 Redfield 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 Redfield equation to read: ∂ρ(t) = [LHS + LB (T ) + LP (∆T )] ρ(t) = Lρ(t), ∂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 finite-size systems We are interested in the t → ∞, stationary state solution ρ (∞) of the above equation, satisfying both Eq. (2.88) and Eq. (2.89), which are rewritten in the following form for convenience, Lρ (∞) = 0, (3.31) L̄ρ (∞) = ν. (3.32) and As we have shown in examples in the previous chapter, L has one nondegenerate zero eigenvalue and all its other (transient) eigenvalues have a negative real part, so ρ (∞) 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 efficient, can be obtained using linear response theory. The first step is to separate the Liouvillian L of Eq. (3.30) into a “large” plus a “small” part. There are two possible choices: ( (1) L0 = LHS + LB (T ) (3.33) ∆L(1) = LP (∆T ) ( or (2) L0 = LHS ∆L(2) = LB (T ) + LP (∆T ) (3.34) n o (1) (1) We begin with the first choice. Assume that L0 has eigenvalues L0,µ and left/right eigenvectors {|Lµ )}, {|Rµ )}. As discussed, the unique (zero (1) order in perturbation theory) steady-state solution of L0 ρ0 = 0 is ρ0 = (1) ρeq (HS ). The deviation δρK due to the perturbation ∆L(1) is obtained like in Eq. (3.13): XZ ∞ (1) (1) δρ (∞) = dteL0,µ t−ηt |Rµ ) (Lµ | ∆L(1) ρ0 =− X |Rµ ) (Lµ | µ (1) L0,µ −η µ 0 ∆L(1) ρ0 = − X |Rµ ) (Lµ | µ>0 (1) L0,µ ∆L(1) ρ0 . (3.35) (1) Note that the only divergent term, due to L0,0 = 0, disappears because (L0 | ∆L(1) ρ0 = (ρ0 | ∆L(1) ρ0 = 0. To see why, we start from Eq. (3.31), 48 3.4. OKF: the Kubo formula for open finite-size systems L(ρ0 + δρ) = 0, project it on (ρ0 | and keep terms only to the first order, (1) (1) to find 0 = (ρ0 | L0 + ∆L(1) (ρ0 + δρ) = (ρ0 | ∆L(1) ρ0 since L0 ρ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 equation [90] instead of the Redfield equation. (1) Eq. (3.35) is difficult to use in practice: finding all eigenstates of L0 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 defined by replacing the first row of the equation Lρ (∞) = 0 by T r (ρ (∞)) = 1, so that ν is a vector whose first 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̄ = (1) (1) L̄0 +∆L̄(1) . Again, the overbar shows that in matrix terms, L̄0 is obtained (1) from L0 by replacing its first row with T r (ρ (∞)) = 1, while ∆L̄(1) is obtained from ∆L(1) by replacing its first row with zeros. We find δ ρ̄(1) (∞) = −[L̄0 ]−1 ∆L̄(1) ρ0 . (1) (3.36) This is much more convenient because inverting the non-singular matrix (1) L̄0 is a much simpler task than finding all the eigenvalues and eigen(1) vectors of L0 . We have verified that both schemes produce identical results on systems on which both of them can be performed. We denote ρ0 + δ ρ̄(1) (∞) = ρ(1) (∞). (2) The second option is to take L0 = LHS and ∆L(2) = LB (T ) + LP (∆T ). (2) In this case, we can still choose the stationary solution associated with L0 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 (2) of HS satisfies L0 ρ0 = 0. Expanding the corresponding analogue of Eq. (3.14) in the eigenbasis of HS , we now find: δρ(2) (∞) = −i X hm|∆L(2) ρ0 |ni n,m m − n − iη |mihn|. (3.37) We call ρ0 + δρ(2) (∞) = ρ(2) (∞). Note that unlike ρ(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 infinite systems. This is not an accident. As discussed, the standard Kubo formula for infinite 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) (∞) or the singular solution ρ(2) (∞) – 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 N −1 X N −1 X c†l cl+1 + c†l+1 cl + V0 c†l+1 cl+1 c†l cl . l=1 (3.38) l=1 coupled to two heat baths, modeled as collections of fermions: X HB = ωk,α b†k,α bk,α , (3.39) k,α=L,R 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: X Vkα c†α bk,α + cα b†k,α , VSB = λ (3.40) k,α where the left (right) bath is coupled to the first (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 ± ∆T 2 . The corresponding Redfield equation reads from Eq. (2.48), −λ2 ∂ρ(t) = −i[HS , ρ(t)] ∂t i o X nh ˆ α ρ(t) + h.c. , c†α , m̂α ρ(t) + cα , m̄ (3.41) α=L,R ˆ α are defined by Eq. (2.49), which is rewritten where operators m̂α and m̄ here for convenience: X m̂α = π |mihn|hm|cα |ni (1 − nα (Ωnm )) , (3.42a) m,n ˆα = π m̄ X |mihn|hm|c†α |ninα (Ωmn ) . (3.42b) m,n 50 3.5. Comparison of the two OKFs 0.08 0.008 160 N (2) 3 0.006 120 0.06 -450 N ex N 80 J D 2.5 0.04 (b) 0 3.5 (a) 0 (2) 0.2 0.4 -900 0.004 D2 0.02 0 ex 40 D1 0 0.2 δT (λ=0.1, V0=0.2) 0.4 0 J (1) J (2) J 0.002 0 0 0.2 δT (λ=0.1, V0=0.2) 0.4 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) (∞) (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 ± ∆T 2 , ∆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) (∞), i = 1, 2 using the norm: sX 2 (3.43) Di = hn|ρex − ρ(i) (∞) |mi . n,m For the proper solution, this difference should be small but still finite, because 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 cutoff η. 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,2) (triangles, squares, respectively circles). Both N (1) and J (1) are very close to the exact values N ex , J ex . However, N (2) is very different from N ex while J (2) is close to J ex . These results confirm that ρ(1) (∞) of Eq. (3.36) is the proper Kubo solution. They also show that ρ(2) (∞) can also be used, but only for quantities A for which hm |A| ni = 0 51 3.5. Comparison of the two OKFs 6 (a) 0.003 D1 D2 5 3.5 -10 -20 N 0.002 ex N (1) -30 -40 2.5 3 J D (2) 3 4 (b) 0 N 0.0025 0 -50 0.05 0.1 0.15 0.2 0.0015 2 ex 0.001 J (1) J (2) J 1 0.0005 0 0 0.05 0.1 λ(δT=0.1, V0=0.2) 0.15 0.2 0 0 0.05 0.1 λ(δT=0.1, V0=0.2) 0.15 0.2 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) (∞) (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 fixed bias ∆T as a function of the strength of the system-bath coupling λ. It confirms again that ρ(1) (∞) is the proper approximation of ρex , but that for the charge current ρ(2) (∞) works well too. This also shows that the non-equilibrium stationary states depend on the coupling strength λ. This is not surprising for a finite-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 finite 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 fixed bias ∆T as a function of the strength of the interaction V0 . It confirms again that ρ(1) (∞) is the proper approximation of ρex . We also find that values of currents J ex are quite different between V0 = 0 and V0 = 1 for example. The charge currents calculated from both ρ(1,2) (∞) indeed capture the major part of such difference. 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 differences among all energy levels. For different values of V0 , such energy differences may either shift continuously or change qualitatively. 52 3.6. Summary and discussion 3 (a) 0.0008 (b) D1 D2 2.5 Nex N1 N2 4 0 2 0.00075 -4 0 0.2 0.4 0.6 0.8 1 J D -8 1.5 0.0007 1 ex J (1) J (2) J 0.5 0.00065 0 0 0.2 0.4 0.6 V0(λ=0.1, δT=0.1) 0.8 1 0 0.2 0.4 0.6 V0(λ=0.1, δT=0.1) 0.8 1 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 steadystate electric current calculated with ρex (triangles) and ρ(1,2) (∞) (squares and circles). The inset shows the corresponding particle numbers. 3.6 Summary and discussion To conclude, we first showed that the standard Kubo formula fails to provide an approximation for non-equilibrium stationary states, instead approximating the equilibrium state of the whole Hamiltonian. Secondly, applying the standard Kubo formula to finite 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 efficient but not always valid. In fact, the above results give us enough grounds to question the validity of the standard Kubo formula, where the effect of the baths are taken care of only via an additional potential. We believe that this is not enough. There are two different types of “driving forces” responsible for a steady-state current flow. The first force is the applied electric field, 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 ρ (∞). 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 effect 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 artificial 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 conceptual 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 efficient ways to solve the Redfield 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 Redfield equation 4.1 Introduction In this chapter we introduce a much more efficient way to solve the Redfield equation, based on a BBGKY-like hierarchy. Its efficiency depends on the level at which the hierarchy is truncated. At first order, where only singleparticle Green’s functions are involved, it leads to a linear system with N 2 unknowns for a N -site spinless fermionic system. This makes it possible to analyze systems larger than those that can be treated with the nonequilibrium 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 noninteracting 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 infinite. 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 Redfield equation. Not surprisingly, for non-interacting central systems we find that the hierarchy is decoupled, just like for systems in equilibrium. In fact, in Ref. [36], the Redfield 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 Redfield 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. Consistent 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 Redfield equation of non-interacting systems. This chapter is organized as follow. In Section §4.2, using a specific example, a BBGKY-like equation hierarchy is derived from the Redfield equation. In Section §4.3, non-interacting systems are discussed and a non-equilibrium Wick’s theorem is proved. In Section §4.4, we introduce the two methods to truncate and solve the coupled hierarchy for interacting systems, and compare them against exact results. We find that both methods are significantly more efficient than the direct methods. The first 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 Redfield equation describing an N -site chain of spinless fermions coupled with two fermionic baths. Our system of interest is described by: HS = −t N −1 X N −1 X c†l cl+1 + c†l+1 cl + V0 c†l+1 cl+1 c†l cl = H0 + VS , l=1 (4.1) l=1 and its Redfield equation, given by Eq. (3.41), is rewritten here for convenience: i o X nh ∂ρ(t) ˆ α ρ(t) + h.c. , = −i[HS , ρ(t)] − λ2 c†α , m̂α ρ(t) + cα , m̄ ∂t α=L,R (4.2) 56 4.2. Derivation of the BBGKY-like hierarchy ˆ are defined as follows: where the operators m̂ and m̄ Z ∞ X m̂α = |Vkα |2 dτ cα (−τ ) e−iωk,α τ h1 − n (ωk,α )i , k ˆα = m̄ 0 X k Z |Vkα |2 ∞ dτ c†α (−τ ) eiωk,α τ hn (ωk,α )i . (4.3a) (4.3b) 0 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) = e−iHS t is known, then so are cα (t) = U † (t) cα U (t) and thereˆ thus requires a full fore the operators m̂. Finding explicit forms for m̂ and m̄ 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 ρ (∞) 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 find: i E X nDh 0 = i h[A, H0 ]i + i h[A, VS ]i + λ2 A, c†α m̂α ˆα − + [A, cα ] m̄ D E m̂†α [A, cα ] α D h iEo ˆ †α A, c†α − m̄ , (4.4) ˆ †α ) is the hermitian conjugate of m̂α (m̄ ˆ α ). where hAi = tr (Aρ (∞)) and 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 equation. For example, the first and the second equation of the hierarchy for respectively single-particle Greens functions: D E G1 m† , n = c†m cn (4.5) and two-particle Greens functions: E D 0 0 G2 m† , n† , m , n = c†m c†n cm0 cn0 (4.6) can be derived by using A = c†m cn , respectively A = c†m c†n cm0 cn0 in Eq(4.4): 57 4.2. Derivation of the BBGKY-like hierarchy E D E D E D E D 0 = it c†m−1 cn + it c†m+1 cn − it c†m cn+1 − it c†m cn−1 D E D E −iV0 c†m c†n−1 cn cn−1 + iV0 c†m+1 c†m cm+1 cn D E D E −iV0 c†n+1 c†m cn+1 cn + iV0 c†m c†m−1 cn cm−1 E XD ˆ α + δnα m̄ ˆ †α c†m − δnα c†m m̂α − δmα m̂†α cn , −λ2 δmα cn m̄ (4.7a) (4.7b) (4.7c) α and D D D E E E 0 = it c†m−1 c†n cm0 cn0 + it c†m+1 c†n cm0 cn0 + it c†m c†n−1 cm0 cn0 D E D D E E +it c†m c†n+1 cm0 cn0 − it c†m c†n cm0 cn0 +1 − it c†m c†n cm0 cn0 −1 D E D E −it c†m c†n cm0 +1 cn0 − it c†m c†n cm0 −1 cn0 D E +iV0 c†m c†n cm0 cn0 δm0 +1,n0 + δm0 −1,n0 − δm+1,n − δm−1,n (4.8a) D E X −iV0 c†l c†m c†n cl cm0 cn0 l=m±1,n±1 +iV0 −λ2 D Xn D X 0 D c†l c†m c†n cl cm0 cn0 E (4.8b) 0 l=m ±1,n ±1 E D E δm0 α c†m c†n cn0 m̂α − δn0 α c†m c†n cm0 m̂α α D E D E ˆ α + δn0 α m̄ ˆ †α c†m c†n cm0 +δmα − δnα c†m cm0 cn0 m̄ D E D E D Eo ˆ †α c†m c†n cn0 + δnα m̂†α c†m cm0 cn0 − δmα m̂†α c†n cm0 cn0 −δm0 α m̄ . (4.8c) ˆα c†n cm0 cn0 m̄ E n o Note that since the set of all polynomials of cl , c†l form a complete basis of the space, operators m̂ are certainly functions of polynomials n operator o † of cl , cl . 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 Redfield 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 different orders can be found in Appendix C. Here we focus on the potential of this BBGKY-like formulation and discuss briefly 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 k1† , k2† , k3 , k4 = G1 k1† , k4 G1 k2† , k3 − G1 k1† , k3 G1 k2† , k4 , (4.9) where G1 and G2 are defined as D E G1 k1† , k2 = c†k1 ck2 , (4.10) respectively D E G2 k1† , k2† , k3 , k4 = c†k1 c†k2 ck3 ck4 . (4.11) In a general non-equilibrium state, one may have hck3 ck4 i 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 specific 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 different. Here it is more convenient to work in the momentum representation than the position representation, therefore we define: N 1 X klπ ck = √ sin cl . N +1 N l=1 (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 infinite-size systems eikl are replaced by the eigenstates sin Nklπ +1 , hence this transformation. In this representation, H0 = N X k c†k ck , (4.13) k=1 where k = −2t cos πk . N +1 (4.14) Starting from Eq(4.4) with the above H0 in momentum space and using † † † A = ck1 ck2 and A = ck1 ck2 ck3 ck4 , we find the equations for G1 k1† , k2 and respectively G2 k1† , k2† , k3 , k4 as follows: 0 = i (k2 − k1 ) G1 k1† , k2 −λ2 k1 πlα k2 πlα 2π X sin sin (n (k1 ) + n (k2 )) N +1 α N +1 N +1 k2 πlα 2π X kπlα sin +λ2 sin G1 k1† , k N +1 N +1 N +1 α,k k1 πlα kπlα + sin sin G 1 k † , k2 N +1 N +1 (4.15) 60 4.3. Non-equilibrium Wick’s theorem and 0 = i (k4 + k3 − k2 − k1 ) G2 k1† , k2† , k3 , k4 k1 πlα 2π X kπlα sin +λ2 sin G2 k † , k2† , k3 , k4 N +1 N +1 N +1 α,k 2π X k2 πlα kπlα +λ2 sin sin G2 k1† , k † , k3 , k4 N +1 N +1 N +1 α,k X k3 πlα kπlα † † 2 2π +λ sin sin G2 k1 , k2 , k, k4 N +1 N +1 N +1 α,k 2π X k4 πlα kπlα +λ2 sin sin G2 k1† , k2† , k3 , k N +1 N +1 N +1 α,k 2π X k2 πlα k4 πlα +λ2 sin sin G1 k1† , k3 (n (k2 ) + n (k4 )) N +1 α N +1 N +1 2π X k3 πlα k2 πlα −λ2 sin G1 k1† , k4 (n (k2 ) + n (k3 )) sin N +1 α N +1 N +1 2π X k1 πlα k4 πlα −λ2 sin sin G1 k2† , k3 (n (k1 ) + n (k4 )) N +1 α N +1 N +1 2π X k1 πlα k3 πlα +λ2 sin sin G1 k2† , k4 (n (k1 ) + n (k3 )) N +1 α N +1 N +1 (4.16) Each of these forms a closed linear system of equations and has a unique solution. Therefore, we only need to find one solution and that must be the unique solution. We first 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 equivalent 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 k2 , k4 together, we will have n G1 k2† , k4 i ( (k3 ) − (k1 )) G1 k1† , k3 k1 πlα k3 πlα 2π X sin sin (n (k1 ) + n (k3 )) −λ2 N +1 α N +1 N +1 X k3 πlα kπlα 2 2π +λ sin sin G1 k1† , k N +1 N +1 N +1 α,k k1 πlα kπlα † + sin sin G 1 k , k3 , (4.17) N +1 N +1 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 N 2 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 functions are coupled to one another. In order to solve Eq(4.7) explicitly, we first n have o to find explicit forms of the operators m̂ in terms of the operators † cl , cl . In Appendix A, we present an exact numerical calculation and a perturbational calculation for these operators. Correspondingly, based on these two methods of finding 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 first 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 efficient method, the nonequilibrium cluster expansion. The latter works only for relatively smaller V0 † but can be applied to much larger systems. The values G 1 m , n obtained Ex † from both methods will be compared against G1 m , 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 sacrificing its efficiency. 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 2N -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 n o of all polynomials of cl , c†l as: m̂α = X dα;l cl + V0 Dα (4.18a) d¯α;l c†l + V0 D̄α , (4.18b) l ˆα = m̄ X l where using the definition of inner product between super-vectors hhA|Bii = tr A† B , we have 1 dα;l = (N −1) tr c†l m̂α , (4.19a) 2 1 ˆα . (4.19b) d¯α;l = (N −1) tr cl m̄ 2 Here 2N −1 is a normalization constant to make dα,l = 1 when m̂α = cl . Also, dα;l and d¯α;l are the expansion coefficients at the linear order in cl and c†l ; operators V0 D and V0 D̄ are the remaining terms in the expansion of the ˆ , respectively. V0 is explicitly factorized out because of operators m̂ and m̄ the fact that when V0 = 0, this remaining part vanishes. With these expressions for m̂, Eq(4.7) becomes, D E D E D E D E 0 = it c†m−1 cn + it c†m+1 cn − it c†m cn+1 − it c†m cn−1 E XD δnα dα;l + d¯∗α;l c†m cl + δmα d¯α;l + d∗α;l c†l cn (4.20a) +λ2 l,α −λ2 X δmα d¯α;n + δnα d¯∗α;m (4.20b) α D E D E −iV0 c†m c†n−1 cn cn−1 + iV0 c†m+1 c†m cm+1 cn D E D E −iV0 c†n+1 c†m cn+1 cn + iV0 c†m c†m−1 cn cm−1 E XD −λ2 V0 δmα cn D̄α + δnα D̄α† c†m − δnα c†m Dα − δmα Dα† cn . (4.20c) (4.20d) α Note that every c0 , c†0 , cN +1 and c†N +1 that appears in the equation should be set to 0. First let us replace all G2 s in Eq(4.20c) by their values at equiE,(0) librium, denoted here as G2 where the superscripts indicate the use of thermal equilibrium (denoted by superscript E) as the zeroth order approximation (denoted by (0)) of the non-equilibrium G2 . Using the first term as an example, E,(0) G2 m† , n = tr c†m c†n−1 cn cn−1 ρeq (HS ) , (4.21) 63 4.4. Truncation of the hierarchy H where ρeq (HS ) = S 1 − kB T Ze . This requires knowing the eigenstates of HS . E,(0) Similarly one can define GD from Eq(4.20d), using the first term as an example, E,(0) GD m† , n = δmα tr cn D̄α ρeq (HS ) . (4.22) E,(1) Next let us calculate G1 from Eq(4.20), where the superscript (1) means that the approximate calculation takes care of the first equation of the hi E,(1) m† , n as a vector, erarchy, Eq(4.20). We organize all G1 h iT E,(1) g1 = G1 1† , 1 , G1 1† , 2 , · · · , G1 N † , N , (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 are N 2 such equations. After substituting E,(0) E,(0) G2 and GD for the exact but unknown G2 and GD , the whole set of E,(1) with dimension N 2 , Eq(4.20) for all m, n becomes a linear system for g1 E,(1) Γ(1) g1 E,(0) = iV0 g2 E,(0) + λ2 ν + λ2 V0 gD , (4.24) E,(1) . where the vector ν comes from ordering Eq(4.20b) in the same way as g1 E,(0) E,(0) The same holds for g2 and gD 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), X δmα d¯α;n + δnα d¯∗α;m , νmN +n = (4.25a) α (1) (4.25b) ΓmN +n,(m−1)N +n = it. E,(0) We can calculate single-particle equilibrium Green’s functions, G1 E,(0) organize them in the same way into a vector denoted as g1 . We define the distance between two vectors A and B as, qP 2 i |Ai − Bi | q dA = . B P 2 i |Bi | , and (4.26) In order to gauge the accuracy, we compare dE,(0) , the difference between E,(0) the zeroth order g1 and the exact solution g1Ex , and dE,(1) , the difference E,(1) between the first order solution above, g1 and the exact solution g1Ex . Here GEx m† , n = tr c†m cn ρex , where ρex is the exact solution from Eq(4.2). 64 4.4. Truncation of the hierarchy Results E,(1) First, we keep V0 = 0.2 as a constant, and check the accuracy of g1 for different 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 E,(1) accuracy of g1 for different 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, defined as: N −1 ie X J= (G (l + 1, l) − G (l, l + 1)) . N −1 (4.27) l=1 E,(0) E,(1) J E,(0) , J E,(1) , J Ex are calculated respectively from g1 , g1 and g1Ex . E,(1) From Fig.4.1(c) and (d) we see that in both cases, J is very close to Ex E,(0) the exact value, J , while J , 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 nonequilibrium 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†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 first 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 ∆T T , both of which 2 are assumed to be small in the following. Thus λ V0 V0 , therefore we know that gD is smaller than the g2 term so we drop it. Similarly, since λ2 ∆T ∆T , we drop the λ2 ∆T term in λ2 ν in Eq(4.24), λ2 ν = λ2 ν0 (T ) + λ2 ∆T ν,T , (4.28) and keep only the large term, λ2 ν0 (T ), which is independent of ∆T . Here ν,T d ν. The general idea is then to denotes formally a derivative of T on ν — dT E,(1) write down equations for g1Ex and g1 (as shown in Eq. (C.2) of Appendix E,(1) E,(1) C), and then compare them to get an estimate for ∆1 = g1 − g1Ex . In 65 4.4. Truncation of the hierarchy 0.02 (b) (a) 0.08 E,(0) d E,(1) d 0.015 0.06 d 0.01 0.04 0.005 0.02 0 0 0.2 0.4 0.6 ∆T/T(V0=0.2) 0.8 1 0 0 0.5 1 1.5 V0(∆T/T=0.4) 2 0.004 E,(0) 0.002 (c) (d) J Ex J E,(1) J 0.003 J 0.0015 0.002 0.001 0.001 0 0.0005 0 0.2 0.4 0.6 ∆T/T(V0=0.2) 0.8 0 0 0.5 1 1.5 V0(∆T/T=0.4) 2 E,(1) Figure 4.1: g1 is compared with g1Ex for interacting systems at nonE,(1) E,(0) equilibrium. dE,(1) (dE,(0) ) is the difference between g1 (g1 ) and g1Ex E,(1) E,(0) and d are plotted vs. respectively defined using Eq. (4.26). Both d ∆T in (a) and V0 in (b). In both cases, dE,(1) is much smaller than dE,(0) . In (c) and (d) J E,(0) and J E,(1) , plotted vs. respectively ∆T and V0 , are compared against J Ex . We see that J E,(0) is zero while J E,(1) is close to J Ex 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 E,(1) E,(0) E,(0) accuracy, we also want to compare ∆1 to ∆1 = g1 − g1Ex , which is estimated in the same way from the difference between the equations E,(0) respectively for g1Ex and g1 , see Appendix C for further details. Here we summarize the results, namely E,(0) (1) −1 (1) Ex (1) −1 E,(0) ∆1 = ∆T Γ0 Γ,T g1 + iV0 Γ0 ∆2 (4.29) and E,(0) E,(0) = (1) −1 Γ0 (2) −1 E,(0) −V02 Γ0 ∆3 (2) −1 (2) Ex (2) −1 E,(0) +iV0 ∆T Γ0 Γ,T g2 + iV0 λ2 Γ0 ∆1 . E,(1) ∆1 E,(1) (4.30) E,(1) Here ∆n = gn − gnEx and ∆n = gn − gnEx for general n-particle Green’s functions gn . We refer readers to Appendix C for definitions of all Γ matrices. Most importantly here we see that (from the last term in Eq. E,(0) (4.30)) ∆1 is multiplied by a small number λ2 V0 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 λ2 V0 t the method is very reasonable. As for the other two terms, they can be regarded as V02 g3Ex + V0 g2Ex ∆T . Therefore, the maximum value of V0 −1 −1 where this method is still accurate is determined by g2Ex or g3Ex 2 . n Ex Such limit could be much larger than 1 since roughly gn = g Ex — 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 diagonalization of HS , it is necessarily not efficient 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 — first level A different way to turn Eq(4.7) into a closed equation is to use the cluster expansion. Taking G2 as an example, this leads to: 0 0 0 0 G2 m† , n† , m , n = −G1 m† , m G1 n† , n 0 0 0 0 +G1 m† , n G1 n† , m + G2 m† , n† , m , n . (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 §4.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 first one. Here, we will study both levels, beginning first with the lowest level approximation, i.e. setting G2 → 0. However, the cluster expansion cannot be applied to the operators D defined in the previous section. Instead, n o we have to expand the operators † m̂ in higher-order polynomials of cl , cl . 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 (0) (1) cl (t) = cl (t) + V0 cl (t) + O V02 , (0) (4.33) (1) and then derive and solve the equations of motion of cl , cl from Heisenberg’s equation. In this way, one avoids the direct diagonalization of HS . This simplifies the calculation but its accuracy depends on the order of V0 at which the expansion stops. Stopping at the linear order of V0 is compatible with the cluster expansion for G2 (first 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 first order, the operators m̂ become X X m̂α = Dα;m cm + V0 Dα;m1 m2 m3 cm1 c†m2 cm3 + O V02 (4.34a) m ˆα = m̄ X m m1 m2 m3 D̄α;m c†m − V0 X Dα;m1 m2 m3 c†m3 cm2 c†m1 + O V02 , (4.34b) m1 m2 m3 where the definitions of Dα;m and Dα;m1 m2 m3 are given in Appendix A. 68 4.4. Truncation of the hierarchy With the first-order cluster expansion and the above expansion of operators m̂ plugged into Eq(4.7), we have 0 = itG1 (m − 1)† , n + itG1 (m + 1)† , n −itG1 m† , n + 1 − itG1 m† , n − 1 Xh δnα Dα;l + D̄∗α;l G1 m† , l +λ2 l,α +λ2 V0 i +δmα D̄α;l + D∗α;l G1 l† , n X (Dα;nm2 m1 − Dα;m1 m2 n ) G1 m†1 , m2 δmα (4.35a) α,m1 ,m2 +λ2 V0 (Dα;mm2 m1 − Dα;m1 m2 m ) G1 m†2 , m1 δnα X α,m1 ,m2 −λ2 +λ2 V0 X X δmα D̄α;n + δnα D̄∗α;m (4.35b) (4.35c) α (Dα;m1 m1 n δmα + Dα;m1 m1 m δnα ) (4.35d) −iV0 G m† , (n − 1)† , n, n − 1 +iV0 G (m + 1)† , m† , (m + 1), n −iV0 G (n + 1)† , m† , n + 1, n +iV0 G m† , (m − 1)† , n, m − 1 . (4.35e) α,m1 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 define a C,(1) vector g1 , where as before superscript C means cluster expansion and (1) symbolizes keeping only the first equation in the hierarchy. For simplicity C,(1) C,(1) we order Eq(4.35e) in the same way and denote it as g2 = Π g1 , where Π refers to the nonlinear function — summation of products — of C,(1) g1 in Eq(4.35e). Then the above equation can be written in matrix form as: (1) (1) C,(1) C,(1) Γ 0 + λ2 V0 Γ D g 1 = λ2 ν0 + λ2 V0 ν1 + iV0 g2 , (4.36) where the five terms are respectively the five 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 [n+1] g1 (1) (1) −1 [n] = Γ 0 + λ 2 V0 Γ D λ2 ν0 + λ2 V0 ν1 + iV0 Π g1 , (4.38) starting from the initial value [0] (1) −1 2 g1 = Γ0 λ ν0 . [0] (4.39) C,(0) Here we start from g1 = g1 , which is the exact solution of Eq(4.36) when C,(1) V0 = 0. Through the iteration defined above we get the solution g1 = [n] [n] [n−1] limn→∞ g1 . In practice we stop at large enough n such that g1 − g1 is small enough. Results Using again Eq. (4.26), we define dC,(0) as the relative distance between C,(0) C,(1) g1 and g1Ex , and dC,(1) as the relative distance between g1 and g1Ex . C,(1) First, we keep V0 = 0.2 as a constant, and check the accuracy of g1 for different 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 C,(1) g1 for different values of V0 . The worst case is d(1) ≈ 2% as shown in Fig.4.2(b). Overall, the difference 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 difference between J C,(1) and J C,(0) grows. We should also note that for larger V0 , J C,(1) starts to deviate from J Ex . 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 C,(1) C,(0) C,(0) = that λ2 and V0 are small. We define ∆n = gn − gnEx and ∆n C,(0) C,(1) C,(1) Ex gn − gn . Again we start from the equations of the three: g1 , g1 and g1Ex (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.04 0.15 (b) (a) C,(0) d C,(1) d 0.12 0.03 d 0.09 0.02 0.06 0.01 0.03 0 0 0.2 0.4 0.6 ∆T/T(V0=0.2) 0.8 0.004 0 0 0.2 0.4 0.6 V0(∆T/T=0.4) 0.8 1 0.0018 C,(0) J 0.003 0.0016 0.002 0.0014 0.001 0.0012 0 0 0.2 J Ex J C,(1) J (d) (c) 0.4 0.6 ∆T/T(V0=0.2) 0.8 0.001 0 0.2 0.4 0.6 V0(∆T/T=0.4) 0.8 1 C,(1) Figure 4.2: g1 is compared with g1Ex 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, dC,(1) is always much smaller than dC,(0) . In (c) and (d) J C,(0) and J C,(1) , plotted vs. respectively ∆T and V0 , are compared against J Ex . 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 J C,(1) are consistent with J Ex . From (d) we find that for relatively larger V0 , J C,(1) agrees much better with J Ex than J C,(0) . 71 4.4. Truncation of the hierarchy are proportional to λ2 V0 , see Appendix C for more details. We arrive at: C,(0) ∆1 2 (1) −1 Ex = −iV0 Γ0 g2 ∼ V0 g1Ex , and C,(1) ∆1 = (1) −1 Γ0 (2) −1 C,(0) ∆1 ∼ V02 g1Ex +λ2 V0 Γ0 (4.40) (2) −1 Ex g3 iV02 Γ0 3 2 + λ2 V02 g1Ex . (4.41) We refer the readers to Appendix C for definitions of all Γ matrices. C,(0) This estimate agrees with the numerical results that ∆1 is proporC,(1) 2 tional to V0 while ∆1 is proportional to V0 . Most importantly, we C,(0) see again that ∆1 is multiplied by a small number λ2 V0 and then ben C,(1) comes a part of ∆1 . Since roughly gnEx = g Ex , the other term, 3 2 C,(0) V02 g3Ex ∼ V02 g1Ex , is also much smaller than ∆1 ∼ V0 g1Ex . 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 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 N 2 . 4.4.3 Method 2: second-level cluster expansion In the previous section, we showed that even truncation at first level leads to quite accurate results, and a very efficient 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 system. While the first-level approximation is equivalent to the Hartree-Fock approximation, this second-level form goes beyond that. The cluster expansion of G3 [94] expresses it in terms of G1 , G2 and G3 (see Appendix B for the explicit expression): X X G m†1 , m†2 , m†3 , m4 , m5 , m6 = (−1)P G1 G1 G1 + (−1)P G1 G2 P +G3 P m†1 , m†2 , m†3 , m4 , m5 , m6 . (4.42) 72 4.4. Truncation of the hierarchy P P Here P (−1) G1 G1 G1 is a short-hand notation for the various ways of pairing {m†1 , m†2 , m†3 , m4 , m5 , m6 } into three groups using G1 and taking care of the anti-commutation relations, for example G1 m†1 , m6 G1 m†2 , m5 G1 m†3 , m4 . (4.43) P Similarly P (−1)P G1 G2 denotes all different ways of pairing these into two groups using G1 and G2 , for example (4.44) G1 m†1 , m6 G m†2 , m†3 , m4 , m5 . 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 V02 . However, including the V02 terms greatly complicates both equations (see Appendix B for details). Therefore, here we use only the expression of operators m̂α truncated at the first order of V0 . In this case, Eq. (4.7) becomes an equation very similar to Eq. (4.35) but with a different 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 = λ2 ν0 + λ2 V0 ν1 + iV0 π (g1 g1 ) + iV0 g2 . (4.45) Note that now we have an additional term iV0 g2 , so this is not a closed equation by itself. Here π (g1 g1 ) is the short-hand notation for the combinations of products between g1 and g1 the cluster expansions in Eq. (4.31). Meanwhile, Eq. (4.8) becomes, Γ(2) g2 = −Γ(2) π (g1 g1 ) + λ2 g1 + λ2 V0 g1 +iV0 π (g1 g1 g1 ) + iV0 π (g1 g2 ) , (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 π (g1 g2 ) is the short-hand notation for the combinations of products between g1 and g2 from Eq. (4.42). Same holds for π (g1 g1 g1 ). The matrices Γ(1) and Γ(2) come from the following terms (1) (1) (1) (4.47a) (2) (2) (2) (4.47b) Γ(1) = ΓH0 + λ2 ΓB0 + λ2 V0 ΓB1 , (2) Γ(2) = ΓH0 + iV0 ΓV + λ2 ΓB0 + λ2 V0 ΓB1 . 73 4.4. Truncation of the hierarchy 0.2 (a) (b) (0) d1 (1) d1 (2) 0.15 d1 0.1 d 0.1 0.05 0.05 0 0 0.2 0.4 0.6 ∆T/T(V0=0.8) 0.8 1 0 0 0.2 0.4 0.6 V0(∆T/T=0.4) 0.8 1 Figure 4.3: G1 calculated with different levels of the approximations are (2) (1) (0) compared. d1 is plotted in addition to the previous shown d1 and d1 . (2) (1) We see that the distance d1 is smaller than, but rather close to d1 , while (0) both are much smaller than d1 . The various Γs are obtained by considering only the term indicated by their (1) subscripts. For example, ΓH0 is the one derived from H0 only, without using (1) (1) the potential V and bath operators m̂α , while λ2 ΓB0 (λ2 ΓB1 ) is related to the component of the bath operators that is of the zeroth (first) 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 finite 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 first-level truncation. (2) First, we define d1 to be the difference between the G1 calculated from this second-order cluster expansion and the exact G1 , similarly as dC,(1) , (2) (1) which is now denoted as d1 . In Fig.4.3 we find that d1 is smaller than (1) but very close to d1 . This means that the second-order cluster expansion has better accuracy, but quantitatively does not lead to a significant improvement of G1 . 74 4.4. Truncation of the hierarchy 0.3 0.4 (b) (a) (0) d2 0.25 (1) d2 (2) 0.3 d2 d 0.2 0.2 0.15 0.1 0.1 0.05 0 0 0.2 0.4 0.6 ∆T/T(V0=0.8) 0.8 1 0 0 0.2 0.4 0.6 V0(∆T/T=0.4) 0.8 1 Figure 4.4: G2 calculated with different levels of the approximations are (2) (2) compared. d2 is plotted v.s. ∆T and V0 . We see that d2 is significantly (1) smaller than d2 . (2) Second, in the same way, we define d2 and compare the differences between G2 s calculated from the various truncations to the exact G2 . From Fig.4.4 we find that the second order truncation leads to a sizable improvement in G2 . To find how significant is this difference for physical quantities, we compare the currents calculated from the approximations and exact calculation. In order to enhance these differences, we increase the coupling constant to λ = 0.5 in the following figures. 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 J Ex is much better than the one between J (1) and J Ex , 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.03 0.08 (a) (b) Ex J (0) J (1) J (2) J 0.06 J 0.025 0.04 0.02 0.02 0 0 0.2 0.4 0.6 ∆T/T(V0=0.8) 0.8 1 0 0.2 0.4 0.6 V0(∆T/T=0.4) 0.8 1 Figure 4.5: Particle currents, J (2) , J (1) , J (0) , calculated respectively from the second-, first- and the zeroth-order truncations are compared against J Ex , which is found from the exact diagonalization method. For all the range of V0 shown in this figure, J (2) is very close to J Ex , while J (1) is accurate only for small V0 . 76 4.4. Truncation of the hierarchy 0.01 (a) (b) 0.025 0.009 Ex Jh (0) Jh (1) 0.02 Jh 0.008 (2) J Jh 0.015 0.007 0.01 0.006 0.005 0.005 0 0 0.2 0.4 0.6 ∆T/T(V0=0.8) 0.8 1 0.004 (2) 0 0.2 (1) (0) 0.4 0.6 V0(∆T/T=0.4) 0.8 1 Figure 4.6: Heat currents, Jh , Jh , Jh are compared against JhEx . For all (2) range of V0 shown in this figure, Jh is close to JhEx . Jh depends strongly on G2 while particle current J in the previous figure relies on G1 . making explicit use of G2 : hJh i = N −2 −2 X Im t2 G1 (l − 1, l + 1) N −2 l=2 +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 = N −2 X B i hB l−1 , hl (4.49) l=2 and † † † † hB l = −tcl cl+1 − tcl+1 cl + V cl cl cl+1 cl+1 . (4.50) (2) Again, we find from Fig.4.6 that the agreement between Jh and JhEx is (1) much better than the one between Jh and JhEx , especially for larger V0 . This indicates that the finite G2 are essential for getting proper values of heat currents. 77 4.5. Conclusion and discussion Keeping orders of V02 in the bath operators m̂α is consistent with truncation 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 different terms (see Appendix A for more details). An investigation of their effect 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 Redfield equation and two systematic approximations are suggested to truncate and solve the hierarchy. Using the first-level of the first method and both the first- and the second-levels of the second method, non-equilibrium stationary 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 find that the first method works for strong interaction V0 . However, because it requires an explicit diagonalization of the Hamiltonian, it is less computationally efficient. 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 efficient and therefore it can be applied to much larger systems, but its accuracy decreases with increasing interaction strength V0 . However, truncation at a higher level is shown to significantly 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 significantly more involved computations. In future work, one can investigate the generalization to higher level truncations. Even more interesting, however, is the application of this method to physical problems. 78 Chapter 5 Solving the Redfield 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 often a quantum master equation for a single photon (boson) mode is solved in the coherent-state representation, where the operator-form quantum master equation becomes a differential equation of c-numbers and a density operator becomes a “distribution” function. The problem of the Hilbert space’s infinite dimension is bypassed. For example, a distribution function over N complex numbers is sufficient 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 stochastic differential equations such as Langevin equations (see Refs. [70, 98] and also Appendix E for the relation between generalized Fokker-Planck equations 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 coherentstate representations). In the P -representation coherent eigenvalues ξ and ξ ∗ are treated as conjugate complex numbers, thus one solves only the equations for ξ. Sometimes it is also useful to work with generalizations of the P -representation [70]. In fact, the generalized P -representation, where coherent eigenvalues ξ and ξ ∗ are treated independently, has been used in simulations 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 efficient method to solve the Redfield equation: using the coherent-state representation, the Redfield equation is mapped onto a generalized Fokker-Planck equation. For noninteracting systems, analytical results are obtained for the non-equilibrium stationary states. For interacting systems, we find that the resulting equations can be solved efficiently 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 simplifies 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 quantum master equation, such as the Redfield 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 Redfield equation efficiently. The fact that there is no analytical expression for non-equilibrium stationary states may simply be due to the technical difficulty in solving the quantum master equation. A more efficient way of solving this equation might lead to different 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 finding such non-equilibrium stationary 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 Redfield equation. The focus of this chapter, however, is to study the efficiency of this approach when applied to interacting systems. The chapter is organized as follows. In section 5.2, we introduce the model system, define its Redfield equation and then convert the operatorform Redfield equation into a differential equation in the coherent-state representation. We find the mapped equation to be in the form of a generalized Fokker-Planck equation. We then proceed to find 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 Redfield equation. This provides further support for the work presented in the previous chapter. In section 5.5, we solve the generalized FokkerPlanck equation via numerical simulations for interacting bosonic systems. 80 5.2. The model and its effective 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 effective equation of motion For concreteness, the system of interest in this chapter is a 1-D BoseHubbard chain [96, 97] coupled to a bath at each of its ends. The left ∆T (right) bath is maintained at temperature TL = T + ∆T 2 (TR = T − 2 ). 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: H S = H 0 + VS = N −1 X l=1 N UX a†l al+1 + a†l+1 al + nl (nl − 1) , 2 l=1 X ωk,α b†k,α bk,α + 1 , HB = k,α HSB = λ X Vkα a†α bk,α + h.c. . (5.1a) (5.1b) (5.1c) k,α=L,R 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 Redfield equation reads, from Eq. (2.48): −λ2 ∂ρ(t) = −i[HS , ρ(t)] ∂t i o X nh ˆ α ρ(t) + h.c. , a†α , m̂α ρ(t) + aα , m̄ (5.2) α=L,R 81 5.2. The model and its effective equation of motion ˆ L (m̄ ˆ R ) is related to a†1 (a†N ), where m̂L (m̂R ) is related to a1 (aN ) while m̄ Z ∞ X α 2 m̂α = |Vk | dτ aα (−τ ) e−iωk,α τ h1 + n (ωk,α )i, (5.3a) 0 k ˆα = m̄ X Here hn (ωk,α )i = e ∞ dτ a†α (−τ ) eiωk,α τ hn (ωk,α )i. (5.3b) 0 k Z |Vkα |2 ω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 |Vkα |2 may depend on . However, for simplicity we set the density of state and also |Vkα |2 to constants and absorb them into the coupling constant λ. Next we map the Redfield equation into a differential equation in the coherent-state representation. A similar approach has been used in earlier works mentioned in Ref. [72]. However, they are based on different equations 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 Redfield equation for finite bath temperatures. From a technical point of view, in the coherent-state representation it is easier to deal with the localoperator Lindblad equation than the Redfield 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 Redfield equation. ˆ 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 X X (n) †,(n) (5.4) al (t) = U n al (t) , a†l (t) = U n al (t) , n n ˆ. 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†l (t)) is a linear combination of all am (a†m ), (0) al (t) = al (t) = πkl πkm −ik t 2 X sin sin e am , N +1 N +1 N +1 (5.5) km 82 5.2. The model and its effective equation of motion where k = 2 cos Nkπ +1 . This is because for finite-size systems with free ends, the Fourier transform is discrete and eikl is replaced by the sin Nπkl +1 wavefunction, see for example Eq. (A.4). For interacting systems, in principle it is possible to perform this perturbational calculation up to a large n, however here we will stop at the first (1) order, al (t), Z t (1) †,(0) (0) (0) al (t) = i dτ al (τ ) al (τ ) al (τ ) . (5.6) 0 After some tedious but straightforward algebra( see Ref. [101] and also Appendix A for details), we find, to the linear order in U , that: X X Dα;m am + U Dα;m1 m2 m3 a†m1 am2 am3 + O U 2 m̂α = (5.7a) m1 m2 m3 m ˆα = m̄ X X D̄α;m a†m + U m Dα;m1 m2 m3 a†m3 a†m2 am1 + O U 2 , (5.7b) m1 m2 m3 where πklα πkm 2 X sin sin [1 + n (k , Tα , µα )] , N +1 N +1 N +1 k πklα πkm 2 X sin sin n (k , Tα , µα ) , D̄α;m = π N +1 N +1 N +1 k 4 2 Dα;m1 m2 m3 = π N +1 X n ( (k2 ) + (k3 ) − (k1 ) , Tα , µα ) − n ( (k) , Tα , µα ) (k2 ) + (k3 ) − (k1 ) − (k) Dα;m = π · (5.8a) (5.8b) k,m,k1 ,k2 ,k2 πklα πkm πk1 m1 πk2 m2 sin sin sin N +1 N +1 N +1 N +1 πk3 m3 πk1 m πk2 m πk3 m sin sin sin . · sin N +1 N +1 N +1 N +1 · sin (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 operators a†l and al become differential operators in the coherent representation, for example: ∂ † ∗ P ξ~ . al ρ ↔ ξl − (5.9) ∂ξl 83 5.2. The model and its effective equation of motion Using the above approximation for the operators m̂α and the transformation from creation/annihilation operators to differential operators, Eq.(5.2) becomes, in the coherent P -representation, a generalized Fokker-Planck equation, ∂P (0) (1) = LH0 + U LV + λ2 LB + λ2 U LB P. ∂t (5.10) where ∂ ∂ ξl + i ξl+1 + c.c. , ∂ξl+1 ∂ξl l=1 N X ∂ ∂ 2 ξl2 ∗ 2 LV = +i ξ ξ + c.c. , −i 2 ∂ξl l l ∂ξl 2 l=1 X ∂ ∂2 = Dα,m − D̄α,m ξm + D̄α,m + c.c. , ∗ ∂ξα ∂ξα ∂ξm α,m X ∂2 (1) ∗ LB = Dα,m1 m2 m3 ξm ξ 1 m3 ∗ ∂ξ ∂ξ α m 2 α,m LH0 = (0) LB N −1 X i (5.11a) (5.11b) (5.11c) i + ∂2 ∗ ∂ξα ∂ξm 3 ∂2 ξm ξm + c.c. ∂ξα ∂ξm1 2 3 ∂3 ξm1 + c.c. . − ∗ ∂ξα ∂ξm2 ∂ξm3 ∗ ξm ξ − 1 m2 (5.11d) (5.11e) Here c.c. refers to the complex conjugate and it is performed formally, e.g. (1) −i ∂ξ∂∗ ξl∗ is the c.c. of i ∂ξ∂l+1 ξl . The last term, Eq. (5.11e), of LB is l+1 separated from the other terms intentionally for reasons that will become apparent soon. If only equilibrium states are of interest, then the first and second term can be discarded since equilibrium states commute with HS . However, for ∂ non-equilibrium stationary statesthoseterms are important. Setting ∂t P = ~ ∞ , which in the following is denoted 0, we get the stationary state P ξ, simply as P ξ~ . Next, we will calculate the stationary of Eq. (5.10). We first 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 (0) For U = 0 only the terms LH0 and LB contribute, and the stochastic equation is a Fokker-Plank equation [98], in terms of complex variables: # " X ∂2 ∂ ∂ Γij ξj + (5.12) P = − Dij ∗ + c.c. P. ∂t ∂ξi ∂ξi ∂ξj∗ i,j For our specific case of a non-interacting system we have: X λ2 πδlα δαm , Γlm = −i (δl,m−1 + δl,m+1 ) − (5.13) α,m and 1 X D̄α,l δmα + D̄α,m δlα . Dlm∗ = λ2 2 α (5.14) Note that we have verified that D † = D. This is in fact a standard Fokker-Planck equation (FPE) for the OrnsteinUhlenbeck process. Such a FPE can be solved analytically [70, 98]. Its stationary solution is found to be a Gaussian function: 1 1 † −1 P ξ~ = e− 2 ξ σ ξ , (5.15) Z where Z is a normalization constant and Γσ + σΓ† + D = 0, (5.16) 1 hν l |D| ν m i |µl ihµm | . κl + κ∗m (5.17) which has the solution, σ=− X l,m Here κ, |µi , hν| are respectively the eigenvalues, right eigenvectors and left eigenvectors of Γ and hµ| , |νi are the complex conjugate of the two corresponding vectors. This expression is correct for both equilibrium and non-equilibrium states. The only difference is in the matrix D, in whether TL is the same as, or is different from TR . The only needed numerical task is to find an eigenvector decomposition of Γ, an N -dimensional matrix. Once σ is known, all correlation functions canbe calculated, in particular the single-particle Green’s function is G l† , m = hξl∗ ξm i = 2σlm . 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 firstly confirm that it leads to the proper thermal equilibrium, and secondly to show the general procedure to find 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 algorithm 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 specific case, −λ2 π −i D11 D12 (5.18) Γ= ,D = , −i −λ2 π D12 D22 where n(1, TL , µL ) + n(−1, TL , µL ) , 2 n(1, TL , µL ) − n(−1, TL , µL ) + n(1, TR , µR ) − n(−1, TR , µR ) = λ2 π , 4 n(1, TR , µR ) + n(−1, TR , µR ) D22 = λ2 π . 2 (5.19) D11 = λ2 π D12 ¿From this we get the stationary distribution, 2 1 D11 − D22 D11 + D22 2D12 λ π −i σ= 2 + . 2D12 D11 + D22 i −λ2 π 4λ π 4 (λ4 π 2 + 1) (5.20) For equilibrium states TL = TR = T , µL = µR = µ and therefore D11 = D22 , so that ∗ +ξ∗ ∗ −ξ∗ ξ1 ξ1 ξ +ξ ξ −ξ 1 1 √ 2 1√ 2 − √ 2 1√ 2 − n(1,T,µ) ~ n(−1,T,µ) 2 2 2 2 , P ξ ∝e (5.21) which is exactly the P -representation form of the Boltzmann distribution 2 in terms of the eigenmodes a1√±a with eigenenergies ±1. Note that this is 2 independent of the coupling constant λ since the second term in σ vanishes and the λ2 in the first term cancel out. This confirms 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,1 − D2,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 previous chapter [101]. Numerical exact solution for non-interacting systems have been reported in Refs. for the Redfield 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,An and Gre,B . 1 1 re,An re,B Similar notations are used for particle currents J and J . We measure the difference between the two correlations, for example the analytical solution and the BBGKY solution, through, r 2 P re,An † (l , m) − Gre,B (l† , m) 1 l,m G1 r (5.22) dre,An = . B 2 P re,B † (l , m) l,m G1 In Fig.5.1, we compare this analytical solution and the BBGKY-like numerical results for the Redfield equation. Since dre ≈ 0, the two methods are indeed equivalent. The inset shows that currents calculated from both methods also agree. This validates the fact that the generalized Fokker-Planck equation approach produces accurate results at least for non-interacting systems. In the same figure, we also plot dre,nu and J re,nu , which are obtained B from a numerical solution that will be discussed in detail in section 5.4. There is a small but finite 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 FokkerPlanck 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 difference equations, which are generalizations of the Langevin equation [98, 102]. We will discuss this technique in the next section §5.5. Before doing that, 87 5.4. Interacting systems: BBGKY-like hierarchy derived from the GFPE 0.6 d d J 0.2 d re,An B re,Nu 0.3 0 0 0.4 ∆T/T B re,B 0.8 J re,An J Re,Nu J 0.1 0 0 0.2 0.4 ∆T/T 0.6 0.8 1 Figure 5.1: Analytical solution of the generalized Fokker-Planck equation is compared against the BBGKY numerical solution for non-interacting systems. The red (circles) curve in the main plot shows the comparison between the above analytical solution and the BBGKY-like numerical solution of the Redfield equation. Both Green’s functions (dre,An , main panel) and particle B currents (J re,An and J re,B , inset) exactly agree with each other. The blue (square) curve is a comparison between a numerical simulation, to be discussed in §5.4, and the BBGKY solution. Here we find 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 mirrors 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 R ∗ † ∗ ~ ~ tions G l ; n , defined as hξl ξn i = dξξl ξn P ξ . Let us then derive a set of equations for them from the generalized Fokker-Planck equation, Eq. (5.10). For example, consider a first-order derivative term in Eq. (5.11b). After performing a partial integration, we have, XZ ∗ 2 ~ ∗ ξn i ∂ ~ ∞ dξξ ξm ξm P ξ, l ∂ξm m Z X ~ ~ ∞ ξ ∗ ξ 2 ∂ (ξ ∗ ξn ) dξP ξ, = −i m m ∂ξm l m Z ~ ~ ∞ ξ∗ξ∗ ξ2 = −i dξP ξ, l n n = −i ξl∗ ξn∗ ξn2 . Therefore, setting ∂ ∗ ∂t hξl ξn i (5.23) to zero, we have ∗ ∗ 0 = −ihξl∗ ξn+1 i − ihξl∗ ξn−1 i + ihξl+1 ξn i + ihξl−1 ξn i X ∗ ξn i] +λ2 D̄α,m − Dα,m [δnα hξl∗ ξm i + δlα hξm α,m +λ2 +λ2 U X X D̄α,n δlα + D̄α,l δnα α +iU hξl∗ ξl∗ ξl ξn i ∗ ξ iδlm3 Dα,m1 m2 m3 δnα hξm 1 m2 − iU hξl∗ ξn∗ ξn ξn i ∗ + hξm ξ iδlm2 1 m3 α,mi ∗ ∗ +δlα hξm ξ iδnm3 + hξm ξ iδnm2 2 m1 3 m1 . (5.24a) (5.24b) (5.24c) (5.24d) (5.24e) This is exactly the same with Eq. (4.7). This confirms that the BBGKYlike 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 finite series of derivatives beyond the second order. This does not, however, rule out the possibility that a generalized Fokker-Planck equation with an infinite series of derivatives may lead to a positively valued distribution function. Our Redfield equation is, in principle, such a generalized FokkerPlanck equation with infinite derivatives, if all terms in the expansion of the m̂ operators are kept. A truncation at a certain finite 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 difference 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 specific kind of stochastic difference 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 finite systems. Another benefit of this approach is that one could also go beyond Markov approximation. On the other hand, one should also note the possibility that dropping such triple derivative terms may be appropriate in certain circumstances [104]. For example, in Eq. (5.24) — the first-order BBGKY equation derived above for two-variable (corresponding to single-particle) correlations, 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 example in the equation of four-variable correlations, which itself appears in the equation of two-variable correlation. Therefore such third order terms affect 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 effects. This is better than not taking the interaction in consideration at all. In this sense, neglecting the third order terms is equivalent to the first order of the cluster expansion, where the four-variable correlations are approximated by products of two-variable ones such that effects 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 effects of their inclusion will be the topic of further investigations. For our specific 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, X ∂ (2) ∂ ∂ (2) − P = Γij ξj − Γ ∗ ∗ ηj ∂t ∂ξi ∂ηi i j i,j ∂2 ∂2 (2) (2) Dij + Di∗ j ∗ ∂ξi ∂ξj ∂ηi ∂ηj 2 ∂ ∂2 (2) (2) + D ∗+ D ∗ P, ∂ξi ∂ηj ij ∂ηi ∂ξj i j + (5.25) (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 diffusion matrix D is always positive definite so that we can work with the P representation. For interacting systems, however, D is not positive definite therefore we have to use the more general representation. Here we use the the positive(2) (2) (2) (2) P representation. We denote the coefficients as Γij , Γi∗ j ∗ and Dij , Di∗ j , (2) (2) (2) Dij ∗ , Di∗ j ∗ . In principle, there could also be terms such as − ∂ξ∂ i Γij ∗ ηj , but for this specific problem such terms do not appear. The explicit definitions of those coefficients are, Γij = Γij + iU ξi ηi δij , Γi∗ j ∗ = (Γij )∗ − iU ξi ηi δij , (2) (2) (5.27) and (2) Dij = −i U 2 λ2 U ξi δij − 2 2 X (Dα,j,m1 ,m2 δα,i + Dα,i,m1 ,m2 δα,j ) ξm1 ξm2 , α,m1 ,m2 (5.28a) 91 5.5. Approximate numerical solution for interacting systems (2) Di∗ j ∗ = i λ2 U U 2 ηi δij − 2 2 X (Dα,j,m1 ,m2 δα,i + Dα,i,m1 ,m2 δα,j ) ηm1 ηm2 , α,m1 ,m2 (5.28b) (2) (2) Dij ∗ = Dj ∗ i = Dij ∗ + λ2 U 2 X [(Dα,m1 ,j,m2 + Dα,m1 ,m2 ,j ) δα,i α,m1 ,m2 + (Dα,m2 ,i,m1 + Dα,m2 ,m1 ,i ) δα,j ] ηm1 ξm2 . (5.28c) Eq. (5.26) can be solved through the following Langevin equation [98] in positive-P representation [70, 72], X (2) √ dξi = Γij ξj dt + 2dwi (t) , (5.29a) j dηi = X (2) Γi∗ j ∗ ηj dt + √ 2dw̄i (t) , (5.29b) j where dw (t) and dw̄ (t) are stochastic processes satisfying, (2) dwi (t) dwj t0 = Dij δtt0 dt, (2) dwi (t) dw̄j t0 = Dij ∗ δtt0 dt, (2) dw̄i (t) dwj t0 = Di∗ j δtt0 dt, (2) dw̄i (t) dw̄j t0 = Di∗ j ∗ δtt0 dt. We can define a 2N -dimensional symmetric matrix D, (2) (2) Dij Dij ∗ D = (2) N ×N (2) N ×N . Di∗ j Di∗ j ∗ N ×N (5.30a) (5.30b) (5.30c) (5.30d) (5.31) N ×N Then these random variables can be simplified further and organized as a linear transformations of 2N independent Wiener noises d~ ω [70], dw ~ = Bd~ ω, (5.32) D = BB T . (5.33) where This requires a diagonalization of the 2N × 2N 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.4 0.3 d J 0.6 re,Nu B 0.3 0 0.2 d 0 0.4 ∆T/T B 0.8 J re,Nu J 0.1 0 0 0.2 0.4 ∆T/T 0.6 0.8 1 Figure 5.2: The numerical simulation of the generalized Fokker-Planck equation approach is compared against the BBGKY-like numerical solution of the Redfield equation, for interacting systems. The calculated current J re,N u is u very close to J re,B while dre,N in the main panel is relatively large. Other B 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 accurate results for non-interacting systems, where analytical solutions are possible. Next, we compare this numerical simulation against the BBGKY-like method for interacting systems, in Fig.5.2. The currents Jre,N u and J re,B are in agreement with each other up to high temperature biases. However, u dre,N — the difference between the two Green’s functions – is relatively B 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 Itō form of the stochastic differential equation. The Euler method is not always stable, especially for large U . In fact, more stable and efficient 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 difficulty in implementing such methods for the Redfield equation: matrices B are found numerically so it is hard to derive the corresponding Stratonovich u could also be dure to the truncation of form. The relatively large dre,N B 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 Redfield equation, some authors prefer to use the localoperator Lindblad equation in the study of transport [58, 59, 99]. In this section, we demonstrate how to find NESSs from the local-operator Lindblad equation using the coherent-state representation. The Redfield equation finds that the operators m̂α for the sites coupled to baths typically involve more than just the operator at the coupling sites. The local-operator Lindblad equation, on the other hand, chooses m̂α ∝ h1 + n (ωk,α )iaα , which results in the much simpler form of the coupling to baths: ∂ρ(t) = −i[HS , ρ(t)] ∂t h i n h i −λ2 π h1 + n (1 , TL , µL )i a†1 , a1 ρ(t) + hn (1 , TL , µL )i a1 , a†1 ρ(t) + h.c. h i h i o +h1 + n (N , TR , µR )i a†N , aN ρ(t) + hn (N , TR , µR )i aN , a†N ρ(t) + h.c. . (5.34) In the case of this specific two-site systems, we have chosen that N = 2 and 1,2 = 0. Note that we get Eq. (5.34) from the Redfield equation Eq. (5.2) by letting Dα;m = πδαm [1 + n (α , Tα , µα )] , D̄α;m = πδαm n (α , Tα , µα ) , Dα;m1 m2 m3 = 0. (5.35) These relations greatly simplify the resulting generalized Fokker-Planck equation. Following the same procedure as for the Redfield equation, we 94 5.6. Solving the LOLE via the GFPE method derive a similar Langevin equation for Eq. (5.34), " # X (n) 2 2 dξn = −i (ξn−1 + ξn+1 ) − iU ξn ηn − λ π (δαn ξα ) dt + B1ν dwν(n) , α " dηn = i (ηn−1 + ηn+1 ) + iU ξn ηn2 − λ2 π X (5.36a) # (n) (δαn ηα ) dt + B2ν dwν(n) , α (5.36b) (n) where the matrix B (n) = Bµν 2×2 B (n) B (n) T = satisfies, 2 2λ2 π P −iU ξn 2 2λ π α n (α , Tα , µα ) δnα P α n (α , Tα , µα ) δnα iU ηn2 , (5.37) n o (n) and dwν 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 interacting bosonic systems. The 2N -dimensional matrices B are decoupled into N 2-dimensional matrices B (n) , which can be diagonalized analytically. Therefore, it is technically easier to solve the local-operator Lindblad equation than the Redfield 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, analytical expressions can be derived, and one finds 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, n (0, TL , µL ) + n (0, TR , µR ) 1 0 lole σ = 0 1 4 2 2 λ π [n (0, TL , µL ) − n (0, TR , µR )] λ π i . + (5.38) −i −λ2 π 4 (λ4 π 2 + 1) In Fig.5.3, we have compared these analytical and numerical results against the BBGKY results. Results for both currents and Green’s functions are found to be consistent for the two approaches. This confirms that our generalized Fokker-Planck equation method works for both the Redfield equation and the local-operator Lindblad equation. 95 5.7. Conclusions 0.4 J d lole,An B lole,Nu B 0.3 lole,Nu B 0.2 0 0.4 ∆T/T 0.8 d d 0.2 0 d 0.4 J d 0.4 0.2 0.2 lole,B 0 J lole,An J lole,Nu J 0.1 0 0.4 ∆T/T 0.8 lole,B J lole,Nu J 0.1 0 0 0.2 0.4 ∆T/T 0.6 0.8 1 0 0 0.2 0.4 ∆T/T 0.6 0.8 1 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 Redfield equation into a generalized 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 infinite (bosons) dimension. For interacting systems, firstly 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 systems, 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 equation. Here the method is even more efficient, since at every time step one needs to solve N 2-dimensional eigenvalue problems instead of a 2N dimensional one, like for the Redfield equation. Results are found to be consistent with those from the BBGKY-like method. Further development of the techniques for solving the generalized FokkerPlanck 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 phenomenological 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 effort [19, 35, 36, 42– 44, 56, 59, 60, 82, 83, 106, 115]. Most of this work [35, 42–44, 60, 82, 115] studied infinite and/or periodic chains, and used the standard Kubo formula [6, 40] where finite/zero Drude weight signals anomalous/normal transport. 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 thermal 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 terminology, 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 effects 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 conflicting 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)14 Cu24 O41 , Sr2 CuO3 and CuGeO3 , [79, 120, 121] seems to validate this conjecture, but Ref. [81] finds normal transport in Sr2 CuO3 at high temperatures. Even if there was no conflict between results based on the local-operator Lindblad equation [56, 59], as discussed in Chapter 2 this approach is less reliable than the Redfield 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 finite spin chains coupled to thermal baths, using the Redfield 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 transport 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 Redfield equation. We study such systems here. We find no direct relation between integrability and anomalous transport. Instead, based on our systematic investigations, 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 diffusion is universally present in interacting 1D systems [42]. It is important to note that, chronologically, this was the first project undertaken once the formalism of the Redfield equation was adopted. This 98 6.2. The model and its Redfield equation explains why the results shown below are derived using (inefficient) direct methods to solve the Redfield 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 inefficiency that drove the subsequent research on more efficient methods, presented in the previous chapters. In future, we plan to use these efficient methods to extend this study to larger systems. The chapter is organized as follows. Section 6.2 defines the system and its Redfield equation. In Sections 6.3 we discuss the numerical methods and in Section 6.4 we define 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 Redfield equation We consider an N -site chain of spins- 21 described by the Hamiltonian: HS = N −1 X N X y y x x z z ~ ~si Jx si si+1 + Jy si si+1 + Jz si si+1 − B · (6.1) i=1 i=1 while the heat baths are collections of bosonic modes: X HB = ωk,α b†k,α bk,α (6.2) k,α where α = R/L indexes the right/left-side baths. The system-baths coupling is taken as: X (α) Vk syα ⊗ b†k,α + bk,α (6.3) V =λ k,α where syL = sy1 and syR = syN , i.e. the left (right) thermal bath is only coupled to the first (last) spin and can induce its spin-flipping. This is because we ~ · ~ey = 0 while |B| ~ is finite, meaning that spins primarily lie in the choose B y x0z plane so that s acts as a spin-flip operator. From Eq. (2.29), the resulting Redfield equation for ρ (t) is: X ∂ρ(t) = −i[HS , ρ(t)] − λ2 ([syα , m̂α ρ(t)] + h.c.) ∂t (6.4) α=L,R where m̂α = syα · Σα . Here, (·) refers to the element-wise product of two matrices, hn|a · b|mi = hn|a|mihn|b|mi. The bath matrices ΣL,R are defined 99 6.3. Numerical methods in terms of the eigenstates of the system’s Hamiltonian HS |ni = En |ni as: Σα = π X h (α) |mihn| Θ (Ωmn ) nα (Ωmn ) Dα (Ωmn ) |Vkmn |2 m,n (α) +Θ (Ωnm ) (1 + nα (Ωnm )) Dα (Ωnm ) |Vknm |2 i where Ωmn = Em − En = −Ωnm and kmn is defined by ωkmn ,α = Ωmn , i.e. a bath mode resonant with this transition. Furthermore, Θ(x) is the Heaviside −1 function, nα (Ω) = eβα Ω − 1 is the Bose-Einstein equilibrium distribution 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 ) |Vkmn |2 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 Redfield 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 → ∞. As already discussed, another approach [58] is to solve directly for the stationary state from: Lρ(∞) = 0, (6.5) i.e. to find 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 efficiency than this eigenvalue problem and also better time efficiency 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 normalization condition tr (ρ) = 1 such that L̄ρ(∞) = ν, (6.6) where ν = [1, 0, · · · ]T and L̄ is found from L by replacing the first row by tr(ρ). Then we use for example the generalized minimal residual method 100 6.4. Definitions of the thermal current and local temperatures (GMRES) [122], which requires only the matrix-vector multiplication rule. This new way of solving the Redfield equation has memory cost of ∼ 22N and time efficiency of a linear system with dimension 22N . It is still a direct method so its efficiency 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 efficient methods can, for example, destroy the integrability of the system. Also, because the commutation relations between spin operators are quite different from those of fermions and bosons, the more efficient methods are not readily generalizable. Developing and applying them to this problem will be the topic of future work. 6.4 Definitions of the thermal current and local temperatures P −1 PN We rewrite HS = N i=1 hi,i+1 + i=1 hi , where hi,i+1 is the exchange between nearest-neighbor spins and hi is the on-site coupling to the magnetic field. We can then define a local site Hamiltonian (S) hi 1 1 = hi−1,i + hi + hi,i+1 2 2 (6.7) with h0,1 = hN,N +1 = 0, and a local bond Hamiltonian (B) hi 1 1 = hi + hi,i+1 + hi+1 2 2 (6.8) such that HS = N X (S) hi i=1 = N −1 X (B) hi . (6.9) i=1 The local bond Hamiltonians can be used to derive the heat current operator from the continuity equation (B) ĵi→i+1 − ĵi−1→i = ∇ĵ = − This results in ∂hi ∂t h i (B) = −i HS , hi . h i (B) (B) ĵi→i+1 = i hi , hi+1 (6.10) (6.11) 101 6.4. Definitions of the thermal current and local temperatures fori = 1, · · · , N − 2. As expected, in the steady state we find hĵi→i+1 i = tr ĵi→i+1 ρ (∞) = J to be independent of i. Similarly, we define the spin polarization hszi i and the spin current for XXZ model, Js = Jxy syi sxi+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 profile 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 finite intrinsic resistivity. In principle, the scaling of the current with the sample size, for a fixed effective 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 effective bias, i.e. the difference 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 profile along the chain, in order to determine the type of transport. This, of course, is the “local temperature”, which is a difficult quantity to define properly. One consistency condition for any definition is that if TL = TR = T , i.e. the system is in thermal equilibrium at T , then all local temperatures should equal T . We define local site temperatures Ti which fulfill this condition in the following way. Since we know all eigenstates of HS , it is straightforward to calculate its P equilibrium density matrix at a given T , ρeq (T ) = Z1 n e−βEn |nihn|, where P (S) (S) Z = n e−βEn . Let then hhi ieq (T ) = tr[ρeq (T ) hi ]. We define Ti to be the solution of the equation: (S) (S) hhi ieq (Ti ) = tr[ρ (∞) hi ]. (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 (B) Of course, we can also use other “local” operators such as hi to calculate a local bond temperature Ti+ 1 . We find that when these definitions are 2 meaningful, the results are in very good agreement no matter what “local” operator is used. This type of definition of Ti is meaningful only if a large magnetic field B (S) is applied. For small B, the expectation values hhi 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 (S) hhi ieq (T ) which varies fast enough with T for values of interest, so that a meaningful Ti can be extracted. Since we could not find a meaningful ~ → 0, we cannot investigate such cases. Note, definition for Ti when B however, that most integrable models remain integrable under addition of ~ = Bêz . an external field B 6.5 Results In all our calculations, we take Bz = 1 and the exchange J ∼ 0.1. Temperatures 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 N B, 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 profiles Ti , Ti+ 1 and (b) local spin polarization profiles Siz . We apply a large bias δ = 2 (TL − TR )/T = 0.4 for clarity, but we find similar results for smaller δ. For these values, it seems that the “contact regions” include only the end spins. The profile of the central part of the chain is consistent with anomalous transport (flat profile) for the XX chain (Jx = Jy , Jz = 0) and shows normal transport (roughly linear profile) in all the other cases. XY chains with Jx 6= Jy behave similarly with the XX chain. We find similar results for ferromagnetic couplings. We also find that the ratio between the effective temperature difference, T2 − TN −1 , and the applied temperature difference TL − TR , is not a constant for different system sizes N . Therefore, in our numerical calculation, it is not possible to keep the effective temperature difference as a constant by applying the same temperature difference 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 Ti S XX XXZ(0.05) XXX XXZ(0.2) XYZ 2.2 z i 0.14 2 0.12 1.8 0.1 1 2 3 4 5 6 i 7 8 9 10 1 2 3 4 5 6 7 8 9 10 i Figure 6.1: Plot of (a) local temperature profile and (b) local spin polarization profile 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 flat Ti and Siz profiles. All other models have almost linearly decreasing profiles of both temperature and S z . 104 6.5. Results 0.002 XX XYZ -0.89 J=0.001*(N-1.9) XXX -0.37 J=0.0017*(N-0.23) J 0.0015 0.001 0.0005 0 4 5 6 7 N 8 9 10 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 effective applied temperature difference, the fitting 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 S z profiles for all available sizes. We normalize all the profiles for various models with different sizes by harvesting only the data in the central regions and converting them into dimensionless values in [0, 1], for example, Ti −TN −1 from (i, Ti ) to ( Ni−2 −3 , T2 −TN −1 ), i = 2, 3, . . . , N − 1. In Fig.6.3, we plot such normalized temperature and S z profiles. We see that curves for all values of N collapse onto one single straight line and furthermore, there is no difference between XY Z and XXX chains. Currently we can only study small systems, but from these normalized profiles, no qualitative difference has been found for different 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 finite Jz leads to nearest-neighbor interactions between 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 1 1 z z -S 2) N-1 z z (S i-S 2)/(S (Ti-TN-1)/(T2-TN-1) N=10,XXX 9 8 7 6 5 4 10, XYZ 0.5 0.5 0 0 0 0.2 0.4 0.6 (i-1)/(N-2) 0.8 1 0 0.2 0.4 0.6 (i-1)/(N-2) 0.8 1 Figure 6.3: Normalized temperature and S z profiles are plotted. Curves for different values of N collapse onto a single straight line and we see no difference between XY Z and XXX chains. normal transport occurs, we plot temperature profiles 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 difference, namely that the contact regions seem larger than just the end spins, as can be seen from the normalized profiles. 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 confidence. As we pointed out before, our definition of local temperature and also the idea of local a spin polarization Siz are only meaningful for sufficiently large Bz . We find 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 field (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 λ ∈ [0.03, 0.2], T ∈ [0.3, 30.0] and δ ≥ 0.01, the spin chain has normal conductivity when Jz ∈ [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 find no difference between thermal transport 106 6.5. Results 2.2 2.2 (a) 0.5 2.1 0 Ti (b) 1 2.1 0 0.4 0.8 2 2 Jz=0.00 Jz=0.01 Jz=0.02 Jz=0.03 Jz=0.04 1.9 1.8 Bz=0.1 Bz=0.3 2 1.9 4 6 i 8 10 1.8 2 4 6 8 10 i Figure 6.4: (a) Temperature profiles for small Jz values are plotted. Here Jx = Jy = 0.1, Bz = 1.0. The effective temperature drop becomes smaller for smaller Jz , but the slope is still finite if Jz ≥ 0.02. The inset plots normalized profiles, and shows no obvious differences for models with Jz ≥ 0.02. The case of Jz = 0.01 shows a finite slope but the normalized profile indicates a slight difference: the contact regions seem to be more extended than in the other cases. (b) Temperature profiles of XXX chains with J = 0.1 and small values of Bz . A linear temperature drop is observed for Bz = 0.3, but a flat profile 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 field Bx . It maps to non-interacting spinless fermions [123] and if we add a Bz field, it becomes interacting. Another closely related model is a system of spinless fermions on a tight-binding chain, with the nearest-neighbor interaction: HS = N X l=1 c†l cl − t N −1 X N −1 X c†l cl+1 + c†l+1 cl + V0 c†l+1 cl+1 c†l cl , l=1 (6.14) l=1 X HB = ωk,α b†k,α bk,α , (6.15) c†α bk,α + cα b†k,α . (6.16) k,α=L,R VSB = λ X k,α The XXZ chains map exactly into this HS model, after using the JordanWigner transformation. However, the coupling to the baths in this fermionic system is different 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 c† 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 influenced by the specific model for the system-bath coupling. In Fig.6.5, we plot local temperature profiles 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 flat temperature profiles. In summary, the first conclusion we draw from these results is that integrability is not sufficient to guarantee anomalous transport: several integrable models show normal heat transport, in agreement with other studies [42, 58, 59]. The second conclusion is that only models that map onto Hamiltonians of non-interacting fermions exhibit anomalous heat transport. This is a reasonable sufficient 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 2.2 (a) (b) V0=0.0 V0=0.2 Isingx Isingxz 2.2 Ti 2 2 1.8 1.8 1.6 1 2 3 4 5 6 i 7 8 9 10 1 2 3 4 5 6 7 8 9 10 i Figure 6.5: Temperature profiles 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 flat temperature profile 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 artifacts due to mainly the boundary effects 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. Diffusion is impossible 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 field 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. Specifically, 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 first manifold of low-energy eigenstates have one spin flipped (single magnon states), followed by states with two spins flipped (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 transport. 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 find 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 differences 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 be110 6.6. Conclusions comes finite. 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 experimentally in compounds such as Sr2 CuO3 [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 Redfield equation, we propose a new conjecture for what determines the appearance of anomalous heat transport at all temperatures in spin chains. Unlike previous 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 finite 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 efficient 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, similar questions had been investigated using either the standard Kubo formula, the non-equilibrium Green’s functions method or the local-operator Lindblad equation. However, the idea of using the Redfield equation for these questions had just been proposed and it had been applied only to noninteracting systems or very small interacting but non-integrable systems. We first developed and applied the direct methods to the Redfield 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 confirm or reject this observation. However, at the time there were no available methods capable of solving the Redfield equation for large interacting systems. Thus, we turned to investigate such methods. We first 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 coupled to baths and instead we had to derive a similar Kubo formula from the Redfield 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 suffer from singularities, like the the standard Kubo formula, however it is not very efficient. A more efficient OKF was also identified, however only for certain types of quantities. Neither approach is efficient 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 approximations are needed to truncate this hierarchy, if there are interactions. We identified two possible approaches. The first 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 efficient 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 first 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 generates correlated two-particle Green’s functions, beyond the Hartree-Fock approximation. In Chapter 5 we proposed another method: solving the Redfield equation in the coherent-state representation. The resulting differential 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 Redfield equation and study also non-equilibrium stationary states. We have confirmed that for non-interacting systems exact results can be found analytically 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 find efficient and stable ways to solve these differential 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 definitely. The second-order form of the BBGKY-like method seems to have a good potential there. However, the study of spin systems requires different kinds of BBGKY-like equations, and also to understand whether the cluster approximation interferes with the integrability of a model. Therefore, although the methods we proposed are in principle applicable to this problem, specific 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 σiz σcz (b) Lc L BL L BR λL λR Figure 7.1: (a) A sketch of exactly solvable decoherence model: a single spin (~σc ) is coupled to a bath of N spins {~σi }. (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 systems of a few coupled quantum dots or molecular devices coupled to reservoirs. 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 decoherence control and classical simulations of quantum computations. In all models studied in this thesis, the baths are assumed to be collections of eigenmodes |ω ν (k)i, connected to the central systems via λVkν (see for example Eq. (2.44)). We then assume Vkν to have certain properties. A more realistic and specific setup is to connect the baths locally to the central system, as sketched in Fig7.1(b). In that case, Vkν is implicitly defined. 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 off-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 off-diagonal elements disappear and if so, which “diagonal” and “off-diagonal” parts are relevant, and how that depends on the specific 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 different approaches: starting from the pure dynamics of the whole composite system, which usually however is only applicable to non-interacting systems (both the central system and the baths); through the Redfield equation, which assumes the Markovian approximation, weak coupling and idealized infinite-size baths; and finally by the non-equilibrium Green’s function method, which is typically applied for configurations similar to the variation in Fig.7.1(b). In fact, there is another approach — the exact master equation or the influence functional method, which relaxes the assumption of the Markovian approximation and the weak coupling limit, but still takes the baths to be infinitely large. However, 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 implement the first 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 temperatures of the baths to be equal then the non-equilibrium Green’s function arrives at the corresponding equilibrium state of the whole composite system, 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 qualitatively different results. If we want to compare them, a proper model, on which the non-equilibrium Green’s function method and the Redfield equation ideally should produce the same result, is required. Comparison between the first two might shed some light on differences between the Markovian and the non-Markovian approximation, the validity of the infinite-size bath assumption, etc. Comparison between the latter two might help us understand them better or even improve both. Lastly, our efficient methods for the Redfield 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 Redfield equation. It should be very straightforward to implement the methods for the localoperator 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 Scientific, 1999. [4] S. Datta. Quantum Transport: Atom to Transistor. Cambridge University Press, 2005. [5] D.K. Ferry and S. M. Goodnick. Transport in Nanostructures. Cambridge 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 lattices. Philosophical Magazine, 21(172):863–867, 1970. [9] R. Landauer. Spatial variation of currents and fields due to localized 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 manychannel 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 calculation of tunneling current. Journal of Physics C-Solid State Physics, 4(8):916–929, 1971. [15] R. Combesco. Direct calculation of tunneling current 3: Effect 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. Redfield. 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 Interactions in Low Dimensions, chapter 12, pages 383–417. Springer, 2005. [20] A. Yacoby. Transport in quatnum wire. In D. Baeriswyl and L. Degiorgi, 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 organic molecules. Advanced Materials, 15(22):1881–1890, 2003. [23] N. J. Tao. Electron transport in molecular junctions. Nature Nanotechnology, 1(3):173–181, 2006. [24] W. Goepel and C. Ziegler. Nanostructures Based on Molecular Materials. 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 INSTITUTES SERIES, SERIES E, APPLIED SCIENCES, pages 105– 214, 1997. [26] D. Natelson. Handbook of Organic Electronics and Photonics. American Scientific 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. Shortsleeves, 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 momentum 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 field. 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 integrable 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 fluids. Physical Review, 112(6):1829–1842, 1958. [40] J. M. Luttinger. Theory of thermal transport coefficients. Physical Review, 135(6):A1505–A1514, 1964. [41] D. S. Fisher and P. A. Lee. Relation between conductivity and transmission matrix. Physical Review B, 23(12):6851–6854, 1981. [42] J. Sirker, R. G. Pereira, and I. Affleck. Diffusion and ballistic transport 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 field. Physical Review B, 71(18):184415, 2005. [44] M. Rigol and B. S. Shastry. Drude weight in systems with open boundary conditions. Physical Review B, 77(16):161101, 2008. [45] C. Toher, A. Filippetti, S. Sanvito, and K. Burke. Self-interaction errors 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 environment. Physical Review E, 77(1):011112, 2008. [51] R. Kubo, M. Toda, and N. Hashitsume. Statistical physics II: nonequilibrium statistical mechanics. Springer, 2nd edition, 1991. [52] A. G. Redfield. On the theory of relaxation processes. IBM Journal of Research and Development, 1(1):19–31, 1957. [53] A. G. Redfield. The theory of relaxation processes. Adv. Magn. Reson., 1:1–32, 1965. [54] W. T. Pollard, A. K. Felts, and R. A. Friesner. The redfield 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 multilevel redfield 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 confirmed 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 nonequilibrium 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 onedimensional quantum system. Europhysics Letters, 61(1):34–40, 2003. [61] W. Huisinga, L. Pesce, R. Kosloff, 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 redfield 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. Efficiency of different numerical methods for solving redfield 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 solution 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 mechanics. Kluwer Acad. Publ., 1995. [69] L. P. Kadanoff 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. Physical 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: Firstprinciples 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 lowdimensional lattices. Physics Reports-Review Section of Physics Letters, 377(1):1–80, 2003. [77] K. Damle and S. Sachdev. Spin dynamics and transport in gapped onedimensional 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. Diffusive 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 antiferromagnet, sr2cuo3, up to 800 k. Physical Review Letters, 87(24):247202, 2001. [82] F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brenig. Zerofrequency 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 Transition. 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. communications in mathematical physics, 48(2):119–130, 1976. [91] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan. Completely positive 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 kubogreenwood conductivity of disordered-systems. Journal of Physics CSolid 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 effects 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 superfluid 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 Letters, 56(3):372–378, 2001. [103] Benoit Palmieri, Yuki Nagata, and Shaul Mukamel. Phase-space algorithm for simulating quantum nonlinear response functions of bosons using stochastic classical trajectories. Physical Review E, 82:046706, 2010. [104] S.E. Hoffmann, J.F. Corney, and P. D. Drummond. Hybrid phasespace simulation method for interacting bose fields. 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 effectivehamiltonians 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 field. 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 fluctuations. 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 levelspacing 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 finite-size effects. 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 field. 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 spin1/2 xxz chain at arbitrary temperature. Journal of Physics AMathematical 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 specific 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 Scientific 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 Letters, 54(17):1879–1882, 1985. [124] M. Schlosshauer. Decoherence, the measurement problem, and interpretations 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 two In this appendix we will derive expressions of operators m̂ and m̄ 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 definition in Eq(4.3), in the representation of eigenvalues {Em } and eigenstates {|mi} of HS , the operator m̂ can be written as = (cα )mn X k Z = (cα )mn π Z |Vkα |2 hm |m̂α | ni = (m̂α )mn ∞ dτ ei(En −Em −ωk,α )τ h1 − nα (ωk,α )i 0 dωD α (ω) |V α (ω) |2 h1 − nα (ω)iδ (Ωmn − ω) = (cα )mn πDα (En − Em ) |V α (En − Em ) |2 h1 − nα (En − Em )i, (A.1) R∞ where we have used 0 dτ eiωτ = πδ (ω)+iP ω1 and neglected the principal value part. We have also assumed that it is possible to perform a change of variable on Vkα such that it becomes V α (knm ), where kmn is defined by ωkmn ,α = Ωmn = Em − En = −Ωnm , i.e. a bath mode resonant with this transition. This limits the possible forms of Vkα and ωk,α . For example, for each given energy Ωmn , there should be a unique value of Vkαnm . In all the work in this thesis, we take Vkα as a constant so this condition is satisfied. 128 A.2. Perturbational expansion starting from non-interacting systems Dα (ω) is the bath’s density of states. We arrive at X m̂α = π |mihn|hm|cα |ni (1 − nα (Ωnm )) Dα (Ωnm ) |Vkαnm |2 , m,n X ˆα = π m̄ |mihn|hm|c†α |ninα (Ωmn ) Dα (Ωmn ) |Vkαmn |2 . (A.2a) (A.2b) m,n We furthermore set Vkαnm Dα (ω) as a constant and absorb it into λ2 . Similar 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 finding such operators m̂ perturbationally. A.2 Perturbational expansion starting from non-interacting systems To present the general idea concretely, let us work on a specific system, the Hamiltonian HS = H0 + VS defined in Eq. (4.1), which we rewrite here for convenience, HS = −t N −1 X N −1 X c†l cl+1 + c†l+1 cl + V0 c†l+1 cl+1 c†l cl = H0 + VS . l=1 (A.3) l=1 Next assuming V0 is small and treating it perturbationally, we want to express the operator m̂α in terms of polynomials of {cm } 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 N 1 X klπ sin ck = √ cl , N +1 N l=1 (A.4) diagonalizes H0 , H0 = N X k c†k ck , (A.5) k=1 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 , (0) cl (t) ≡ cl (t) = 2 X πkl πkm −ik t sin sin e cm . N +1 N +1 N +1 (A.7) km 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. (0) Let us call cl (t) when V0 = 0 as cl (t). Treating this as the zeroth order form for the full dynamical cl (t), and expanding X (n) cl (t) = V0n cl (t) , (A.8) n (n) we derive a perturbation equation of cl (t): (n) (n) (n) (n−1) ċl = it cl−1 + cl+1 − idl , (A.9) where we define, for non-negative integers n, n1 , n2 , n3 , o X n (n ) †,(n ) (n ) (n) (n ) †,(n ) (n ) dl = cl 1 cl−1 2 cl−13 + cl 1 cl+1 2 cl+13 . (A.10) n1 ,n2 ,n3 P i ni =n The solution of the above equation can be written as Z t 2 X πkl πkm −ik (t−τ ) (n−1) (n) dτ sin cl (t) = −i sin e dm (τ ) . N +1 N +1 N +1 0 km (A.11) Here the initial condition c(n) (0) = 0(∀n ≥ 1) is used. As the expression of (0) (1) cl is given in Eq. (A.7), explicitly the expression of cl is (1) cl (t) = 2 N +1 4 X i((k2 )−(k1 )−(k3 ))t e − ei(k)t (k1 ) + (k3 ) − (k2 ) − (k) km ki mi πkl πk1 l πkm πk1 m1 πk2 m2 πk3 m3 · sin sin sin sin sin sin N + 1 N + 1 N + 1 N + 1 N + 1 N + 1 πk3 (l − 1) πk2 (l + 1) πk3 (l + 1) πk2 (l − 1) sin + sin sin · sin N +1 N +1 N +1 N +1 ·cm1 c†m2 cm3 . (A.12) 130 A.2. Perturbational expansion starting from non-interacting systems (2) Similarly we can derive cl , which will be shown explicitly later when we discuss terms proportional to V02 in bath operators m̂. Plugging this general solution into Eq(4.3), and after straightforward but tedious algebra, ˆ α in Eq(4.34) with expansion we arrive at the decomposition of m̂α and m̄ coefficients defined as follows, 2 X πklα πkm sin sin [1 − n (k , Tα )] , N +1 N +1 N +1 k 2 X πklα πkm D̄α;m = π sin sin n (k , Tα ) , N +1 N +1 N +1 Dα;m = π (A.13a) (A.13b) k for the zeroth order expansion, Dα;m1 m2 m3 = π 2 N +1 4 X k,m,k1 ,k2 ,k2 n ( (k) , Tα ) − n ( (k1 ) + (k3 ) − (k2 ) , Tα ) (k1 ) + (k3 ) − (k2 ) − (k) k1 πm1 k2 πm2 k3 πm3 kπm k1 πm kπlα sin sin sin sin sin · sin N +1 N +1 N +1 N +1 N +1 N +1 k2 π (m + 1) k3 π (m + 1) k2 π (m − 1) k3 π (m − 1) · sin sin + sin sin , N +1 N +1 N +1 N +1 (A.14) · 131 A.2. Perturbational expansion starting from non-interacting systems for the first order expansion. For the second order expansion we make use (2) of cl and we have Dα;m1 m2 m3 m4 m5 = π 7 X k,m,k0 ,m0 ,k1 ,k2 ,k3 ,k4 ,k5 πki mi πklα πk1 lα 5 Πi=0 sin sin sin N +1 N +1 N +1 n ( (k)) − n ( (k5 ) − (k0 ) − (k4 )) 1 (k1 ) + (k3 ) − (k2 ) − (k0 ) (k0 ) + (k4 ) − (k5 ) + (k) n ( (k)) − n ( (k1 ) + (k3 ) + (k5 ) − (k2 ) − (k4 )) − (k2 ) + (k4 ) − (k1 ) − (k3 ) − (k5 ) + (k) πk0 lα X πk2 (lα + δ) πk3 (lα + δ) · sin sin sin N +1 N +1 N +1 δ=±1 X πk4 (lα + δ) πk5 (lα + δ) · sin sin N +1 N +1 δ=±1 1 n ( (k)) − n ( (k0 ) + (k1 ) + (k5 )) + (k2 ) + (k4 ) − (k3 ) − (k0 ) (k) − (k0 ) − (k1 ) − (k5 ) n ( (k)) − n ( (k1 ) + (k3 ) + (k5 ) − (k2 ) − (k4 )) − (k2 ) + (k4 ) − (k1 ) − (k3 ) − (k5 ) + (k) X πk0 (lα + δ) πk4 (lα + δ) πk5 (lα + δ) sin sin · N +1 N +1 N +1 δ=±1 πk2 lα πk3 lα πk2 (lα + 2δ) πk3 (lα + 2δ) · sin sin + sin sin N +1 N +1 N +1 N +1 n ( (k)) − n ( (k1 ) − (k0 ) − (k2 )) 1 + (k3 ) + (k5 ) − (k4 ) − (k0 ) (k) + (k0 ) − (k1 ) + (k2 ) n ( (k)) − n ( (k1 ) + (k3 ) + (k5 ) − (k2 ) − (k4 )) − (k2 ) + (k4 ) − (k1 ) − (k3 ) − (k5 ) + (k) X πk0 (lα + δ) πk2 (lα + δ) πk3 (lα + δ) · sin sin sin N +1 N +1 N +1 δ=±1 πk4 lα πk5 lα πk4 (lα + 2δ) πk5 (lα + 2δ) · sin sin + sin sin . N +1 N +1 N +1 N +1 (A.15) πkm · sin N +1 2 N +1 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 fit into the page width. Therefore up the the order of V02 , the bath operators can be written down as X X m̂α = Dα;m cm + V0 Dα;m1 m2 m3 cm1 c†m2 cm3 +V02 X m1 m2 m3 m Dα;m1 m2 m3 m4 m5 cm1 c†m2 cm3 c†m4 cm5 + O V03 m1 m2 m3 m4 m5 ˆα = m̄ −V02 X X m D̄α;m c†m − V0 X (A.16a) Dα;m1 m2 m3 c†m3 cm2 c†m1 m1 m2 m3 Dα;m1 m2 m3 m4 m5 c†m5 cm4 c†m3 cm2 c†m1 + O V03 . (A.16b) m1 m2 m3 m4 m5 The same procedure is applicable for bosonic systems and similar expressions can be derived. This has been used in Chapter 5 to calculate the ˆ α , and find non-equilibrium stationary states for bosonic operators m̂α and m̄ 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 expansion 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, 0 0 0 0 G2 m† , n† , m , n = −G1 m† , m G1 n† , n 0 0 0 0 +G1 m† , n G1 n† , m + G2 m† , n† , m , n . (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 m†1 , m†2 , m†3 , m4 , m5 , m6 = G m†1 , m6 G m†2 , m5 G m†3 , m4 + G m†1 , m5 G m†2 , m4 G m†3 , m6 + G m†1 , m4 G m†2 , m6 G m†3 , m5 − G m†1 , m6 G m†2 , m4 G m†3 , m5 − G m†1 , m5 G m†2 , m6 G m†3 , m4 − G m†1 , m4 G m†2 , m5 G m†3 , m6 + G m†1 , m6 G m†2 , m†3 , m4 , m5 − G m†1 , m5 G m†2 , m†3 , m4 , m6 + G m†1 , m4 G m†2 , m†3 , m5 , m6 + G m†2 , m5 G m†1 , m†3 , m4 , m6 − G m†2 , m4 G m†1 , m†3 , m5 , m6 − G m†2 , m6 G m†1 , m†3 , m4 , m5 + G m†3 , m4 G m†1 , m†2 , m5 , m6 − G m†3 , m5 G m†1 , m†2 , m4 , m6 + G m†3 , m6 G m†1 , m†2 , m4 , m5 + G m†1 , m†2 , m†3 , m4 , m5 , m6 . (B.2a) (B.2b) (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 first equation of the equation hierarchy becomes 0 = itG (m − 1)† , n + itG (m + 1)† , n −itG m† , n + 1 − itG m† , n − 1 Xh +λ2 δnα Dα;l + D̄∗α;l G m† , l l,α +λ2 V0 i +δmα D̄α;l + D∗α;l G l† , n X (Dα;nm2 m1 − Dα;m1 m2 n ) G m†1 , m2 δmα α,m1 ,m2 +λ2 V0 (Dα;mm1 m2 − Dα;m2 m1 m ) G m†1 , m2 δnα X α,m1 ,m2 +λ2 V02 X [Dα;m3 m2 m1 m3 n + Dα;m3 m3 nm2 m1 α,m1 ,m2 ,m3 +Dα;nm2 m3 m3 m1 − Dα;m1 m2 m3 m3 n −Dα;m3 m3 m1 m2 n − Dα;m3 m2 nm3 m1 ] G m†1 , m2 δmα X +λ2 V02 [Dα;mm1 m3 m3 m2 + Dα;m3 m3 mm1 m2 α,m1 ,m2 ,m3 +Dα;m3 m1 m2 m3 m − Dα;m3 m1 mm3 m2 −Dα;m3 m3 m2 m1 m − Dα;m2 m1 m3 m3 m ] G m†1 , m2 δnα (B.3a) X −λ2 δmα D̄α;n + δnα D̄∗α;m 2 +λ V0 +λ2 V02 X X α (Dα;m1 m1 n δmα + Dα;m1 m1 m δnα ) α,m1 (Dα;m1 m1 m2 m2 n δmα + Dα;m1 m1 m2 m2 m δnα ) (B.3b) α,m1 ,m2 −iV0 G m† , (n − 1)† , n, n − 1 + iV0 G (m + 1)† , m† , (m + 1), n −iV0 G (n + 1)† , m† , n + 1, n + iV0 G m† , (m − 1)† , n, m − 1 n X +λ2 V02 G m†1 , m†2 , m3 , m4 α,m1 ,m2 ,m3 ,m4 · [δmα (Dα;m1 m3 nm4 m2 + Dα;nm3 m2 m4 m1 − Dα;m2 m4 m1 m3 n ) + δnα (Dα;m3 m1 mm2 m4 − Dα;mm1 m3 m2 m4 − Dα;m4 m2 m3 m1 m )]} (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) (1) (1) (1) ΓH0 + λ2 ΓB0 + λ2 V0 ΓB1 + λ2 V0 ΓB2 g1 = λ2 ν0 + λ2 V0 ν1 + λ2 V02 ν2 + iV0 + λ2 V02 π (g1 g1 ) + iV0 + λ2 V02 g2 . (B.4) Similarly we can derive an equation for g2 truncated at the order of V02 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) (2) (2) (2) ΓH0 + λ2 ΓB0 + λ2 V0 ΓB1 + λ2 V0 ΓB2 (g2 + π (g1 g1 )) = λ2 g1 + λ2 V0 g1 + λ2 V02 g1 + iV0 + λ2 V02 π (g1 g1 g1 ) + iV0 + λ2 V02 π (g1 g2 ) . (B.5) In the calculation illustrated in the main text, we truncated the bath operators at the order of V0 instead of V02 . In this case, the second equation of 137 B.2. Second-order equations of G1 and G2 the hierarchy becomes, 0 0 0 0 0 = itG m† , (n − 1)† , m , n + itG m† , (n + 1)† , m , n 0 0 0 0 +itG (m − 1)† , n† , m , n + itG (m + 1)† , n† , m , n 0 0 0 0 −itG m† , n† , m , n + 1 − itG m† , n† , m , n − 1 0 0 0 0 −itG m† , n† , m + 1 , n − itG m† , n† , m − 1 , n 0 0 +iV0 G m† , n† , m , n δm0 +1,n0 + δm0 −1,n0 0 0 −iV0 G m† , n† , m , n (δm+1,n + δm−1,n ) X 0 0 −iV0 G l† , m† , n† , l, m , n l=m±1,n±1 +iV0 −λ2 Xn X 0 0 0 G l† , m† , n† , l, m , n (B.6a) (B.6b) 0 l=m ±1,n ±1 h i 0 0 δmα G n† , m D̄α,n0 − G n† , n D̄α,m0 α i h 0 0 −δnα G m† , m D̄α,n0 − G m† , n D̄α,m0 i h 0 0 +δm0 α G m† , n D̄α,n − G n† , n D̄α,m h io 0 0 −δn0 α G m† , m D̄α,n − G n† , m D̄α,m X +λ2 Dα,l + D̄α,l h (B.6c) α,l 0 0 0 0 · δmα G l† , n† , m , n + δnα G m† , l† , m , n i 0 0 +δm0 α G m† , n† , l, n + δn0 α G m† , n† , m , l (B.6d) · · · to be continued 138 B.2. Second-order equations of G1 and G2 continued h Xn 0 0 +V0 λ2 δm0 α G m† , n Dα,m1 m1 n − G n† , n Dα,m1 m1 m α,m1 i 0 +G m†1 , n (Dα,nm1 m − Dα,mm1 n ) h 0 0 −δn0 α G m† , m Dα,m1 m1 n − G n† , m Dα,m1 m1 m i 0 +G m†1 , m (Dα,nm1 m − Dα,mm1 n ) h 0 0 +δmα G n† , m Dα,m1 m1 n0 − G n† , n Dα,m1 m1 m0 i +G n† , m1 Dα,n0 m1 m0 − Dα,m0 m1 n0 h 0 0 −δnα G m† , m Dα,m1 m1 n0 − G m† , n Dα,m1 m1 m0 io +G m† , m1 Dα,n0 m1 m0 − Dα,m0 m1 n0 h X n +V0 λ2 δm0 α Dα,m1 n0 m2 G m† , n† , m1 , m2 (B.6e) α,m1 ,m2 0 † † − (Dα,m1 m2 m − Dα,mm2 m1 ) G m2 , n , m1 , n i 0 + (Dα,m1 m2 n − Dα,nm2 m1 ) G m†2 , m† , m1 , n h − δn0 α Dα,m1 m0 m2 G m† , n† , m1 , m2 0 − (Dα,m1 m2 m − Dα,mm2 m1 ) G m†2 , n† , m1 , m i 0 + (Dα,m1 m2 n − Dα,nm2 m1 ) G m†2 , m† , m1 , m h 0 0 + δmα Dα,m1 nm2 G m†1 , m†2 , m , n 0 + Dα,m1 m2 n0 − Dα,n0 m2 m1 G n† , m†1 , m , m2 i 0 − Dα,m1 m2 m0 − Dα,m0 m2 m1 G n† , m†1 , n , m2 h 0 0 − δnα Dα,m1 mm2 G m†1 , m†2 , m , n 0 + Dα,m1 m2 n0 − Dα,n0 m2 m1 G m† , m†1 , m , m2 io 0 − Dα,m1 m2 m0 − Dα,m0 m2 m1 G m† , m†1 , n , m2 (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 V02 . 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 E,(1) C,(1) E,(1) as ∆1 and ∆1 , defined in Chapter 4. We will see that ∆1 in fact E,(0) involves ∆2 , which in turn needs Eq. (4.8) — the equation of G2 . A C,(1) similar equation is needed for estimation of ∆1 . C.1 E,(1) ∆1 from method 1 In order to estimate the accuracy of the first-order form of this approximation and also to get an estimate of accuracy for higher orders, we study the leading order of residues in terms of λ2 and ∆T T , both of which are assumed to be small in the following. Thus λ2 V0 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 λ2 ∆T ∆T , we drop the λ2 ∆T term in λ2 ν in Eq(4.24), λ2 ν = λ2 ν0 (T ) + λ2 ∆T ν,T , (C.1) and keep only the large term, λ2 ν0 (T ), which is independent of ∆T . Here d ν,T denotes formally a derivative of T on ν — dT ν. The general idea is then E,(1) to write down respectively equations for g1Ex and g1 , and then compare E,(1) E,(1) Ex the two to estimate ∆1 = g1 − g1 . In order to understand how an approximation to the next order improves the accuracy, we also want to E,(1) E,(0) E,(0) compare ∆1 to ∆1 = g1 − g1Ex , which is estimated in the same way E,(0) from the difference between the equations respectively for g1Ex and g1 . E,(1) E,(1) E,(0) E,(0) Ex Ex = gn − gn for More generally we define ∆n = gn − gn and ∆n n-particle Green’s functions gn . 141 E,(1) C.1. ∆1 from method 1 After dropping the gD term and the term which is proportional to λ2 ∆T , E,(0) E,(1) and keeping only up to the linear order of ∆T , g1Ex , g1 and g1 respectively satisfy (1) (1) Γ0 + Γ,T ∆T g1Ex = iV0 g2Ex + λ2 ν0 , (C.2a) (1) E,(0) E,(0) Γ0 g1 = iV0 g2 + λ2 ν 0 , (1) (1) E,(1) E,(0) Γ0 + Γ,T ∆T g1 = iV0 g2 + λ2 ν 0 , (1) (C.2b) (C.2c) (1) where Γ0 +Γ,T ∆T is the zeroth and first order in ∆T from Γ(1) of Eq(4.24). (1) Γ,T denotes formally a derivative of T on Γ(1) — one may use Eq(C.2a) and Eq(C.2b), E,(0) ∆1 E,(1) ∆1 d (1) dT Γ . E,(0) To consider ∆1 (1) −1 (1) Ex (1) −1 E,(0) = ∆T Γ0 Γ,T g1 + iV0 Γ0 ∆2 . , (C.3) can be estimated from Eq(C.2a) and Eq(C.2c), E,(1) ∆1 −1 (1) (1) E,(0) = iV0 Γ0 + Γ,T ∆T ∆2 , E,(0) (C.4) E,(1) where ∆2 is required. We find that roughly speaking ∆1 takes the E,(0) second term of ∆1 but drops the first term. Therefore, next we only need (1) −1 E,(0) to show that the second term, iV0 Γ0 ∆2 , is much smaller than the E,(0) first term, or equivalently, smaller than the whole ∆1 . E,(0) Estimation of ∆2 involves the second equation of the hierarchy, i.e. equation of G2 , which can be derived from substituting Eq(4.18) into Eq(4.8): 142 E,(1) C.1. ∆1 from method 1 D E D E 0 = it c†m c†n cm0 cn0 +1 + it c†m c†n cm0 cn0 −1 D E D E +it c†m c†n cm0 +1 cn0 + it c†m c†n cm0 −1 cn0 D E D E −it c†m c†n−1 cm0 cn0 − it c†m c†n+1 cm0 cn0 D D E E −it c†m−1 c†n cm0 cn0 − it c†m+1 c†n cm0 cn0 D E D E Xn +λ2 δn0 α dα;l c†m c†n cm0 cl + δm0 α dα;l c†m c†n cl cn0 l,α D E D E +δnα d¯α;l c†m c†l cm0 cn0 + δmα d¯α;l c†l c†n cm0 cn0 D E D E +δnα d∗α;l c†m c†l cm0 cn0 + δmα d∗α;l c†l c†n cm0 cn0 D E D Eo +δn0 α d¯∗α;l c†m c†n cm0 cl + δm0 α d¯∗α;l c†m c†n cl cn0 D E +iV0 c†m c†n cm0 cn0 δm0 +1,n0 + δm0 −1,n0 − δm+1,n − δm−1,n (C.5a) D E X −iV0 c†l c†m c†n cl cm0 cn0 l=m±1,n±1 +iV0 −λ 2 Xh X 0 D c†l c†m c†n cl cm0 cn0 E (C.5b) 0 l=m ±1,n ±1 D E D E δnα d¯α;m0 c†m cn0 + δmα d¯α;n0 c†n cm0 α D E D E +δm0 α d¯∗α;n c†m cn0 + δn0 α d¯∗α;m c†n cm0 D E D E −δnα d¯α;n0 c†m cm0 − δmα d¯α;m0 c†n cn0 D E D Ei −δn0 α d¯∗α;n c†m cm0 − δm0 α d¯∗α;m c†n cn0 (C.5c) D E D E Xn −λ2 V0 δm0 α c†m c†n cn0 Dα − δn0 α c†m c†n cm0 Dα D α D E D E +δmα − δnα c†m cm0 cn0 D̄α + δn0 α D̄α† c†m c†n cm0 D E D E D Eo −δm0 α D̄α† c†m c†n cn0 + δnα Dα† c†m cm0 cn0 − δmα Dα† c†n cm0 cn0 . (C.5d) c†n cm0 cn0 D̄α E This can be written in a compact matrix form as Ex Γ(2) g2Ex = iV0 g3Ex + λ2 g1Ex + λ2 V0 gD3 , (C.6) 143 E,(1) C.1. ∆1 from method 1 where vector g2Ex is defined as follows, h iT g2Ex = G2 1† , 1† , 1, 1 , G2 1† , 1† , 1, 2 , · · · , G2 N † , N † , N, N . (C.7) Note that some of the elements of g2Ex are naturally zero but we still keep Ex come from ordering respectively them in this vector. g3Ex , g1Ex and gD3 Eq(C.5b), Eq(C.5c) and Eq(C.5d) in the same way as g2Ex . The matrix Γ(2) 0 0 can be read off from Eq(C.5a), for example assuming m, n, m , n are all different and not equal to 1 or N , we have (2) ΓmN 3 +nN 2 +m0 N +n0 ,mN 3 +nN 2 +m0 N +n0 +1 = it. (C.8) Ex term in the following estimation of the Since λ2 V0 V0 we drop the gD3 accuracy. Therefore, the exact solution and the equilibrium solution satisfy respectively, (2) (2) Γ0 + Γ,T ∆T g2Ex = iV0 g3Ex + λ2 g1Ex , (C.9a) (2) E,(0) Γ0 g2 (2) where Γ0 E,(0) = iV0 g3 E,(0) + λ2 g 1 , (C.9b) (2) and Γ,T ∆T stand for the zeroth and first order in ∆T in Γ(2) . E,(0) Now ∆2 E,(0) ∆2 can be analyzed: (2) −1 2 E,(0) (2) −1 (2) (2) −1 E,(0) λ ∆1 . Γ,T ∆T g2Ex + Γ0 ∆3 + Γ0 = iV0 Γ0 (C.10) Here, in fact Γ(1) and Γ(2) have different dimensions. However, in this estimation of order of magnitudes, we ignore this difference and furthermore the matrices are regarded as constants with order 1. For the moment, E,(0) let us focus on the last term of ∆2 . Recall that we want to compare −1 (1) E,(0) E,(0) iV0 Γ0 ∆2 against ∆1 . Focusing only on the last term, we have E,(0) iV0 ∆2 E,(0) ≈ iV0 λ2 ∆1 , (C.11) E,(0) which is much smaller than ∆1 as long as V0 λ2 1. All together we have arrived at: E,(0) (1) −1 (1) Ex (1) −1 E,(0) ∆1 = ∆T Γ0 Γ,T g1 + iV0 Γ0 ∆2 (C.12) 144 C,(1) C.2. ∆1 while E,(1) ∆1 from method 2 (1) −1 (2) −1 E,(0) (2) −1 (2) Ex = Γ0 −V02 Γ0 ∆3 + iV0 ∆T Γ0 Γ,T g2 (2) −1 E,(0) 2 ∆1 . +iV0 λ Γ0 (C.13) E,(0) Most importantly here we see that ∆1 is multiplied by a small number E,(1) 2 λ V0 and then becomes a part of ∆1 . Judging from this it follows that, as long as λ2 V0 is a small number compared with t the method is very reasonable. As for the other two additional terms, they can be regarded as V02 g3Ex + V0 g2Ex ∆T . Therefore, the limit of maximum value of V0 , where −1 −1 this method is still applicable, is determined by g2Ex or g3Ex 2 . Such n limit could be much larger than 1 since roughly gnEx = g Ex — smaller for larger n. This explains why this method is applicable even for V0 larger than t, as confirmed 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 C,(0) C,(0) C,(1) C,(1) λ2 and V0 are small. We define ∆n = gn − gnEx and ∆n = gn − C,(0) C,(1) Ex Ex gn . Again we start from the equations of the three: g1 , g1 and g1 , and then compare the three equations while ignoring certain higher-order terms such as terms which are proportional to λ2 V0 . (1) In this case, after dropping the ΓD term and terms which are proporC,(1) tional to λ2 V0 , we find that g1Ex , g1C,0 and respectively g1 satisfy the following equations: (1) Γ0 g1Ex = iV0 g2Ex + λ2 ν0 , (1) C,(0) = λ2 ν 0 , (C.14b) C,(1) + λ2 ν 0 . (C.14c) (1) −1 Ex = −iV0 Γ0 g2 . (C.15) Γ0 g1 (1) C,(1) Γ0 g1 (C.14a) = iV0 g2 From Eq(C.14b) and Eq(C.14a) we find C,(0) ∆1 145 C,(1) C.2. ∆1 from method 2 Comparing Eq(C.14c) and Eq(C.14a), we get C,(1) ∆1 (1) −1 C,(1) = iV0 Γ0 g2 − g2Ex . (C.16) C,(1) C,(1) Note that the magnitude of ∆2 = g2 − g2Ex is in fact smaller than C,(0) C,(0) Ex the magnitude of ∆2 = g2 − g2 . 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 C,(0) are respectively the solutions of to λ2 V0 , it follows that g2Ex and g2 (2) Γ0 g2Ex = iV0 g3Ex + λ2 g1Ex , (2) C,(0) Γ0 g2 C,(0) = λ2 g1 (C.17a) . (C.17b) Comparing these two equations, we find that C,(0) ∆2 (2) −1 Ex (2) −1 C,(0) = −iV0 Γ0 g 3 − λ2 Γ 0 ∆1 . (C.18) Focusing only on the last term, we have C,(0) iV0 ∆2 C,(0) which is much smaller than ∆1 We arrived at: C,(0) ∆1 and C,(1) ∆1 C,(0) ≈ −iV0 λ2 ∆1 , (C.19) as long as V0 λ2 1. 2 (1) −1 Ex = −iV0 Γ0 g2 ∼ V0 g1Ex , (C.20) (1) −1 (2) −1 Ex (2) −1 C,(0) 2 2 = Γ0 iV0 Γ0 g3 + λ V0 Γ0 ∆1 ∼ V02 g1Ex 3 2 + λ2 V02 g1Ex . (C.21) C,(0) This agrees with the numerical results that ∆1 is proportional to V0 while C,(1) C,(0) 2 ∆1 is proportional to V0 . Most importantly, we see again that ∆1 is C,(1) 2 multiplied by a small number λ V0 and then becomes a part of ∆1 . Since 146 C,(1) C.2. ∆1 n from method 2 3 roughly gnEx = g Ex , the other term, V02 g3Ex ∼ V02 g1Ex , is also much C,(0) 2 smaller than ∆1 ∼ V0 g1Ex . 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 briefly 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 defined similarly with the usual complex numbers replaced by Grassman numbers[71]. In the §D.1, coherent states for single-mode bosonic systems is introduced and the simplest P -representation is discussed. In the §D.2, we discuss correspondingly the coherent-state representations for multi-mode bosonic systems. We also briefly review the more general positive P -representation. D.1 Coherent states for systems with a single mode A single-mode bosonic system is defined by a Hamiltonian H a, a† , where a, a† are respectively annihilation and creation operators of a particle with the particular mode. Its Hilbert space, expanded in occupation-number representation, is {|ni}, where n = 0, 1, 2, · · · . A coherent state |ξi is defined to be an eigenstate of annihilation operator, a |ξi = ξ |ξi (D.1) with eigenvalue ξ. Operator a is not a hermitian operator so ξ is not a real but a complex number (c-number). Let us define a displacement operator D (ξ) = eξa † −ξ ∗ a (D.2) , which equals to D (ξ) = e− |ξ|2 2 † ∗ eξa e−ξ a , (D.3) 148 D.1. Coherent states for systems with a single mode utilizing the Baker-Campbell-Hausdorff formula[70]. From there one can verify that |ξi = D (ξ) |0i , (D.4) where |0i is the vacuum state. This can be seen in the following. Let us express D (ξ) |0i in terms n (a† ) occupation-number states |ni and note |ni = √n! |0i, |ξ|2 † D (ξ) |0i = e− 2 eξa |0i ∞ n |ξ|2 X ξ n − 2 =e a† |0i n! n=0 − =e |ξ|2 2 ∞ X ξn √ |ni . n! n=0 (D.5) Therefore aD (ξ) |0i = e− = ξe− |ξ|2 2 ∞ X m=0 |ξ|2 2 ∞ X n=0 ξm p (m)! p ξn (n − 1)! |n − 1i |mi = ξD (ξ) |0i . (D.6) Inner product between two coherent states |ξi and |ηi is, hξ | ηi = eξ 2 2 ∗ η− |ξ| − |η| 2 2 , (D.7) so 2 |hξ | ηi|2 = e−|ξ−η| , (D.8) which is close to zero when ξ is very different 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 different. 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 Z 1 I= d2 ξ |ξi hξ| . (D.9) π 149 D.1. Coherent states for systems with a single mode This can be verified by considering arbitrary |ni such that Z 1 |ni = d2 ξ |ξi hξ |ni , π which becomes X1Z |ni = d2 ξ |mi hm |ξi hξ |ni π m Z X 2 ξ m (ξ ∗ )n 1 p d2 ξe−|ξ| p = |mi π (m)! (n)! m X = |mi δmn = |ni . (D.10) (D.11) m Here we have used 1 π Z 2 ξ m (ξ ∗ )n p d2 ξe−|ξ| p = δmn , (m)! (n)! (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] ZZ |η|2 |ξ|2 1 ρ= 2 d2 ξd2 ηe− 2 − 2 R (ξ ∗ , η) |ξi hη| . (D.13) π where R (ξ ∗ , η) = e |η|2 |ξ|2 + 2 2 hξ| ρ |ηi . (D.14) This distribution function R (ξ ∗ , η) holds generally for all ρ and it is called the R representation. However, since the overcompleteness and the nonuniqueness of the coherent-state representation, there are other expansions, which hold not for all but a large class of density matrices[70]. For example, the P representation and the Positive P representation are respectively defined as Z ρ = d2 ξP (ξ ∗ , ξ) |ξi hξ| , (D.15) and ZZ ρ= d2 ξd2 ηP (ξ, η) e |ξ|2 |η|2 + 2 −ξη 2 |ξi hη ∗ | . (D.16) 150 D.2. Coherent states for multi-mode systems Next we convert operator a and a† to be operators acting on such distribution function P , in for example, the P representation. For aρ, Z aρ = d2 ξP (ξ ∗ , ξ) a |ξi hξ| Z = d2 ξξP (ξ ∗ , ξ) |ξi hξ| , (D.17) so a → ξ. For a† ρ, noticing a† |ξi = † Z 2 d dξ |ξi ∗ d |ξi hξ| dξ a ρ = d ξP (ξ , ξ) Z d 2 ∗ = d ξ ξ − P (ξ ∗ , ξ) |ξi hξ| , dξ (D.18) d so a† → ξ ∗ − dξ . Similarly we need to consider also ρa and ρa† . Summarize all those operators, we have aρ ↔ ξP (ξ ∗ , ξ) , ρa† ↔ ξ ∗ P (ξ ∗ , ξ) , ∂ † ∗ P (ξ ∗ , ξ) , a ρ↔ ξ − ∂ξ ∂ ρa ↔ ξ − ∗ P (ξ ∗ , ξ) . ∂ξ D.2 (D.19) Coherent states for multi-mode systems For a multi-mode (countably many, denoted as N ) system, a coherent basis can be defined similarly as E |ξ~|2 ~ † ~ † ~∗ ξ~ = D ξ~ |0i = eξ·~a −ξ ·~a |0i = e− 2 eξ·~a |0i , (D.20) T where ξ~ = (ξ1 , ξ2 , · · · , ξN )T and ~a† = a†1 , a†2 , · · · a†N . a†l is the creation operator at site l, ξl is a c number and |0i is the vacuum state. Coherent states are over complete. One convenient choice of representation is the P -representation[70], which decomposes the identity as Z ED 1 I= d2 ξ~ ξ~ ξ~ , (D.21) π 151 D.2. Coherent states for multi-mode systems and it maps a density matrix ρ (a physical quantity A) to a function P ξ~∗ , ξ~ (A ξ~∗ , ξ~ ) such that Z ED ~ ρ = d2 ξP ξ~∗ , ξ~ ξ~ ξ~ , Z ~ hAi = d2 ξP ξ~∗ , ξ~ A ξ~∗ , ξ~ . (D.22) (D.23) If A is in normal order of al , a†l . † In this representation, operators a, a become differential operators on P ξ~∗ , ξ~ , al ρ ↔ ξl P ξ~∗ , ξ~ , ρa†l ↔ ξl∗ P ξ~∗ , ξ~ , ∂ † ∗ P ξ~∗ , ξ~ , al ρ ↔ ξl − ∂ξl ∂ ρal ↔ ξl − ∗ P ξ~∗ , ξ~ . ∂ξl (D.24) In this way, the operator-form QME can be mapped to a differential equation. 152 Appendix E Fokker-Planck equation and Langevin equation A standard Fokker-Planck equation with time-independent coefficients reads, 2 X X ∂ ∂ 1 ∂ ~ t = − ~ t . (E.1) P ξ, Aj ξ~ + Dij ξ~ P ξ, ∂t ∂ξj 2 ∂ξi ∂ξj j i,j When there are higher-order derivative terms besides the above first two terms, it is called a generalized Fokker-Planck equation. ThePawula theorem[98] ~ t ) non-negative, states that in order to keep the distribution function (P ξ, there has to be either only 2 terms or infinity terms. In the later case, unless the equation is exactly solvable usually the generalized Fokker-Planck equation is truncated first 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 FokkerPlanck 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 n o equation besides brute-force numerical solution in the whole space of ξ~ . Converting 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 coefficients Aj ξ~ and Dij ξ~ , 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 define a correspondence between the standard Fokker-Planck equation and the Langevin equation. Assuming Dij ξ~ is positive semidefinite such that there is a matrix G satisfying D ξ~ = G ξ~ G† ξ~ , (E.4) then gij ξ~ = Gij ξ~ , hi ξ~ = Ai ξ~ . (E.5) (E.6) In numerical simulation, one works with directly dξi = hi ξ~ dt + gij ξ~ dWj (t) , where hdWi (t)i = 0, hdWi (t) dWj t0 i = δij dtδt,t0 . (E.7) (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, hxl i and hxl xm i etc.. Let us work out the first derivative term for the above one-variable average, X ∂ d hξl i = −hξl Ai i dt ∂ξj j X ∂ X ∂ = −h ξl Ai i + hAi ξl i = hAl i. ∂ξj ∂ξj j (E.9) j Here we have used the following property: for arbitrary F ξ~ , as long as limξ −→∞ F ξ~ does not diverge, we have j Z ∂ ∂ h ~ ~i ~ h F ξ i = d2 ξ~ F ξ P ξ = 0. ∂ξj ∂ξj (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 diffusion matrix D ~ We further require Γ is positive definite in or~ = Γξ. and a linear drift term A der to have a non-divergent distribution function. Although time-dependent solution is also possible, here we focus on stationary solution, P = (2π)− 2 (Detσ)− 2 e− 2 ξ N where σ= X l,m 1 1 ~T σ −1 ξ~ , D E 1 l m v |D| v |ul ihum | , κl + κ∗m and Γ= X κl ul ihv l . (E.11) (E.12) (E.13) l The last expression defines v l (|ul i) 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 equations. E.3 Stochastic difference equation for truncated generalized Fokker-Planck equations A generalized Fokker-Planck equation with time-independent coefficients truncated at the third-order derivatives reads, 1 X ∂2 X ∂ ∂ P ξ~ = − Aj ξ~ + Dij ξ~ ∂t ∂ξj 2 ∂ξi ∂ξj j i,j 3 X 1 ∂ ~ + Dijk ξ P ξ~ . (E.14) 6 ∂ξi ∂ξj ∂ξk i,j,k Now let us solutions of this generalized Fokker-Planck equation. For such generalized Fokker-Planck equation an equivalent stochastic difference 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 difference equation for truncated GFPEs used. The basic idea is to use proper random variables to mimic the effect of the third order term just like what was done for the second derivative terms. For example, ∂ Ai → h∆ξi i = Ai ∆T, ∂ξi (E.15) ∂2 Dij → h∆ξi ∆ξj i = Dij ∆T, ∂ξi ∂ξj (E.16) ∂3 Dijk → h∆ξi ∆ξj ∆ξk i = Dijk ∆T. ∂ξi ∂ξj ∂ξk (E.17) 1 Therefore, there will be terms like ηijk (∆T ) 3 appearing in the stochastic difference equation to make sure the above average holds, 1 1 ∆ξi = hi ∆t + gi (∆t) 2 + ki (∆T ) 3 . (E.18) just like gi happens to be gij wj such that Eq. (E.16) is satisfied, ki has to be defined properly to obey Eq. (E.17). Generally relation between ki and Dijk is still missing, but specific forms for example equations can be found in Ref.[102]. 156 Appendix F Solving the Redfield 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 Redfield equation, Eq. (4.2) into equation of motion of Green’s Dfunctions, which is defined by singleE † † particle Green’s function G l ; n = al an ρ (∞) (G1 for short), two-article D E Green’s function G i† , j † ; l, n = a†i a†j al an ρ (∞) (G2 for short) and so on. Here again the density matrix ρ (∞) refers to a stationary state and it is defined by Lρ (∞) = 0. Then generally for interacting systems the stationary equation becomes an equation hierarchy, in the first 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 difference between the Redfield 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 coefficients Dγ,m should then be substituted properly for the local-operator Lindblad equation. The operator-form Redfield equation reads ∂ρ = LH0 ρ + U LV ρ + λ2 L0B ρ + λ2 U L1B ρ, ∂t (F.1) 157 F.1. Exact solution on non-interacting systems where LH0 ρ = −i N −1 h X a†l al+1 + a†l+1 al , ρ i (F.2a) l=1 LV ρ = N −i X h i a†l al a†l al − 1 , ρ 2 l=1 h i h i o n X (0) Dγ,m a†γ , am ρ + D̄γ,m aγ , a†m ρ + h.c. LB ρ = − γ=L,R nh X (1) LB ρ = − Dγ,m1 m2 m3 γ,m1 ,m2 ,m3 a†γ , a†m1 am2 am3 ρ (F.2b) (F.2c) i i o h + aγ , a†m3 a†m2 am1 ρ + h.c. . (F.2d) Let us now transform the above Redfield equation in the operator form into equations of Green’s functions by multiplying a†l an at the RHS of Eq.(3.30) and then performing a trace, 0 = −iha†l an+1 i − iha†l an−1 i + iha†l+1 an i + iha†l−1 an i i X h +λ2 D̄γ,m − Dγ,m δnγ ha†l am i + δlγ ha†m an i γ,m +λ2 X D̄γ,n δlγ + D̄γ,l δnγ (F.3a) (F.3b) (F.3c) γ +λ2 U X γ,mi F.1 h Dγ,m1 m2 m3 δnγ +iU ha†l a†l al an i − iU ha†l a†n an an i ha†m1 am2 iδlm3 + ha†m1 am3 iδlm2 i +δlγ ha†m2 am1 iδnm3 + ha†m3 am1 iδnm2 . (F.3d) (F.3e) 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 (0) case. This is an N 2 -dimension linear system for unknown G1 , where (0) refers to U = 0, (0) (0) ΓH0 g1 + λ2 ΓB0 g1 = λ2 ν (0) . (F.4) Here we arrange all unknown G1 into a vector g1 = (G1 (1, 1) , G1 (1, 2) , · · · , )T . Elements of the two matrices and the vector can be read off from Eq. (4.7). 158 F.2. Approximate solution on interacting systems For example, (ΓH0 )l∗N +n,l∗N +n+1 = −i, (F.5) and when l = 1 (or l = N , or n = 1, N ), (0) ν1∗N +n = λ2 D̄1,n , (F.6) otherwise, it is zero. Solving this N 2 -dimension linear system we find 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 a†i a†j al an and then performing a trace. In fact, if one do so, one will find 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 first order by using the cluster expansion[101], ha†l a†l al an i = 2ha†l al iha†l an i + 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, iU 2ha†l al iha†l an i − iU 2ha†l an iha†n an i. (F.8) Therefore, Eq. (F.3) becomes a closed equation, (1) (1) (1) (1) ΓH0 g1 + λ2 ΓB0 g1 + λ2 U ΓB1 g1 = λ2 ν (0) + iU π g1 , where (1) π g1 l∗N +n (1) (1) = 2 g1 g1 l∗N +n n∗N +n (1) (1) −2 g1 g1 . l∗N +l l∗N +n (F.9) (F.10) 159 F.2. Approximate solution on interacting systems This equation can be solved iteratively, [n+1] g1 h i (0) (1) −1 [n] λ2 ν (0) + iU Π g1 , = ΓH0 + λ2 ΓD + λ2 U ΓD (F.11) (0) with initial solution g [0] = g1 for the corresponding non-interacting system. Then the solution is (1) [n] g1 = lim g1 , (F.12) n→∞ [n] [n−1] which in practice stops at large enough n such that g1 − g1 enough. is small 160