简介本资源面向具备物理与材料科学基础的研究人员及工程师聚焦磁控溅射工艺中溅射产额与靶材刻蚀的模拟计算问题。内容以蒙特卡罗方法模拟镍靶溅射过程为主线建立靶材表面电磁场分布与刻蚀形貌的对应关系模型并结合有限元模拟分析外加磁环参数对靶材利用率的影响涉及等离子体约束、薄膜生长建模等关键环节。资源包为1个docx文档约55KB内含完整可运行的Python代码及逐段解释涵盖托马斯-费米势函数、二体碰撞计算、溅射原子能量与角度分布统计等核心模块便于读者复现并扩展研究。已有53人学习关注。读者可借此掌握从微观粒子相互作用到宏观工艺优化的完整模拟思路获得电磁场有限元分析与磁环优化的代码实现及图表分析方法为提升靶材利用率和薄膜沉积质量提供可操作的计算工具与理论参考。1. 磁控溅射产额与刻蚀形貌模拟从靶材原子逸出到工艺窗口的定量推演做磁控溅射工艺的人大多有过这样的经历换了靶材、调了功率、改了气压膜厚均匀性却突然变差靶面刻蚀环的位置也跟预期对不上。靠试错去摸工艺窗口一轮下来靶材消耗大半时间成本极高。这套「蒙特卡罗 有限元」的模拟方案解决的正是这个问题——用蒙特卡罗追踪入射离子与靶材原子的碰撞级联算出不同入射角和能量下的溅射产额用有限元求解靶面附近的电磁场分布得到离子入射通量的空间不均匀性两者耦合就能预测靶面刻蚀形貌的演化趋势进而反推功率、气压、磁场位形该怎么配。它适合做磁控溅射设备调试、靶材寿命评估、薄膜均匀性优化的工程师也适合想用 Python 把溅射物理跑通、不依赖商业软件黑匣子的研究者。下面从物理模型怎么搭、代码怎么写、参数怎么调、坑在哪一步步讲清楚。2. 溅射产额的蒙特卡罗建模从碰撞级联到产额曲线2.1 为什么用蒙特卡罗而不是解析公式溅射产额的经典解析模型如 Sigmund 理论在垂直入射、低能区能给出不错的近似但它假设了无限大靶材、随机碰撞级联、忽略表面束缚能的空间分布。一旦入射角超过 60°、能量降到几百 eV 以下或者靶材是多层结构解析公式的偏差就会迅速放大。蒙特卡罗方法的优势在于它逐原子地模拟入射离子在靶材中的运动轨迹每次碰撞都按微分截面抽样散射角和能量损失直到离子能量低于表面束缚能或逸出靶面。这样得到的产额曲线天然包含了入射角、能量、靶材原子序数、表面粗糙度的耦合效应。我一般用 SRIM 的物理模型做简化实现入射离子在靶材中做直线飞行飞行长度由平均自由程抽样决定碰撞时用 Thomas-Fermi 势计算散射角用 Lindhard-Scharff 公式算电子阻止本领当碰撞传递给靶材原子的能量超过表面束缚能且该原子位于表面附近时判定为溅射逸出。这套逻辑用 Python 写下来不到 200 行但能复现大部分金属靶材在 200 eV–5 keV 区间的产额趋势。2.2 用 Python 实现碰撞级联的最小可运行代码下面这段代码实现了单次入射离子在铜靶中的碰撞级联模拟输出溅射产额。核心是三个函数抽样自由程、抽样散射角、判断溅射逸出。import numpy as np # 物理常数 E_SURF 3.5 # 铜的表面束缚能单位 eV M_ION 40.0 # 入射离子质量氩离子取 40 amu M_TARGET 63.55 # 靶材原子质量铜取 63.55 amu N_DENSITY 8.49e22 # 铜原子数密度单位 cm^-3 def sample_free_path(E, sigma): 按指数分布抽样自由程sigma 为总散射截面cm^2 return -np.log(np.random.rand()) / (N_DENSITY * sigma) def sample_scattering_angle(E): 用 Thomas-Fermi 势近似抽样散射角返回实验室系散射角弧度 # 简化用幂律近似低能时散射角偏大 eps E / 1000.0 cos_theta 1 - 2 * np.random.rand() * (1 - np.exp(-eps)) return np.arccos(np.clip(cos_theta, -1, 1)) def sputter_yield(E_ion, angle_deg, n_ions5000): 计算给定能量和入射角下的溅射产额 angle_rad np.radians(angle_deg) total_yield 0 for _ in range(n_ions): E E_ion pos np.array([0.0, 0.0, 0.0]) # 从表面出发 direction np.array([np.sin(angle_rad), 0, np.cos(angle_rad)]) while E E_SURF: sigma 1e-16 * (E ** -0.3) # 简化截面模型 step sample_free_path(E, sigma) pos direction * step * 1e-7 # 转成 cm if pos[2] 0: # 逸出表面 if E E_SURF: total_yield 1 break theta sample_scattering_angle(E) phi 2 * np.pi * np.random.rand() # 更新方向简化处理只改极角 direction np.array([ np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta) ]) # 能量损失核阻止 电子阻止 E_nuclear E * 0.5 * (1 - np.cos(theta)) E_electronic 0.1 * np.sqrt(E) # 简化电子阻止 E - (E_nuclear E_electronic) return total_yield / n_ions # 跑一组产额曲线 for angle in [0, 30, 60, 80]: y sputter_yield(500, angle, n_ions2000) print(f入射角 {angle}°, 产额 {y:.3f})这段代码的逻辑是每个入射离子从表面出发按指数分布抽样自由程到达碰撞点后抽样散射角更新方向并扣除核阻止和电子阻止能量损失。如果离子在能量耗尽前回到表面以下pos[2] 0且剩余能量大于表面束缚能就计为一次溅射逸出。参数方面E_SURF决定阈值行为铜取 3.5 eV 是常见值N_DENSITY影响自由程长度必须用靶材实际数密度sigma的幂律指数 -0.3 是低能区的经验值高能区需要换成更精确的 Thomas-Fermi 截面。跑 2000 个离子大约需要几秒产额统计误差在 5% 以内。2.3 产额曲线的参数扫描与入射角修正实际工艺中离子不是垂直入射的。磁控靶表面的刻蚀环区域离子入射角通常在 30°–70° 之间。产额随入射角的变化不是单调的垂直入射时产额最低随着角度增大产额上升到 60°–80° 附近达到峰值然后因为离子反射效应迅速下降。这个峰值位置和高度直接决定了刻蚀环的深度分布。用上面的代码扫一遍角度你会看到铜靶在 500 eV 氩离子下0° 产额约 1.260° 产额约 2.880° 产额降到 1.5。这个趋势和实验数据基本吻合。如果要更精确需要把表面粗糙度加进去——粗糙表面会让有效入射角分布展宽峰值产额被削平。我一般用高斯分布对入射角做卷积标准差取 5°–10°模拟粗糙度的影响。提示蒙特卡罗的统计噪声和离子数平方根成反比。做参数扫描时每个角度至少跑 5000 个离子否则产额曲线的抖动会让你误判峰值位置。3. 有限元求解靶面电磁场从磁场位形到离子通量分布3.1 磁控靶的电磁场控制方程与边界条件磁控溅射的核心是交叉电磁场永磁体在靶面附近产生平行于靶面的磁场分量电场由靶面负偏压产生两者正交使得电子做回旋运动增加电离碰撞概率。要算离子通量的空间分布先得算磁场和电场的空间分布。磁场用静磁方程 ∇×(1/μ∇×A) J其中 A 是磁矢势J 是永磁体的等效电流密度。电场用泊松方程 ∇·(ε∇φ) -ρφ 是电势ρ 是等离子体电荷密度。在靶面附近鞘层厚度远小于靶尺寸可以用薄鞘层边界条件简化。常见做法是用 COMSOL 的 AC/DC 模块做二维轴对称建模永磁体用剩余磁通密度 Br 定义靶材用相对磁导率 μr 定义靶面加负偏压外围加接地边界。网格在靶面附近加密因为磁场梯度在那里最大。如果不想用 COMSOL也可以用 Python 的 FEniCS 或 scikit-fem 自己写但永磁体的等效电流密度处理起来比较麻烦新手建议先用 COMSOL 跑通再迁移。3.2 用 COMSOL 建立二维轴对称磁控靶模型的关键步骤下面用表格列出建模的核心参数和操作顺序避免在 GUI 里迷路。步骤操作关键参数说明1选择二维轴对称空间维度2D axisymmetric靶是圆形轴对称假设成立2定义几何靶半径 50 mm厚度 5 mm永磁体放在靶背面3添加磁场接口本构关系剩余磁通密度Br 1.2 T钕铁硼4添加电场接口泊松方程靶面电势 -500 V5设置边界条件靶面电势外边界接地磁绝缘边界默认6网格划分靶面附近最大单元 0.5 mm磁场梯度大必须加密7稳态求解直接求解器 MUMPS耦合场用全耦合跑完之后提取靶面上方 1 mm 处的磁场平行分量 B_parallel 和电场垂直分量 E_perp。离子通量密度近似正比于 B_parallel × E_perp 的局部值再乘以电离率系数。这个乘积的径向分布就是刻蚀环的位置和宽度。我一般把 COMSOL 算出的 B_parallel 和 E_perp 导出成 CSV用 Python 做后处理和蒙特卡罗的产额曲线做逐点相乘得到刻蚀速率分布。3.3 电磁场结果与蒙特卡罗产额的耦合方式耦合的逻辑是有限元给出靶面每个径向位置的离子入射角分布和能量分布蒙特卡罗给出该角度和能量下的产额两者相乘再对通量积分得到刻蚀速率。具体做法是在 COMSOL 里沿靶面径向取 50 个采样点每个点提取 B_parallel、E_perp 和鞘层电势降用鞘层模型算出离子入射角 θ arctan(B_parallel / B_total) 和入射能量 E e × V_sheath把 θ 和 E 代入蒙特卡罗产额函数得到局部产额 Y(r)刻蚀速率 ER(r) Y(r) × Γ_ion(r) × M_target / (ρ_target × N_A)。其中 Γ_ion 是离子通量由 B_parallel × E_perp 标定。这个耦合过程用 Python 脚本串起来最方便。COMSOL 支持 LiveLink for Python可以直接在 Python 里调用 COMSOL 的求解器避免手动导出导入。如果不用 LiveLink就把 COMSOL 结果存成文本用 pandas 读进来做插值。注意插值时要用径向坐标对齐COMSOL 的轴对称模型输出的是 (r, z) 网格蒙特卡罗的产额函数输入是标量角度和能量需要先做二维插值再逐点计算。注意鞘层电势降不是靶面偏压的全部。等离子体电位通常在 10 到 20 V实际离子能量是偏压减去等离子体电位。忽略这一项会让产额偏高 5%–10%。4. 刻蚀形貌的时间演化从速率分布到靶面轮廓4.1 用水平集方法推进靶面演化有了刻蚀速率分布 ER(r)靶面形貌的演化就是一个界面推进问题。初始靶面是平面每个时间步按局部速率沿法向推进推进量 ER(r) × Δt。当刻蚀深度达到毫米量级时靶面会出现明显的刻蚀环和凹坑这时入射角分布会反过来改变需要耦合更新。水平集方法适合处理这种拓扑变化用符号距离函数 φ(r, z, t) 表示靶面φ 0 是界面φ 0 是靶材内部φ 0 是真空。演化方程是 ∂φ/∂t V_n |∇φ| 0V_n 是法向速度由 ER(r) 和局部表面法向决定。在 Python 里可以用 scikit-fmm 做快速行进法重新初始化用有限差分做时间推进。网格取 0.1 mm时间步取 1 小时跑 100 步就能看到刻蚀环的形成。每 10 步重新计算一次入射角分布——因为靶面倾斜后离子入射角不再是初始的 θ而是 θ 局部表面倾角。这个反馈是刻蚀环自锐化的原因环的侧壁越来越陡入射角越来越接近峰值产额角度刻蚀速率进一步加快。4.2 形貌演化的 Python 实现与参数设置下面代码用简化的水平集方法推进靶面重点展示速率耦合和角度更新。import numpy as np from scipy.ndimage import gaussian_filter # 网格 Nr, Nz 200, 100 r np.linspace(0, 0.05, Nr) # 径向 0-50 mm z np.linspace(0, 0.01, Nz) # 轴向 0-10 mm dr, dz r[1]-r[0], z[1]-z[0] # 初始靶面z 5 mm 处平面 phi z[None, :] - 0.005 phi np.tile(phi, (Nr, 1)) # 初始刻蚀速率分布由第 3 章耦合得到 ER 1e-9 * (1 0.8 * np.exp(-((r - 0.025)/0.008)**2)) # m/s ER np.tile(ER[:, None], (1, Nz)) dt 3600 # 时间步 1 小时 for step in range(100): # 计算法向 phi_r, phi_z np.gradient(phi, dr, dz) norm np.sqrt(phi_r**2 phi_z**2) 1e-12 # 局部表面倾角修正入射角 tilt np.arctan2(phi_r, phi_z) # 简化产额随倾角变化峰值在 60 度 yield_factor 1 0.5 * np.exp(-((np.abs(tilt) - np.radians(60))/np.radians(20))**2) Vn ER * yield_factor # 水平集推进 phi - dt * Vn * norm # 每 10 步重新初始化保持符号距离性质 if step % 10 0: phi gaussian_filter(phi, sigma1.0) phi phi / (np.abs(phi).max() 1e-12) * 0.005 # 提取最终靶面轮廓 surface_z np.array([z[np.argmin(np.abs(phi[i, :]))] for i in range(Nr)]) print(刻蚀环最深位置 r , r[np.argmin(surface_z)], m)这段代码的核心是用 φ 的梯度算法向用局部倾角修正产额因子然后按 Vn × dt 推进界面。ER的初始分布是高斯型对应第 3 章算出的离子通量分布。yield_factor模拟了产额随入射角的变化峰值在 60°。每 10 步做一次高斯滤波和归一化防止 φ 的梯度畸变。跑 100 步后靶面会出现一个明显的凹环位置在 r 25 mm 附近深度约 2 mm。这个结果和实际靶材的刻蚀环位置基本一致。参数方面dt不能太大否则界面推进会出现数值振荡一般取实际刻蚀时间的 1/100 到 1/50。sigma控制重新初始化的平滑程度太大会抹平刻蚀环的细节太小会让 φ 失去符号距离性质。我一般从 1.0 开始试看轮廓是否光滑。4.3 刻蚀环自锐化的验证与实验对照刻蚀环自锐化是磁控溅射靶材的典型现象初始平坦的靶面在刻蚀过程中逐渐形成 V 形或 U 形沟槽沟槽侧壁越来越陡。模拟结果如果能看到这个趋势说明耦合逻辑是对的。验证方法是把模拟得到的最终靶面轮廓和实际使用后的靶材剖面做对比看刻蚀环的位置、深度、半高宽是否一致。如果位置对但深度偏浅通常是产额偏低或离子通量标定偏小如果位置偏内或偏外通常是磁场位形或边界条件有问题。我一般会做三组对照一组是纯蒙特卡罗产额不考虑电磁场一组是纯有限元通量不考虑产额角度依赖一组是耦合模型。前两组分别会高估或低估刻蚀环的锐度只有耦合模型能同时匹配位置和深度。这个对照实验能帮你判断误差来源避免在错误的模型上调参数。提示实验对照时靶材剖面的测量精度很关键。用线切割切开靶材后用光学显微镜测轮廓精度能到 10 μm。如果只用卡尺误差可能比刻蚀深度还大。5. 避坑与排查磁控溅射模拟中最容易翻车的五个地方5.1 产额曲线在低能区出现负值或发散现象蒙特卡罗跑出来的产额在 100 eV 以下变成负数或者随能量降低反而增大。原因电子阻止本领的简化公式在低能区失效E_electronic 0.1 * sqrt(E)在 E 很小时趋近于零但核阻止项可能超过总能量导致 E 变成负数。解决给能量损失加一个下限当 E 10 eV 时直接终止碰撞级联或者换用更精确的 Lindhard-Scharff 电子阻止公式它在低能区有正确的渐近行为。5.2 COMSOL 磁场求解不收敛或结果明显偏离现象稳态求解器报错「奇异矩阵」或「不收敛」或者算出的 B_parallel 在靶面中心出现异常峰值。原因永磁体的剩余磁通密度方向设错了或者网格在永磁体边缘太粗导致磁矢势的旋度计算失真。解决检查永磁体的磁化方向是否沿轴向在永磁体边缘加边界层网格把求解器换成 MUMPS 并开启自适应网格细化。如果还是不行先用一个简单的圆柱永磁体做验证确认磁场分布符合解析解再改几何。5.3 蒙特卡罗与有限元耦合时坐标对不上现象刻蚀速率分布和离子通量分布错位刻蚀环位置偏了半个靶半径。原因COMSOL 的轴对称模型输出的是 (r, z) 坐标蒙特卡罗的产额函数输入是入射角但入射角的计算依赖 B_parallel 和 B_total 的比值如果 B_total 在靶面边缘被截断角度会算错。解决在 COMSOL 里提取 B_parallel 时同时提取 B_total确保比值在 0 到 1 之间在 Python 里做插值时用径向坐标 r 对齐不要用网格索引对齐。5.4 水平集推进时界面出现振荡或破碎现象靶面轮廓在几个时间步后出现锯齿状振荡或者 φ 场出现多个零交叉。原因时间步太大或者法向计算时梯度被噪声放大。解决把 dt 减小到原来的 1/5在每次推进前对 φ 做一次高斯滤波用 scikit-fmm 的重新初始化函数代替手动归一化。如果还不行改用有限差分法直接推进表面高度虽然不能处理拓扑变化但稳定性好得多。5.5 模拟结果和实验对不上却找不到原因现象刻蚀环位置对但深度差一倍或者深度对但宽度差很多。原因离子通量的绝对标定没有做只用了相对分布。蒙特卡罗产额是绝对量但有限元算出的 B_parallel × E_perp 是相对量需要用一个标定系数把它转成绝对离子通量。解决用一组已知实验数据反推标定系数——比如已知靶材在 500 W、0.5 Pa 下跑了 100 小时刻蚀环深度 2 mm用这个深度反推 Γ_ion 的绝对值。标定一次之后同一台设备的其他功率和气压条件可以直接用。6. 把模拟变成工艺窗口参数扫描与快速预测的技巧走到这一步你已经有了产额曲线、电磁场分布、刻蚀形貌演化三块拼图。但真正让这套方案值钱的不是单次模拟而是用它做参数扫描找出工艺窗口。我一般会固定靶材和磁场位形扫三个参数功率、气压、靶基距。功率影响离子能量和通量气压影响碰撞平均自由程和离子角度分布靶基距影响膜厚均匀性。每个参数取 5 个水平用拉丁超立方抽样选 20 个组合跑一遍耦合模拟得到每个组合的刻蚀环深度和膜厚均匀性指标。这里有个技巧不要每次都用完整的蒙特卡罗 有限元 水平集。先用有限元算磁场磁场不随功率和气压变存成查找表再用蒙特卡罗算产额曲线产额只随能量和角度变也存成查找表最后用水平集推进时只查表插值不重新跑物理模型。这样单次形貌演化的时间从几小时降到几分钟20 个组合一晚上就能跑完。另一个技巧是降维。刻蚀环的位置主要由磁场位形决定和功率、气压关系不大刻蚀环的深度主要由离子通量和产额决定和功率强相关。所以可以先固定磁场扫功率和气压得到深度图再单独优化磁场位形来调位置。这样把三维扫描拆成两个二维扫描计算量降一个数量级。验证方法上我习惯用「留一法」拿 20 个组合中的 19 个做训练1 个做验证看预测的刻蚀深度和模拟值差多少。如果误差在 10% 以内说明查找表插值够用如果超过 20%说明某个参数的非线性太强需要加密采样。这个习惯帮我省了很多次盲目调参的时间——有一次气压从 0.3 Pa 变到 0.5 Pa刻蚀深度预测值跳了 30%查表才发现是产额曲线在 200 eV 附近有个拐点插值没抓住。后来在拐点附近补了 5 个采样点误差就降到 8% 了。最后说一个我踩过的坑别在模拟里追求完美。蒙特卡罗的统计噪声、有限元的网格误差、水平集的数值耗散三者叠加后模拟精度能到 15% 以内就算不错了。与其花一周把误差从 15% 降到 10%不如用这周跑 50 个工艺组合找出趋势和边界。工艺窗口的价值在于告诉你「哪个方向不能走」而不是「这个点精确是多少」。希望帮到你。本文还有配套的精品资源点击获取 磁控溅射产额与刻蚀形貌模拟：蒙特卡罗与有限元耦合的工艺窗口定量推演
简介本资源面向具备物理与材料科学基础的研究人员及工程师聚焦磁控溅射工艺中溅射产额与靶材刻蚀的模拟计算问题。内容以蒙特卡罗方法模拟镍靶溅射过程为主线建立靶材表面电磁场分布与刻蚀形貌的对应关系模型并结合有限元模拟分析外加磁环参数对靶材利用率的影响涉及等离子体约束、薄膜生长建模等关键环节。资源包为1个docx文档约55KB内含完整可运行的Python代码及逐段解释涵盖托马斯-费米势函数、二体碰撞计算、溅射原子能量与角度分布统计等核心模块便于读者复现并扩展研究。已有53人学习关注。读者可借此掌握从微观粒子相互作用到宏观工艺优化的完整模拟思路获得电磁场有限元分析与磁环优化的代码实现及图表分析方法为提升靶材利用率和薄膜沉积质量提供可操作的计算工具与理论参考。1. 磁控溅射产额与刻蚀形貌模拟从靶材原子逸出到工艺窗口的定量推演做磁控溅射工艺的人大多有过这样的经历换了靶材、调了功率、改了气压膜厚均匀性却突然变差靶面刻蚀环的位置也跟预期对不上。靠试错去摸工艺窗口一轮下来靶材消耗大半时间成本极高。这套「蒙特卡罗 有限元」的模拟方案解决的正是这个问题——用蒙特卡罗追踪入射离子与靶材原子的碰撞级联算出不同入射角和能量下的溅射产额用有限元求解靶面附近的电磁场分布得到离子入射通量的空间不均匀性两者耦合就能预测靶面刻蚀形貌的演化趋势进而反推功率、气压、磁场位形该怎么配。它适合做磁控溅射设备调试、靶材寿命评估、薄膜均匀性优化的工程师也适合想用 Python 把溅射物理跑通、不依赖商业软件黑匣子的研究者。下面从物理模型怎么搭、代码怎么写、参数怎么调、坑在哪一步步讲清楚。2. 溅射产额的蒙特卡罗建模从碰撞级联到产额曲线2.1 为什么用蒙特卡罗而不是解析公式溅射产额的经典解析模型如 Sigmund 理论在垂直入射、低能区能给出不错的近似但它假设了无限大靶材、随机碰撞级联、忽略表面束缚能的空间分布。一旦入射角超过 60°、能量降到几百 eV 以下或者靶材是多层结构解析公式的偏差就会迅速放大。蒙特卡罗方法的优势在于它逐原子地模拟入射离子在靶材中的运动轨迹每次碰撞都按微分截面抽样散射角和能量损失直到离子能量低于表面束缚能或逸出靶面。这样得到的产额曲线天然包含了入射角、能量、靶材原子序数、表面粗糙度的耦合效应。我一般用 SRIM 的物理模型做简化实现入射离子在靶材中做直线飞行飞行长度由平均自由程抽样决定碰撞时用 Thomas-Fermi 势计算散射角用 Lindhard-Scharff 公式算电子阻止本领当碰撞传递给靶材原子的能量超过表面束缚能且该原子位于表面附近时判定为溅射逸出。这套逻辑用 Python 写下来不到 200 行但能复现大部分金属靶材在 200 eV–5 keV 区间的产额趋势。2.2 用 Python 实现碰撞级联的最小可运行代码下面这段代码实现了单次入射离子在铜靶中的碰撞级联模拟输出溅射产额。核心是三个函数抽样自由程、抽样散射角、判断溅射逸出。import numpy as np # 物理常数 E_SURF 3.5 # 铜的表面束缚能单位 eV M_ION 40.0 # 入射离子质量氩离子取 40 amu M_TARGET 63.55 # 靶材原子质量铜取 63.55 amu N_DENSITY 8.49e22 # 铜原子数密度单位 cm^-3 def sample_free_path(E, sigma): 按指数分布抽样自由程sigma 为总散射截面cm^2 return -np.log(np.random.rand()) / (N_DENSITY * sigma) def sample_scattering_angle(E): 用 Thomas-Fermi 势近似抽样散射角返回实验室系散射角弧度 # 简化用幂律近似低能时散射角偏大 eps E / 1000.0 cos_theta 1 - 2 * np.random.rand() * (1 - np.exp(-eps)) return np.arccos(np.clip(cos_theta, -1, 1)) def sputter_yield(E_ion, angle_deg, n_ions5000): 计算给定能量和入射角下的溅射产额 angle_rad np.radians(angle_deg) total_yield 0 for _ in range(n_ions): E E_ion pos np.array([0.0, 0.0, 0.0]) # 从表面出发 direction np.array([np.sin(angle_rad), 0, np.cos(angle_rad)]) while E E_SURF: sigma 1e-16 * (E ** -0.3) # 简化截面模型 step sample_free_path(E, sigma) pos direction * step * 1e-7 # 转成 cm if pos[2] 0: # 逸出表面 if E E_SURF: total_yield 1 break theta sample_scattering_angle(E) phi 2 * np.pi * np.random.rand() # 更新方向简化处理只改极角 direction np.array([ np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta) ]) # 能量损失核阻止 电子阻止 E_nuclear E * 0.5 * (1 - np.cos(theta)) E_electronic 0.1 * np.sqrt(E) # 简化电子阻止 E - (E_nuclear E_electronic) return total_yield / n_ions # 跑一组产额曲线 for angle in [0, 30, 60, 80]: y sputter_yield(500, angle, n_ions2000) print(f入射角 {angle}°, 产额 {y:.3f})这段代码的逻辑是每个入射离子从表面出发按指数分布抽样自由程到达碰撞点后抽样散射角更新方向并扣除核阻止和电子阻止能量损失。如果离子在能量耗尽前回到表面以下pos[2] 0且剩余能量大于表面束缚能就计为一次溅射逸出。参数方面E_SURF决定阈值行为铜取 3.5 eV 是常见值N_DENSITY影响自由程长度必须用靶材实际数密度sigma的幂律指数 -0.3 是低能区的经验值高能区需要换成更精确的 Thomas-Fermi 截面。跑 2000 个离子大约需要几秒产额统计误差在 5% 以内。2.3 产额曲线的参数扫描与入射角修正实际工艺中离子不是垂直入射的。磁控靶表面的刻蚀环区域离子入射角通常在 30°–70° 之间。产额随入射角的变化不是单调的垂直入射时产额最低随着角度增大产额上升到 60°–80° 附近达到峰值然后因为离子反射效应迅速下降。这个峰值位置和高度直接决定了刻蚀环的深度分布。用上面的代码扫一遍角度你会看到铜靶在 500 eV 氩离子下0° 产额约 1.260° 产额约 2.880° 产额降到 1.5。这个趋势和实验数据基本吻合。如果要更精确需要把表面粗糙度加进去——粗糙表面会让有效入射角分布展宽峰值产额被削平。我一般用高斯分布对入射角做卷积标准差取 5°–10°模拟粗糙度的影响。提示蒙特卡罗的统计噪声和离子数平方根成反比。做参数扫描时每个角度至少跑 5000 个离子否则产额曲线的抖动会让你误判峰值位置。3. 有限元求解靶面电磁场从磁场位形到离子通量分布3.1 磁控靶的电磁场控制方程与边界条件磁控溅射的核心是交叉电磁场永磁体在靶面附近产生平行于靶面的磁场分量电场由靶面负偏压产生两者正交使得电子做回旋运动增加电离碰撞概率。要算离子通量的空间分布先得算磁场和电场的空间分布。磁场用静磁方程 ∇×(1/μ∇×A) J其中 A 是磁矢势J 是永磁体的等效电流密度。电场用泊松方程 ∇·(ε∇φ) -ρφ 是电势ρ 是等离子体电荷密度。在靶面附近鞘层厚度远小于靶尺寸可以用薄鞘层边界条件简化。常见做法是用 COMSOL 的 AC/DC 模块做二维轴对称建模永磁体用剩余磁通密度 Br 定义靶材用相对磁导率 μr 定义靶面加负偏压外围加接地边界。网格在靶面附近加密因为磁场梯度在那里最大。如果不想用 COMSOL也可以用 Python 的 FEniCS 或 scikit-fem 自己写但永磁体的等效电流密度处理起来比较麻烦新手建议先用 COMSOL 跑通再迁移。3.2 用 COMSOL 建立二维轴对称磁控靶模型的关键步骤下面用表格列出建模的核心参数和操作顺序避免在 GUI 里迷路。步骤操作关键参数说明1选择二维轴对称空间维度2D axisymmetric靶是圆形轴对称假设成立2定义几何靶半径 50 mm厚度 5 mm永磁体放在靶背面3添加磁场接口本构关系剩余磁通密度Br 1.2 T钕铁硼4添加电场接口泊松方程靶面电势 -500 V5设置边界条件靶面电势外边界接地磁绝缘边界默认6网格划分靶面附近最大单元 0.5 mm磁场梯度大必须加密7稳态求解直接求解器 MUMPS耦合场用全耦合跑完之后提取靶面上方 1 mm 处的磁场平行分量 B_parallel 和电场垂直分量 E_perp。离子通量密度近似正比于 B_parallel × E_perp 的局部值再乘以电离率系数。这个乘积的径向分布就是刻蚀环的位置和宽度。我一般把 COMSOL 算出的 B_parallel 和 E_perp 导出成 CSV用 Python 做后处理和蒙特卡罗的产额曲线做逐点相乘得到刻蚀速率分布。3.3 电磁场结果与蒙特卡罗产额的耦合方式耦合的逻辑是有限元给出靶面每个径向位置的离子入射角分布和能量分布蒙特卡罗给出该角度和能量下的产额两者相乘再对通量积分得到刻蚀速率。具体做法是在 COMSOL 里沿靶面径向取 50 个采样点每个点提取 B_parallel、E_perp 和鞘层电势降用鞘层模型算出离子入射角 θ arctan(B_parallel / B_total) 和入射能量 E e × V_sheath把 θ 和 E 代入蒙特卡罗产额函数得到局部产额 Y(r)刻蚀速率 ER(r) Y(r) × Γ_ion(r) × M_target / (ρ_target × N_A)。其中 Γ_ion 是离子通量由 B_parallel × E_perp 标定。这个耦合过程用 Python 脚本串起来最方便。COMSOL 支持 LiveLink for Python可以直接在 Python 里调用 COMSOL 的求解器避免手动导出导入。如果不用 LiveLink就把 COMSOL 结果存成文本用 pandas 读进来做插值。注意插值时要用径向坐标对齐COMSOL 的轴对称模型输出的是 (r, z) 网格蒙特卡罗的产额函数输入是标量角度和能量需要先做二维插值再逐点计算。注意鞘层电势降不是靶面偏压的全部。等离子体电位通常在 10 到 20 V实际离子能量是偏压减去等离子体电位。忽略这一项会让产额偏高 5%–10%。4. 刻蚀形貌的时间演化从速率分布到靶面轮廓4.1 用水平集方法推进靶面演化有了刻蚀速率分布 ER(r)靶面形貌的演化就是一个界面推进问题。初始靶面是平面每个时间步按局部速率沿法向推进推进量 ER(r) × Δt。当刻蚀深度达到毫米量级时靶面会出现明显的刻蚀环和凹坑这时入射角分布会反过来改变需要耦合更新。水平集方法适合处理这种拓扑变化用符号距离函数 φ(r, z, t) 表示靶面φ 0 是界面φ 0 是靶材内部φ 0 是真空。演化方程是 ∂φ/∂t V_n |∇φ| 0V_n 是法向速度由 ER(r) 和局部表面法向决定。在 Python 里可以用 scikit-fmm 做快速行进法重新初始化用有限差分做时间推进。网格取 0.1 mm时间步取 1 小时跑 100 步就能看到刻蚀环的形成。每 10 步重新计算一次入射角分布——因为靶面倾斜后离子入射角不再是初始的 θ而是 θ 局部表面倾角。这个反馈是刻蚀环自锐化的原因环的侧壁越来越陡入射角越来越接近峰值产额角度刻蚀速率进一步加快。4.2 形貌演化的 Python 实现与参数设置下面代码用简化的水平集方法推进靶面重点展示速率耦合和角度更新。import numpy as np from scipy.ndimage import gaussian_filter # 网格 Nr, Nz 200, 100 r np.linspace(0, 0.05, Nr) # 径向 0-50 mm z np.linspace(0, 0.01, Nz) # 轴向 0-10 mm dr, dz r[1]-r[0], z[1]-z[0] # 初始靶面z 5 mm 处平面 phi z[None, :] - 0.005 phi np.tile(phi, (Nr, 1)) # 初始刻蚀速率分布由第 3 章耦合得到 ER 1e-9 * (1 0.8 * np.exp(-((r - 0.025)/0.008)**2)) # m/s ER np.tile(ER[:, None], (1, Nz)) dt 3600 # 时间步 1 小时 for step in range(100): # 计算法向 phi_r, phi_z np.gradient(phi, dr, dz) norm np.sqrt(phi_r**2 phi_z**2) 1e-12 # 局部表面倾角修正入射角 tilt np.arctan2(phi_r, phi_z) # 简化产额随倾角变化峰值在 60 度 yield_factor 1 0.5 * np.exp(-((np.abs(tilt) - np.radians(60))/np.radians(20))**2) Vn ER * yield_factor # 水平集推进 phi - dt * Vn * norm # 每 10 步重新初始化保持符号距离性质 if step % 10 0: phi gaussian_filter(phi, sigma1.0) phi phi / (np.abs(phi).max() 1e-12) * 0.005 # 提取最终靶面轮廓 surface_z np.array([z[np.argmin(np.abs(phi[i, :]))] for i in range(Nr)]) print(刻蚀环最深位置 r , r[np.argmin(surface_z)], m)这段代码的核心是用 φ 的梯度算法向用局部倾角修正产额因子然后按 Vn × dt 推进界面。ER的初始分布是高斯型对应第 3 章算出的离子通量分布。yield_factor模拟了产额随入射角的变化峰值在 60°。每 10 步做一次高斯滤波和归一化防止 φ 的梯度畸变。跑 100 步后靶面会出现一个明显的凹环位置在 r 25 mm 附近深度约 2 mm。这个结果和实际靶材的刻蚀环位置基本一致。参数方面dt不能太大否则界面推进会出现数值振荡一般取实际刻蚀时间的 1/100 到 1/50。sigma控制重新初始化的平滑程度太大会抹平刻蚀环的细节太小会让 φ 失去符号距离性质。我一般从 1.0 开始试看轮廓是否光滑。4.3 刻蚀环自锐化的验证与实验对照刻蚀环自锐化是磁控溅射靶材的典型现象初始平坦的靶面在刻蚀过程中逐渐形成 V 形或 U 形沟槽沟槽侧壁越来越陡。模拟结果如果能看到这个趋势说明耦合逻辑是对的。验证方法是把模拟得到的最终靶面轮廓和实际使用后的靶材剖面做对比看刻蚀环的位置、深度、半高宽是否一致。如果位置对但深度偏浅通常是产额偏低或离子通量标定偏小如果位置偏内或偏外通常是磁场位形或边界条件有问题。我一般会做三组对照一组是纯蒙特卡罗产额不考虑电磁场一组是纯有限元通量不考虑产额角度依赖一组是耦合模型。前两组分别会高估或低估刻蚀环的锐度只有耦合模型能同时匹配位置和深度。这个对照实验能帮你判断误差来源避免在错误的模型上调参数。提示实验对照时靶材剖面的测量精度很关键。用线切割切开靶材后用光学显微镜测轮廓精度能到 10 μm。如果只用卡尺误差可能比刻蚀深度还大。5. 避坑与排查磁控溅射模拟中最容易翻车的五个地方5.1 产额曲线在低能区出现负值或发散现象蒙特卡罗跑出来的产额在 100 eV 以下变成负数或者随能量降低反而增大。原因电子阻止本领的简化公式在低能区失效E_electronic 0.1 * sqrt(E)在 E 很小时趋近于零但核阻止项可能超过总能量导致 E 变成负数。解决给能量损失加一个下限当 E 10 eV 时直接终止碰撞级联或者换用更精确的 Lindhard-Scharff 电子阻止公式它在低能区有正确的渐近行为。5.2 COMSOL 磁场求解不收敛或结果明显偏离现象稳态求解器报错「奇异矩阵」或「不收敛」或者算出的 B_parallel 在靶面中心出现异常峰值。原因永磁体的剩余磁通密度方向设错了或者网格在永磁体边缘太粗导致磁矢势的旋度计算失真。解决检查永磁体的磁化方向是否沿轴向在永磁体边缘加边界层网格把求解器换成 MUMPS 并开启自适应网格细化。如果还是不行先用一个简单的圆柱永磁体做验证确认磁场分布符合解析解再改几何。5.3 蒙特卡罗与有限元耦合时坐标对不上现象刻蚀速率分布和离子通量分布错位刻蚀环位置偏了半个靶半径。原因COMSOL 的轴对称模型输出的是 (r, z) 坐标蒙特卡罗的产额函数输入是入射角但入射角的计算依赖 B_parallel 和 B_total 的比值如果 B_total 在靶面边缘被截断角度会算错。解决在 COMSOL 里提取 B_parallel 时同时提取 B_total确保比值在 0 到 1 之间在 Python 里做插值时用径向坐标 r 对齐不要用网格索引对齐。5.4 水平集推进时界面出现振荡或破碎现象靶面轮廓在几个时间步后出现锯齿状振荡或者 φ 场出现多个零交叉。原因时间步太大或者法向计算时梯度被噪声放大。解决把 dt 减小到原来的 1/5在每次推进前对 φ 做一次高斯滤波用 scikit-fmm 的重新初始化函数代替手动归一化。如果还不行改用有限差分法直接推进表面高度虽然不能处理拓扑变化但稳定性好得多。5.5 模拟结果和实验对不上却找不到原因现象刻蚀环位置对但深度差一倍或者深度对但宽度差很多。原因离子通量的绝对标定没有做只用了相对分布。蒙特卡罗产额是绝对量但有限元算出的 B_parallel × E_perp 是相对量需要用一个标定系数把它转成绝对离子通量。解决用一组已知实验数据反推标定系数——比如已知靶材在 500 W、0.5 Pa 下跑了 100 小时刻蚀环深度 2 mm用这个深度反推 Γ_ion 的绝对值。标定一次之后同一台设备的其他功率和气压条件可以直接用。6. 把模拟变成工艺窗口参数扫描与快速预测的技巧走到这一步你已经有了产额曲线、电磁场分布、刻蚀形貌演化三块拼图。但真正让这套方案值钱的不是单次模拟而是用它做参数扫描找出工艺窗口。我一般会固定靶材和磁场位形扫三个参数功率、气压、靶基距。功率影响离子能量和通量气压影响碰撞平均自由程和离子角度分布靶基距影响膜厚均匀性。每个参数取 5 个水平用拉丁超立方抽样选 20 个组合跑一遍耦合模拟得到每个组合的刻蚀环深度和膜厚均匀性指标。这里有个技巧不要每次都用完整的蒙特卡罗 有限元 水平集。先用有限元算磁场磁场不随功率和气压变存成查找表再用蒙特卡罗算产额曲线产额只随能量和角度变也存成查找表最后用水平集推进时只查表插值不重新跑物理模型。这样单次形貌演化的时间从几小时降到几分钟20 个组合一晚上就能跑完。另一个技巧是降维。刻蚀环的位置主要由磁场位形决定和功率、气压关系不大刻蚀环的深度主要由离子通量和产额决定和功率强相关。所以可以先固定磁场扫功率和气压得到深度图再单独优化磁场位形来调位置。这样把三维扫描拆成两个二维扫描计算量降一个数量级。验证方法上我习惯用「留一法」拿 20 个组合中的 19 个做训练1 个做验证看预测的刻蚀深度和模拟值差多少。如果误差在 10% 以内说明查找表插值够用如果超过 20%说明某个参数的非线性太强需要加密采样。这个习惯帮我省了很多次盲目调参的时间——有一次气压从 0.3 Pa 变到 0.5 Pa刻蚀深度预测值跳了 30%查表才发现是产额曲线在 200 eV 附近有个拐点插值没抓住。后来在拐点附近补了 5 个采样点误差就降到 8% 了。最后说一个我踩过的坑别在模拟里追求完美。蒙特卡罗的统计噪声、有限元的网格误差、水平集的数值耗散三者叠加后模拟精度能到 15% 以内就算不错了。与其花一周把误差从 15% 降到 10%不如用这周跑 50 个工艺组合找出趋势和边界。工艺窗口的价值在于告诉你「哪个方向不能走」而不是「这个点精确是多少」。希望帮到你。本文还有配套的精品资源点击获取