郭 爽，李建清

(電子科技大學(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)法和球諧波展開(SHE)法。MC 法從微觀的角度把半導(dǎo)體器件中的載流子看成分立的個(gè)體，主要研究載流子的個(gè)體行為，通過模擬一個(gè)隨機(jī)過程來求解BTE。MC 法非常靈活，而且可分析復(fù)雜的能帶結(jié)構(gòu)和散射過程，不需要做很多近似，因此廣泛應(yīng)用于器件模擬［7］。但由于MC 法的隨機(jī)特性，該方法具有很多缺陷，比如，計(jì)算量大，很難模擬小事件和小信號分析，小電流下需要很大的CPU 時(shí)間等。另一種直接求解BTE 的方法是SHE 法，SHE 法對BTE 進(jìn)行球諧波展開，離散后獲得代數(shù)矩陣方程，從而求得BTE 的確定解。SHE 法能夠提供整個(gè)器件的分布函數(shù)，不依靠遷移率模型，不易受統(tǒng)計(jì)噪聲的影響，避免了MC 法的隨機(jī)性缺點(diǎn)，并且能夠保存載流子能量分布的信息，而且計(jì)算量比較小，尤其在小信號分析和噪聲分析中，SHE 法比MC 法更有優(yōu)勢［1－6］。長時(shí)間以來，由于SHE 法需要很大的內(nèi)存，導(dǎo)致其在器件模擬方面受到限制。最近幾十年，隨著計(jì)算機(jī)內(nèi)存和速度的指數(shù)級增長，該方法的可行性也大大增加。早在二十世紀(jì)九十年代，一階展開的SHE 法就應(yīng)用于一維器件的模擬，后來逐漸擴(kuò)展到任意階展開［8－10］和二維實(shí)空間的模擬［11－14］，并且在物理模型上也有了很大的改善，包括一些全帶效應(yīng)的低階展開 ，全帶結(jié)構(gòu)的價(jià)帶模擬［1，16］，考慮磁場作用的體硅模擬［15］，相關(guān)散射機(jī)制的MOSFET 模擬［17，19］，采用低階展開的二維器件模擬［12，17］以及高階展開的二維模擬［13］。為了克服高階展開的巨大內(nèi)存需要，尤其是對于二維、三維器件，K. Rupp 等提出了壓縮矩陣存儲定理［3］。通過結(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)展開為一系列球諧波加權(quán)求和:

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

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

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

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

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

3 Galerkin 方法

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

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

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

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

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

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

4 投影法

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

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

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

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

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

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

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

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

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

5 總結(jié)

用于直接求解BTE 的項(xiàng)匹配法、Galerkin 法以及投影法的共同點(diǎn)是將分布函數(shù)展開為一系列球諧波函數(shù)之和。不同之處在于前兩種方法是將分布函數(shù)展開為關(guān)于波矢模的展開系數(shù)與球諧波函數(shù)乘積之和，而投影法是展開為關(guān)于能量的展開系數(shù)。而且項(xiàng)匹配法和Galerkin 法得到的是以球諧波展開系數(shù)為待求解的代數(shù)矩陣方程，但是投影法是將最初的連續(xù)方程轉(zhuǎn)化為網(wǎng)格點(diǎn)上以球諧波展開系數(shù)與一般狀態(tài)密度之積為待求解的離散的代數(shù)矩陣方程。目前，投影法作為求解BTE 的一種SHE 方法得到了廣泛的應(yīng)用，主要是由于高階展開能夠準(zhǔn)確的描述器件的彈道效應(yī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］ 王祖軍，劉書煥，唐本奇，等.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)，韓雁，董樹榮，等.Mos 器件二次擊穿行為的電路級宏模塊建模［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.