統計物理學基礎研究新進展

王 馨， 張 江， 吳金閃

（北京師范大學管理學院，北京100875）

1 系統理論與統計物理學

系統理論與統計物理之間聯系緊密，首先統計物理學所關注的系統一直是系統理論的重要研究對象，其統計物理學是系統理論研究方法的重要來源.系統理論是研究相互作用的大量個體構成的系統科學，研究這樣的系統的描述方法、性質、結構與層次，功能以及干預調控.其中一個核心的研究方向就是構建從微觀描述到宏觀性質的橋梁，以及反過來從宏觀現象探索微觀機制.常常人們還額外關注系統宏觀性質發生定性變化的現象.這樣的現象稱為涌現現象或者相變.在具體系統的研究中，這樣的現象往往具有重要的實際意義.這些問題正好就是統計物理學的研究內容，只不過傳統統計物理學更注重物理系統中這些問題的研究，而系統科學則研究更廣泛的各個具體學科中的這樣的系統.因此，統計物理學是系統科學，或者說是復雜性研究的研究方法的主要來源.另一方面，包含相互作用的大量個體的物理系統，自然也是系統科學的研究對象之一.

從系統科學的角度研究統計物理學更多地是從非平衡統計的角度展開.平衡態統計物理學給出一個這樣的大圖景（盡管作者個人不認為這樣的圖景能夠從平衡態統計物理學得到，但是在此還是用通常的說法——“熱寂說”）：在低溫條件下，系統會產生結構；而在高溫條件下，系統趨向于沒有結構的各向同性狀態.然而，結構的涌現是很多系統的普遍現象，尤其是在生物、社會、經濟系統之中.那么，什么系統，在什么條件下會產生結構，自然就成了系統理論研究的核心問題.既然平衡態統計物理學給出了不一樣的圖景，那么Prigogine等就認為與外界存在能量流和粒子流，也就是開放性和非平衡性，是結構產生的根源.這個問題與從可逆動力學到不可逆熱力學的問題也是聯系在一起的.因此，非平衡與不可逆性，從系統科學創立階段開始，一直是其最核心和最基礎的研究問題.

2 統計物理學的基本問題及其研究進展的粗略圖景

平衡態統計物理學的基本理論是系綜理論.以哈密頓量H所描述的處于溫度T的正則系綜為例，平衡態統計物理學認為，其狀態由Boltzmann分布

描述，其中，β＝1／kBT（kB為參數），歸一化常數Z被稱為配分函數.一個正則系綜通常被認為是處于一個大熱浴中的一個小系統，與熱浴存在著能量交換但沒有粒子交換.平衡態還表示，這樣的一個小系統經過與大熱浴的長時間接觸，無論其初始狀態如何，都會“演化”到一個穩定態上.這樣的統計力學演化與基于第一原理（牛頓方程或者薛定諤方程）的動力學演化是不一致的.前者通常不可逆，存在穩定態；后者是可逆的，往往不能給出一個穩態.這就是著名的從微觀可逆性到宏觀不可逆性的問題.這就是統計物理學基本問題之一：平衡態如何從第一原理推導出來.基于第一原理的演化，或者說動力學演化遵循Newton方程或Schro¨dinger方程，統一記為

其中，ρ稱為密度分布函數（經典力學）或者密度矩陣（量子力學），LH稱為Liouville算符.它完全由系統哈密頓量H決定.

通常認為，這樣的正則分布是可以從微正則分布——孤立系統的狀態符合等能面上的等幾率分布推導出來，因而，上面的問題就成了如何從第一原理得到微正則分布.由此，如果能夠證明微正則分布或者找到系統出現微正則分布的充分必要條件，則正則分布也得到了證明.研究微正則系統比正則系統要簡單得多，因為前者是孤立系統，而后者必須考慮熱浴.因此，很多的工作都從動力學演化的角度來探討微正則分布出現的充分必要條件.微正則分布是指與外界不存在能量與物質交換的孤立系統，其狀態是給定能量所對應的等能面上的均勻分布.這是一個很可以理解的假設：孤立系統的能量自然是常數（其它的動力學守恒量也是常數，這會把等能面分成幾個獨立的部分，在此為簡化問題假設沒有這樣的守恒量），因此系統只能處于給定等能面上；等能面上的各個點沒什么不同，在不知道確切初始條件的情況下，也只能假設它們都等價了.在這方面，最近的研究表明，不用微正則分布，只要考慮了系統與熱浴構成的大系統（有時候稱為宇宙），大系統甚至處于純態而不用處于微正則分布，對于絕大多數的純態，子系統的約化密度矩陣都處于正則分布，這被稱為正則典型性［1］.這個結果相當于用“絕大多數”純態代替了微正則分布.這方面的研究看起來為統計物理學基本問題的研究開辟了一個新的方向，引起了很多研究者的興趣.

然而，最近的研究注意到這樣一個事實［2－3］：不增加任何別的假設，正則分布不能通過微正則分布（或者以上的絕大多數態）推導出來.通常教科書式的標準做法是把子系統和熱浴一起看成一個孤立系統，然后從這個孤立系統的微正則分布，加上熱浴的態密度是隨著能量指數增長的假設，導出子系統的正則分布.注意，這個推導沒有要求熱浴自身處于熱平衡態，只要求熱浴和子系統構成的整體系統是孤立系統，符合微正則分布.盡管這個關于熱浴態密度的假設也很可以理解，但是并不是所有的熱浴都具有這個性質.是否只有這樣的熱浴才能驅動系統演化到熱平衡，還是一個問題.對于考慮了正則典型性以后的表述，也就是說整體系統處于絕大多數純態都可以（在同樣的態密度假設下）得到子系統的正則分布，所要求的條件看起來弱了很多.但是，目前尚無法確定，正則典型性是否可以保證系統在演化過程中保持正則分布.

基于以上兩方面的原因，于是有研究者就試圖直接從系統與熱浴耦合的整體系統出發，來討論系統演化到定態是正則分布的充分必要條件.鑒于正則分布在描述實際系統上很高的適用性，這樣的條件應該不是非常特殊.研究者們希望證明類似這樣的一個定理：只要系統哈密頓量滿足某些條件，熱浴哈密頓量滿足某些條件，相互作用形式滿足某些條件，熱浴的大小（粒子數、自由度數目）滿足某些條件，熱浴的初始狀態滿足某些條件，系統就會在整體系統的動力學演化下趨向平衡.當然，如果有必要，系統的初態也可以做約束.所有的這些條件越容易滿足越好，而且最好也是必要條件.

回答這個問題是幾代物理學家的夢想.20世紀60年代發展起來的投影算符方法，可以從整體系統的動力學方程得到子系統的有效動力學方程，通常具有如下形式，也被稱為動理學方程

它在第一原理動力學方程的基礎上增加了一項來自于熱浴的算符LB，同時在某些條件下也有可能改變原來的子系統本身的算符.

從這個方程出發，假設熱浴很大，而且處于熱平衡態，以至于改變子系統狀態時其自身狀態不變，那么從得到的動理學方程可以計算出定態解，而且這個定態解正好就是熱平衡Boltzmann分布.也就是說，至少有一個充分條件：處于熱平衡的理想無限大熱浴，弱耦合（耦合強度小，Markovian近似適用，耦合項的二階微擾適用），則子系統符合熱平衡分布.當然，這個結果并沒有回答前述統計物理學的根本問題——從第一原理推導出熱平衡分布.實際上，只是從一個熱平衡（熱浴的）導出了另一個熱平衡（子系統的）.甚至這個結論本身，也還有很多欠缺之處：有限但是非常大的熱浴是否可行，Markovian近似和二階微擾是都必要等.更一般地說，熱浴是否必須處于平衡態，處于別的某種狀態是否也可以導致相應的定態（當然此時就不一定是熱平衡態），還是說熱平衡態具有某種特殊的穩定性，只有它才能從熱浴“傳到”子系統上去，這些問題都還沒有回答.

最近有研究者宣稱他們已經回答了這個問題：只要滿足非常簡單的條件，熱浴的自由度大于子系統的自由度［4］，子系統就會趨向定態，而且在滿足進一步的條件下，這個定態就是熱平衡態.然而，這些作者判斷子系統趨向定態的方法是定義子系統密度矩陣的長時平均

其中，ρS＝trB（ρSB）為熱浴部分自由度求跡以后的密度矩陣.文獻［4］的作者們發表了很多相關的文章，有興趣的讀者可以做一個系統地跟蹤閱讀.如果ρS存在長時穩態極限，那么上面定義的極限確實是其長時極限.不過，反之就不正確，也就是說，就算如上定義的極限存在，也不能肯定ρS存在長時穩態極限.有關這方面的反例以及相關工作正在展開.

統計物理學的第二個基本問題是系綜理論為什么能夠描述熱力學測量.系綜理論是一個關于狀態分布函數的理論，也就是一系列系統構成的整體系綜的理論，熱力學測量原則上是在測量一個演化的系統或者一個處于定態的系統.一個由概率論所描述的隨機個體，某個概率分布描述了這個個體，含義是從對大量個體的獨立測量構成的整體來看，測量的統計特征與分布函數相符.但是，在熱力學測量中，通常只在一個系統上做測量，最多重復很少的次數，就把測量量與基于分布函數的計算量相比，而且符合得很好，為什么？通常統計物理學教科書的回答是：對于處于熱平衡的系統，其長時平均等于系綜平均.然而，通過對具體系統的計算，發現實現長時平均等于系綜平均所需要的時間（以后稱為回復時間，相當于遍歷時間的尺度）遠遠大于實際測量時間［3］.這使得這個曾經認為已經得到回答的問題成了新問題.本文初步研究發現，如果只考慮系統自身的動力學演化，回復時間確實非常長，但是考慮子系統與熱浴耦合之后的動理學方程，這個回復時間就小得多.因此，可以認為，熱力學測量就是對分布函數的多次抽樣，在每一次抽樣之后，系統都要經過一個回復時間來重新回到成為平衡分布中的一個樣本點，而這個回復時間不是完全由系統自身決定的遍歷時間而是由系統和熱浴共同決定的回復時間.目前，具體系統上的計算還在進行中.

統計物理學的第三個基本問題是如何處理非平衡定態，這也是本綜述的重點.非平衡定態的計算是研究輸運問題的基礎.考慮例如熱輸運問題，一個準一維系統連接到兩個不同溫度的熱浴上，經過一段時間，系統到達定態，其上有熱流通過.希望能夠計算熱流與溫度、溫度差、系統材料之間的關系，還希望能處理電流、粒子流、自旋流等輸運現象.更廣義地說，結構的產生，也是系統的非平衡定態的性質.因此，對非平衡定態的描述方法的探索是系統科學基本理論研究的一個非常重要的方面.

如果有第一個問題的完整答案，這個問題也就差不多解決了.既然知道某種條件下耦合到一個大熱浴的子系統會演化到平衡態上去，只要讓這個子系統耦合到兩個大熱浴，那么可見子系統就會演化到一個新的定態上去，而這個定態就是要找的非平衡定態.但是，沒有這樣的一個充分必要條件，因此，也不能從根本上解決這個問題.換一個角度，沒有解決平衡態的由來這一根本問題，這并不妨礙用系綜理論來計算平衡態物理量.同樣，對于非平衡定態，也希望先得到一套可計算的框架，然后有可能的時候再回去解決根本問題.

在這樣一個指導思想下，發現有可能推廣投影算符方法到非平衡態的計算.從耦合到單個熱浴的子系統的動理學方程，可以類似地得到耦合到多個熱浴的子系統的動理學方程，然后求解這樣的動理學方程的定態.這個想法首先由日本學者Saito在耦合到諧振子浴的一維諧振子鏈上實現［5］.在這個工作中，Saito確實發現了非平衡定態，得到了半解析半數值的解，注意到一維諧振子鏈其實可以轉化為無相互作用系統，這使得Saito能夠求解這一動理學方程.對于相互作用系統費米子系統，這個關于密度矩陣的方程有4N個變量；對于相互作用系統波色子系統，這個方程有無窮多個變量.求解這樣的方程是一件非常困難的工作.后來的工作在運用這個方法討論具體系統和具體物理問題的同時，都把大量的工作用在了尋找求解這一方程的高效方法上.目前來說，除了Runge－Kutta方法、直接對角化等直接方法之外，還發展了基于局域算符展開的密度矩陣重整化方法［6］、Monte Carlo方法［7］、Green函數方法［8－9］、相干態量子隨機微分方程方法［9］等一系列間接方法.本文將在下面展開論述其中一部分方法.這些方法的對比，以及在具體系統上的應用等問題都還沒有得到深入的研究，在這些方面還有大量的工作要開展.如圖1所示的是連接到兩個熱浴的準一維系統.其中TL，TR為左右熱浴的溫度，jh為流密度.宏觀輸運定律假設系統上存在著溫度或者電壓的梯度，而從微觀角度希望能夠從第一原理出發，計算系統的穩態以及各個物理量.

圖1 連接到兩個熱浴的準一維系統Fig.1 Sketch of a typical setup of heat conduction

在具體物理問題的研究方面，正常導熱性的充分必要條件，也就是研究在什么樣的系統上熱傳導滿足傅里葉定律，一直是非平衡統計研究的一個重要方向.傅里葉定律是一個適用于宏觀客體的實驗定律，其經典力學基礎和量子力學基礎都是不明確的，適用范圍也不明確.費米等在經典FPU模型上的工作，以及后來大量的在FPU上的工作，都企圖從動力學演化的角度回答這個問題.但是迄今為止沒有找到正常熱導的充要條件，或者從動力學演化到熱平衡分布的充要條件.當然也得到了一定的結果，例如完全可積系統不能實現熱平衡，沒有正常熱導.其物理圖像可以作以下理解：完全可積系統存在非常多的動力學守恒量，也就是說系統可以看成獨立的很多個本征運動模式，每一個運動模上的能量是獨立的，不能相互混合傳遞，因此不能實現熱平衡.在僅考慮動力學演化的范圍內，這個圖像似乎是有道理的.然而，所有的數學形式以外的物理學理解，也就是所謂的物理學圖像，都可能是欺騙性的.例如量子力學是不是要有現實性，以及什么樣的現實性，等等這樣的問題，都可能是用人們用習慣了經典力學圖像的大腦來理解量子力學世界的結果.在這里，唯一可靠的就是Hilbert空間的結構和薛定諤方程本身，除了數學結構，沒有獨立的有意義的非欺騙性的“物理圖像”.人們能夠問的問題是為什么量子力學的實驗就要求有這么奇怪的數學形式，而不是形成所謂的量子力學的獨立于數學形式之外的物理圖像.回過頭來看這個關于可積性的物理圖像，如果系統存在與外界的耦合，而這個往往是討論趨向熱平衡或者熱傳導必須的，這些不能直接混合傳遞的本征模之間，完全可以通過與外界交換能量來實現間接地混合與傳遞.因此，不考慮熱浴，不考慮動理學演化的關于熱傳導或者趨向熱平衡的討論是不完整的.從這個角度來說，有了這個非平衡定態研究的框架，才能夠正確地討論正常熱導的充分必要條件.在經典低維系統的正常導熱性方面，基于具體系統的研究結果，中國學者趙鴻等提出了動量守恒是反常熱傳導的根源的假說，以及最近關于非對稱勢導致正常熱傳導的理論.在此，主要關注量子系統，對經典系統有興趣的讀者請參考相關的綜述文獻［10］.

另外一方面，本文還可以用同樣的框架研究電流與自旋流的問題.目前電子輸運的實驗研究已經表明在微觀尺度，歐姆定律不成立，電阻與導體長度之間沒有正比關系.這就要求有一個基于第一原理的新的處理電子輸運過程的框架.除了動理學方程這一方法之外，通常（或者說更經常地）可用另外一套半原理半唯象的處理方法：Landauer－Buttikker公式［11］以及非平衡格林函數理論（NEGF）［12］或者Kubo（久保）公式［13］.Kubo公式只適用于無限大系統，近平衡，而且是外勢驅動的輸運問題［14］.在此本文不展開討論.Landauer－Buttikker公式的物理圖像是零溫無相互作用系統的輸運，可以用兩個有偏分布函數來描述.假設左端化學勢μL比右端μR高，那么右行粒子（來自于左端電極）占滿到μL的能級，而左行粒子占滿到μR的能級.自然右行粒子多于左行粒子，算出這個差別和相應的電荷，就得到了電流.這個方法的基本計算是先找到所有的散射波函數，然后按照各自的分布把左行和右行散射波函數的電流加起來，即

式中，q為電荷量；?為常數；T（E）為能量是E的右行散射波函數的透射系數.通常在計算這樣的散射波函數過程中，還假設這個準一維系統是無窮長的.在這個框架下面，如果要問子系統的狀態（由約化密度矩陣描述）是什么，那么它可以認為是如下整體系統的密度矩陣的約化：整體密度矩陣是左行部分和右行部分的直積，其中左行（右行）部分占滿到μR（μL）的能級.在這個理論框架中，這樣的密度矩陣是一個假設，盡管有道理也可以理解，但不是從第一原理推導出來的.實際應用中，人們發現對大量的系統，這個理論與實際測量符合得比較好.這是針對無相互作用子系統的理論，考慮到相互作用時，這個理論就不適用了.一個代替它的理論是NEGF.NEGF是一個微擾計算理論，沒有非微擾形式.從不相連的兩個半無窮長鏈（每一側的鏈都有自己的溫度和化學勢，因而可以算出自己的格林函數）出發，把連接部分當成微擾項加入系統，然后計算這個微擾對系統格林函數的影響.由于系統在微擾前后的定態上有定性區別，必須要用閉路格林函數來做微擾計算，而不是通常的場論微擾展開方法.這個理論也是經驗性的，盡管在無相互作用系統上它被認為與Landauer－Buttikker公式等價.考慮一個電壓線性下降的導體，不知道經過微擾計算，能不能從階躍函數型的電壓降得到線性下降的電壓降.當然，這個方法的優點是，可以利用格林函數的方法直接研究相互作用的效果.實際計算過程中，對于相互作用往往要加入額外的近似，例如Hartree－Fock近似、集團展開等.這對此方法的精度帶來了額外的影響.目前，在某些系統上，這個方法的表現還可以，但是在另一些系統上計算結果存在數量級上的差別［15－16］.用這個非平衡定態計算的一般框架，就可以直接處理相互作用系統，而且避免產生兩重微擾的問題.再加上這個框架基本上是第一原理性的，不是通過微擾形式定義的，也可以為改進和檢驗以上方法提供參考.

以上簡略地介紹了統計物理學基礎研究的大圖景以及輸運問題的基本處理方法.然而，綜述以非平衡定態的框架來研究輸運問題是本文的重點，因此，后面將集中在這個方面，而不再涉及處理輸運問題的其它方法.另外，在介紹完這個基本框架之后，在最后部分本文也會提到這個計算框架，也有助于解決統計物理學前面的兩個基本問題.統計物理學，乃至整個物理學的中心問題之一——如何處理相互作用的多體系統，也會是貫穿整個綜述的另一條線.

3 非平衡定態計算的基本框架

首先，簡單地介紹一下投影算符技術和動理學方程方法，然后推廣到應用于輸運問題的情形；其次，舉幾個研究實例，展示這個基本框架；最后，將在簡略地介紹幾個計算方法的基礎上，詳細地闡述其中幾個高效的計算方法.

3.1 投影算符技術與動理學方程

投影算符技術的細節可以參考Kubo的經典統計物理學教材［17］.其基本思想是從子系統與熱浴的整體動理學方程出發，積分掉熱浴的自由度，得到子系統的約化密度矩陣的有效運動方程.考慮整體系統

式中，HS，HB為子系統和熱浴的哈密頓量；HSB為子系統與熱浴之間的耦合項.整體系統的密度矩陣ρT滿足方程式（2），其中，用到的Liuville算法是整體系統

定義子系統約化密度矩陣

式中，trB為在熱浴B的狀態空間作求跡運算，對于定義在子系統上的物理量AS，有

因此，只要能夠求解ρ而不用整個ρT，就可以知道所有子系統上的觀測量.定義投影算法

可以檢驗P確實是投影算法P2＝P，定義Q＝IP，其中I是單位算符.可由式（7）得到PρT和QρT的方程

第二個方程可以形式上求解，然后代入第一個方程，結合初始條件

得到

這是一個關于ρ的方程，當然其中還有整體系統的Liuville算符項L以及算符Q.這個方程描述了依賴于歷史的非Markovian過程：方程中包含PρT（T）的積分.如果只關心定態而不是暫態的問題，可引入Markovian近似（粗略地說是讓t→∞，ρ（t－T）→ρ（t））來大大簡化這個方程.另外，指數上的L算符也很難處理，將其拆分成幾個部分，例如

e（t－T）QL＝e（t－T）（1－P）L＝e（t－T）L（1＋O（P L））（14）然后在合理的條件下只保留第一項.對于方程，這相當于保留到耦合項的二階微擾.細節請參考文獻［9，17］.考慮這兩個假設之后，這個方程可以得到進一步簡化.下面以二次量子化形式下產生湮滅算符形式的哈密頓量為例，寫出這一方程簡化以后的具體形式，同時將方程推廣到多個熱浴的情形.

考慮準一維子系統與多個熱浴通過第ν個格點與第ν個熱浴耦合，其中，第ν個熱浴的參數為溫度Tν，化學勢μν，耦合常數Vk，ν和λ，得

在這里aν（T）＝eiHSt aνe－iHSt是算符aν的Heisenberg繪景形式.這里已經作了Markovian近似，并只考慮了耦合項的二級微擾.這個方程有時候也被稱為Redfield方程（RE）［18－20］.有關非平衡定態的計算可以直接從這個方程開始.這里nk，ν是第ν個熱浴的第k個模上的平均粒子數

其中，±分別對應著費米子和波色子熱浴.為簡單起見，以后在不引起混淆的情況下〈nk，ν〉直接記為nk，ν.與孤立系統的動力學方程相比，這個考慮了熱浴的動理學方程僅僅多了后面的有關算符和的項，而且所有熱浴的信息都包含在這些算副符中，因此，稱這些算符為熱浴算符.下面對具體系統給出幾個熱浴算符的例子，并作簡單計算來闡述一個需要注意的細節.

3.2 舉例

從簡單到復雜，從平衡到非平衡作幾個具體系統的計算.一方面，用實際系統實現了上面的一般理論的計算；另一方面，從這些例子中將得到一些定性的結論.這些結論對于將投影算符方法推廣到非平衡系統是有意義的.

3.2.1 單熱浴中的單格點費米子系統：直接求解

考慮由哈密頓量描述的如圖2所示的單格點費米子模型

對此得到動理學方程

和熱浴算符

這個組合表達式J由熱浴的態密度、子系統本征模的能量以及耦合常數決定.

圖2 單格點系統與單熱浴耦合草圖Fig.2 Sketch of a one－site system coupled to a heat bath

在單格點單本征模系統上，為簡單計，設J＝1.定義密度矩陣

并把它表示成為一個4維向量

湮滅算符與熱浴算符可以表示成

把這些算符以及密度矩陣都代入式（20），可以得到

其中，Γ是一個4×4的矩陣，或者更一般地說是一個4N×4N維的矩陣.對于一個N格點的無自旋費米子子系統，這個方程的演化完全由Γ的譜決定，其定態P0由零本征值對應的本征向量給定，即

當且僅當Γ只有一個零本征值，而且其它本征值小于零時，方程具有唯一定態解.對于單熱浴中的單格點費米子系統可以驗證Γ滿足這一條件.計算出來的平衡態也正好是巨正則系綜對應的熱平衡態同時，注意到λ與定態沒有關系，以后會看到，非平衡

態的細節與λ是有關系的.還要注意到對于平衡定態，系統自身動力學LHs項不直接進入動理學方程（它決定采用什么基矢——它們是HS的本征態）.它實際上只影響HS表象下密度矩陣的非對角項，而平衡態這樣的非對角元素為零.這一點，在非平衡定態中，也會發生變化，對系統自身動力學LHs將有重要影響.由于空間維數的問題，波色子系統的相應計算會更復雜，但是整體框架和物理圖景是一致的.

單格點系統的Redfield方程式（20）可以拋開熱浴算符和，寫成

這一形式被稱為Lindblad方程.有的研究者也用這個方程而不是Redfield方程來研究非平衡輸運現象.可以看到對于單格點系統，兩者是等價的.但是，以后會看到對于多格點系統，Lindblad方程實際上是在Redfield方程的基礎上引入了進一步的近似：稱之為弱連接條件或者局域哈密頓量條件.稱這個推廣以后的對于多格點系統的Lindblad方程為局域算符Lindblad方程（LOLE）.在弱連接條件不滿足的情況下（格點之間的hopping強度大于或者類似于格點上粒子的能量），Lindblad方程甚至不能給出正確的平衡態.

3.2.2 單熱浴中的雙格點費米子系統：選擇RE而不是LOLE

考慮如圖3所示的單熱浴中的雙格點緊束縛模型

其中，格點1與熱浴連接

圖3 單熱浴中的雙格點緊束縛費米子系統草圖Fig.3 Sketch of a two－site system coupled to a heat bath

按照方程式（13），有

其中，熱浴算符為

有的研究工作對于多格點系統，還是簡單地直接推廣單格點系統的Lindblad方程式（31），即

這里平均粒子數由ν格點的參數（化學勢、溫度、onsite能量）決定：n＝n（εν；Tν，νν）.可以看到這時Redfield方程與Lindblad方程并不相同.后者是前者在弱連接（T?ε）下的近似.再看兩個方程（具體來說這里就是比較式（34）與式（36）的特例——ν＝1的情況）給出的定態是不是相同，并把定態與熱平衡分布做比較.

取上述系統HS的本征模

其中

利用這些本征模的粒子數表象｜n＋n－〉，具體來說｜00〉，｜10〉，｜01〉，｜11〉相應本征值E1＝0，E2＝ε－T，E3＝ε＋T，E4＝2ε.定義這個表象下的密度矩陣

并把它表示成為一個16維向量

寫出這個表象下的產生湮滅算符與熱浴算符，代入式（34），可以得到形如式（28）的線性方程.求解這個方程可發現，定態P0也正好是熱平衡

其中，Z＝1＋e－β（ε－T－μ）＋e－β（ε＋T－μ）＋e－β（2ε－2μ）

對于Lindblad方程式（36），也計算了其穩態并與平衡分布做了對比.兩個分布函數ρA與ρB之間的距離定義為

可以看到，只有在T?ε時Lindblad方程的解比較接近Redfield方程的解，而后者就是熱平衡分布.隨著t的增大，其距離也逐漸增大.基于這個結果，作者認為，在非平衡態的研究中，應該優先考慮Refield方程.兩者在平衡態下的對比詳見圖4，其中，Redfield方程給出了正確的平衡態.

圖4 Redfield方程與局域Lindblad方程在平衡態下的對比Fig.4 Stationary states from the local－operator Lindblad equation compared against those from the Redfield equation

3.2.3 雙熱浴中的雙格點費米子系統：非零非對角元

同樣的思路可以處理多個熱浴中的子系統，研究其定態.事實上式（16）就可以直接用于多熱浴的計算，只要注意熱浴算符和，用第ν個熱浴的溫度、化學勢等參量.以下計算雙熱浴中的雙格點費米子系統.

如圖5一個雙格點系統耦合到兩個熱浴，Redfield方程為

其中，用到新的熱浴算符

如果兩個熱浴的溫度與化學勢相同，重復上面的計算，可發現相同的熱力學平衡態.

圖5 雙熱浴中的雙格點費米子系統草圖Fig.5 Sketch of a two－site system coupled to two heat baths

那么一個自然的問題就是，如果兩個熱浴有不同的參數，還能不能得到穩態；如果得到了穩態，這個穩態是否就是要尋找的非平衡定態.對此，沒有一般的理論來支持這個猜想.這個動理學方程也不是第一原理性的.在第一原理動力學方程之上，還加入了諸如無限大熱浴、熱浴處于熱平衡、Markovian近似、耦合項的二級微擾等很多近似.因此，不能保證這個動理學方程的處理方法就能給出正確的非平衡定態.然而，包括前面提到的Landauer－Buttikker公式，NEGF等同樣非第一原理性的方法，都不能寫出非平衡分布的一般形式.既然這個動理學方程方法在合理的條件下能夠給出正確的平衡態，有理由推廣，在相應的條件下，其穩態也可以是非平衡定態.真實情況如何就只能依賴于實驗檢驗了.

現在可以讓兩個熱浴具有不同的溫度TL和TR，重復以上的計算.在此，不再詳細寫出計算過程.從計算結果可以看到密度矩陣的非對角元，此時非零，這一點與平衡態是不一樣的.圖6所示Redfield方程的解與平衡態對比.當δT＝0，也就是平衡態，兩者完全一致.對于非平衡態的非對角元非零，對角元也與平衡態的值不同.其中平衡態取與平均溫度相對應的Boltzmann分布.

圖6 Redfield方程解與平衡狀態對比Fig.6 Comparison between the solution of Redfield equation and equilibrium

下面將這樣的一個框架用來計算多個熱浴中的N－格點系統.

3.3 基本框架的計算問題

總結求解Redfield方程的方法：

a.求得aν，a?ν的動力學演化，并進而得到熱浴算符的表達式.這一步通常需要HS的譜展開.對于N－qubit系統，這是一個2N維本征值問題.

b.得到矩陣Γ的表達式.這是一個4N維的矩陣.

c.求解矩陣Γ的零本征值對應的本征向量.這是求解一個4N維矩陣的本征值.

可以看到，直接求解方法要比求解Schro¨dinger方程困難得多.Schro¨dinger方程的求解已經需要用到各種各樣的近似與間接方法，更何況Redfield方程.因此，發展和運用這一基本框架的重要問題，就是研究出高效率的求解上述方程的方法.

在平衡態統計力學中，有量子Monte Carlo方法、以無相互作用系統上的精確解為基礎的微擾論（Feynman圖展開、BBGKY鏈）、密度矩陣重整化、密度泛函理論等方法來處理相互作用系統.那么，對于非平衡Redfeild方程，能否也研究一套方法，先簡略地總結一下他人提出的方法，然后再重點介紹一下自己的工作.

4 相互作用系統非平衡定態的計算技術

4.1 現有方法小結

量子主方程，不管是Redfield方程還是Lindblad方程，都是一階微分方程，記為

可以用直接積分的方法求解.這樣的積分求解方法包含對一般系統都是用的Runge－Kutta方法（RK）、short－iterative－Arnoldi（SIA）傳播子、短時Chebyshev多項式（CP）傳播子、Newtonian多項式傳播子、以及辛積分子（symplectic integrator）方法.對于這些方法，文獻［21］做了比較和綜述，在此僅做簡單介紹.

Euler方法就是直接計算數值積分，例如ρ（t＋δt）＝ρ（t）＋Lρ（t）δt.它有正比于（δt）2的誤差.RK就是用更合理的多次迭代方式，在不增加太多計算量、不用改變步長的前提下，減小誤差，例如將一步的迭代改成兩步、四步或者更多步.這樣的計算方法對于很多的演化方程都是適用的，但是它并沒有利用動理學方程是線性的這一特征.這個方程的計算復雜度主要取決于步長，每一步的計算需要處理2N矩陣的相乘運算.

線性微分方程存在一般解

所以，只要知道算符（有時候稱為超算符，因為算符通常指作用在波函數上的線性映射，而L作用在密度矩陣，也就是通常所說的算符上）L的本征向量與本征值（當說超算符的本征態和本征向量時，把超算符看在作一個大矩陣，例如方程式（28）中的Γ.在這個意義上，把Γ與L等同，P與ρ等同），就可以知道一切

特別地，如果只研究定態，則

其中要求

只有

所以，問題轉化為求4N維矩陣的零本征值問題.考慮歸一化條件

可以用tr（ρ）替換L的任意一行，不妨說第一行，然后方程的右邊也需要作相應的變化，也就是第一個元素為1，得到

這里diag是指把所操作的向量放在對角元上形成一個矩陣.直接求解4N維線性系統的零本征向量或者線性方程的解，盡管數值意義上精確，但仍然是一個困難的問題.下面一類近似解法的基本想法是用更簡單高效的方法來近似代表演化算符eLt.細節請參見文獻［21］以及其中相應的文獻.下面以SIA為例，介紹這些方法的思想.從一個t時刻的狀態ρ（t）出發，構造一個Krylov空間，就是做n次重復計算Lnρ（t），記錄這些矢量組成的空間.在n維這個子空間內，算符L擁有的形式為

其中，矩陣V可以由Lanczos方法在構造以上Krylov空間的過程中得到.下一步，在子空間內對角化算符l.記對角化以后的算符為L，變換的矩陣為S，則

這樣，就可以得到密度矩陣演化過程的近似解.在這個方法中，還是要計算很多次4N維矩陣的代數運算.

另一個經常使用的方法是隨機波函數方法［7］.其基本思想是利用一系列波函數來等價地代表密度矩陣

然后用這一系列波函數的演化來代替密度矩陣的演化.在這些波函數的演化中，往往會出現跳躍項等需要隨機過程來模擬的動力學，而不再是厄米算符所描述的幺正演化.這個方法的算法復雜性相當于求解Schro¨dinger方程再加上隨機跳躍過程，也就是說大概相當于求解2N維的線性系統.

以上這些方法，對一般的量子主方程都是適用的，沒有要求方程具有任何特殊結構.最近，Prosen等提出了針對密度矩陣的tDMRG方法［6］，但是只能用于局域算符Lindblad方程.其中用到了演化算符的局域展開，這一點在Redfield方程上不能實現.這個方法對于無相互作用系統，相當于把2N的線性系統轉化為N2維的問題.對于相互作用系統，這個方法的效率Prosen等還在進一步研究中.

下面主要介紹本文提出的兩個方法：基于Green函數的非平衡BBGKY鏈方法和基于量子相干態表象的量子隨機微分方程方法.

4.2 類BBGKY方法

從處理平衡態的多體系統的技術［22］知道BBGKY鏈方法與Feynman圖微擾展開或者Dyson方程等是等價的［23－24］.對于無相互作用系統，Wick定理表明單粒子Green函數的方程式封閉，而且只要知道單粒子Green函數，多粒子的Green函數就可以計算出來.因此，只需求解單粒子Green函數.記n粒子Green函數為Gn.對于相互作用系統，G1的方程中會出現G2，G2的方程中會出現G3等.對于這樣的一個方程鏈，必須引入截斷才能處理.其中一種截斷方式是集團展開［24－25］.現在就要將這樣的方法應用于非平衡動理學方程.

為了方便和具體，用V－t2模型（緊束縛模型加上近鄰相互作用）來表達一般理論.子系統定義為

耦合到左右兩端的兩個熱浴上，其Redfield方程可以按照式（16）寫出.

同樣，只對定態感興趣ρ（∞），并且用這個定態來計算各種物理量.對于任何定義在子系統的物理量A，從定態Redfield方程（式（16）取左邊等于零）得到A滿足的方程

這是中心方程.其它方程都將從這個方程中推導出來.根據多體理論，任何物理量都可以表示成為產生湮滅算符的函數.如果知道了所有的各階Green函數，就能夠計算任何物理量.下面計算Green函數，從單粒子

和雙粒子開始

就像在平衡態BBGKY鏈與集團展開中一樣，先從無相互作用的子系統開始.這就是令式（60）中V0＝0.發現這時熱浴算符也只有產生湮滅算符的線性項.合起來，發現此時G1的方程是封閉的.還可以證明此時G2的方程與G1的方程等價，也就是說G2完全由G1給定.這正好就是平衡態Wick定理的內容

以它為基礎，就可以來構造微擾展開.對于相互作用系統，方程右側還有額外的項，記為

通常也稱前者為Hartree－Fock項，后者為關聯項.更高階的Green函數也有類似的展開關系，加上額外的關聯項［9，25］.集團展開的基本思想就是在某一級忽略這樣的關聯項.例如取G2＝0就得到了G1的封閉方程.這相當于是一級集團展開.還可以在G3中保留G1和G2而取G3＝0.這就相當于二級集團展開.同時對于相互作用系統，熱浴算符也不僅僅包含產生湮滅算符的一階項.通過微擾計算可以得到熱浴算符的逐級展開形式，例如

其中，Dα；m與Dα；m1m2m3的定義請參見文獻［8－9］.

把集團展開與熱浴算符的展開形式代入G1的方程，可以得到G1的封閉方程.其主要部分是一個線性方程，同時包含非線性的微擾項.這樣的方程可以通過迭代的方式求解.在小系統上，將用這一階近似求得的近似解gC′，（1），不考慮相互作用的零階近似求得的近似解gC′，（0）和數值精確解gEx作對此，分別計算了密度矩陣的距離與粒子流，見下頁圖7［8］.圖7（a）和圖7（b），展示了一階解的距離dC，（1）和零階解的dC，（0）隨著ΔT和V0的變化，即dC，（1）都要比dC，（0）小很多.在圖7（c）和圖7（d）中，展示了JC，（0）和JC，（1）隨著ΔT和V0的變化，并與數值精確解JEx做了比較.從圖7（c）中可見，在相互作用比較小的時候（V0＝0.2），JC，（0）與JC，（1）都離JEx不遠.在圖7（d）中，看到當相互作用比較大時，JC，（1）要遠遠好于JC，（0）.由此可初步得出結論：在相互作用比較小的條件下，兩種結果符合的很好；從零階近似（也就是忽略相互作用）到一階近似確實有很大的提高.用二階近似作了同樣的計算（計算結果還沒有正式發表）.初步的處理可以在http：∥ccr.bnu.edu.cn／wujinshan上看到.其結果比一階近似有了更大的提高.

圖7 N＝4系統上類BBGKY一階近似解gC，（1）和零階解gC，（0）與精確數gEx對比圖Fig.7 The first order approximation solution of BBGKY－like compared with the exact solution for interacting systems at non－equilibrium

4.3 量子隨機微分方程方法

用確定性動理學方程描述系統是統計力學的一種思路，用隨機過程表述系統是統計力學的另一種思路.現在將上面這個確定性動理學方程轉化為隨機微分方程，然后用Langevin方程、Monte Carlo模擬等方法來求解這樣的方程定態.

將產生湮滅算法轉化為導數算符，密度矩陣轉化為密度分布函數的方法是利用相干態表象.在量子光學中，光子的量子主方程經常通過在相干態表象中轉化為隨機微分方程來求解［27］：密度矩陣成為“概率密度分布函數”（形式上，實際有時取負值），產生湮滅算符成為導數算符.在量子光學中所討論的問題經常是動力學的或者是平衡態的.在這個意義上的工作就是將這樣的技術推廣到非平衡態的研究中去.

對于一個N格點系統，通過這樣一個變換得到的隨機微分方程有2N個復隨機變量.下面以N格點Bose－Hubbard［28－29］為例簡單介紹這個方法，得到的方程有時可以解析求解，有時要通過Monte Carlo模擬來求解隨機微分方程對應的Langevin方程［27，30］.

相干態表象有很多種，這里只介紹其中的一種.一個常用的相干態表象是P表象［27］.在P表象下，第i個格點所對應的隨機變量ξi和ξ＊i看成相互共軛的變量，因此，只考慮N個隨機變量.對于很多系統，兩者是要分開考慮的.一個關于各種不同表象的綜述可以在文獻［27，31］中找到.已有的研究表明這種技術可以處理約15×104個原子在106個動量本征態上的BEC動力學［32］.本文僅利用這種方法給出無相互作用系統的非平衡態的解析表達式，對于相互作用系統，數值工作還不很成熟.

為了討論更具體，以雙熱浴中的一維Bose－Hubbard模型［28－29］為例，介紹相干態表象下的量子隨機微分方程方法.左右熱浴的溫度分別為TL＝和子系統、熱浴、耦合項分別定義為

約化密度矩陣的Redfield方程為式（13）

下一步，利用P表象，將a?l和al替換成導數算符，例如

這樣就得到一個廣義Fokker－Planck方程.對于無相互作用系統，這一方程退化成為標準Fokker－Planck方程［30］.在此，只給出這個無相互作用系統的方程.

其中，Z是歸一化常數，矩陣σ是下述方程的解

圖8為這個解析解與類BBGKY數值解的對比Green函數d和粒子流J，參數取N＝2，λ＝0.5，T＝2.0，μ＝－2.0，紅色（藍色）曲線是解析解（數值模擬）與類BBGKY數值解的對比，結果符合得較好.在考慮相互作用的系統中，必須通過數值模擬來求解廣義Fokker－Planck方程.因此，在圖8中，還將無相互作用系統的數值解與類BBGKY數值解做了對比，用來對數值方法做一個檢驗，也得到了符合較好的結果.

對于包含相互作用的中心系統，得到的廣義Fokker－Planck方程在標準Fokker－Planck的基礎上加入了高階導數項.有的學者認為這樣的廣義Fokker－Planck方程仍然能夠轉化為廣義Langevin方程求解［30，33］.而有的學者認為，必須發明新的方法或者想辦法去掉這樣的高階導數項，例如通過使用Gaussian相空間的方法［31］或者把整體系統而不僅僅是中心系統的演化納入考慮［34］.兩者都有可能成功，尤其是后者，并認為可以得到更加可靠的結果，應該被更加廣泛的應用；而且這樣的處理能夠很容易地求解非Markovian過程.目前，正在開展沿著這個思路的工作.

圖8 無相互作用系統上相干態表象下非平衡態的精確解與類BBGKY數值解的對比Fig.8 Analytical solution compared against the BBGKY numerical solution

本文利用相干態表象找到了無相互作用系統的Redfield方程的解析解，對相互作用系統也大致闡述了求解的思路與技術.實際上，這一技術也能夠用來處理Lindblad方程，而且更簡單.還可以證明這一方法也能給出前面類BBGKY方法中的級聯方程.但是，注意到完全發揮這一方法的威力還需要解決含有三階導數項的廣義Fokker－Planck方程的求解問題，這一問題還有待研究.

5 結論與展望

主要介紹處理量子系統非平衡定態的量子主方程方法的一般框架、計算技術，以及最新提出的兩種技術.這兩種技術不僅可以處理Redfield方程，也可以處理Lindblad方程等其它量子主方程.其方法不僅可以處理非平衡態也可以處理平衡態與動力學.在處理平衡態的時候，不需要從Boltzmann分布出發來計算，主方程演化的穩態自然就是Boltzmann分布.另外，相干態表象下量子隨機微分方程的技術能夠處理非常大的系統，有可能用來研究熱浴和子系統構成的整體系統的動力學.這樣，就可以研究統計物理學的第一個基本問題：熱平衡分布的充分必要條件.可以讓系統處于不同的初始狀態，可以將二級微擾、Markovian近似、無限大熱浴這些假設都去掉或者換成別的假設，來探討子系統的演化末狀態.

另外，統計物理學的第二個基本問題：系綜分布與單一宏觀系統的測量，也可以通過量子主方程與穩態分布的框架來討論.前面提到，系統自身回復時間的尺度很大，通常大于宏觀測量的時間尺度，因此長時平均等于系綜平均的說法實際上不合理：系統在宏觀測量下沒時間作長時平均.實際上與宏觀測量時間對比的不是系統自身的回復時間，而應該是系統在環境耦合下的回復時間.這樣做的物理圖景相當于宏觀測量擾動了系統，使其偏離了平衡態（或者說得到了平衡態分布的一個抽樣），接著環境迫使系統回復到平衡態，然后下一次宏觀測量又擾動了系統，如此反復，直到宏觀測量結束.所以，宏觀測量在環境的幫助下，相當于給定分布函數的抽樣過程，自然其平均等于概率平均，實際模型上的計算過程還在進行.

總之，這個一般框架加上計算方法不僅可以提供一個用來計算非平衡態的一般道路（統計物理學的第三個基本問題），還可以解決第二個基本問題，還有助于解決第一個基本問題.

本文首先從統計物理學本身對非平衡定態研究做了總結，現在再從系統科學的角度來看這個方向的意義.科學可以從研究對象來劃分，也可以從研究深入的層次來劃分為基礎科學與應用科學，從研究內容上還可以分為研究處理研究對象方法的還是研究具體系統性質與規律的.整個物理學基本上是研究物理系統的學問，研究具體系統的性質和規律，同時也研究處理物理系統的方法.具體到統計物理學，它更多地研究處理相互作用的多體系統的處理方法，當然，也研究典型物理系統、典型物理模型.系統科學研究更加一般的（物理的或者非物理的）系統的處理方法.系統科學沒有特定的研究對象，它關注的具體學科領域有很多，包括經濟、工程、社會、物理、地理、生物、心理等各個領域.也許從研究對象來看它是一門交叉科學，但是系統科學是方法論性的學科，不是以研究對象來安身立命的學科.從這個意義上說，它不是交叉科學，僅僅是研究對象交叉起來的科學，其研究方法（至少希望）是一致的統一的，成系統的，有內在邏輯關系的.它從具體的學科來，到一般的方法中去，然后再從一般的方法到具體的學科中去.也許前后這兩個具體的系統屬于不同的學科或者同一個學科.系統科學注重事物（來源于相同或者不同學科的具體系統）之間的聯系，發現其共同與不同的地方，借鑒和發展某一具體系統的研究與處理方法，應用于解決另一具體系統的問題.在這里，關鍵詞就是：聯系、共同點、方法、解決問題.尤其是聯系，這意味著從事物之間，或者事物內的部分之間的聯系，也就是從相互作用出發，來考慮問題.這是系統科學非常重要的考慮問題的角度，有時候稱它為“系統的思想”.有的人把它稱為整體與部分，用整體的觀點來考慮問題，這也有一定道理.但是，也應盡量避免這個提法，因為這看起來好像部分的研究不重要了，整體才是重要的，更好的理解方式應該是說，部分的研究是不能完整地描述整體的，還需要從整體來考慮，有的時候相互作用系統的整體行為與系統的個體的直接行為之間是有差別的，這也被稱為涌現性問題.由于系統科學主要關注相互作用（有時候甚至是多種相互作用的混合）的多體系統，并考察多個層次上的涌現性問題，在這個意義上系統科學也常常被稱為復雜性研究.

對于系統科學這樣一個研究方法，解決多個具體系統的問題的一門學問，它的主要方法又從哪里來，物理學是研究所有自然現象最基礎的學科，它既研究具體系統，也研究方法.物理學探討對自然現象的描述，進而研究現象出現和演化的推動力來理解現象，或者有的時候還要預測和控制現象.物理學的精神就是：尋找和發現客觀現實，探尋其原由，以至（希望）所有自然現象都能夠用最簡單的統一的方式來描述和理解.物理學的發展擴大了其研究對象，一方面，物理學的方法，尤其是統計物理學的方法，還被應用于社會現象的研究，而不僅局限于自然現象.當然，這還沒有成為物理學的主流.物理學這種很強的侵略性使得一部分研究非物理系統的學者對其愛恨交加，也使得最傳統的物理學研究者恨愛交加.另一方面物理學對復雜系統的研究使得某些研究問題和某些研究方法，開始從最傳統的物理學稍微地分離開來.例如，統計物理學和熱力學展示物理系統通常具有趨向熱平衡的特性，可是結構和花樣的出現也是一個普遍事實，例如Poiseuille－Bernard流中花樣的形成、Belousov－Zhabotinsky反應中花樣的出現，甚至生物的進化.在一定條件下，系統可以從各項同性的沒有結構的狀態演化成有結構的狀態.這些問題的研究有其自身的共性，雖然還可以看作是物理學、化學或者生物學的研究，但是把它們聯合起來從原來的學科中獨立出來還是有意義的.這也促使系統科學的研究從傳統物理學中區別開來.

這樣的現象盡管主要出現在物理學研究中，在其它學科，例如生物學、生態學、經濟學、社會學中也有類似風格的工作：發展和利用其它學科或者本學科的方法，尋找不同具體系統之間的共性和差異，然后解決本學科或者其它學科的問題.所以，系統科學方法的一個重要來源就是物理學，尤其是統計物理學.當然，廣義上說，數學本身就是這樣一門學問：從實際問題中（主動地或者被動地）提煉出抽象結構，研究和發展這樣的結構，然后用于實際問題的研究.從這個意義上，系統科學也可以看成是數學的一個分支學科.不過，能看成數學分支的學科太多了，例如理論物理學、理論生物學等，這些都是廣義的數學模型方面的學科.因此，把系統科學獨立出來看，還是有一定意義的.

系統科學就是在這樣的一個背景下得到了比較大的發展：類似風格的研究者和研究工作聚集在一起，用“系統的思想”來研究一個或多個具體系統.基于這樣的特征，人們給這些人和工作取了個名字：系統科學.而且發現這樣一個處理相互作用的多體系統的思想與方法與統計物理學是緊密結合的.因此，希望這個關于統計物理學基本問題研究的綜述不僅僅對物理學工作者有意義，還能夠幫助系統科學的研究者建立一個粗略的圖景.大量來自于平衡態統計物理學的模型與方法，例如Ising模型、相變與臨界現象理論、動力學系統的行為與底層幾何結構的關系、有限大小標度、系綜理論、Metropolis模擬方法及Green函數方法等，已經進入了系統科學的研究.可以預見，在不久的將來，非平衡態的模型與方法，包括處理非平衡系統的一般框架、處理相互作用的Green函數方法與集團展開方法、量子隨機方程的技術等等，也將成為系統理論基礎的一部分.

［1］ Goldstein S，Lebowitz J L，Tumulka R，et al.Canonical typicality［J］.Physical Review Letters，2006，96（5）：050403.

［2］ Dubey S，Silvestri L，Finn J，et al.The approach to typicality in many－body quantum systems［J］.Phys Rev E，2011，85（1）：011141.

［3］ Singh N.How and why does statistical mechanics work［EB／OL］.［2011－03－21］.http：∥arxiv.org／abs／1103.4003.

［4］ Linden N，Popescu S，Short A J，et al.Quantum mechanical evolution towards thermal equilibrium［J］.Phys Rev E，2009，79（6）：061103.

［5］ Saito K，Takesue S，Miyashita S.Energy transport in the integrable system in contact with various types of phonon reservoirs［J］.Physical Review E，2000，61（3）：2397－2409.

［6］ Prosen T，Znidaric M.Matrix product simulations of non－equilibrium steady states of quantum spin chains［J］.Journal of Statistical Mechanics：Theory and Experiment，2009，2009：p02035.

［7］ Wichterich H，Henrich M J，Breuer H P，et al.Modeling heat transport through completely positive maps［J］.Physical Review E，2007，76（3）：031115.

［8］ Wu J.Non－equilibrium stationary states from the equation of motion of open systems［J］.New Journal of Physics，2010，12（8）：083042.

［9］ Wu J.Quantum transport through open systems［D］.Vancouver：The University of British Columbia，2011：79－96.

［10］ Lepri S，Livi R，Politi A.Thermal conduction in classical low－dimensional lattices［J］.Physics Reports－Review Section of Physics Letters，2003，377（1）：1－80.

［11］ Landauer R.Electrical resistance of disordered onedimensional lattices［J］.Philosophical Magazine，1970，21（172）：863－867.

［12］ Cresti A，Farchioni R，Grosso G，et al.Keldysh－green function formalism for current profiles in mesoscopic systems［J］.Physical Review B，2003，68（7）：075306.

［13］ Kubo R.Statistical－mechanical theory of irreversible processes 1：general theory and simple applications to magnetic and conduction problems［J］.Journal of the Physical Society of Japan，1957，12（6）：570－586.

［14］ Wu J，Berciu M.Kubo formula for open finite－size systems［J］.Europhysics Letters，2010，92（3）：30003.

［15］ Toher C，Filippetti A，Sanvito S，et al.Self－interaction errors in density－functional calculations of electronic transport［J］.Phys Rev lett，2005，95（14）：146402.

［16］ Evers F，Weigend F，Koentopp M.Conductance of molecular wires and transport calculations based on density－functional theory［J］.Phys Rev B，2004，69（23）：235411.

［17］ Kubo R，Toda M，Hashitsume N.Statistical physics II：nonequilib－rium statistical mechanics［M］.Berlin：Springer，1991.

［18］ Blum K.Density matrix theory and applications：physics of atoms and molecules［M］.New York：Springer，1996.

［19］ Redfield A G.On the theory of relaxation processes［J］.IBM Journal of Research and Development，1957，1（1）：19－31.

［20］ Redfield A G.The theory of relaxation processes［J］.Adv Magn Reson，1965，1：1－32.

［21］ Kondov I，Kleinekathofer U，Schreiber M.Efficiency of different numerical methods for solving redfield equations［J］.Journal of Chemical Physics，2001，114（4）：1497－1504.

［22］ Mahan G D.Many－particle physics［M］.New York：Plenum Press，1990.

［23］ Petrina D Y.Mathematical foundations of quantum statistical mechanics［M］.Dutch：Kluwer Acad Publ，1995.

［24］ Kadano L P，Martin P C.Theory of many－particle systems II：superconductivity［J］.Physical Review，1961，124（3）：670－697.

［25］ Kira M，Koch S W.Many－body correlations and excitonic effects in semiconductor spectroscopy［J］.Progress in Quantum Electronics，2006，30（5）：155－296.

［26］ Martin P C，Schwinger J.Theory of many－particle systemsi［J］.Physical Review，1959，115（6）：1342－1373.

［27］ Gardiner C W，Zoller P.Quantum noise［M］.Berlin：Springer，2000.

［28］ Jaksch D，Bruder C，Cirac J I，et al.Cold bosonic atoms in optical lattices［J］.Physical Review Letters，1998，81（15）：3108－3111.

［29］ Greiner M，Mandel O，Esslinger T，et al.Quantum phase transition from a superfluid to a mott insulator in agas of ultracold atoms［J］.Nature，2002，415（6867）：39－44.

［30］ Risken H.The Fokker－planck equation：methods of solutions and appli－cations［M］.Berlin：Springer，1989.

［31］ Drummond P D，Deuar P，Corney J F.Quantum many－body simula tions using gaussian phase－space representations［J］.Optics and Spectroscopy，2007，103（1）：7－16.

［32］ Deuar P，Drummond PD.Correlations in a bec collision：first－principles quantum dynamics with 150 000atoms［J］.Physical Review Letters，2007，98（12）：120402.

［33］ Plimak L I，Olsen M K，Fleischhauer M，et al.Beyond the fokker－planck equation：stochastic simulation of complete wigner representation for the optical parametric oscillator［J］.Europhysics Letters，2001，56（3）：372－378.

［34］ Palmieri B，Nagata Y，Mukamel S.Phase－space algorithm for simulating quantum nonlinear response functions of bosons using stochastic classical trajectories［J］.Physical Review E，2010，82（4）：046706.