球諧波展開(kāi)法求解玻爾茲曼傳輸方程綜述

郭 爽，李建清

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

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

1 玻爾茲曼方程(BTE)

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

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

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

2 項(xiàng)匹配法

將分布函數(shù)在實(shí)空間中每一點(diǎn)展開(kāi)為一系列球諧波加權(quán)求和:

考慮一維實(shí)空間 ，以極坐標(biāo)表示的二維動(dòng)量空間(k，θ)，由于分布函數(shù)關(guān)于電場(chǎng)的對(duì)稱(chēng)性，球諧波函數(shù)退化為勒讓德多項(xiàng)式，與φ 無(wú)關(guān)，即:

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

文獻(xiàn)［11，24］對(duì)分布函數(shù)展開(kāi)到一階得到的以上兩個(gè)偏導(dǎo)方程式進(jìn)行了離散。圖1 為文獻(xiàn)［20］中采用項(xiàng)匹配法求解BTE，模擬硅中電子傳輸?shù)玫降膬蓚€(gè)結(jié)果圖，左右兩個(gè)圖分別是在電場(chǎng)為70 kV/cm 和500 kV/cm 時(shí)，將分布函數(shù)展開(kāi)到一階球諧波得到的分布函數(shù)隨能量的變化曲線(xiàn)，并與MC 法對(duì)比的結(jié)果。從圖中可以看出，包括二階諧波展開(kāi)系數(shù)能夠得到與MC 法非常接近的結(jié)果，并且在較高的能量下，MC 法引入了數(shù)值噪聲(如圖1(a))，而采用SHE 法能得到穩(wěn)定的結(jié)果。

項(xiàng)匹配法存在的問(wèn)題是，偏導(dǎo)方程組由展開(kāi)階數(shù)決定，展開(kāi)階數(shù)必須提前確定，而且如果要求高階解，則必須重新對(duì)方程進(jìn)行展開(kāi)和離散，過(guò)程繁瑣，很難擴(kuò)展到任意階展 另一方面，高階展開(kāi)對(duì)器件傳輸特性的模擬是非常重要的，如果忽略一些高階展開(kāi)將會(huì)影響結(jié)果的準(zhǔn)確性。

圖1 MC 方法和項(xiàng)匹配法得到的分布函數(shù)隨能量的變化

3 Galerkin 方法

Galerkin 方法的思路是將分布函數(shù)展開(kāi)為球諧波函數(shù)之和，并將球諧波展開(kāi)插入到BTE 中，然后對(duì)BTE 的每一項(xiàng)同乘以球諧波函數(shù)的共軛，并在整個(gè)動(dòng)量空間進(jìn)行單位球積分，最后采用有限差分法在能量和實(shí)空間對(duì)展開(kāi)系數(shù)進(jìn)行離散，得到以諧波展開(kāi)系數(shù)為未知量的代數(shù)矩陣方程［8，21］:

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

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

圖2 單邊和迎風(fēng)離散得到的f0 對(duì)比圖

圖3 一階和三階展開(kāi)得到的分布函數(shù)

Galerkin 方法的優(yōu)點(diǎn)在于更高階展開(kāi)只增加代數(shù)方程組產(chǎn)生的系數(shù)矩陣的大小，而對(duì)于項(xiàng)匹配法，高階展開(kāi)不僅會(huì)擴(kuò)大耦合偏導(dǎo)方程組而且需要對(duì)其進(jìn)行重新離散。Galerkin 方法允許在不同的空間區(qū)域采用不同的展開(kāi)階數(shù)，而且將球諧波展開(kāi)階數(shù)作為程序的一個(gè)輸入?yún)?shù)，便于檢驗(yàn)高階展開(kāi)對(duì)于器件模擬結(jié)果的影響［8］。

4 投影法

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

為了在等能面上展開(kāi)分布函數(shù)，需要在能量空間與波矢空間建立一種映射關(guān)系，根據(jù)這種轉(zhuǎn)換關(guān)系，利用微觀量投影表達(dá)式［1，22］，將BTE 在等能面上投影為關(guān)于展開(kāi)系數(shù)與一般狀態(tài)密度乘積的球諧波方程，詳細(xì)的推導(dǎo)過(guò)程可參見(jiàn)文獻(xiàn)［2］，得到球諧波方程:

為了得到穩(wěn)定的SHE 解，文獻(xiàn)［1－2］采用MEDS，引入熵函數(shù)

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

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

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

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

圖4 展開(kāi)到一階和九階時(shí)的電流和相對(duì)誤差

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

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

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

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

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

5 總結(jié)

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

［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］ 王祖軍，劉書(shū)煥，唐本奇，等.SiGe HBT 器件特征參數(shù)的數(shù)值模擬與分析［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］ 崔強(qiáng)，韓雁，董樹(shù)榮，等.Mos 器件二次擊穿行為的電路級(jí)宏模塊建模［J］.傳感技術(shù)學(xué)報(bào)，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.