多脉冲EUV诱导氢等离子体的PIC-漂移扩散混合模型数值模拟

Numerical simulation of PIC-drift-diffusion hybrid model for multi-pulse EUV-induced hydrogen plasma

摘要

基于粒子模拟（PIC）-漂移扩散混合模型，开发了适用于极紫外（EUV）诱导等离子体的二维柱坐标混合模型计算程序。该混合模型充分结合EUV诱导等离子体的特点，采用PIC方法跟踪离子的运动过程，利用漂移扩散模型处理快速进入准平衡态的电子。基于该混合模型，模拟多脉冲EUV诱导氢等离子体的长时间动态演化过程。结果表明，该混合模型不仅可以准确地描述等离子体的动力学行为，还可以显著提高计算效率，扩大可模拟的时间尺度。等离子体的演化呈现出一定的累积效应，随着EUV脉冲数量的增加，氢等离子体的平均密度逐渐上升，电子的平均动能逐渐下降，到达壁面的H离子通量逐渐增加，所有这些参数在脉冲数达到一定数值后趋于稳定。此外，背景气体的压力对于多脉冲累积效应具有显著影响，随着背景气体压力的增加，等离子体参数需要更多的脉冲达到稳定状态。

Abstract

Based on the PIC-drift-diffusion hybrid model， a two-dimensional cylindrical-coordinate hybrid model computational procedure suitable for EUV-induced plasma has been developed， which fully takes into account the characteristics of EUV-induced plasma， and adopts the PIC method to track the motion process of ions， and adopts the drift-diffusion model to deal with electrons rapidly entering the quasi-equilibrium state. Based on this hybrid model， the long-time dynamic evolution of multi-pulse EUV-induced H plasma is simulated. The results show that the hybrid model can accurately describe the dynamical behavior of the plasma， and also significantly improve the computational efficiency and expand the time scale. The plasma evolution shows a certain cumulative effect， with the increase of the number of EUV pulses， the average density of the H plasma gradually rises， the average kinetic energy of the electrons gradually decreases， and the flux of H ions reaching the wall gradually increases， all of these parameters tend to stabilize after the number of pulses reaches a certain value. In addition， the pressure of the background gas has a significant effect on the multi-pulse cumulative effect， and the plasma parameters require more pulses to reach a steady state as the background gas pressure increases.

Graphical abstract

关键词

极紫外光刻 / EUV诱导等离子体 / 多脉冲 / 混合模型

Key words

extreme ultraviolet lithography / EUV-induced plasma / multi-pulse / hybrid model

在当代半导体工程领域，极紫外（extreme ultraviolet， EUV）光刻技术已经成为推进芯片制程极限的先进工具，借助波长仅为13.5 nm的EUV光，可以进一步缩小晶体管的特征尺寸^［1］。然而，随着技术的深入应用，EUV光刻也面临着一系列挑战。特别是，EUV光的使用要求光学系统中的镜面具有极高的反射率。然而，残余气体分子可能会吸附在镜面上，导致高能量的EUV光子引发分子解离和化学反应污染镜面^［2］。因此要保持高效的光束传输和反射，就必须在光刻机的腔室中维持清洁环境，通过特定的背景气体减少镜面的污染和腐蚀^［3］。目前使用氢气作为低吸收率的背景气体，一方面氢分子对EUV光的吸收率较低，能够保证光束传输效率^［4］；另一方面，氢气在EUV辐射下会电离形成等离子体，其所包含的活性粒子可以与吸附在镜面上的污染物发生反应，起到清洁作用，但同时也不可避免地会造成一定程度的侵蚀^［5］。为此，必须对这种等离子体的行为进行深入的了解与分析，从而采取相应的控制与优化策略，延长光学元件的使用寿命。

目前对EUV诱导等离子体的研究主要集中在单脉冲阶段。最早的研究是van der Velden等^［6］进行的，他们研究了氩离子对反射镜的溅射损伤，模拟时间被限制为1 µs。Astakhov等^［7］采用粒子模拟（particle-in-cell，PIC） 方法对微波谐振腔中EUV诱导等离子体的演化过程进行模拟，并与实验结果进行了对比，但受限于PIC方法的计算量，模拟时间被限制在10 µs的范围内。到目前为止，大多数实验使用的EUV光源频率较低（如500 Hz），在此重复率下，通过EUV脉冲产生的等离子体在施加下一个EUV脉冲时将完全消失^［8］。目前的EUV光刻系统通常使用基于锡液滴的激光等离子体光源，具有高度瞬态性，脉冲持续时间通常为100 ns，脉冲频率为50 kHz，对应的脉冲周期为20 µs^［9］。提高重复率会显著改变EUV诱导等离子体的动力学，在更高的重复率下，在2个连续的照射脉冲之间会残留更多的背景等离子体^［10］。因此在使用现代激光等离子体光源的EUV光刻系统中，高频率的EUV脉冲导致脉冲间隔时间小于等离子体衰减时间，这意味着下一道EUV脉冲是建立在之前留下的冷等离子体背景上的。随着脉冲的累积，等离子体会产生积聚，直至达到动态平衡，从而形成冷等离子体背景。在多脉冲等离子体模拟的长时间尺度下，完全使用PIC模拟是不现实的。van de Kerkhof等^［11］先前已经对掩模板附近的EUV诱导等离子体特性进行了研究，他们使用PIC与漂移扩散模型形成的混合模型，实现了长时间尺度的多脉冲等离子体模拟。但目前对多脉冲等离子体的研究相对匮乏，亟需进一步深入。

基于此，本文通过结合PIC和漂移扩散模型，开发了适用于模拟EUV诱导等离子体的二维柱坐标混合模型计算程序。该混合模型充分结合EUV诱导等离子体的特点，采用PIC方法跟踪离子的运动过程，利用漂移扩散模型处理快速进入准平衡态的电子。基于该混合模型，模拟多脉冲EUV诱导氢（H）等离子体的长时间动态演化过程，分析等离子体参数的多脉冲累积效应。研究结果有助于更好地理解EUV等离子体的长时间演化规律。

1 模型与方法

1.1 计算模型

本文采用PIC-漂移扩散混合模型模拟多脉冲EUV诱导H等离子体的动态演化过程，模型的设置参考了van de Ven等^［12］的实验装置，如图1所示。其中模拟区域为圆柱形空腔，长度为100 mm，半径为50 mm，空腔内充满H₂气体。EUV光束由左侧进入腔室，考虑到光束的自然发散现象，光束在左侧20 mm的长度范围内半径为2 mm，随后线性扩大，到达右侧壁面时半径为16 mm，整个光束区域如图1中紫色区域所示。模拟中假设所有壁面均导电并接地（电势为0），且当离子和电子到达壁面时将直接被壁面吸收（从模拟区域移除）。为了模拟真实EUV光刻机中的光学元件，假设右侧壁面具有钌涂层，在EUV光的辐照下因光电效应发射二次电子，二次电子的产率为2.1%^［13］，即平均入射1个光子发射出0.021个电子。其能量分布由下式^［14］描述

S\left(E_{\mathrm{e}},W\right)=\frac{\mathrm{6}{W}^{\mathrm{2}}{E}_{\mathrm{e}}}{{\left(E_{\mathrm{e}}+W\right)}^{\mathrm{4}}}.（1）

其中，S是二次电子能量的概率密度（eV^-1），E_e是二次电子的能量；W为功函数，对于钌，W=6.74 eV。

1.2 PIC模拟

当EUV穿过充满H₂气体的空腔时，部分光子会被H₂吸收形成H等离子体，此即所谓的EUV诱导等离子体。在该过程的初期，EUV电离H₂产生的等离子体处于高度非平衡状态。为更好地重现此阶段等离子体的动态演化过程，我们采用PIC方法进行模拟。模拟所用的程序为本文作者采用C++编程语言基于经典PIC算法^［15］独立开发。模拟过程所采用的参数如表1所示。

考虑到EUV光电离是一个复杂的物理过程，本文的模拟中并没有直接考虑光子与H₂分子之间的碰撞反应，而是通过Lambert-Beer 吸收定律^［6］计算光电离事件的发生频率，具体公式如下

N_{\mathrm{p}\mathrm{h}}=\frac{I}{h\nu }\times n\times \sigma \times L.（2）

其中，N_{\mathrm{p}\mathrm{h}}为光电离事件的发生总次数，I为EUV脉冲能量，h为普朗克常数，\nu为光子频率，n为背景气体的密度，\sigma是相应光电离过程的反应截面，L=100 mm是光束通过的距离。

本文共考虑了3种光电离反应，即非解离光电离、单解离光电离和双解离光电离^［16］。具体的反应式如下

{\mathrm{H}}_{\mathrm{2}}+h\nu \to {\mathrm{H}}_{\mathrm{2}}^{+}+{\mathrm{e}}^{-},（3）

{\mathrm{H}}_{\mathrm{2}}+h\nu \to {\mathrm{H}}^{+}+\mathrm{H}+{\mathrm{e}}^{-},（4）

{\mathrm{H}}_{\mathrm{2}}+h\nu \to \mathrm{2}{\mathrm{H}}^{+}+\mathrm{2}{\mathrm{e}}^{-}.（5）

对应的反应截面分别为4.6×10^-24、1.4×10^-24和0.27×10^-24 m²，能量阈值分别为15.4、18.1和48.9 eV^［16-17］。由于电子的质量远小于各种H离子的质量，在真实的光电离过程中，电子将获得剩余光子能量的绝大部分。在模拟中进行简化处理，假定光子的剩余能量全部由电子吸收，离子的动能则根据室温（300 K）下的麦克斯韦分布进行随机取样，速度方向随机。例如对于式（3）的光电离反应，EUV光子的能量为92 eV，反应阈值为15.4 eV，因此光电离产生的电子能量为92-15.4=76.6 eV，{\mathrm{H}}_{\mathrm{2}}^{+}离子能量则遵循室温下的麦克斯韦分布。对于反应（4），电子的能量为92-18.1= 73.9 eV，{\mathrm{H}}^{+}离子和H原子的能量均根据温度设定。反应（5）的剩余能量为92-48.9=43.1 eV，由于产物中有2个电子，因此需要对能量进行分配，本文采用随机分配方案^［18］，即生成一个随机数R\in (\mathrm{0,1})，其中一个电子的能量为R\times \mathrm{43.1}\mathrm{e}\mathrm{V}，另一个电子的能量为\left(\mathrm{1}-R\right)\times \mathrm{43.1}\mathrm{e}\mathrm{V}。

由前面的描述可知，在EUV光电离过程中会产生大量的高能电子，大部分电子的能量超过H₂的电离能（15.4 eV），因此会进一步电离背景气体。在此，考虑如下几种电子电离反应

{\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}},（6）

{\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}}^{\mathrm{*}},（7）

{\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{2}}^{+}+\mathrm{2}{\mathrm{e}}^{-},（8）

{\mathrm{e}}^{-}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}^{+}+\mathrm{H}+\mathrm{2}{\mathrm{e}}^{-}.（9）

反应（6）表示弹性碰撞，碰撞后电子的运动方向会发生改变。反应（7）为激发碰撞，包括振动激发、转动激发、b^{\mathrm{3}}{\Sigma }_u^{+}和B^{\mathrm{1}}{\Sigma }_u^{+}电子激发，对应的能量阈值分别为0.516、0.044、7.93和12.75 eV。发生激发碰撞后，电子会失去相应阈值的能量。反应（8）对应非解离电离碰撞，能量阈值为15.4 eV。反应（9）是解离电离碰撞，能量阈值为30.6 eV。这些反应的截面数据取自文献［19］。

除电子与背景气体分子的碰撞外，本文还考虑了离子与背景气体之间的碰撞反应，具体反应过程如下

{\mathrm{H}}^{+}\mathrm{动量}\mathrm{转移}:{\mathrm{H}}^{+}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}^{+}+{\mathrm{H}}_{\mathrm{2}},（10）

{\mathrm{H}}^{+}\mathrm{电荷}\mathrm{交换}:{\mathrm{H}}^{+}+{\mathrm{H}}_{\mathrm{2}}\to \mathrm{H}+{\mathrm{H}}_{\mathrm{2}}^{+},（11）

{\mathrm{H}}_{\mathrm{2}}^{+}\mathrm{动量}\mathrm{转移}:{\mathrm{H}}_{\mathrm{2}}^{+}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{2}}^{+}+{\mathrm{H}}_{\mathrm{2}},（12）

{\mathrm{H}}_{\mathrm{2}}^{+}\mathrm{质子}\mathrm{转移}:{\mathrm{H}}_{\mathrm{2}}^{+}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{3}}^{+}+\mathrm{H},（13）

{\mathrm{H}}_{\mathrm{3}}^{+}\mathrm{动量}\mathrm{转移}:{\mathrm{H}}_{\mathrm{3}}^{+}+{\mathrm{H}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{3}}^{+}+{\mathrm{H}}_{\mathrm{2}}.（14）

相关的碰撞截面数据取自文献［20］。

上述电子和离子与背景气体之间的碰撞反应全部采用Monte Carlo碰撞理论^［21］进行处理。考虑到等离子体密度较低，本文的模拟中并未考虑电子和离子之间的复合反应。

1.3 漂移扩散模型

在PIC模拟中，由于需要追踪等离子体中每一个带电粒子的运动轨迹，因此计算量很大。虽然增加宏粒子的权重可以在一定程度上降低计算量，但是可模拟的时间尺度仍然在几个μs的量级，无法达到真实多脉冲EUV的时间尺度。实际上，电子的运动速度很快，在EUV脉冲结束后能够快速充满整个腔体，达到近似平衡状态，因此可以采用基于连续介质理论的漂移扩散模型对电子的运动进行描述^［22］。具体来说，电子满足如下连续性方程

\frac{\partial n_{\mathrm{e}}}{\partial t}+\nabla ·{\mathit{\Gamma }}_{\mathrm{e}}=R_{\mathrm{e}}.（15）

其中，n_{\mathrm{e}}为电子密度，{\mathit{\Gamma }}_{\mathrm{e}}为电子通量，R_{\mathrm{e}}为产生电子的源项。在本文的模拟中，切换到漂移扩散模型时所有电子的能量已经低于H₂的电离碰撞阈值，不会再发生碰撞反应产生新的电子，因此可忽略源项R_{\mathrm{e}}。电子通量{\mathit{\Gamma }}_{\mathrm{e}}可以通过漂移扩散方程进行计算

{\mathit{\Gamma }}_{\mathrm{e}}=-{\mu }_{\mathrm{e}}{n}_{\mathrm{e}}\mathit{E}-D_{\mathrm{e}}\nabla n_{\mathrm{e}}.（16）

其中，{\mu }_{\mathrm{e}}为电子迁移率，\mathit{E}为电场强度矢量，D_{\mathrm{e}}为电子扩散系数。此外，电子还应满足能量守恒方程

\begin{array}{l}\frac{\partial {\omega }_{\mathrm{e}}}{\partial t}+\nabla ·{\mathit{\Gamma }}_{\epsilon }=-e{\mathit{\Gamma }}_{\mathrm{e}}·\mathit{E}+\\ \mathrm{3}{\nu }_{m\mathrm{e}}\frac{m_{\mathrm{e}}}{M}{n}_{\mathrm{e}}{k}_{\mathrm{B}}\left(T_{\mathrm{e}}-T_{\mathrm{b}}\right)+n_{\mathrm{e}}{n}_{\mathrm{b}}\sum _j{\epsilon }_j{k}_{\mathrm{e}j},\end{array}（17）

{\mathit{\Gamma }}_{\epsilon }=-{\mu }_{\epsilon }{\omega }_{\mathrm{e}}\mathit{E}-D_{\epsilon }\nabla {\omega }_{\mathrm{e}}.（18）

其中，电子能量密度 {\omega }_{\mathrm{e}}=\overline{\epsilon }{n}_{\mathrm{e}}，\overline{\epsilon }为平均电子能量。{\mathit{\Gamma }}_{\epsilon }是能量通量，{\nu }_{m\mathrm{e}}是动量转移碰撞频率，m_{\mathrm{e}}和M分别代表电子和H₂分子质量，k_{\mathrm{B}}为玻尔兹曼常数，T_{\mathrm{e}}和T_{\mathrm{b}}分别为电子和背景气体温度，n_{\mathrm{b}}为背景气体密度，k_{\mathrm{e}j}是电子与背景气体之间发生的第j种非弹性碰撞的速率系数，{\epsilon }_j为第j种碰撞的能量阈值。{\mu }_{\epsilon }是电子能量迁移率，D_{\epsilon }是能量扩散系数。

在基于漂移扩散模型的模拟中，边界条件设置与PIC模拟保持一致。由于电子运动至边界时被壁面吸收，因此可直接将边界处的电子密度和能量密度设置为0。PIC模型与漂移扩散模型之间通过泊松方程进行耦合

{\nabla }^{\mathrm{2}}\phi =-\frac{\rho }{{\epsilon }_{\mathrm{0}}},（19）

\mathit{E}=-\nabla \phi .（20）

其中，\phi为电势，\rho为电荷密度。在求解以上方程组之前，需要确定方程中的部分参数，包括电子迁移率{\mu }_{\mathrm{e}}、扩散系数D_{\mathrm{e}}、能量迁移率{\mu }_{\epsilon }、能量扩散系数D_{\epsilon }和反应速率系数k_{\mathrm{e}j}。本文采用局部平均能量近似^［23］计算这些系数，即假设这些系数是平均电子能量\overline{\epsilon }的函数。然后使用电子Boltzmann方程求解软件BOLSIG+^［24］根据PIC模型中的碰撞截面进行计算，结果如图2所示。

1.4 模型验证

为了验证混合模型的准确性，分别采用纯PIC方法和混合模型模拟计算单个EUV脉冲期间H等离子体的演化过程。在混合模型模拟过程中，前2 μs采用PIC方法，从2 μs开始采用漂移扩散模型模拟电子的运动过程。图3所示为离子通量随时间的变化，除PIC和混合模型计算结果外，图3同时展示了van de Ven等^［12］的实验结果。从图中可以看出，模拟结果和实验结果吻合得很好，证明了本文所采用数值模型的准确性。此外，通过对比PIC和混合模型的结果可以看出，两种结果之间的差异很小。而且，对于相同时间（2~15 μs）的模拟，采用纯PIC方法需要72 h，而采用混合模型只需要3.5 h。这些结果表明，从2 μs开始采用混合模型对等离子体的动态演化进行模拟是完全可行的，不仅计算精度完全在可接受的范围内，同时可以大大节省模拟时间。

在EUV诱导等离子体的研究领域，多脉冲等离子体的演化时间往往处于ms量级甚至更长。但此前的研究模型普遍使用PIC模型，很少有文献采用混合模型进行多脉冲长时间尺度的模拟。而若完全采用PIC方法，所需的计算资源是非常庞大的。通过在等离子体趋于平衡后切换为漂移扩散模型，可以突破PIC方法的时间尺度限制，使得多脉冲、长时间的模拟成为可能。这使得我们能够更加全面和深入地认识EUV诱导等离子体的演化规律，优化光刻机的工作参数。特别是随着EUV光源脉冲重复频率的提高，脉冲间的累积效应和长时间演化行为的重要性日益凸显。本文通过对比分析证实了混合模型在处理多脉冲EUV诱导等离子体这一特定问题时的有效性，为后续研究提供了参考。

2 结果与讨论

2.1 单个脉冲内等离子体的演化过程

首先采用混合模型模拟在单个EUV脉冲期间H等离子体的动态演化过程，其中背景气体H₂的压力为5 Pa，EUV脉冲的频率为50 kHz，对应的脉冲周期为20 μs，脉冲宽度为0.1 μs。在脉冲开始的前2 μs内，采用PIC方法模拟离子和电子的运动，此后采用漂移扩散模型模拟电子的运动，离子的动力学行为全程采用PIC方法模拟。

图4展示了模拟区域内电势随时间的演化过程。从图4（a）可以看出，在脉冲初期（0.01 μs），右侧壁面附近出现明显的负电势区，这正是二次电子发射的结果。根据式（2）可以计算得出，对于5 Pa的背景气体压力，EUV光电离H₂产生电子的比例为0.075%，即平均1个EUV光子产生7.5×10^-4个电子，而金属钌在EUV光照射下的二次电子产率为2.1%，即平均入射1个EUV光子发射2.1×10^-2个电子。显然二次电子的产率高得多，也正因如此，在EUV脉冲初期，右壁面附近因电子的聚集而形成负电势区。随着该区域电势的不断降低，部分电子在局部电场的作用下重新回到壁面，并逐渐与二次电子的发射形成动态平衡。与此同时，随着EUV持续电离背景H₂，光束内的等离子体密度逐渐升高，在双极扩散的驱动下电子迅速沿径向向外围扩散，逐渐形成中心高、四周低的电势分布，如图4（b）所示。在轴线处，从左到右随着坐标Z的增加，电势首先迅速升高到最大值，然后逐渐降低，并在右侧壁面附近出现明显的负电势区。这是由于左侧壁面同样会吸收靠近的电子，从而形成鞘层，但是只有右侧壁面发射二次电子，因此电势呈现左高右低的分布规律。随着EUV脉冲的结束（0.1 μs），右侧壁面附近聚集的电子迅速被壁面吸收，负电势区消失，如图4（c）所示。此后，在双极扩散的驱动下，等离子体逐渐向周围扩散，当t=0.2 μs时，整个空间内的负电势区已完全消失（见图4（d）），但仍然可以看出轴线区域的电势明显较高。随着等离子体的进一步演化，电势的空间分布进一步均匀化（见图4（e）和4（f）），当t=10 μs时，脉冲周期过半，此时空间电势值已经非常低，但电势梯度依旧存在，由于图4所有子图共用同一图例，因此图4（f）中的电势梯度看起来并不明显。需要注意的是，虽然在单个EUV脉冲周期结束时电势值非常低，但是等离子体并未消失，根据我们计算的结果，当第1个脉冲周期结束时，等离子体平均密度约为0.3×10¹⁴ m^-3。

2.2 多脉冲的累积效应

为进一步考察多个EUV脉冲产生的累积效应，采用混合模型模拟了等离子体在10个脉冲周期内的演化过程。图5所示为空间平均电子密度随时间的变化曲线。由于EUV脉冲驱动下的光电离以及二次电子发射，等离子体密度迅速上升至接近1.4×10¹⁴ m^-3的峰值，随后由于等离子体在双极扩散的作用下运动至壁面而不断被吸收，等离子体密度迅速下降。Astakhov等^［7］模拟微波谐振腔中的单脉冲EUV诱导氢等离子体的演化过程，将计算得到的平均电子密度和实验结果进行了对比。我们计算得到的单脉冲内电子密度变化趋势与其类似，但由于他们使用的是完全PIC模型，因此只模拟了10 μs的时间尺度。而使用混合模型，可以将模拟时间尺度拓展到200 μs，这进一步说明了混合模型的优势，不仅能够准确捕捉EUV诱导等离子体的动力学行为，而且大大拓展了模拟的时间尺度，为多脉冲效应的研究提供了有力的工具。

随着脉冲数的增加，背景等离子体密度呈现出逐步累积的效应。这是由于连续脉冲作用下，前一脉冲产生的等离子体在下一脉冲到来前并未完全被壁面吸收，新的脉冲在此基础之上产生新的等离子体，二者叠加使得整体密度水平逐步抬升。但密度的累积不会无限增长，随着密度的增加，单位时间内扩散至壁面的等离子体数量增加，即等离子体的损失速率加快，由于混合模型中没有考虑电子与离子之间的复合反应，平衡条件完全取决于等离子体被壁面吸收的速率，当单个脉冲周期内等离子体产生与损失的量相同时，腔体内的等离子体密度达到动态平衡。在当前模拟条件下，等离子体密度约在第6个脉冲达到稳定状态。

图6展示了平均电子能量随时间的变化曲线。可以看出，在第1个脉冲期间，平均电子能量远低于光电离产生的电子能量。根据我们的计算结果，单纯光电离产生的平均电子能量约为70 eV。但是由于EUV脉冲期间右侧壁面会发射大量的低能二次电子，使得全场平均电子能量仅约27.5 eV。随着高能电子与背景H₂气体分子的碰撞反应，电子能量迅速下降至1 eV以下，此后电子电离过程逐渐结束，等离子体通过双极扩散达到平衡。当下一道EUV脉冲开始时，由于模拟区域内已经存在一定密度的低能电子，光电离产生的高能电子占比降低，因此平均电子能量的峰值相比前一个脉冲更低，此即电子能量的多脉冲累积效应。与图5中所示的电子密度相同，平均电子能量约在第6个脉冲达到平衡。

EUV诱导H等离子体的一个重要用途就是对镜面的清洗作用，而清洗的效果取决于H等离子体的通量和能量。图7所示为到达壁面的离子通量随时间的变化曲线。与图5中平均电子密度的变化趋势类似，在单脉冲周期内，到达壁面的离子通量先迅速上升，随后逐渐下降，随着脉冲数的增加，离子通量呈现出明显的累积效应，约在第6个脉冲达到稳定。需要注意的是，图7中离子通量的峰值随脉冲数增加的变化并不明显，相反每个脉冲结束时的通量极小值变化更加显著。这是由于图7中的纵轴采用了对数坐标，离子通量峰值与通量极小值均表现出明显多脉冲累积效应。

2.3 背景气体压力的影响

根据式（2），EUV光电离事件的发生频率与背景气体压力成正比，因此增加或减小背景气体压力会成比例地增加和减小通过光电离产生的等离子体密度。图8所示为不同背景气体压力下H等离子体平均电子密度随时间的变化曲线。从图中可以看出，当压力为1 Pa时，等离子体密度从第2个脉冲开始几乎不随EUV脉冲的增加而变化，且整体变化幅值不明显。根据我们的计算结果，第1个脉冲的空间平均等离子体密度峰值为0.25×10¹⁴ m^-3，第10个EUV脉冲期间的峰值为0.27×10¹⁴ m^-3。而当压力增加到10 Pa时，等离子体密度随EUV脉冲的增加而增加，表现出明显的累积效应。计算结果显示，第1个脉冲时的等离子体密度峰值为3.0×10¹⁴ m^-3，第10个脉冲时的峰值为5.5×10¹⁴ m^-3。由此可以得出结论：背景气体的压力越大，多脉冲产生的累积效应越大。

图9所示为不同背景气体压力下平均电子能量随时间的变化曲线。从图中可以清晰地看出，当背景气体压力为1 Pa时，平均电子能量仅在第1、2脉冲之间存在明显的差异，从第2个脉冲开始，平均电子能量不再发生明显变化，与图8中电子密度的变化规律一致。而当压力为10 Pa时，平均电子能量的累积效应更加明显，下降的幅值也更大，计算结果显示第1道脉冲期间的平均电子能量峰值约为20.5 eV，第10道脉冲期间的平均电子能量峰值仅为2.4 eV。需要注意的是，当背景气体压力增加时，由式（2）可知光电离产生的高能电子数量增加，考虑到二次电子的产额不变，因此单个电子的平均能量应该有所上升。但本文中计算的是等效到网格点的平均电子能量，在数值上等效于空间体积平均，不仅与高能电子的占比有关，还与电子的数量有关。当平均能量峰值出现时，电子已经经历了一系列电离碰撞，在高压力的气体环境中，电子的碰撞频率增加，冷却速度加快，整体的平均能量反而降低。因此，即使在第1个脉冲期间，10 Pa条件下的平均电子能量也低于1 Pa的工况，随着EUV脉冲数的增加，这种差距进一步放大。

图10所示为不同背景气体压力下到达壁面的离子通量随时间的变化曲线。对于1 Pa的背景气体压力，第1道脉冲的离子通量峰值为2.8×10¹⁶ m^-2·s^-1，第10道脉冲的离子通量峰值约为3.8×10¹⁶ m^-2·s^-1。随着脉冲数的增加，离子通量变化不大，这一点与等离子体平均密度的变化类似。当背景气体压力增加到10 Pa，第1道脉冲期间的离子通量峰值为1.6×10¹⁷ m^-2·s^-1，第10道脉冲期间的离子通量峰值为3.8×10¹⁷ m^-2·s^-1。离子通量增加至约2.4倍，累积效应更加明显。这是由于增加背景气体压力一方面会使得到达壁面的离子通量增加，从而需要更多的脉冲使得新产生的离子与壁面吸收的离子达到平衡状态；另一方面，双极扩散时间与背景气体的压力成正比^［25］，背景气体压力越大，离子到达壁面的时间越长，因此越难达到平衡。同时，双极扩散时间还与电子温度成反比^［25］，增加气体压力会增加电子与背景气体分子之间的碰撞频率，导致电子温度下降更快，从而降低双极扩散的速度，增强多脉冲的累积效应。

3 结论

本文通过结合PIC和漂移扩散模型，开发了适用于模拟EUV诱导等离子体的二维柱坐标混合模型。与纯PIC模拟和实验结果的对比，验证了混合模型的正确性和准确性。同时也证实了采用混合模型可以显著提高计算效率，扩大可模拟的时间尺度。

基于该混合模型，模拟了多脉冲EUV诱导H等离子体的长时间动态演化过程。结果表明，在单个EUV脉冲结束时，虽然电势梯度几乎消失，但是等离子体仍有大量残留，因此随着EUV脉冲的增加，等离子体的演化呈现出一定的累积效应。具体表现为：随着EUV脉冲数量的增加，H等离子体的平均密度逐渐上升，电子的平均动能逐渐下降，到达壁面的H离子通量逐渐增加，所有这些参数在脉冲数达到一定数值后趋于稳定。此外，背景气体的压力对于多脉冲累积效应具有显著影响，随着背景气体压力的增加，等离子体参数需要更多的脉冲达到稳定状态，其根本原因是背景气体压力的增加直接影响EUV诱导等离子体的密度，而电子的能量变化也与背景气体的密度密切相关。

参考文献

[1] van de Kerkhof M A， Benschop J P H， Banine V Y. Lithography for now and the future［J］. Solid-State Electronics， 2019， 155： 20-26. DOI： 10.1016/j.sse.2019.03.006 .

[2] Morgan C G， Naulleau P P， Rekawa S B， et al. Removal of surface contamination from EUV mirrors using low-power downstream plasma cleaning［C］// Extreme Ultraviolet （EUV） Lithography. SPIE， 2010， 7636：595-604. DOI： 10.1117/12.846386 .

[3] Dolgov A， Lopaev D， Rachimova T， et al. Comparison of H₂ and He carbon cleaning mechanisms in extreme ultraviolet induced and surface wave discharge plasmas［J］. Journal of Physics D： Applied Physics， 2014， 47（6）： 065205. DOI： 10.1088/0022-3727/47/6/065205 .

[4] Braginsky O V， Kovalev A S， Lopaev D V， et al. Removal of amorphous C and Sn on Mo： Si multilayer mirror surface in hydrogen plasma and afterglow［J］. Journal of Applied Physics， 2012， 111（9）： 093304. DOI： 10.1063/1.4709408 .

[5] Beckers J， van de Ven T， van de Horst R， et al. EUV-induced plasma： a peculiar phenomenon of a modern lithographic technology［J］. Applied Sciences， 2019， 9（14）： 2827. DOI： 10.3390/app9142827 .

[6] van der Velden M H L， Brok W J M， van der Mullen J J A M， et al. Particle-in-cell Monte Carlo simulations of an extreme ultraviolet radiation driven plasma［J］. Physical Review E， 2006， 73（3）： 036406. DOI： 10.1103/physreve.73.036406 .

[7] Astakhov D I， Goedheer W J， Lee C J， et al. Exploring the electron density in plasma induced by EUV radiation： Ⅱ. Numerical studies in argon and hydrogen［J］. Journal of Physics D： Applied Physics， 2016， 49（29）： 295204. DOI： 10.1088/0022-3727/49/29/295204 .

[8] van der Horst R M， Beckers J， Osorio E A， et al. Exploring the electron density in plasma induced by EUV radiation： Ⅰ. Experimental study in hydrogen［J］. Journal of Physics D： Applied Physics， 2016， 49（14）： 145203. DOI： 10.1088/0022-3727/49/14/145203 .

[9] Brandt D C， Fomenkov I V， Farrar N R， et al. LPP EUV source readiness for NXE 3300B［C］// Extreme Ultraviolet （EUV） Lithography V. SPIE， 2014， 9048： 69-76. DOI： 10.1117/12.2048184 .

[10] van der Horst R M. Electron dynamics in EUV-induced plasmas［D］. Eindhoven University of Technology： Eindhoven， The Netherlands， 2015.

[11] van de Kerkhof M， Yakunin A M， Astakhov D， et al. EUV-induced hydrogen plasma： pulsed mode operation and confinement in scanner［J］. Journal of Micro/Nanopatterning， Materials， and Metrology， 2021， 20（3）： 033801. DOI： 10.1117/1.jmm.20.3.033801 .

[12] van de Ven T H M. Ion fluxes towards surfaces exposed to EUV-induced plasmas［D］. Eindhoven University of Technology： Eindhoven， The Netherlands， 2018.

[13] Yakshinskiy B V， Wasielewski R， Loginova E， et al. Carbon accumulation and mitigation processes， and secondary electron yields of ruthenium surfaces［C］// Emerging Lithographic Technologies XI. SPIE， 2007， 6517：875-885. DOI： 10.1117/12.711785 .

[14] Henke B L， Smith J A， Attwood D T. 0.1-10 keV x-ray-induced electron emissions from solids： models and secondary electron measurements［J］. Journal of Applied Physics， 1977， 48（5）： 1852-1866. DOI： 10.1063/1.323938 .

[15] Birdsall C K， Langdon A B. Plasma physics via computer simulation［M］. Boca Raton， Florida， USA： CRC Press， 2018. DOI： 10.1201/9781315275048 .

[16] Chung Y M， Lee E M， Masuoka T， et al. Dissociative photoionization of H₂ from 18 to 124 eV［J］. The Journal of Chemical Physics， 1993， 99（2）： 885-889. DOI： 10.1063/1.465352 .

[17] Kossmann H， Schwarzkopf O， Kammerling B， et al. Photoionisation cross section of H₂ ［J］. Journal of Physics B： Atomic， Molecular and Optical Physics， 1989， 22（14）： L411-L414. DOI： 10.1088/0953-4075/22/14/004 .

[18] van der Velden M H L， Brok W J M， van der Mullen J J A M， et al. Kinetic simulation of an extreme ultraviolet radiation driven plasma near a multilayer mirror［J］. Journal of Applied Physics， 2006， 100（7）： 073303. DOI： 10.1063/1.2356085 .

[19] Yoon J S， Song M Y， Han J M， et al. Cross sections for electron collisions with hydrogen molecules［J］. Journal of Physical and Chemical Reference Data， 2008， 37（2）： 913-931. DOI： 10.1063/1.2838023 .

[20] Tabata T， Shirai T. Analytic cross sections for collisions of H⁺， H₂ ⁺， H₃ ⁺， H， H₂， and H^- with hydrogen molecules［J］. Atomic Data and Nuclear Data Tables， 2000， 76（1）： 1-25. DOI： 10.1006/adnd.2000.0835 .

[21] Nanbu K. Probability theory of electron-molecule， ion-molecule， molecule-molecule， and Coulomb collisions for particle modeling of materials processing plasmas and cases［J］. IEEE Transactions on Plasma Science， 2000， 28（3）： 971-990. DOI： 10.1109/27.887765 .

[22] Economou D J. Hybrid simulation of low temperature plasmas： a brief tutorial［J］. Plasma Processes and Polymers， 2017， 14（1/2）： 1600152. DOI： 10.1002/ppap.201600152 .

[23] Hagelaar G J M， Kroesen G M W. Speeding up fluid models for gas discharges by implicit treatment of the electron energy source term［J］. Journal of Computational Physics， 2000， 159（1）： 1-12. DOI： 10.1006/jcph.2000.6445 .

[24] Hagelaar G J M， Pitchford L C. Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models［J］. Plasma Sources Science Technology， 2005， 14（4）： 722-733. DOI： 10.1088/0963-0252/14/4/011 .

[25] Ellis H W， Pai R Y， McDaniel E W， et al. Transport properties of gaseous ions over a wide energy range［J］. Atomic Data and Nuclear Data Tables， 1976， 17（3）： 177-210. DOI： 10.1016/0092-640X（76）90001-2 .

基金资助

国家自然科学基金(12172356)

中国科学院青年交叉团队项目资助