射频(RF)容性耦合等离子体(capacitively coupled plasmas, CCPs)作为一种典型的低温等离子体源, 因其结构简单被广泛应用于微电子领域^[1–5], 例如薄膜沉积^[6,7]、刻蚀^[8–11]等. 在这些应用中, 电子功率吸收模式在很大程度上决定了外部电源(射频)产生的电磁场能量如何耗散到等离子体中, 从而影响电子驱动的碰撞过程与反应性自由基的生成速率, 最终决定材料表面的加工质量.

在不同外部控制参数下, 等离子体中电子获取能量的动理学机制可发生显著的变化, 并与放电运行的多种模式密切相关. 在典型的低气压电正性放电(如氩气和氦气放电)中, 实验与模拟均观察到两种电子功率吸收模式, 具体呈现取决于放电条件. 在α模式下, 电子在每个射频周期内鞘层扩张期间被电场加速, 因此电离率在扩张中鞘层边缘附近达到峰值^[12,13]. 而在γ模式下, 电离主要依靠电极受离子轰击后发射的二次电子来维持. 这些电子在鞘层电压较高时被强电场加速, 并在鞘层内部通过碰撞倍增. 因此, 大部分电离发生在完全扩张后的鞘层边缘附近^[14,15].

氧气在工业过程中具有重要应用价值, 例如硅的氧化、聚合物薄膜的表面改性以及半导体制造中的刻蚀工艺. 作为一种典型的电负性气体, 氧气中诸如解离附着等过程会在体区形成大量负离子并显著耗尽电子数密度. 随之而来的体区电导率下降将促使系统建立较强的漂移电场以加速电子, 从而维持电流连续性. 同时, 由于电子密度空间分布发生显著的变化, 其在鞘层边缘处出现局域极大值(即电正性边缘), 鞘层边缘附近亦形成较强的双极电场. 与电正性气体放电中鞘层边缘所形成的双极扩散场不同, 该电场加速电子向极板运动. 此时, 电子功率吸收主要发生在体区的漂移电场和鞘层边缘的双极扩散场中, 因而名为漂移-双极模式(DA-mode)^[16]. 当驱动频率等于或低于高电负性等离子体的本征频率时, 正、负离子可在体区内部响应射频交变电场, 从而形成一种自组织的条纹结构, 即条纹(STR)模式. 自首次报道不同模式存在的研究以来, 已在多种气体中观测到这些截然不同的能量吸收模式之间的转变, 其触发条件包括气压、射频电压及驱动频率的变化. α模式与γ模式之间的转变已有大量的实验^[17]及数值模拟^[18,19]研究; 而DA模式与α模式之间的转变则已被证实可由驱动频率^[20]、气压^[21,22]、电极间距^[21,23]以及外加磁场^[24]等参数调控.

针对不同放电模式下的等离子体参数分布, 如时空分辨的电场结构、电子密度演化以及功率吸收特性已有深入且详尽的研究. 然而, 现有研究总体上主要聚焦在电子功率吸收过程, 这在本质上仅描述了电磁能向电子体系输入能量的瞬时位置与初始机制; 相对而言, 电子在获取能量后如何在空间输运并最终热化耗散的动力学过程仍缺乏系统刻画. 事实上, 这些内部输运和耗散过程很大程度上决定了等离子体的能量分布和相关反应速率, 并影响着工业应用中的工艺表现.

在宏观层面, 上述能量演化过程通常可借助流体模型进行描述. 但传统流体模型往往依赖闭合假设, 并且其方程体系通常难以实现总能量与流体动能、流体热能的严格解耦, 从而不利于从能量守恒链条上清晰辨识定向加速与能量转化这两类物理过程.

基于电子动量平衡方程、并辅以粒子模拟/蒙特卡罗碰撞(particle in cell/Monte Carlo collision, PIC/MCC)数值模拟的电子功率吸收统一而精确的处理方法, 最早由Surendra与Dalvie^[25]提出. 其核心在于使用PIC/MCC方法计算流体方程中的各项, 从而可以在动理学精度上解析总电子功率吸收的构成项. 该方法能够自洽地、在时空上分辨地描述容性耦合射频放电(CCRF)中的电子功率吸收过程, 被称为玻尔兹曼项分析法^[26,27]. 现有研究表明, 在绝大多数工况下, “欧姆加热”和与压力梯度相关的“压力加热”是总电子功率吸收的主要贡献项, 仅在特定条件下, 惯性功率吸收的时间平均有非零贡献.

但从更严格的动理学视角来看, 电子被电场加速获得能量导致其速度分布变得各向异性, 与随后通过碰撞或无碰撞机制重新变得各向同性, 在物理上是两种不同过程, 分别对应电子功率吸收与电子加热^[28]. 但玻尔兹曼项分析基于动量守恒方程(一阶矩方程), 其不仅难以将上述两个过程在理论上解耦, 也无法进一步识别能量注入之后的输运、转化与损失过程.

延续玻尔兹曼项分析的思想, 将矩方程扩充至能量方程(二阶矩方程)并解耦了宏观定向运动和微观无序热运动, 建立自洽完整的动理学矩框架, 从而对电子能量“吸收-转换-输运-耗散”动力学过程进行完整描述^[29]. 对于电负性气体而言, 电负性的增强和电子的剧烈损耗不仅会触发放电从α模式向DA模式的转变, 也会显著地影响后续能量的输运方向、动能向热能转化通道以及能量耗散的空间分布. 本工作系统对比了氧气CCP中α模式和DA模式下电子能量输运的动力学过程. 基于动理学矩框架, 详细阐述了外部电源(射频)产生的电磁场耦合给等离子体系统的能量-流体动能-流体热能的能量转换过程, 为电负性CCP中的电子能量动理学演化提供新的见解和理论研究工具.

输运过程是指系统受扰动后形成的非平衡态向平衡态转化的不可逆过程, 当宏观系统不处于热平衡时会自然地出现守恒量的输运现象. 射频等离子体系统在周期性外部激励下会偏离热力学平衡, 质量、动量、能量等内部守恒量将沿自发形成的梯度发生空间输运, 从而在空间上重新分布粒子和能量. 通过不可逆过程耗散多余的粒子和能量, 系统在每一瞬间皆向稳态弛豫. 在等离子体中, 所有宏观输运现象的根源在于携带守恒量的微观粒子运动. 由于等离子体是一个包含大量相互作用的带电粒子的系统, 直接通过追踪每个粒子的运动轨迹来分析等离子体的宏观行为是不切实际的. 流体模拟将等离子体中的各种粒子视为流体介质并求解由质量、动量、能量守恒方程组成的流体力学方程组以获得等离子体中宏观量的空间分布. 其优点是计算速度快, 但存在大量物理假设且忽略了微观粒子的动理学细节, 无法捕捉粒子的分布函数和能量分布, 难以处理非局域效应和复杂的化学反应, 并且需要假设粒子处于平衡态, 电子能量需服从麦克斯韦分布. 然而在动理学框架下, 等离子体在任意时刻的完整状态可由物种s在相空间中的分布函数$ {f}_{s}(\boldsymbol{r}, \boldsymbol{v}, t) $唯一确定, 据此我们可以推导出物种s的所有宏观物理量. 具体而言, $ {f}_{s}(\boldsymbol{r}, \boldsymbol{v}, t) $表示在t时刻, 相空间坐标($ \boldsymbol{r}, \boldsymbol{v}) $处, 单位体积相空间中s粒子的代表点密度, 因此相空间体积元内的粒子数为

等离子体中的每个粒子都具有一些微观属性$ \chi (\boldsymbol{r}, \boldsymbol{v}, t) $, 一般来说, 它是粒子位置$ \boldsymbol{r} $、速度$ \boldsymbol{v} $和时间t的函数. 将s粒子属性$ \chi (\boldsymbol{r}, \boldsymbol{v}, t) $的平均值定义为

其中$ \langle \cdot {\rangle }_{s} $表示s粒子对速度空间的平均值, 平均值总是与$ \boldsymbol{v} $无关, 只是$ \boldsymbol{r} $和t的函数. 宏观量, 如数密度、通量、动理学压力、热能通量等, 可被视为物理量的平均值, 涉及大量粒子的集体行为. 这些宏观量与分布函数的各阶矩相关. 分布函数$ {f}_{s} $的第k阶矩写为

它是一个k阶张量. 集合$ \left\{{\mathcal{M}}_{k}, \, k=0{, }1, 2, \cdots \right\} $可以看作是分布函数的另一种表达形式, 当分布函数足够光滑时, 这个集合可以唯一确定该分布函数^[30,31]. 低阶矩都有物理含义. 其中f_s的零阶矩表示粒子数密度,

一阶矩表示粒子通量密度,

二阶矩表示动量流, 也被称为应力张量或并矢,

三阶矩表示能流密度:

粒子分布函数对独立变量$ \boldsymbol{r}, \boldsymbol{v} $和t的依赖性受玻尔兹曼方程所支配. 首先考虑其一般形式,

其中, $ \boldsymbol{a} $为外力引起的加速度, 方程右侧$ {\left(\delta f/\delta t\right)}_{\text{coll}} $代表分布函数由于碰撞而产生的变化率. 在(8)式中, 我们省略了表示物种的下标s, 且只考虑了电磁力引起的加速度. 将(8)式每项乘以物理量并在速度空间积分, 得到一般形式的输运方程:

(9)式第一项$ \partial (n\langle \chi \rangle )/\partial t $为输运量在空间某一点的时间变化率, 对处于周期性稳态的容性耦合射频等离子体放电, 该项在一个射频周期内的平均值为零. 第二项$ \nabla \cdot (n\langle \chi \boldsymbol{v}\rangle ) $为输运通量的散度, 描述由于该通量在空间上的不均匀导致的局域密度变化, 通量散度项为正值表示流出该区域, 负值表示流入该区域. 第三项$ -n\langle \boldsymbol{a}\cdot {\nabla }_{v}\chi \rangle $为外部作用力对输运量的影响; 第四项[$ {\delta (n\langle \chi \rangle )/\delta t]}_{\text{coll}} $为碰撞导致的输运量的时间变化率. 分别取输运量为粒子质量、动量和能量, 即$ \chi =m,\; \chi =m\boldsymbol{v} $, $ \chi =m{v}^{2}/2 $, 可以得到相应的宏观流体方程组, 即流体的质量、动量和能量输运方程, 这三个方程也称为玻尔兹曼方程的前三阶速度矩方程:

(10)式中第一项分别为流体质量、动量和总能量密度, 总能量密度又可分为流体动能密度和流体热能密度. 第二项散度项中的$ \boldsymbol{{\varGamma }} $或$ {\boldsymbol{\mathcal{T}}} $代表输运量的通量, 即单位时间、单位面积内通过的输运量的值. 方程右边代表由于碰撞导致的质量、动量和总能量的时间变化率, 即碰撞导致的质量、动量、能量密度的损失率. 从量纲分析的角度来看, 通量通常可以表示为相应物理量与流体速度的乘积, 然而需要注意的是, 通量定义本身同时包含了流体的宏观整体运动与微观热运动. 因此, 通量一般可进一步分解为两部分: 由整体漂移运动引起的对流输运, 以及由微观随机热运动导致的非对流输运, 后者在宏观上表现为应力和热传导.

高阶矩的物理意义在流体共动系中会变得更加清晰. 在这一参考系下, 粒子相对于流体的本动速度可写为$ \boldsymbol{w}=\boldsymbol{v}-\boldsymbol{u} $, 此时动量通量张量(应力张量)退化为压力张量, 而能量通量密度则相应地成为热流矢量:

其中压力张量$ {\boldsymbol{\mathcal{P}}}$代表随机热运动对动量通量的贡献, 热通量$ \boldsymbol{q} $代表由于粒子微观热运动导致的能量输运. 通过伽利略变换, 在这两个不同参考系下(实验室系与共动系)中测量的矩是相关的. 通过直接代入, 很容易验证这一点:

流体的总能量密度可进一步分解为与宏观漂移运动相关的流体动能密度, 以及与微观随机热运动相关的流体热能密度. 具体来说, 流体动能密度$ {\rho }_{\varepsilon , \mathrm{k}}=\dfrac{1}{2}{\rho }_{\mathrm{m}}{u}^{2} $(也常称为机械能密度或流体动压)取决于流体的宏观运动速度; 而流体热能密度$ {\rho }_{\varepsilon , \text{th}}= \dfrac{1}{2}m\displaystyle\int{w}^{2}f(\boldsymbol{x}, \boldsymbol{v}, t)\mathrm{d}\boldsymbol{v}=\dfrac{3}{2}p $(也常称为内能密度或流体静压)反映了流体内分子无规则热运动的强度. 其中, $ p=\dfrac{1}{3}\text{Tr}({\boldsymbol{\mathcal{P}}}) $为标量形式的各向同性压力, 即压力张量的迹.

我们将能流具体表达式$ {\boldsymbol{{\varGamma }}}_{\varepsilon } = \dfrac{1}{2}{\rho }_{m}\langle {v}^{2}\boldsymbol{v}\rangle = {\rho }_{\varepsilon }\boldsymbol{u}+ {\boldsymbol{\mathcal{P}}}\cdot \boldsymbol{u}+\boldsymbol{q} $代入总能量输运方程, 并将总能量密度分解为流体动能与流体热能$ {\rho }_{\varepsilon }={\rho }_{\varepsilon , \mathrm{k}}+{\rho }_{\varepsilon , \text{th}} $; 之后通过动量方程与漂移速度的点积, 结合连续性方程, 构造动能输运方程; 最后利用总能量守恒的性质, 用总能量输运方程减去上述得到的流体动能输运方程, 从而得到流体动能输运方程与热能输运方程如下:

本文采用自主研发的1d3v (一维位形空间、三维速度空间)粒子网格/蒙特卡罗(PIC/MCC)模拟程序ASTRA^[32–35]对平行板容性耦合氧气放电进行数值模拟. 图1示意性地展示了本模拟所描述的放电结构. 射频电极位于x = 0 cm处, 电极间距为L_gap = 25 mm, 施加一个射频电源, 电压波形为V_rf(t) = 200 cos(2πft), 单位为V. 选取两个频率, 分别为10和30 MHz. 气体压力为p_g = 150 mTorr (1 mTorr = 0.133 Pa), 气体温度固定为T_g = 300 K, 时间步长为∆t = T_rf /1600, 空间步长为∆x = d/256. 模拟中追踪的带电粒子包括电子、$ \text{O}_{2}^{+} $和O^–. 模拟中涵盖的碰撞反应集来自于Derzsi等^[36]报道. 其中, 电子与O₂的碰撞包括弹性散射、转动/振动激发、电离、解离及解离附着/激发. 针对氧负离子O^–, 模型考虑了其与O₂的弹性散射, 与电子和O₂碰撞的脱附反应、与$ \text{O}_{2}^{+} $的中和反应, 以及与O₂(a¹∆_g)的脱附反应; 对于$ \text{O}_{2}^{+} $, 考虑了$ \text{O}_{2}^{+} $与O₂的弹性碰撞. 特别需要指出的是, 亚稳态分子O₂(a¹∆_g)在氧气放电中具有重要作用. 鉴于中性粒子输运的时间尺度远大于射频周期, 直接模拟其动态平衡极耗机时. 在我们的模拟中, 假设O₂(a¹∆_g)粒子在腔室内均匀分布, 其密度采用与Derzsi等^[36]相同的方法确定, 即通过平衡方程计算得出. 亚稳态粒子的表面猝灭概率α是模拟中的重要参数, 它直接影响O₂(a¹∆_g)的密度, 进而显著改变带电粒子的分布. 本模拟采用α = 6×10^–3的设定, 该数值与Derzsi等^[36]报道一致, 且经多项实验验证具有合理性. 为了简化计算, 同时为了在物理上解耦复杂的加热机制, 从而在纯净的背景下孤立的研究α和DA模式下的能量输运特性, 本研究在忽略了粒子在极板处的反射及二次电子发射.

在PIC/MCC方法中, 各种场量通常都是在网格点上定义的, 包括电场、粒子密度、电荷密度等. 通过精准统计时空上每一点的分布函数再重新计算各阶矩量在计算上是重复且低效的, 我们通过追踪足够多数量的超粒子来直接统计分布函数的矩. 在每一个时间步, 将粒子的微观信息(位置r、速度v)映射到空间网格上. 例如, 压力张量$ {\boldsymbol{\mathcal{P}}}$是通过统计网格内的所有粒子的随机热运动速度w = v – u计算得到的二阶矩: $ {\boldsymbol{\mathcal{P}}}({\boldsymbol{r}}_{i})=\displaystyle\sum\nolimits_{p}S({\boldsymbol{r}}_{i}, {\boldsymbol{r}}_{p}){m}_{p}{\boldsymbol{w}}_{p}\times {\boldsymbol{w}}_{p} $, 其中S为形状因子. 一旦在网格上获得了宏观物理量(如流体速度场u和压力张量场$ {\boldsymbol{\mathcal{P}}} $), 涉及空间导数的项则可以采用中心差分法在网格上直接计算. 类似地, 输运方程中的碰撞项可以通过蒙特卡罗碰撞(MCC)模块前后, 精准统计参与碰撞粒子的动量和能量变化量, 并映射到网格上积累得到.

通过上述方法, 输运方程中的各项可以在无需引入任何经验性闭合参数的情况下, 实现独立诊断. 各方程在不同放电模式下, 左右两端各项之和在时间平均和空间平均上严格相等, 从而实现了动理学层面上的“自洽诊断”. 在达到周期稳态后, 跟踪的每种粒子超粒子数都约为1×10⁵, 对约5000个射频周期的数据进行平均以获得高精度数据. 模拟中采用了静电近似, 模拟中每个输运项均被独立诊断, 通过比较输运方程等号两边各项之和, 可以检验该方法在不同研究场景中的一致性与完整性.

值得说明的是, 电负性放电的模式受多个参数影响. 我们选取驱动频率作为控制变量以调节放电的电负性程度, 并聚焦于分析α模式和DA模式主导下的电子能量动理学行为.

当容性耦合等离子体(CCPs)达到稳态时, 其等离子体参数会在射频周期内呈现快速周期性变化, 这意味着相邻射频周期中参数的时空演化规律保持一致. 此时, 粒子的各项守恒量均处于动态平衡状态, 若要揭示其时间演化与空间输运的物理本质, 就需要对总能量输运方程及其拆分后的动能输运方程、热能输运方程中的每一项以及相关物理量进行时空分辨分析. 同时, 在射频周期内对各项分别进行时间与空间平均处理, 可清晰地呈现其空间分布特征与时间演化规律, 进而为电子能量平衡的动理学机制提供自洽且完整的理解.

图2展示了放电到达稳态后, 电子数密度、电场以及电子碰撞电离率在一个射频周期内的时空分布, 横坐标为一个射频周期内的时间, 纵坐标表示电极间距归一化的空间位置. 由图2(a1)可见, 在频率为30 MHz电压波形驱动时, 表现出“类电正性”特性^[21]. 射频电压为负时, 电子迅速响应外加电场并被排斥到主体区, 电子在射频周期的不同阶段被拉入或推出鞘层区域, 在体区保持恒定, 其电子数密度时空分布与电正性气体类似. 在电正性气体放电中, 电子密度随与电极距离的增加而单调递增. 从图2(b1)可以看出, 在鞘层塌缩阶段, 鞘层内会出现与其自身电场方向相反的电场, 即双极扩散场. 与鞘层内电场不同, 双极扩散场能够加速电子向极板运动, 从而确保电子通量的守恒. 同时, 从图3(a1)可以看出, 该电场还会在鞘层塌缩阶段引发局部区域内的电子功率吸收现象. 图2(c1)中, 电离峰主要集中在扩张鞘层的边缘, 表明在30 MHz驱动频率下, 放电表现出显著的鞘层边缘功率吸收特征, 表明其维持在α模式.

在频率为10 MHz时, 电子数密度分布呈现截然不同的特征, 这是由于低频率下, 氧气放电具有高电负性所致. 高电负性使得鞘层边缘最大值处形成“电正性边缘”(见图2(a2)), 导致该位置出现电子密度的局部最大值. 由于负离子的存在会导致准中性引起的电子密度耗尽, 因此电子密度将在该“负离子区域”的边缘达到最大值. 相较于电正性气体放电(其中不存在“负离子核心”), “负离子核心”与“电正性边缘区域”将从根本上改变电负性放电中电子功率吸收机制的特性^[27]. 如图2(b2), 等离子体体区电子密度的降低导致导电性下降, 因此在给定电流下会产生较高的电场, 即漂移电场. 同时, 由于正电边缘区域的形成, 在每个电极靠近最大鞘层宽度位置的最大电子密度梯度位置处产生强双极性电场. 如图2(a2)所示, 10 MHz时, 由于电子密度从体区中心向鞘层边缘逐渐降低, 为了维持电流连续性, 漂移电场必然会随密度变化而呈现空间梯度, 因此电场向扩张鞘层逐渐减小, 电子被这些漂移电场和双极性电场加速. 这样在塌缩阶段, 电离率的最大值出现在靠近鞘层边缘的最强双极电场位置, 如图2(c2)所示. 在驱动频率为10 MHz时, 放电在漂移-双极(drift-ambipolar, DA)模式下进行.

为了对电负性的氧气放电进行能量输运分析, 对玻尔兹曼方程做二阶速度矩, 得到标准形式的能量输运方程. 对其进行整理, 将能量源项置于等式左边, 得到如下形式:

其中$ {\rho }_{\varepsilon } $是总能量密度, $ {\boldsymbol{{\varGamma }}}_{\varepsilon } $是能流, $ \boldsymbol{J} $和$ \boldsymbol{E} $分别是带电粒子的电流密度和电场强度. 将能量输运方程划分为以下四项:

1) $ \boldsymbol{J}\cdot \boldsymbol{E} $表征洛伦兹力做功, 对应带电粒子从电磁场吸收的功率密度, 它作为能量输运方程中的源项.

2) $ \partial {\rho }_{\varepsilon }/\partial t $表示能量密度的局域时间变化率, 在周期性稳态过程中该项平均值为零.

3) $ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $代表能流的散度, 反映了能量输运引起的能量密度空间变化. 该项的正值表示能量从该区域净流出, 而负值则表示能量净流入.

4) $ -{\left(\delta {\rho }_{\varepsilon }/\delta t\right)}_{\text{coll}} $反映了由于碰撞导致的能量密度变化率.

各项对能量输运的贡献反映了等离子体内部复杂的能量交换情况. 在稳态CCPs中, 粒子的各守恒量均处于周期性动态平衡. 这意味着净源区(吸收超过耗散)必然对应能量净流出, 而净汇区必然对应能量净流入. 图3显示了电子在一个射频周期内的能量密度以及能量输运方程中各项的时空分布. 图4进一步提供了各能量输运项的时间和空间平均结果. 图4表明, 无论是时间平均还是空间平均, (14)式左侧等于右侧各项之和, 即第一项等于第2至4项之和$ \left(\displaystyle\sum\nolimits_{i=2}^{4}i\right) $, 验证了模拟与诊断的自洽性与准确性.

带电粒子的能量来源于电磁场对电子、离子流体所做的功. 实际上, 受射频功率驱动, CCPs像一个电容伴随着充电和放电过程不断地储存和释放能量. 射频功率产生的能量以电场能流的形式流入或流出放电室, 大部分输入能量通过电场能量形式存储, 并通过电流从电极流出的能量释放过程返回外部电路. 射频功率输出的一小部分能量通过粒子受电场力在电场中的加速和减速过程实现了场与粒子之间的能量转换, 被粒子直接吸收. 我们这里仅关注带电粒子获得能量后的输运过程, 电磁能的输运过程可参考Wu等^[37]的报道.

当放电处于α模式下, 如图3(a0)—(a4), 与电正性气体放电类似, 电子的能量主要在鞘层附近获得. 不同的是, 体区微弱的漂移电场也会使电子在体区出现部分能量吸收. 电子能量吸收主要发生在鞘层扩张阶段, 在鞘层塌缩过程中, 由于鞘层消退位置上出现电场反转, 导致鞘层塌缩阶段出现显著的电子能量吸收现象.

从图4(a1)的时间平均结果来看, $ {\boldsymbol{J}} \cdot {\boldsymbol{E}} $作为源项, 在鞘层边缘出现极大值. 这说明电子的能量吸收主要通过鞘层电场获得, 少量通过微弱的漂移电场获得; 时间变化项$ \partial {\rho }_{\varepsilon }/\partial t $经时间平均后在空间任一点净贡献均为零, 这归因于CCPs达到稳态时, 等离子体密度或温度在周期上都不会改变. 这意味着粒子从电磁场获得的能量将全部通过边界流出或与背景中性粒子碰撞耗散, 分别对应于第三项$ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $和第四项$ -{\left(\delta {\rho }_{\varepsilon }/\delta t\right)}_{\text{coll}} $. 第三项能流散度项$ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $在体区为负值, 在鞘层区域为正值, 且变化趋势与源项相似, 说明其将鞘层区域吸收的能量输运到体区并通过碰撞耗散. 在α模式下, 电子能量获得位于鞘层、能量损失位于体区, 稳态能量方程必然要求能量由鞘层向体区输运, 表现为时间平均能流指向体区. 从图4(b1)空间平均结果来看, 碰撞项在整个射频周期都为正值, 第三项能流散度项 $ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $仅在鞘层塌缩区域为正值, 表明电子在鞘层塌缩阶段到达电极并携其能量流出系统. 该过程产生的能量损失较少, 说明能量主要通过碰撞耗散, 少量通过系统流出. 能流散度项$ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $主要负责电子能量在空间位置上的再分配.

当放电处于DA模式下, 如图3(b0)—(b4)所示, 从时空分布可以看出, 在鞘层扩张和主体区域均存在显著的电子功率吸收现象, 这是DA模式的典型现象. 对于电负性的氧气放电, 由于体区漂移电场的作用, 体区出现了明显的电子功率吸收; 在最大鞘层宽度位置, 由于形成了电正性边缘区域, 导致电子密度梯度最大位置处产生了强烈的双极电场. 这些漂移和双极电场加速了电子. 通过这种方式, 在鞘层塌缩阶段, 电离峰出现在靠近鞘层边缘的最强双极电场的位置, 如图2(c2)所示.

从图4(a2)的时间平均结果来看, 与α模式相同的是, 时间变化项经时间平均后在空间任一点净贡献均为零, 且粒子从电磁场获得的能量仍将全部通过边界流出或与背景中性粒子碰撞耗散. 然而DA模式下, 源项不再集中于鞘层区域, 而是主要分布在体区; 第三项能流散度项$ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $在体区为正值, 在鞘层区域为负值, 说明其将能量由体区输运到鞘层区域. 在DA模式下主要功率吸收区转移至体区, 而主要的能量损失位于鞘层塌缩阶段, 因此形成从体区指向鞘层的能流. 碰撞项极值与能量散度项的负极值同时出现在鞘层边缘, 表明体区吸收的能量部分通过碰撞耗散, 部分通过能流散度项输运到鞘层区域再通过碰撞耗散.

从图4(b2)空间平均结果来看, 碰撞项具有较明显的周期性变化, 其趋势类似于源项但存在相位滞后. 这是因为电子部分能量需先从体区输运至鞘层区域, 再经历与中性粒子的碰撞过程完成能量耗散. 该相位差反映了能量输运与沉积所需的时间延迟. 第三项能流散度项$ \nabla \cdot {\boldsymbol{{\varGamma }}}_{\varepsilon } $在鞘层塌缩阶段携带电子能量流出系统导致的能量损失仍然较少, 其主要功能仍是负责电子能量在空间位置上的再分配.

为清晰地呈现流体动能(机械能)的输运过程, 将所有能够改变能量形式的项, 包括动能源项(表征电磁能向流体动能的转换)和汇项(表征流体动能向流体热能的转换), 放在方程左端, 而将动能本身的时空演化项置于右侧. 流体动能(机械能)输运方程可改写为

1) $ \boldsymbol{J}\cdot \boldsymbol{E} $ 表示电场对电子流体做功, 导致电子沿电场方向的速度增加, 反映电磁能向带电粒子有序动能的转化过程, 是带电粒子动能的源项.

2) $ ({\boldsymbol{\mathcal{P}}}\cdot \nabla )\cdot \boldsymbol{u} $ 表示压力张量$ {\boldsymbol{\mathcal{P}}} $对粒子流体速度$ \boldsymbol{u} $的空间变化的作用效果, 可理解为流体系统中因流场不均匀而产生的应力导致流体元体积和形状的改变, 从而引发动能向热能的转化. 此项也称为“压力-应变相互作用”, 是宏观流体动能被耗散为热能的通道之一^[38–40].

3) $ -\dfrac{1}{2}{\left(\dfrac{\delta {\rho }_{\mathrm{m}}}{\delta t}\right)}_{\text{coll}}{u}^{2} $ 表示由于质量变化引起的动能变化率, 当碰撞使得局部带电粒子数发生变化时, 等效地使系统的有序动能($ \dfrac{1}{2}{\rho }_{\mathrm{m}}{u}^{2} $)成比例的发生改变. 当带电粒子通过电离产生 ($ {\left(\delta {\rho }_{\mathrm{m}}/\delta t\right)}_{\text{coll}}> 0 $) 或通过复合/附着消失 ($ {\left(\delta {\rho }_{\mathrm{m}}/\delta t\right)}_{\text{coll}}< 0 $) 时, 其携带的宏观动能 ($ \dfrac{1}{2}{\rho }_{\mathrm{m}}{u}^{2} $)会被加入或移除系统. 此项的负号表示: 如果产生粒子 ($ {\left(\delta {\rho }_{\mathrm{m}}/\delta t\right)}_{\text{coll}}> 0 $), 需要能量将新产生的粒子加速到当地漂移速度, 因此这是一个动能的“汇”项.

4) $ \boldsymbol{u}\cdot {\left(\delta {{\rho }}_{\mathrm{p}}/\delta t\right)}_{\text{coll}} $ 表示由于碰撞动能传递(拖曳力)导致的流体动能变化率, 该项描述了定向的流体动能因与中性粒子碰撞(动量散射)而不可逆地转化为无序热运动(热能)的过程, 是主要的动能耗散机制之一.

5) $ \partial {\rho }_{\varepsilon , \mathrm{k}}/\partial t $是动能密度随时间的变化项, 对于稳态放电而言其时空平均为零.

6) $ \nabla \cdot ({\rho }_{\varepsilon , \mathrm{k}}\boldsymbol{u}) $为动能密度的对流输运项, 描述流体携带动能宏观运动导致动能密度的空间变化.

7) $ \nabla \cdot ({\boldsymbol{\mathcal{P}}}\cdot \boldsymbol{u}) $为空间不均匀压力导致的动能空间输运. 当流体流经一个存在压力梯度的区域时, 压力对其所做的功不仅改变了局部能量密度, 还通过流动的方式调节不同区域之间动能密度的分布.

图5给出了不同放电模式下, 带电粒子机械能密度和机械能输运方程各项的时空分布. 图6进一步给出了带电粒子机械能输运方程各项时间平均的空间分布, 以及空间平均的时间变化.

我们首先讨论放电处于α模式下的电子流体动能输运过程. $ {\boldsymbol{J}}\cdot {\boldsymbol{E}} $是带电粒子获得能量的唯一机制, 代表其在电场作用下的加速过程. 在此过程中, 外部电磁场被转化成流体的有序动能, 对应于流体单一方向漂移速度的提升. 如图5(a1)所示, 在α模式下, 该过程主要发生在鞘层区域, 这是因为此模式下电负性较低, 电场主要集中于鞘层区域, 包括鞘层内部强电场以及鞘层最大宽度处由于电子密度和离子密度的不同梯度形成的双极场, 而体区电场很弱. 由图5(a2)和图6(a1)—(b1)可见, 压力-应变相互作用(第二项)不论是在空间平均的时间演化还是时间平均的空间分布上, 其变化趋势基本与动能源项相反, 且同样集中在鞘层区域. 进一步观察图5(a1)和图5(a2), 可见在鞘层扩张阶段, 电子动能增加, 而压力-应变相互作用在同一时空区域内几乎“镜像地”表现为强负值, 对应电子流体动能被转换为热能; 反之在鞘层塌缩时, 压力-应变相互作用使电子热能反向转换为动能并进一步退还回电磁能. 这表明压力-应变相互作用是电子动能和热能转换的双向通道. 其中压力-应变相互作用$ ({\boldsymbol{\mathcal{P}}}\cdot \nabla )\cdot \boldsymbol{u} $并不直接依赖于碰撞过程, 其数学作用类似于流体动力学中的宏观体积功$ \boldsymbol{p}\nabla \cdot \boldsymbol{u} $. 它的双向(可逆)特性严格由高频振荡鞘层驱动的电子流体可压缩性所支配. 在α模式(30 MHz)下, 电子动理学受控于高频振荡的鞘层. 当鞘层扩张时, 它向体区推挤电子流体, 导致局部电子被压缩$ \nabla \cdot \boldsymbol{u}< 0 $. 此时, 压力-应变项-$ ({\boldsymbol{\mathcal{P}}} \cdot \nabla )\cdot \boldsymbol{u} $表现为正的能量转化, 将流体宏观动能转化为随机热能. 当鞘层塌缩时, 电子流体回流并膨胀$ \nabla \cdot \boldsymbol{u}> 0 $, 部分热能通过压力做功重新转化为定向动能. 动能输运方程(15)左侧的第三、第四两个碰撞项共同描述了电子在与中性粒子发生碰撞时, 其大尺度、有序流动动能因动量与质量的重新分布而被转化为微观、无序的热能, 其中碰撞阻力项$ \boldsymbol{u}\cdot {\left(\delta {{\rho }}_{\mathrm{p}}/\delta t\right)}_{\text{coll}} $导致的动能向热能的转化, 是主要的碰撞耗散机制.

方程(15)右端的各项描述了电子流体动能在时间上的调制以及在空间上的重新分布, 除边界损失外整体不涉及动能总量的增减. 结合图5所示的各项时空分布可见, 电子动能的各输运项主要集中在鞘层区域. 电子在鞘层区域获得的动能大部分通过压力-应变相互作用和碰撞摩擦转化为热能, 实际保留下的动能仅占总能量(图3(a0))的极小一部分, 并于鞘层塌缩阶段以能流形式流出系统. 由于电子流体动能的输运大部分发生在鞘层区域, 表明电子流体动能对整个等离子体总能量的空间再分配贡献较小. 图6(a1)所示的电子动能对流输运项, 经时间平均后在各空间位置几乎为零, 表明其对电子动能输运贡献较弱; 相比之下, 压力梯度作用引起的动能输运在电子动能的空间再分配过程中起主导作用. 如图6所示, 我们对方程(15)左端和右端分别求和, 可见整体动能守恒关系得到了良好满足, 进一步验证了动能输运方程在当前模拟条件下的数值准确性与物理一致性.

放电处于DA模式时, 能量吸收过程主要发生在等离子体体区, 这是因为在DA模式下, 电负性较高, 电场主要集中在体区. 从动能的时空分布图5(b0)可见, 动能主要集中在体区. 碰撞项导致的动能向热能的转换, 可以分为两部分, 即方程(15)的第三和第四项. 两个碰撞项共同描述了电子在与其他粒子发生碰撞时, 其大尺度、有序流动动能因质量与动量的重新分布而被转化为微观、无序的热能. 由图6(a2)空间分布可以看出, 第三项用于加速新产生的电子到其漂移速度, 这部分很小, 可以忽略; 最主要的部分是第四项, 碰撞摩擦力导致的动能向热能的转化. 在当前的研究情况下, 电子速度远大于中性粒子且弹性碰撞是主导的碰撞过程, 碰撞导致的电子流体在电场方向的动量变化率可以表示为$ -{\left(\delta {\rho }_{\mathrm{ex}, \mathrm{p}}/\delta t\right)}_{\text{coll}}={m}_{\mathrm{e}}{n}_{\mathrm{g}}\displaystyle\int{\sigma }_{\mathrm{m}}|{\boldsymbol{v}}_{\mathrm{e}}|{v}_{\text{ex}}{f}_{\mathrm{e}}{\mathrm{d}}^{3}v $, 其中m_e为电子质量, $ {n}_{\text{g}} $为气体密度, $ {\sigma }_{\text{m}} $为有效电子动量转移碰撞截面, v_e为电子速率, v_ex为电子在x方向的速度, f_e为电子分布函数. 在宏观上, 它表现为一种摩擦力, 因此其作用方向始终与粒子的宏观漂移方向相反. 在不考虑非弹性碰撞的情况下, 电子与中性粒子的弹性碰撞导致的能量(或速率$ |v| $)损失极小, 但其速度$ v $的方向会发生显著的随机偏转. 正是这种在微观层面持续发生的速度方向随机化(或称各向同性化), 导致由漂移速度$ \boldsymbol{u} $体现的粒子群整体的定向动能被转化为无规热能. 从时间分布图6(b2)来看, 第二项与第四项的时间变化趋势基本与机械能源项相反, 并且第四项的值远大于第二项, 即动能通过压力做功向热能的转换要小于碰撞项的贡献. 在DA模式(10 MHz)下, 放电维持主要依赖于体区的漂移电场和双极电场, 这本质上是一种“电阻性”或“欧姆”加热机制. 电子在体区电场驱动下产生定向漂移运动, 但由于该区域远离鞘层边界且中性气体密度较高, 电子获得的定向动能迅速通过与中性粒子的频繁碰撞被“打散”为无序热运动. 这一过程由碰撞摩擦项$ \boldsymbol{u}\cdot {\left(\delta {{\rho }}_{\mathrm{p}}/\delta t\right)}_{\text{coll}} $描述. 这类似于经典欧姆定律中的电阻耗散, 是一个单向的、不可逆的熵增过程. 因此, 在体区主导的DA模式中, 碰撞摩擦成为了连接动能与热能的主要通道.

方程(15)右端的各项表示电子流体动能在时间上被调制、在空间上得到重新分配. 结合图5(b0)—(b7)所示的各项时空分布可见, 电子动能的各项高度局限于主体区域. 与电子动能的对流输运项相比, 压力梯度作用引起的动能输运在电子动能的空间再分配过程中起主导作用.

进一步, 我们讨论“压力-应变相互作用”与“碰撞摩擦”在不同放电模式下的竞争机制. 在α模式(30 MHz)下, 电子主要通过与鞘层的周期性相互作用获得能量, 这被称为随机加热(Stochastic heating)或鞘层膨胀加热. 高驱动频率意味着极快的鞘层运动速度, 电子流体经历了剧烈的周期性形变, 使得压力-应变相互作用在动能-热能转换中占据主导.

当驱动频率降至10 MHz 时, 氧气放电表现出高电负性, 导致体电子密度严重耗尽(从$ \text{30 MHz} $ 时的约$ \text{2×1}{0}^{15}\;{\text{m}}^{-3} $骤降至约$ \text{7×1}{0}^{14}\;{\text{m}}^{-3} $), 形成低电导率的负离子核心. 为了维持外部电路所要求的射频电流($ J\approx e{n}_{\mathrm{e}}{u}_{\mathrm{e}} $)连续性, 在电子密度极低条件下, 电子必须获得更大的漂移速度. 这就迫使体区自洽地建立起较强的漂移电场($ {E}_{\text{drift}} $)来持续加速电子. 在占据绝大部分空间的体区, 电子流体的宏观运动可近似视为准不可压缩响应(速度散度$ \nabla \cdot \boldsymbol{u}\approx 0 $), 因此压力-应变相互作用自然趋近于零, 这与图5(b2)中该项在体区幅值极小的结果一致. Drude近似下, 电子电导率$ \sigma \sim {n}_{\mathrm{e}}{e}^{2}/(m{\nu }_{\mathrm{m}} $), 当电负性增强导致$ {n}_{\text{e}} $降低时, 为维持相同电流密度需满足$ E\propto {\nu }_{\mathrm{m}}/{n}_{\mathrm{e}} $, 从而漂移速度与定向动能增强. 因此碰撞阻力项近似$ u\cdot \left(\dfrac{\delta {\rho }_{\mathrm{p}}}{\delta t}\right)_{\text{coll}}\propto m{n}_{\mathrm{e}}{\nu }_{\mathrm{m}}{u}^{2} $, 可等效写为$ \propto {\nu }_{\mathrm{m}}{j}^{2}/{n}_{\rm e} $. 可见, 电子密度降低会显著放大碰撞阻力项的相对贡献. 与此同时, 驱动频率从30 MHz降至10 MHz也会使相对碰撞尺度$ {\nu }_{\mathrm{m}}/\omega $增大, 从时间尺度上进一步强化碰撞阻力耗散.

为了定量地比较不同放电模式下两者的贡献, 我们对其进行了计算, 在30 MHz (α模式)下, “压力-应变相互作用”与“碰撞摩擦”导致的能量损失分别为1089.87 W/m³与8148.81 W/m³, 两者的比值约为13%. 在10 MHz (DA模式)下, 比值约为7%.

最后, 我们对流体热能(内能)的输运机制进行分析. 同样地, 我们将热能的源项都集中到方程左端, 流体热能输运方程可重写为

1)$ -\left({\boldsymbol{\mathcal{P}}} \cdot \nabla \right) \cdot \boldsymbol{u} + \dfrac{1}{2}{\left(\dfrac{\delta {\rho }_{\mathrm{m}}}{\delta t}\right)}_{\text{coll}}{u}^{2} - \boldsymbol{u} \cdot {\left(\dfrac{\delta {{\rho }}_{\mathrm{p}}}{\delta t}\right)}_{\text{coll}} $这一项描述带电粒子在受到电场加速后, 其定向动能通过压力-应变相互作用和碰撞过程转化为热能. 对应于机械能输运方程(15)中第二、三、四项之和的负值, 是热能源项.

2) $ \partial {\rho }_{\varepsilon , \text{th}}/\partial t $ 表示单位体积热能随时间的局部变化率, 反映热能的瞬时演化.

3) $ \nabla \cdot \left({\rho }_{\varepsilon , \text{th}}\boldsymbol{u}\right) $表示由流体整体运动所携带的热能输运, 反映单位体积内因对流造成的热能净流出或流入.

4) $ \nabla \cdot \boldsymbol{q} $是电子热流矢量的散度, 描述由粒子微观热运动引起的非对流热传导.

5) $ -{\left(\delta {\rho }_{\varepsilon }/\delta t\right)}_{\text{coll}} $表示电子与中性粒子碰撞造成的热能密度变化率, 在不考虑超弹性碰撞的本系统中为正值, 是真正导致热能(总能量)减少的项.

图7给出了电子热能密度和热能输运方程各项的时空分布. 图8进一步给出方程各项时间平均的空间分布, 以及空间平均的时间变化. 如图7(a0)—(b0)所示, 在30 MHz (α模式)下, 电子动能的平均密度约为8.38×10^–6 J/m³, 而电子热能的频率密度约为6.65×10^–4 J/m³, 两者的比值约为1.26%. 在10 MHz (DA 模式)下, 约为2.6%. 在电负性的氧气CCPs放电中, 无论处于α还是DA模式, 电子流体热能皆占据了电子流体总能量的绝大部分, 其热能时空分布与总能量时空分布高度一致.

我们首先讨论放电处于α模式下的电子热能输运过程. 从电子流体动能转换而来的热能可视为热能源项, 其主要位于鞘层区域, 少量位于体区. 热能的增长可由两种机制导致: 一是电子温度升高, 即实际加热; 二是电子密度升高, 即等离子体压缩, 两者可共同作用. 电子热能的对流输运项$ \nabla \cdot \left({\rho }_{\varepsilon , \text{th}}\boldsymbol{u}\right) $几乎局限于鞘层区域, 该项以源项峰值位置为分界, 靠近电极一侧为正值、靠近体区一侧为负值, 表明该项将热能从鞘层靠近电极区域向靠近体区方向输运. 在整个能量再分配过程中, 真正负责将热能从鞘层传递到体区的, 是热流矢量散度项$ \nabla \cdot \boldsymbol{q} $. 如图7(a4), 与其他输运项不同, 它的作用范围覆盖整个放电区域, 且是唯一能够将能量输送至体区的通道. 总体来看, 该项本身不改变系统热能总量, 而是在空间上实现热能的再分配. 电子热能, 或者说电子总能量的真正损耗是由非弹性碰撞引起, 对应热能输运方程中的碰撞损耗项$ -{\left(\delta {\rho }_{\varepsilon }/\delta t\right)}_{\text{coll}} $. 由微观热能通量从鞘层输送的热能和体区的热能最终通过碰撞耗散.

放电处于DA模式时, 由图8(a2)可以看出, 第三项热能的对流项$ \nabla \cdot \left({\rho }_{\varepsilon , \text{th}}\boldsymbol{u}\right) $在鞘层区域为负值、体区为正值, 表明该项将部分热能从体区向鞘层区域输运, 与α模式下的输运方向相反. 热流矢量散度项$ \nabla \cdot \boldsymbol{q} $在体区几乎为零, 这是由于主导热化机制的转变直接决定了能量在空间上的输运方式. 在DA模式下, 由于压力-应变作用在体区显著减弱, 电子在体区获得的宏观定向动能并不会在当地立刻完全转换为热能. 因此与α模式相比, 电子更容易同时携带动能和通过摩擦阻力转化的热能完成跨区输运. 并且体区热能的对流输运项超越了微观热流矢量, 成为主导DA模式下热能空间再分配的核心物理机制. 从图8(b2)时间分布可以看出, 在鞘层塌缩时会有部分热能以离开系统的形式损失. 最后, 与动能分析一致, 我们对热能输运方程(16)的左右两端在时空上分别求和, 其数值结果吻合得很好, 表明整个热能输运方程在数值上满足能量守恒.

本研究构建了一套基于PIC/MCC模拟和玻尔兹曼方程前三阶速度矩的动理学分析框架, 从能流方向和电子流体动能向热能转化机制的视角, 给出了α与DA模式下电子能量吸收-输运-耗散过程的统一对比框架与清晰物理图像, 系统性地揭示了电负性气体O₂在α和DA两种放电模式下, 电子能量从吸收到耗散的完整输运过程. 在α放电模式下, 电子能量吸收主要集中在鞘层区域. 电子在该区域受到电场的定向加速, 获得沿电场方向的动能, 随后通过压力-应变相互作用和碰撞摩擦转化为热能, 鞘层区域的电子流体热能部分经历非弹性碰撞耗散, 部分以热流矢量的形式向体区输运, 并最终在体区经历非弹性碰撞损失. 在DA放电模式下, 电子能量吸收主要集中在体区, 通过体区漂移/双极扩散电场获得沿电场方向的动能, 随后主要通过以弹性碰撞为主的碰撞过程在局部区域实现电子流体动能向热能的转化. 体区的热能部分经历非弹性碰撞损失, 部分以对流输运向鞘层区域输运, 最终通过非弹性碰撞损失. 两种模式下, 电子能量输运的核心差异在于电子能流方向相反, 其根本原因可以归结为电子能量“源-汇”空间分离结构的重构且由于电子流体动能占总能量比例很小, 总能流与热流的输运方向一致.