球諧波展開法求解玻爾茲曼傳輸方程綜述

郭 爽，李建清

(電子科技大學物理電子學院，成都610054)

在半導體產業發展初期，采用宏觀的漂移擴散模型或流體力學模型能夠很好地模擬器件的性能。最近幾十年，由于半導體器件特征尺寸的不斷縮小，漂移擴散模型或流體力學模型已經不能夠準確地描述其傳輸特性。如果忽略量子力學效應，半經典的玻爾茲曼方程(BTE)是半導體器件中載流子傳輸特性及其分布最好的數學描述，它能夠比較準確地描述電子的微觀行為［1－4］。常用的直接求解BTE 的方法主要有蒙特卡洛(MC)法和球諧波展開(SHE)法。MC 法從微觀的角度把半導體器件中的載流子看成分立的個體，主要研究載流子的個體行為，通過模擬一個隨機過程來求解BTE。MC 法非常靈活，而且可分析復雜的能帶結構和散射過程，不需要做很多近似，因此廣泛應用于器件模擬［7］。但由于MC 法的隨機特性，該方法具有很多缺陷，比如，計算量大，很難模擬小事件和小信號分析，小電流下需要很大的CPU 時間等。另一種直接求解BTE 的方法是SHE 法，SHE 法對BTE 進行球諧波展開，離散后獲得代數矩陣方程，從而求得BTE 的確定解。SHE 法能夠提供整個器件的分布函數，不依靠遷移率模型，不易受統計噪聲的影響，避免了MC 法的隨機性缺點，并且能夠保存載流子能量分布的信息，而且計算量比較小，尤其在小信號分析和噪聲分析中，SHE 法比MC 法更有優勢［1－6］。長時間以來，由于SHE 法需要很大的內存，導致其在器件模擬方面受到限制。最近幾十年，隨著計算機內存和速度的指數級增長，該方法的可行性也大大增加。早在二十世紀九十年代，一階展開的SHE 法就應用于一維器件的模擬，后來逐漸擴展到任意階展開［8－10］和二維實空間的模擬［11－14］，并且在物理模型上也有了很大的改善，包括一些全帶效應的低階展開 ，全帶結構的價帶模擬［1，16］，考慮磁場作用的體硅模擬［15］，相關散射機制的MOSFET 模擬［17，19］，采用低階展開的二維器件模擬［12，17］以及高階展開的二維模擬［13］。為了克服高階展開的巨大內存需要，尤其是對于二維、三維器件，K. Rupp 等提出了壓縮矩陣存儲定理［3］。通過結合有限體積法數值離散方法以及迎風差分、H－轉換以及最大熵消散(MEDS)等數值穩定技術［1，3－4，8，23］，大大地改善了高階SHE 解的穩定性。

1 玻爾茲曼方程(BTE)

BTE 用于研究半導體中非平衡狀態下電子的分布函數，可得到不同條件下電子的分布函數f(→，→，t)，進而能夠求出電子的各種輸運參量。BTE 形式如下:

其中Q 為散射算符，由各種散射機制的散射率s確定:

利用SHE 法求解BTE 主要有3 種方法:項匹配法、Galerkin 方法和投影法。下面分別介紹這3 種方法。

2 項匹配法

將分布函數在實空間中每一點展開為一系列球諧波加權求和:

考慮一維實空間 ，以極坐標表示的二維動量空間(k，θ)，由于分布函數關于電場的對稱性，球諧波函數退化為勒讓德多項式，與φ 無關，即:

將分布函數插入到BTE 中，通過匹配相同諧波階數的展開項，得到一系列以展開系數為未知量的耦合偏導方程［11－12］，最后離散偏導數方程組得到關于展開系數fn(x，→)的代數方程組。如果將分布函數展開到一階，即n 分別取0 和1 則得到［8］:

文獻［11，24］對分布函數展開到一階得到的以上兩個偏導方程式進行了離散。圖1 為文獻［20］中采用項匹配法求解BTE，模擬硅中電子傳輸得到的兩個結果圖，左右兩個圖分別是在電場為70 kV/cm 和500 kV/cm 時，將分布函數展開到一階球諧波得到的分布函數隨能量的變化曲線，并與MC 法對比的結果。從圖中可以看出，包括二階諧波展開系數能夠得到與MC 法非常接近的結果，并且在較高的能量下，MC 法引入了數值噪聲(如圖1(a))，而采用SHE 法能得到穩定的結果。

項匹配法存在的問題是，偏導方程組由展開階數決定，展開階數必須提前確定，而且如果要求高階解，則必須重新對方程進行展開和離散，過程繁瑣，很難擴展到任意階展 另一方面，高階展開對器件傳輸特性的模擬是非常重要的，如果忽略一些高階展開將會影響結果的準確性。

圖1 MC 方法和項匹配法得到的分布函數隨能量的變化

3 Galerkin 方法

Galerkin 方法的思路是將分布函數展開為球諧波函數之和，并將球諧波展開插入到BTE 中，然后對BTE 的每一項同乘以球諧波函數的共軛，并在整個動量空間進行單位球積分，最后采用有限差分法在能量和實空間對展開系數進行離散，得到以諧波展開系數為未知量的代數矩陣方程［8，21］:

以上矩陣方程中，f 為能量－空間網格中任意給定一點的未知球諧波系數矢量。矩陣方程的右端取決于散射機制和邊界條件。對分布函數展開系數對能量和實空間方向的導數采用不同的差分格式離散，在不考慮數值穩定性的情況下，僅僅影響非零矩陣塊在矩陣中的位置。

Khalid Rahmat 等作者在文獻［21］中利用Galerkin 方法模擬一維N+NN+管，采用單邊兩點近似，對于偶數諧波展開系數對能量和空間方向的求導采用前向差分，而對于奇數諧波展開系數對能量和空間方向的求導采用后向差分。盡管該差分定理在大多數情況下是穩定的，但是在N+N 結處出現了不穩定。針對單邊差分格式存在的不穩定性，Khalid Rahmat 等作者在文獻［8］中提出了一種解決方法，即采用迎風格式的差分來離散展開系數對能量的導數。當電場為正時，奇數諧波系數對能量的求導采用前向差分，偶數對能量的求導采用后向差分;當電場為負時，奇數諧波系數對能量的求導采用后向差分，偶數對能量的求導采用前向差分。圖2為分別采用迎風離散和單邊離散得到的關于能量的一階諧波系數對比圖。圖3 為在N+NN+管的實空間方向兩點處，分別展開到一階和三階得到的分布函數。由圖2 和圖3 可以看出，相比低階展開，高階展開能夠得到更加圓滑的分布函數，同時，迎風格式的差分離散大大提高了SHE 解的穩定性。

圖2 單邊和迎風離散得到的f0 對比圖

圖3 一階和三階展開得到的分布函數

Galerkin 方法的優點在于更高階展開只增加代數方程組產生的系數矩陣的大小，而對于項匹配法，高階展開不僅會擴大耦合偏導方程組而且需要對其進行重新離散。Galerkin 方法允許在不同的空間區域采用不同的展開階數，而且將球諧波展開階數作為程序的一個輸入參數，便于檢驗高階展開對于器件模擬結果的影響［8］。

4 投影法

投影法的思路是利用球諧波函數將BTE 投影為以分布函數展開系數與一般狀態密度乘積為待求解的平衡方程:即將BTE 轉化為球諧波方程，并在能量空間和實空間中利用有限體積方法對球諧波方程離散。投影法將分布函數展開為關于能量的函數，其優點為由于熱平衡時的分布函數在等能面上是各向同性的，在許多散射模型中，散射率是能量的函數而且散射期間能量轉換也是常數，從而散射積分在等能面上可以很容易的計算出來。因此分布函數在等能面上的展開比起關于波矢的模的展開有許多優點［1］。

為了在等能面上展開分布函數，需要在能量空間與波矢空間建立一種映射關系，根據這種轉換關系，利用微觀量投影表達式［1，22］，將BTE 在等能面上投影為關于展開系數與一般狀態密度乘積的球諧波方程，詳細的推導過程可參見文獻［2］，得到球諧波方程:

為了得到穩定的SHE 解，文獻［1－2］采用MEDS，引入熵函數

在球諧波方程(10)兩邊同時乘以熵函數(11)，得到:

從而將球諧波方程(10)劃分為奇偶兩部分。對穩定后的平衡方程(12)采用有限體積法［2］在能量和實空間中進行離散，保證BTE 的粒子數守恒特性。

由于數值穩定性問題主要與空間偏導與其他變量導數之間的耦合有關，而平衡方程中既有對空間變量的偏導，又有對能量的偏導，從而會導致數值不穩定。文獻［4］采用H－轉換，即將其中，H=ε－ψ(r)表示電子總能量，ψ(r)為電子靜電勢。該方法相對于MEDS 的優勢在于即使在彈道限制下也能準確的處理自由流算符，但是存在的不足是引入了與電勢有關的能量網格。

C.Jungemann 等作者在文獻［9］中采用投影法以及MEDS 和H－轉換模擬了一維N+NN+。圖4 為將球諧波展開到一階和九階時得到的電流隨外加偏壓的變化以及相對誤差。由圖可以看出，電流作為外加偏壓的函數，在小偏壓下一階展開就已經出現比較大的相對誤差，主要原因是N+NN+內部存在的內建電場，導致了準彈道傳輸 因此， 展開的最大階數對于器件傳輸特性具有較大的影響。

圖4 展開到一階和九階時的電流和相對誤差

目前存在5 種能帶模型可以用于SHE 法，分別是:Modena 模型［27］、Vecchi 模型［6］、各向異性模型［25］、擴展的Vecchi［26］模型以及全帶模型。S.－M.Hong 等作者在文獻［22］中采用投影法模擬N+NN+管，圖5 ～圖7 為采用Modena 模型、擴展的Vecchi模型以及全帶模型得到的電流、電子速度和電子能量的分布圖，并與全帶結構的蒙特卡洛(FBMC)法比較結果。從圖中可以看出，擴展的Vecchi 模型和FBMC 法取得了基本一致的模擬結果，但是仍然存在略微的差異，主要原因就是對能帶結構做出的近似。文獻［25］中采用各項異性能帶結構的投影法模擬N+NN+管，該能帶結構模型能夠準確的模擬碰撞電離等高能效應，并且可以得到與FBMC 法基本吻合的模擬結果。

圖5 不同能帶模型以及MC 方法得到的終端電流

圖6 不同能帶模型以及MC 方法得到的電子速度圖

圖7 不同能帶模型以及MC 方法得到的電子能量

由此可以看出，各向異性模型、擴展的Vecchi 模型以及全帶模型都能取得與FBMC 法相一致的模擬結果，全帶模型得到的結果更加接近FBMC 法。Modena 模型不能準確的描述高場傳輸和碰撞電離。Vecchi 模型由于受限于最低階的展開不適合于短溝道器件終端電流的計算。各向同性的擴展的Vecchi模型在計算熱載流子特性和終端特性方面比Modena模型準確，并且不需要額外的內存和CPU 時間。從數值效率上比較，各向異性模型存在的不足是增加了耦合項的數目，從而導致需要的內存增大。全帶模型在取得與最簡單的Modena 模型相同的CPU 效率的同時，能夠得到比其他能帶模型更準確的模擬結果。因此，從物理準確度和數值效率兩方面考慮，全帶模型是最好的選擇，擴展的Vecchi 模型次之［22］。

5 總結

用于直接求解BTE 的項匹配法、Galerkin 法以及投影法的共同點是將分布函數展開為一系列球諧波函數之和。不同之處在于前兩種方法是將分布函數展開為關于波矢模的展開系數與球諧波函數乘積之和，而投影法是展開為關于能量的展開系數。而且項匹配法和Galerkin 法得到的是以球諧波展開系數為待求解的代數矩陣方程，但是投影法是將最初的連續方程轉化為網格點上以球諧波展開系數與一般狀態密度之積為待求解的離散的代數矩陣方程。目前，投影法作為求解BTE 的一種SHE 方法得到了廣泛的應用，主要是由于高階展開能夠準確的描述器件的彈道效應和高場傳輸特性，并且該方法能夠結合MEDS 和H－轉換穩定技術，從而大大地改善SHE 解的數值穩定性。

［1］ Jungemann C，Pham A T，Meinerzhagen B.Stable Discretization of the Boltzmann Equation Based on Spherical Harmonics，Box Integration，and a Maximum Entropy Dissipation Principle［J］.Journal of Applied Physics，2006，100(2):024 502－024502－13.

［2］ Rupp K. Numerical Solution of the Boltzmann Transport Equation using Spherical Harmonics Expansions［D］.Dissertation of Master，Vienna University of Technology，2011.

［3］ Rupp K，Jüngel A，Grasser T. Matrix Compression for Spherical Harmonics Expansions of the Boltzmann Transport Equation for Semiconductors［J］. Journal of Computational Physics，2010，229(23):8750－8765.

［4］ Hong S M，Jungemann C，Bollh?fer M.A Fully Coupled Scheme for a Boltzmann－Poisson Equation Solver Based on a Spherical Harmonics Expansion［J］.Journal of Computational Electronics，2009，8(3):225－441.

［5］ Jungemann C，Meinerzhagen B.Analysis of the Stochastic Error of Stationary Monte Carlo Device Simulations［J］. IEEE Trans.Electron Devices，2001，48(5):985－992.

［6］ Vecchi M C，Rudan M.Modeling Electron and Hole Transport with Full－Band Structure Effects by Means of the Spherical－Harmonics Expansion of the BTE［J］.IEEE Trans.Electron Devices，1998，45(1):230－238.

［7］ 王祖軍，劉書煥，唐本奇，等.SiGe HBT 器件特征參數的數值模擬與分析［J］.電子器件，2009，32(3):534－537.

［8］ Rahmat K，White J，Antoniadis D A. Simulation of Semiconductor Devices Using a Galerkin/Spherical Harmonic Expansion Approach to Solving the Coupled Poisson－Boltzmann System［J］. IEEE Trans.Computer－Aided Design of Integrated.Circuits and Systems，1996，15(10):1181－1196.

［9］ Jungemann C，Hong S M，Matz G.High－Order Spherical Harmonics Solution of the Boltzmann Equation and Noise Modeling［C］//International Workshop on Computational Electronics (IWCE)，2010:1－6.

［10］ Wu Y J Hennacy K，Goldsman N，Mayergoyz I. Two Dimensional Submicron MOSFET Simulation Using Generalized Expansion Method and Fixed Point Iteration Technique to the Boltzmann Transport Equation［C］//VLSI Technology，Systems，and Applications，Proceedings of Technical Papers，1995:122－126.

［11］ Ventura D，Gnudi A，Baccarani G.Multidimensional Spherical Harmonics Expansion of Boltzmann Equation for Transport in Semiconductors［J］.Applied Mathematics Letters，1992，5(3):85－90.

［12］ Gnudi A，Ventura D，Baccarani G.Two－Dimensional MOSFET Simulation by Means of a Multidimensional Spherical Harmonics Expansion of the Boltzmann Transport Equation［J］.Solid State Electron.，1993，36(4):575－581.

［13］ Hong S M，Jungemann C，Bollh?fer M. A Deterministic Boltzmann Equation Solver for Two－Dimensional Semiconductor Devices［C］//Proc.SISPAD，2008:293－296.

［14］ Vecchi M C，Mohring J，Rudan M. An Efficient Solution Scheme for the Spherical Harmonics Expansion of the Boltzmann Transport Equation［J］.IEEE Trans.on CAD of Integrated Circuits and Systems，1997，16(4):353－361.

［15］ Smirnov S，Jungemann C. A Full Band Deterministic Model for Semiclassical Carrier Transport in Semiconductors［J］. Journal of Applied Physics，2006，99(6):063707－063707－11.

［16］ Pham A C，Meinerzhagen B.A Full－Band Spherical Harmonics Expansion of the Valence Bands up to High Energies［C］//Proceedings of SISPAD，2006:361－364.

［17］ Liang W，Goldsman N，Mayergoyz I.2－D MOSFET Modeling Including Surface Effects and Impact Ionization by Self－Consistent Solution of the Boltzmann，Poisson，and Hole－Continuity Equation［J］.IEEE Trans.Electron Devices，1997，44(2):257－267.

［18］ Reggiani S，Vecchi M C，Rudan M. Investigation on Electron and Hole Transport Properties Using the Full－Band Spherical－Harmonics Expansion Method［J］. IEEE Trans. Electron Devices，1998，45(9):2010－2017.

［19］ 崔強，韓雁，董樹榮，等.Mos 器件二次擊穿行為的電路級宏模塊建模［J］.傳感技術學報，2008，21(2):361－364.

［20］ Ventura D，Gnudi A，Baccarani G.An Efficient Method for Evaluating the Energy Distribution of Electrons in Semiconductors Based on the Spherical Harmonics Expansion［C］//International Workshop on VLSI Process and Device Modeling，Oiso，Japan，1991:27－29.

［21］ Khalid Rahmat，Jacob White，Dimitri A Antoniadis. A Galerkin Method for the Arbitrary Order Expansion in Momentum Space of the Boltzmann Equation Using Spherical Harmonics［C］//Proc.NUPAD V，1994:133－136.

［22］ Hong Sung－Min，Pham Anh－Tuan，Christoph Jungemann.Deterministic Solvers for the Boltzmann Transport Equation ［M］. Springer Press，2011.

［23］ Ringhofer C.Numerical Methods for the Semiconductor Boltzmann Equation Based on Spherical Harmonics Expansions and Entropy Discretization［J］.Transp.Theory Stat.Phys.，2002，31(4):431－452.

［24］ Gnudi A，Ventura D，Baccarani G.One－Dimensional Simulation of a Bipolar Transistor by Means of Spherical Harmonic Expansion of the Boltzmann Transport Equation［C］//Proceedings of SISPED，1991，4:205－213.

［25］ Hong S M，Matz G，Jungemann C A.Deterministic Boltzmann Equation Solver Based on a Higher Order Spherical Harmonics Expansion With Full－Band Effects［J］.IEEE Trans.Electron Devices，2010，57(10):2390－2397.

［26］ Jin S，Hong S M，Jungemann C.An Efficient Approach to Include Full－Band Effects in Deterministic Boltzmann Equation Solver Based on High－Order Spherical Harmonics Expansion［J］. IEEE Trans.Electron Devices，2011，58:1287－1294.

［27］ Brunetti R，Jacoboni C，Nava F，et al.Diffusion Coeffcient of Electrons in Silicon［J］. Journal of Applied Physics，1981，52(11):6713－6722.