搞电磁仿真的朋友应该都遇到过这个尴尬在CST或者HFSS里建模一个石墨烯器件翻遍材料库找不到石墨烯这个选项。即便你手动填一个三维材料的介电常数算出来的结果往往也是离谱的。石墨烯是单原子层的二维材料厚度只有一个碳原子级别正确的建模方式是在边界条件里给它赋一个复表面电导率或者等效的表面阻抗。而这个表面电导率最权威的计算路径就是Kubo公式——它从量子线性响应理论出发把带内intraband和带间interband两种电子跃迁机制统一到一个表达式里既能描述太赫兹频段的Drude型输运也能覆盖红外和光频段的量子跃迁吸收。这篇文章把整个计算流程用MATLAB完整实现一遍包括公式怎么拆、奇异点怎么处理、单位怎么统一、参数怎么扫描以及我实际踩过的几个坑。全文给出可以直接复制的完整代码同时对每个关键步骤解释为什么这么做让不同基础的读者都能拿起来用。适合正在做石墨烯太赫兹调制器、超表面吸波器、石墨烯等离激元波导或者任何需要把二维材料放进电磁仿真软件里的朋友参考。1. 为什么算石墨烯电导率要上Kubo公式1.1 石墨烯器件仿真的痛点石墨烯本质上是一个二维电子气系统载流子被限制在单原子层平面内运动。常规电磁仿真软件里给的材料参数都是三维体材料介电常数ε、体电导率σ的单位是S/m但石墨烯的厚度根本没有明确定义——你总不能拿0.334纳米去当真厚度填进去。有文献喜欢给石墨烯设一个等效厚度比如1纳米然后用体电导率σ_3Dσ_2D/d来近似这在某些频段误差还能接受但在红外以上频段由于带间跃迁的出现等效厚度法就很难还原真实的频散关系。更常见的做法是用边界条件把石墨烯建模成一个无限薄的导电面用复表面电导率σ_s单位S即西门子而不是S/m作为输入。CST里可以用Impedance Sheet边界HFSS里可以用Thin Layer或Lumped RLC等形式。这个表面电导率怎么给就成为仿真精度的决定性因素。如果你只是随手填一个常数电导率那么所有关于频率色散、费米能级调控、温度依赖的物理就全丢了。1.2 Kubo公式相比Drude模型强在哪很多人一开始会想到Drude模型毕竟它简单σ σ_DC/(1-iωτ)。Drude模型在描述金属薄膜、掺杂半导体时非常好用但在石墨烯的场景下至少有两个问题。第一Drude模型只考虑自由载流子的带内散射过程。当光子能量ℏω超过2E_FE_F为费米能级时价带电子可以直接跃迁到导带这个带间吸收通道会突然打开在电导率谱上表现为一个明显的台阶色散结构。Drude模型完全抓不到这个特征。第二石墨烯的能带是线性的狄拉克锥态密度随能量线性变化这跟传统抛物带半导体完全不同。所以带内项的权重因子那个2k_BT ln[2cosh(E_F/2k_BT)]系数也不能简单用一个固定的载流子浓度来套。Kubo公式从线性响应理论出发自然地把这些量子效应都包含进来了。从应用角度说太赫兹波段0.1~10 THz正好是石墨烯表面等离激元的强束缚频段也是带内项主导的区域而中红外到可见光带间跃迁对吸收的贡献变得关键。如果你要设计一个从太赫兹到红外的宽频可调器件只用Drude模型会漏掉一半的物理。1.3 Kubo公式的适用边界Kubo公式也不是万能的。它建立在单层石墨烯、无衬底或弱衬底耦合的假设上。如果石墨烯是嵌入在多层异质结构里跟高介电常数衬底强耦合或者存在显著的局域场增强效应那么还需要考虑衬底介电环境对库仑相互作用的屏蔽修正比如动态介电常数ε(q,ω)进入积分。另外Kubo公式本身是局域近似准经典RPA线性响应对于波长远大于石墨烯晶格常数的情况太赫兹到近红外都满足这个近似是足够精确的。做仿真时还要注意Kubo公式默认石墨烯是均匀、无缺陷的理想体系实际样品的掺杂不均匀性、皱褶、晶界散射都会让测量值和理论值有偏差。所以仿真结果更适合用来做趋势预测和相对对比而不是精确预测某一颗具体器件的绝对性能。2. 公式拆解带内、带间与表面阻抗2.1 带内项的物理意义与解析式石墨烯电导率拆成两项σ σ_intra σ_inter带内项描述的是导带内或价带内电子在外场驱动下的运动物理上对应载流子的漂移-扩散输运频谱特征就是典型的Drude色散。完整的Kubo带内项表达式为σ_intra -i (e²k_B T)/(πℏ²(ω - i/τ)) × ln(2cosh(E_F/(2k_BT)))这里e是电子电荷k_B是玻尔兹曼常数ℏ是约化普朗克常数τ是散射时间E_F是费米能级化学势T是温度。注意分母里的虚部符号跟时间因子约定有关我这里用的是e^{-iωt}约定。对数项ln(2cosh(x))这个形式很有意思。在低温极限k_BT E_F下它可以化简ln(2cosh(E_F/2k_BT)) ≈ E_F/(2k_BT) ln(1e^{-E_F/k_BT})于是带内项趋近于σ_intra ≈ e²E_F/(πℏ²(ω - i/τ)) × (-i) ≈ e²E_Fτ/(πℏ²(1-iωτ))这就是标准的Drude模型DC电导率σ_DC e²E_Fτ/πℏ²。随便代入几个数E_F0.2eV、τ0.5psσ_DC≈11.8mS对应表面阻抗约85Ω/sq跟实验上测得的高质量石墨烯直流方块电阻量级完全一致。这个验证很重要说明公式和代码没写错。在室温300K下E_F0.2eV时k_BT≈0.0259eVE_F/(2k_BT)≈3.86ln(2cosh3.86)≈3.61跟低温极限的3.86有一定差距说明室温热激发对带内电导率有不可忽略的修正。这也解释了为什么在做室温柔性电子器件仿真时不要随便用低温近似公式。2.2 带间项的主值积分带间项来自价带电子吸收光子跃迁到导带的过程。它的表达式比较复杂σ_inter e²/(4ℏ) × [ G(ℏω/2) i(4ℏω/π) ∫₀^∞ (G(ε) - G(ℏω/2))/((ℏω)² - 4ε²) dε ]辅助函数G(ε)定义为G(ε) sinh(ε/k_BT) / (cosh(E_F/k_BT) cosh(ε/k_BT))这个G函数可以理解为有限温度下的跃迁相位空间因子。在T→0极限下G(ε)会退化为一个台阶函数当ε E_F时G0当ε E_F时G1。这时带间项可以解析积分出来σ_inter e²/(4ℏ) × [ θ(ℏω-2E_F) (i/π) ln |(ℏω2E_F)/(ℏω-2E_F)| ]θ是单位阶跃函数。这个公式是大家最常用的低温解析近似非常方便手算验证代码。它清楚说明了一个物理事实当ℏω 2E_F时带间吸收被Pauli阻塞抑制实部很小当ℏω 2E_F时实部突变为e²/4ℏ≈60.8μS这是石墨烯在红外波段的统一透明率对应的电导率量子值。但做室温仿真时T300K下热展宽会把这条突变边抹平成平滑过渡G函数不再是一个简单的台阶所以还是老老实实用完整Kubo积分更稳妥。积分中有一个关键难点当ε ℏω/2时分子分母同时为零被积函数存在一个有限极限。用洛必达法则可以求出极限值lim -G(ℏω/2)/(4ℏω)这个极限在处理数值积分时必须显式处理否则integral函数会在奇异点处报错或给出粗糙结果。2.3 表面阻抗换算与量纲陷阱表面电导率和表面阻抗的关系非常简单Z_s 1/σ这是二维导电面的定义面电流密度J_s单位A/m跟切向电场E_t单位V/m之间满足J_s σE_t所以σ的单位是A/V S西门子Z_s的单位是欧姆。但为了跟薄膜方块电阻区分一般写成Ω/sq。这个sq是square的意思不是单位平方很多新手会以为它跟面积有关其实它表示的是无论取多大面积的正方形面电阻都一样这是二维导电膜的本征属性。跟三维体电导率换算时最常见的错误就是单位搞混。如果你在CST里用体电导率建石墨烯模型要把二维电导率除以等效厚度。比如σ_2D0.1mS等效厚度取1nm则体电导率应为σ_3D0.1mS/1nm1e5 S/m。这个数值看起来很大很多人觉得不靠谱实际上是对的——因为石墨烯把电流压缩在原子厚度里体电导率必然异常高。表面阻抗的实部对应损耗电阻虚部则反映了电感的性质。在太赫兹低频端石墨烯表面阻抗虚部通常为负容性还是感性实部随频率下降。随着频率逼近带间跃迁阈值表面阻抗谱会出现明显拐点。这些特征在器件设计时直接决定了匹配网络的形式。3. MATLAB实现完整代码与逐段讲解3.1 参数设置与物理常量先定义物理常数全部采用国际单位制SI。我见过太多人代码跑出来差十几个数量级最后发现是把eV和J搞混了。%% 物理常数SI e_charge 1.602176634e-19; % 元电荷 [C] hbar 1.054571817e-34; % 约化普朗克常数 [J*s] kB 1.380649e-23; % 玻尔兹曼常数 [J/K] %% 模型参数 T 300; % 温度 [K] E_F 0.2 * e_charge; % 费米能级 [J]即 0.2 eV tau 0.5e-12; % 散射时间 [s]0.5 ps费米能级用E_F 0.2×e_charge是因为我们习惯用eV表达能级但公式里所有能量都是焦耳。换算关系1 eV 1.602e-19 J。如果你从文献里看到E_F0.2eV一定要乘e_charge再代入公式否则结果会差1.6e-19倍。温度参数T默认300K对应室温。如果是低温测量模拟可以设成77K或4.2K但要注意温度越低带间跃迁边越陡峭数值积分对奇异点处理的要求越高delta截断就要选得更小。3.2 频率向量与带内项实现频率范围根据你的目标频段设置。这里覆盖0.1~10THz基本覆盖了整个太赫兹窗口。采样点数Nf设500个够画平滑曲线计算速度快结果也不受积分误差影响。f linspace(0.1e12, 10e12, 500); % 频率 [Hz] omega 2*pi*f; % 角频率 [rad/s] % 带内电导率Drude-like ln_term abs(E_F)/(2*kB*T) log(1 exp(-abs(E_F)/(kB*T))); sigma_intra -1i * e_charge^2 * kB * T ./ (pi * hbar^2 * (omega - 1i/tau)) .* ln_term;注意我用了abs(E_F)的写法。这个对数项本质上只依赖|E_F|因为石墨烯电子-空穴对称P型掺杂和N型掺杂的带内输运特性相同。如果你直接写成E_F当E_F为负值时可能会在log里产生假想成分这是新手常犯的错误。带内项的向量化写法省去了循环。MATLAB中直接用点除./和点乘.*速度比for循环快两个数量级。这里没有奇异点问题因为分母(ω - i/τ)在实轴ω上永远不为零τ有限时也不会遇到奇异性。3.3 带间项数值积分奇异点处理是关键带间项是整段代码最需要小心的部分。直接照抄积分表达式跑一遍integral函数会在ε→ℏω/2的位置警告Minimum step size reached甚至返回NaN。%% 辅助函数 G(ε) Gfun (eps) sinh(eps/(kB*T)) ./ (cosh(E_F/(kB*T)) cosh(eps/(kB*T))); %% 带间电导率 sigma_inter zeros(size(omega)); for idx 1:length(omega) w omega(idx); Ep hbar * w; % 光子能量 E_h Ep/2; % ε₀ ℏω/2 Ghalf Gfun(E_h); integrand (eps) (Gfun(eps) - Ghalf) ./ (Ep^2 - 4*eps.^2); % 奇异点邻域截断 delta max(E_h*1e-9, 1e-30); % 左区间 [0, E_h - delta] if E_h - delta 0 I_left integral(integrand, 0, E_h - delta, RelTol, 1e-9, AbsTol, 1e-13); else I_left 0; end % 右区间 [E_h delta, Emax] Emax max([20*abs(E_F), 20*kB*T, 10*Ep]); % 截断上限 I_right integral(integrand, E_h delta, Emax, RelTol, 1e-9, AbsTol, 1e-13); PV I_left I_right; sigma_inter(idx) (e_charge^2/(4*hbar)) * (Ghalf 1i*(4*Ep/pi)*PV); end关于奇异点这里有个细节值得展开。被积函数在εE_h处分子分母同时为零极限值有限但数值积分算法在奇点附近会不断细分区间导致速度变慢甚至失败。我的处理方式是直接挖掉一个delta邻域。delta太小积分结果对delta敏感delta太大会牺牲准确性。实测下来取E_h的1e-9倍是一个平衡点结果对delta的依赖低于积分容差。还有一点值得注意积分上限直接取Inf会让integral函数在无穷区间的尾部反复自适应有时会告警不收敛。更稳的做法是像我这样取一个物理上足够大的截断值Emax。当ε远远超过max(E_F, k_BT, E_h)时G(ε)趋近于1被积函数近似为(1-Ghalf)/(Ep²-4ε²)在ε→∞时按1/ε²衰减所以截断到几十倍特征能量之后剩余尾部积分可以忽略。如果你需要更高精度可以把奇异点邻域单独提出来用解析极限替换。被积函数在εE_h附近的极限为-G(E_h)/(4Ep)而G(ε)可以解析求出G(ε) [cosh(x)cosh(a)1] / (k_BT × (cosh(a)cosh(x))²)其中xε/(k_BT)aE_F/(k_BT)。在编码中可以把这个极限值赋给奇异点处的被积函数值然后integral函数在含内部奇异点的区间上也能正常工作。3.4 可视化与结果验证总电导率把两项加起来然后直接取倒数得到表面阻抗sigma_total sigma_intra sigma_inter; Z_s 1 ./ sigma_total; figure; subplot(1,2,1); plot(f/1e12, real(sigma_total), b-, LineWidth, 1.5); hold on; plot(f/1e12, imag(sigma_total), r--, LineWidth, 1.5); grid on; xlabel(频率 (THz)); ylabel(电导率 (S)); legend(Re σ,Im σ,Location,best); title(石墨烯复电导率 (E_F0.2eV, T300K)); subplot(1,2,2); plot(f/1e12, real(Z_s), b-, LineWidth, 1.5); hold on; plot(f/1e12, imag(Z_s), r--, LineWidth, 1.5); grid on; xlabel(频率 (THz)); ylabel(表面阻抗 (Ω/sq)); legend(Re Z_s,Im Z_s,Location,best); title(石墨烯表面阻抗);跑完代码后先别急着做其他事情。我强烈建议你做两个验证一是DC验证。在f0.1THz附近近似静态极限看σ_total的实部是否接近理论值。E_F0.2eV、τ0.5ps、T300K时σ_intra在低频段实部应该在几十mS量级表面阻抗实部在几十到上百Ω/sq量级。二是带间边验证。E_F0.2eV时带间跃迁阈值在ℏω0.4eV对应频率f2E_F/h≈96.8THz。所以在我们0.1~10THz的频率窗口内理论上看不到带间吸收的突变。但你可以把最高频率扩展到200THz扫一扫会在约96.8THz处看到σ_inter实部突然上升。如果这个特征位置不对说明能量单位换算或E_F设置有问题。4. 参数扫描物理规律与设计参考4.1 费米能级如何调控表面阻抗费米能级是石墨烯器件最关键的调控旋钮因为在仿真中你可以通过改变栅极电压来连续调控E_F。我把E_F分别设为0.1eV、0.2eV、0.5eV其他参数固定T300Kτ0.5ps扫描频率0.1~10THz。结果非常直观E_F越高带内电导率实部越大表面阻抗越小。在1THz处E_F0.1eV时表面电阻可能接近几百Ω/sqE_F0.5eV时可以降到十几Ω/sq。这跟DC方块电阻的E_F成正比关系是一致的。更关键的是虚部。低频段石墨烯表面阻抗虚部为负呈电容性还是感性这个需要按你用的时间约定来判断。在我的符号约定e^{-iωt}下表面阻抗虚部为负对应容抗表现为电场能量存储。随着频率升高虚部的绝对值先增大后减小。这说明石墨烯在太赫兹频段的行为像一个损失较大的RC分布式电路而不是简单的电阻层。做天线匹配时这个虚部必须纳入考量否则你设计的阻抗匹配点会整体偏移。我还发现一个规律E_F越高带间阈值频率越高。阈值频率f_inter 2E_F/hE_F0.5eV时阈值在241THz已经远超出太赫兹窗口。所以在太赫兹频段做设计时如果费米能级调得比较高带间项对电导率的贡献确实可以忽略Drude近似就够用。但如果你做的是中红外30~50THz传感器那么带间项就必须保留了。4.2 温度与散射时间的作用把温度从10K扫到500K观察电导率谱的变化会发现两个效应。第一个是带内项的幅度整体变化温度越高ln(2cosh(E_F/2k_BT))越偏离低温极限因为热激发让更多载流子参与输运DC电导率轻微上升。第二个是带间跃迁边被展宽低温时跃迁边非常锐利接近阶跃函数室温300K下这个边被热展宽抹成了平滑过渡500K时则进一步降低跃迁边的陡峭度。这种现象对红外探测器的光谱选择性设计有直接影响——低温器件更容易实现锐利的吸收边。散射时间τ对电导率的影响集中在带内项。τ从0.1ps到1ps变化时DC电导率线性增大σ_DC∝τ但频谱的Drude拐点频率f_D 1/(2πτ)会左移。τ0.1ps时拐点在1.6THzτ1ps时拐点降到160GHz。如果你的工作频段在几个THzτ1ps的样品在太赫兹高端已经处于欧姆损耗主导的区域表面阻抗实部随频率上升明显而τ0.1ps时整个太赫兹窗口仍在惯性主导区域电导率随频率变化较平缓。这个趋势对判断器件是损耗受限还是惯量受限很重要。散射时间的取值很有讲究。高迁移率悬浮石墨烯的τ可以到几皮秒而跟衬底强耦合的CVD石墨烯通常只有0.1~0.3ps。做仿真时尽量根据自己样品的霍尔迁移率来估算τ μE_F/(e v_F²)其中v_F≈1e6m/s。迁移率10000cm²/(V·s)的样品E_F0.2eV时τ≈0.63ps迁移率5000cm²/(V·s)时τ≈0.31ps。不要闭着眼睛填τ。4.3 从复电导率到等效介电常数有些仿真场景不提供阻抗边界条件选项比如某些时域有限差分软件这时可以把石墨烯等效为厚度为d的介质层用等效相对介电常数来建模ε_eq(ω) 1 iσ(ω)/(ωε₀d)推导过程很简单把石墨烯看作厚度d的介质面电流J_s σE体电流J_3D σ_3D_E而σ_3D σ/d。另一方面介质中的电流J σ_3D E -iωε₀(ε_r-1)E所以上面那个公式就出现了。实现时取d1nm或者0.5nm取决于你想把多少额外厚度强加给模型。在太赫兹频段ε_eq的虚部会非常大因为σ大、频率低这会让时域仿真中的Courant稳定性条件变得苛刻。如果你用FDTD做石墨烯器件仿真强烈建议还是用阻抗边界条件而不是等效介质层否则网格需要做得很细计算资源会被吃掉一大块。等效介质法有一个坑厚度d是一个人为拟合参数换一个d等效介电常数就变了。不同文献里d取值从0.3nm到10nm都有这导致了文献间结果对不齐的乱象。如果你需要跟某篇文献对比仿真数据一定先看它用的d是多少。5. 常见问题与排查实录5.1 积分告警与不收敛最常遇到的报错是integral函数提示Maximum function count exceeded或者Minimum step size reached通常出现在频率接近带间阈值、温度较低、或者奇异点delta取太小的时候。我总结了一套排查顺序先看奇异点delta。如果delta取得过小比如E_h的1e-12倍左右两个积分区间几乎贴着奇异点自适应积分会在边界处疯狂细化导致函数调用次数爆炸。调到E_h的1e-8~1e-9相对稳妥。再看截断上限Emax。如果Emax取得不够大比如只取5kBT而E_F远大于k_BT那么G(ε)在截断点附近还没有完全饱和到1被积函数尾部截断误差会传进主值积分。建议Emax至少取max(20|E_F|, 20k_BT, 10Ep)。最后看频率采样。当你扫频时如果带间阈值恰好落在你的频率网格中间你会看到电导率实部出现异常的毛刺这不是物理问题而是数值分辨率不够。加密采样点或者把频率扫描设成对数间隔可以改善。5.2 低频段表面阻抗虚部跳变我在最初写代码时发现σ_total的实部在非常低的频率比如0.01THz以下会轻微震荡虚部则乱七八糟。排查后发现是带间积分在低频极限下出了问题。当ω→0时E_h→0奇异点ε0恰好落在积分左端点。此时被积函数在ε→0时本身有定义可以用洛必达法则算但integral函数不知道这个内部奇异点来自端点自适应过程容易被误导。解决方法是把左区间积分拆成[0, 某个极小值]和[极小值, E_h-delta]两段或者干脆在低频段比如ℏω0.01k_BT用解析的低频展开式代替数值积分。实际上低频段远小于带间阈值带间项贡献本身就很小实部趋近于零虚部有一个有限的负值。很多应用场景中在ω远小于2E_F/ℏ时直接忽略带间项也不会差太多。我通常设一个判断如果f 0.1×2E_F/h就跳过带间积分的精细计算只保留带内项。这样既加速又避免了奇异积分。5.3 量纲与单位自查清单这个领域最常见的翻车现场就是单位制混乱。我列一个自查清单每次跑代码前过一遍E_F和k_BT必须同为焦耳或同为eV。如果E_F用eVk_BT也用eV300K对应约0.0259eV那么ln项无量纲没问题。但交给MATLAB时e_charge²、hbar²、kB这些常数全是SI单位你混用eV就会在前面系数上出错。ℏ和h别搞混。Kubo公式里用的是约化普朗克常数ℏ不是普朗克常数h。两者差2π倍误差巨大。类似地频率转能量时用Ehf但如果是角频率ω则Eℏω。我的代码里用的是omega2πf所以光子能量Ephbar*omega没错。电导率单位是S不是S/m。如果你从某篇论文里看到σ1e-3单位标注S这是二维电导率可以直接代入Z_s1/σ。如果看到S/m那是三维体电导率必须搞清楚他们用的等效厚度是多少才能转换。表面阻抗的虚部符号依赖时间约定。CST、HFSS、FDTD方法、频域求解器各自的约定可能不同有的用e^{iωt}有的用e^{-iωt}对应的电导率虚部符号会翻转。在把MATLAB计算出的σ赋给仿真软件前务必确认软件的约定。通常CST的Drude模型选项会带你设置参数的imaginary part sign别默认瞎填。5.4 性能优化与代码加速上面的代码在Nf500时大概几秒钟跑完完全够用。如果你要做大规模参数扫描比如E_F扫100个点、τ扫50个点500×100×50250万个频率点多次积分就不太妙了。这时候有几个提速手段。第一个是并行循环。MATLAB的parfor可以直接把for loop改成并行前提是装有Parallel Computing Toolbox。带间项积分之间相互独立天然适合并行。第二个是减少重复计算。Gfun里cosh(E_F/(kB*T))是一个常数可以提到循环外面算一次。同理带内项的ln_term也是一个常量向量没必要在每次频率循环里重复调用。第三个是缓存积分结果或使用查表插值。如果不要求每一点的绝对精度可以先用稀疏频率网格比如50个点算完积分再用spline插值得到500个点的值误差通常能控制在0.5%以内。这在设计迭代中绰绰有余速度能提升一个数量级以上。第四个如果你用的是低温解析近似θ对数形式带间项根本不需要积分整个σ可以写成全解析表达式向量化直接跑速度最快。但记得只在T远小于E_F/k_B时使用否则结果跟完整Kubo对不齐。6. 写在最后的实操体会代码分享到这里最后说几点我自己的体会。第一石墨烯Kubo公式虽然看起来像一团数学符号但拆开之后物理非常清晰带内项就是载流子的摩擦输运带间项就是光生电子-空穴对的量子跃迁两者在频域上互补低频靠带内高频靠带间中间的过渡区域由费米能级和温度共同决定。理解了这个图像你在调试代码时至少不会走错方向。第二无论数值结果怎样一定要先过一个常识检查DC电导率是否跟E_Fτ成正比表面阻抗是否在合理的几十到几百Ω/sq范围带间吸收边是否落在2E_F的位置这三个检查能在五分钟内帮你抓住90%的错误。第三真实项目中石墨烯的射频/微波特性非常依赖衬底。SiO₂衬底上的石墨烯跟悬浮石墨烯的散射时间可以差一个数量级从而让表面阻抗谱完全不同。仿真时把τ当可调参数跟实验S参数做拟合反推样品实际的τ和E_F这是一个非常实用的标定思路。碰到测不准的样品参数就用这个办法来校准。如果用这套代码复现出问题建议先跑一遍T300K、E_F0.2eV、τ0.5ps的默认参数把得到的结果跟我上面描述的量级对照一下。确认基线正确后再去扫描其他参数这样就不会一头扎进数值迷雾里出不来。祝各位在石墨烯器件仿真的路上一次跑通不用像我当初那样在奇点积分里泡了整整两天。
搞电磁仿真的朋友应该都遇到过这个尴尬在CST或者HFSS里建模一个石墨烯器件翻遍材料库找不到石墨烯这个选项。即便你手动填一个三维材料的介电常数算出来的结果往往也是离谱的。石墨烯是单原子层的二维材料厚度只有一个碳原子级别正确的建模方式是在边界条件里给它赋一个复表面电导率或者等效的表面阻抗。而这个表面电导率最权威的计算路径就是Kubo公式——它从量子线性响应理论出发把带内intraband和带间interband两种电子跃迁机制统一到一个表达式里既能描述太赫兹频段的Drude型输运也能覆盖红外和光频段的量子跃迁吸收。这篇文章把整个计算流程用MATLAB完整实现一遍包括公式怎么拆、奇异点怎么处理、单位怎么统一、参数怎么扫描以及我实际踩过的几个坑。全文给出可以直接复制的完整代码同时对每个关键步骤解释为什么这么做让不同基础的读者都能拿起来用。适合正在做石墨烯太赫兹调制器、超表面吸波器、石墨烯等离激元波导或者任何需要把二维材料放进电磁仿真软件里的朋友参考。1. 为什么算石墨烯电导率要上Kubo公式1.1 石墨烯器件仿真的痛点石墨烯本质上是一个二维电子气系统载流子被限制在单原子层平面内运动。常规电磁仿真软件里给的材料参数都是三维体材料介电常数ε、体电导率σ的单位是S/m但石墨烯的厚度根本没有明确定义——你总不能拿0.334纳米去当真厚度填进去。有文献喜欢给石墨烯设一个等效厚度比如1纳米然后用体电导率σ_3Dσ_2D/d来近似这在某些频段误差还能接受但在红外以上频段由于带间跃迁的出现等效厚度法就很难还原真实的频散关系。更常见的做法是用边界条件把石墨烯建模成一个无限薄的导电面用复表面电导率σ_s单位S即西门子而不是S/m作为输入。CST里可以用Impedance Sheet边界HFSS里可以用Thin Layer或Lumped RLC等形式。这个表面电导率怎么给就成为仿真精度的决定性因素。如果你只是随手填一个常数电导率那么所有关于频率色散、费米能级调控、温度依赖的物理就全丢了。1.2 Kubo公式相比Drude模型强在哪很多人一开始会想到Drude模型毕竟它简单σ σ_DC/(1-iωτ)。Drude模型在描述金属薄膜、掺杂半导体时非常好用但在石墨烯的场景下至少有两个问题。第一Drude模型只考虑自由载流子的带内散射过程。当光子能量ℏω超过2E_FE_F为费米能级时价带电子可以直接跃迁到导带这个带间吸收通道会突然打开在电导率谱上表现为一个明显的台阶色散结构。Drude模型完全抓不到这个特征。第二石墨烯的能带是线性的狄拉克锥态密度随能量线性变化这跟传统抛物带半导体完全不同。所以带内项的权重因子那个2k_BT ln[2cosh(E_F/2k_BT)]系数也不能简单用一个固定的载流子浓度来套。Kubo公式从线性响应理论出发自然地把这些量子效应都包含进来了。从应用角度说太赫兹波段0.1~10 THz正好是石墨烯表面等离激元的强束缚频段也是带内项主导的区域而中红外到可见光带间跃迁对吸收的贡献变得关键。如果你要设计一个从太赫兹到红外的宽频可调器件只用Drude模型会漏掉一半的物理。1.3 Kubo公式的适用边界Kubo公式也不是万能的。它建立在单层石墨烯、无衬底或弱衬底耦合的假设上。如果石墨烯是嵌入在多层异质结构里跟高介电常数衬底强耦合或者存在显著的局域场增强效应那么还需要考虑衬底介电环境对库仑相互作用的屏蔽修正比如动态介电常数ε(q,ω)进入积分。另外Kubo公式本身是局域近似准经典RPA线性响应对于波长远大于石墨烯晶格常数的情况太赫兹到近红外都满足这个近似是足够精确的。做仿真时还要注意Kubo公式默认石墨烯是均匀、无缺陷的理想体系实际样品的掺杂不均匀性、皱褶、晶界散射都会让测量值和理论值有偏差。所以仿真结果更适合用来做趋势预测和相对对比而不是精确预测某一颗具体器件的绝对性能。2. 公式拆解带内、带间与表面阻抗2.1 带内项的物理意义与解析式石墨烯电导率拆成两项σ σ_intra σ_inter带内项描述的是导带内或价带内电子在外场驱动下的运动物理上对应载流子的漂移-扩散输运频谱特征就是典型的Drude色散。完整的Kubo带内项表达式为σ_intra -i (e²k_B T)/(πℏ²(ω - i/τ)) × ln(2cosh(E_F/(2k_BT)))这里e是电子电荷k_B是玻尔兹曼常数ℏ是约化普朗克常数τ是散射时间E_F是费米能级化学势T是温度。注意分母里的虚部符号跟时间因子约定有关我这里用的是e^{-iωt}约定。对数项ln(2cosh(x))这个形式很有意思。在低温极限k_BT E_F下它可以化简ln(2cosh(E_F/2k_BT)) ≈ E_F/(2k_BT) ln(1e^{-E_F/k_BT})于是带内项趋近于σ_intra ≈ e²E_F/(πℏ²(ω - i/τ)) × (-i) ≈ e²E_Fτ/(πℏ²(1-iωτ))这就是标准的Drude模型DC电导率σ_DC e²E_Fτ/πℏ²。随便代入几个数E_F0.2eV、τ0.5psσ_DC≈11.8mS对应表面阻抗约85Ω/sq跟实验上测得的高质量石墨烯直流方块电阻量级完全一致。这个验证很重要说明公式和代码没写错。在室温300K下E_F0.2eV时k_BT≈0.0259eVE_F/(2k_BT)≈3.86ln(2cosh3.86)≈3.61跟低温极限的3.86有一定差距说明室温热激发对带内电导率有不可忽略的修正。这也解释了为什么在做室温柔性电子器件仿真时不要随便用低温近似公式。2.2 带间项的主值积分带间项来自价带电子吸收光子跃迁到导带的过程。它的表达式比较复杂σ_inter e²/(4ℏ) × [ G(ℏω/2) i(4ℏω/π) ∫₀^∞ (G(ε) - G(ℏω/2))/((ℏω)² - 4ε²) dε ]辅助函数G(ε)定义为G(ε) sinh(ε/k_BT) / (cosh(E_F/k_BT) cosh(ε/k_BT))这个G函数可以理解为有限温度下的跃迁相位空间因子。在T→0极限下G(ε)会退化为一个台阶函数当ε E_F时G0当ε E_F时G1。这时带间项可以解析积分出来σ_inter e²/(4ℏ) × [ θ(ℏω-2E_F) (i/π) ln |(ℏω2E_F)/(ℏω-2E_F)| ]θ是单位阶跃函数。这个公式是大家最常用的低温解析近似非常方便手算验证代码。它清楚说明了一个物理事实当ℏω 2E_F时带间吸收被Pauli阻塞抑制实部很小当ℏω 2E_F时实部突变为e²/4ℏ≈60.8μS这是石墨烯在红外波段的统一透明率对应的电导率量子值。但做室温仿真时T300K下热展宽会把这条突变边抹平成平滑过渡G函数不再是一个简单的台阶所以还是老老实实用完整Kubo积分更稳妥。积分中有一个关键难点当ε ℏω/2时分子分母同时为零被积函数存在一个有限极限。用洛必达法则可以求出极限值lim -G(ℏω/2)/(4ℏω)这个极限在处理数值积分时必须显式处理否则integral函数会在奇异点处报错或给出粗糙结果。2.3 表面阻抗换算与量纲陷阱表面电导率和表面阻抗的关系非常简单Z_s 1/σ这是二维导电面的定义面电流密度J_s单位A/m跟切向电场E_t单位V/m之间满足J_s σE_t所以σ的单位是A/V S西门子Z_s的单位是欧姆。但为了跟薄膜方块电阻区分一般写成Ω/sq。这个sq是square的意思不是单位平方很多新手会以为它跟面积有关其实它表示的是无论取多大面积的正方形面电阻都一样这是二维导电膜的本征属性。跟三维体电导率换算时最常见的错误就是单位搞混。如果你在CST里用体电导率建石墨烯模型要把二维电导率除以等效厚度。比如σ_2D0.1mS等效厚度取1nm则体电导率应为σ_3D0.1mS/1nm1e5 S/m。这个数值看起来很大很多人觉得不靠谱实际上是对的——因为石墨烯把电流压缩在原子厚度里体电导率必然异常高。表面阻抗的实部对应损耗电阻虚部则反映了电感的性质。在太赫兹低频端石墨烯表面阻抗虚部通常为负容性还是感性实部随频率下降。随着频率逼近带间跃迁阈值表面阻抗谱会出现明显拐点。这些特征在器件设计时直接决定了匹配网络的形式。3. MATLAB实现完整代码与逐段讲解3.1 参数设置与物理常量先定义物理常数全部采用国际单位制SI。我见过太多人代码跑出来差十几个数量级最后发现是把eV和J搞混了。%% 物理常数SI e_charge 1.602176634e-19; % 元电荷 [C] hbar 1.054571817e-34; % 约化普朗克常数 [J*s] kB 1.380649e-23; % 玻尔兹曼常数 [J/K] %% 模型参数 T 300; % 温度 [K] E_F 0.2 * e_charge; % 费米能级 [J]即 0.2 eV tau 0.5e-12; % 散射时间 [s]0.5 ps费米能级用E_F 0.2×e_charge是因为我们习惯用eV表达能级但公式里所有能量都是焦耳。换算关系1 eV 1.602e-19 J。如果你从文献里看到E_F0.2eV一定要乘e_charge再代入公式否则结果会差1.6e-19倍。温度参数T默认300K对应室温。如果是低温测量模拟可以设成77K或4.2K但要注意温度越低带间跃迁边越陡峭数值积分对奇异点处理的要求越高delta截断就要选得更小。3.2 频率向量与带内项实现频率范围根据你的目标频段设置。这里覆盖0.1~10THz基本覆盖了整个太赫兹窗口。采样点数Nf设500个够画平滑曲线计算速度快结果也不受积分误差影响。f linspace(0.1e12, 10e12, 500); % 频率 [Hz] omega 2*pi*f; % 角频率 [rad/s] % 带内电导率Drude-like ln_term abs(E_F)/(2*kB*T) log(1 exp(-abs(E_F)/(kB*T))); sigma_intra -1i * e_charge^2 * kB * T ./ (pi * hbar^2 * (omega - 1i/tau)) .* ln_term;注意我用了abs(E_F)的写法。这个对数项本质上只依赖|E_F|因为石墨烯电子-空穴对称P型掺杂和N型掺杂的带内输运特性相同。如果你直接写成E_F当E_F为负值时可能会在log里产生假想成分这是新手常犯的错误。带内项的向量化写法省去了循环。MATLAB中直接用点除./和点乘.*速度比for循环快两个数量级。这里没有奇异点问题因为分母(ω - i/τ)在实轴ω上永远不为零τ有限时也不会遇到奇异性。3.3 带间项数值积分奇异点处理是关键带间项是整段代码最需要小心的部分。直接照抄积分表达式跑一遍integral函数会在ε→ℏω/2的位置警告Minimum step size reached甚至返回NaN。%% 辅助函数 G(ε) Gfun (eps) sinh(eps/(kB*T)) ./ (cosh(E_F/(kB*T)) cosh(eps/(kB*T))); %% 带间电导率 sigma_inter zeros(size(omega)); for idx 1:length(omega) w omega(idx); Ep hbar * w; % 光子能量 E_h Ep/2; % ε₀ ℏω/2 Ghalf Gfun(E_h); integrand (eps) (Gfun(eps) - Ghalf) ./ (Ep^2 - 4*eps.^2); % 奇异点邻域截断 delta max(E_h*1e-9, 1e-30); % 左区间 [0, E_h - delta] if E_h - delta 0 I_left integral(integrand, 0, E_h - delta, RelTol, 1e-9, AbsTol, 1e-13); else I_left 0; end % 右区间 [E_h delta, Emax] Emax max([20*abs(E_F), 20*kB*T, 10*Ep]); % 截断上限 I_right integral(integrand, E_h delta, Emax, RelTol, 1e-9, AbsTol, 1e-13); PV I_left I_right; sigma_inter(idx) (e_charge^2/(4*hbar)) * (Ghalf 1i*(4*Ep/pi)*PV); end关于奇异点这里有个细节值得展开。被积函数在εE_h处分子分母同时为零极限值有限但数值积分算法在奇点附近会不断细分区间导致速度变慢甚至失败。我的处理方式是直接挖掉一个delta邻域。delta太小积分结果对delta敏感delta太大会牺牲准确性。实测下来取E_h的1e-9倍是一个平衡点结果对delta的依赖低于积分容差。还有一点值得注意积分上限直接取Inf会让integral函数在无穷区间的尾部反复自适应有时会告警不收敛。更稳的做法是像我这样取一个物理上足够大的截断值Emax。当ε远远超过max(E_F, k_BT, E_h)时G(ε)趋近于1被积函数近似为(1-Ghalf)/(Ep²-4ε²)在ε→∞时按1/ε²衰减所以截断到几十倍特征能量之后剩余尾部积分可以忽略。如果你需要更高精度可以把奇异点邻域单独提出来用解析极限替换。被积函数在εE_h附近的极限为-G(E_h)/(4Ep)而G(ε)可以解析求出G(ε) [cosh(x)cosh(a)1] / (k_BT × (cosh(a)cosh(x))²)其中xε/(k_BT)aE_F/(k_BT)。在编码中可以把这个极限值赋给奇异点处的被积函数值然后integral函数在含内部奇异点的区间上也能正常工作。3.4 可视化与结果验证总电导率把两项加起来然后直接取倒数得到表面阻抗sigma_total sigma_intra sigma_inter; Z_s 1 ./ sigma_total; figure; subplot(1,2,1); plot(f/1e12, real(sigma_total), b-, LineWidth, 1.5); hold on; plot(f/1e12, imag(sigma_total), r--, LineWidth, 1.5); grid on; xlabel(频率 (THz)); ylabel(电导率 (S)); legend(Re σ,Im σ,Location,best); title(石墨烯复电导率 (E_F0.2eV, T300K)); subplot(1,2,2); plot(f/1e12, real(Z_s), b-, LineWidth, 1.5); hold on; plot(f/1e12, imag(Z_s), r--, LineWidth, 1.5); grid on; xlabel(频率 (THz)); ylabel(表面阻抗 (Ω/sq)); legend(Re Z_s,Im Z_s,Location,best); title(石墨烯表面阻抗);跑完代码后先别急着做其他事情。我强烈建议你做两个验证一是DC验证。在f0.1THz附近近似静态极限看σ_total的实部是否接近理论值。E_F0.2eV、τ0.5ps、T300K时σ_intra在低频段实部应该在几十mS量级表面阻抗实部在几十到上百Ω/sq量级。二是带间边验证。E_F0.2eV时带间跃迁阈值在ℏω0.4eV对应频率f2E_F/h≈96.8THz。所以在我们0.1~10THz的频率窗口内理论上看不到带间吸收的突变。但你可以把最高频率扩展到200THz扫一扫会在约96.8THz处看到σ_inter实部突然上升。如果这个特征位置不对说明能量单位换算或E_F设置有问题。4. 参数扫描物理规律与设计参考4.1 费米能级如何调控表面阻抗费米能级是石墨烯器件最关键的调控旋钮因为在仿真中你可以通过改变栅极电压来连续调控E_F。我把E_F分别设为0.1eV、0.2eV、0.5eV其他参数固定T300Kτ0.5ps扫描频率0.1~10THz。结果非常直观E_F越高带内电导率实部越大表面阻抗越小。在1THz处E_F0.1eV时表面电阻可能接近几百Ω/sqE_F0.5eV时可以降到十几Ω/sq。这跟DC方块电阻的E_F成正比关系是一致的。更关键的是虚部。低频段石墨烯表面阻抗虚部为负呈电容性还是感性这个需要按你用的时间约定来判断。在我的符号约定e^{-iωt}下表面阻抗虚部为负对应容抗表现为电场能量存储。随着频率升高虚部的绝对值先增大后减小。这说明石墨烯在太赫兹频段的行为像一个损失较大的RC分布式电路而不是简单的电阻层。做天线匹配时这个虚部必须纳入考量否则你设计的阻抗匹配点会整体偏移。我还发现一个规律E_F越高带间阈值频率越高。阈值频率f_inter 2E_F/hE_F0.5eV时阈值在241THz已经远超出太赫兹窗口。所以在太赫兹频段做设计时如果费米能级调得比较高带间项对电导率的贡献确实可以忽略Drude近似就够用。但如果你做的是中红外30~50THz传感器那么带间项就必须保留了。4.2 温度与散射时间的作用把温度从10K扫到500K观察电导率谱的变化会发现两个效应。第一个是带内项的幅度整体变化温度越高ln(2cosh(E_F/2k_BT))越偏离低温极限因为热激发让更多载流子参与输运DC电导率轻微上升。第二个是带间跃迁边被展宽低温时跃迁边非常锐利接近阶跃函数室温300K下这个边被热展宽抹成了平滑过渡500K时则进一步降低跃迁边的陡峭度。这种现象对红外探测器的光谱选择性设计有直接影响——低温器件更容易实现锐利的吸收边。散射时间τ对电导率的影响集中在带内项。τ从0.1ps到1ps变化时DC电导率线性增大σ_DC∝τ但频谱的Drude拐点频率f_D 1/(2πτ)会左移。τ0.1ps时拐点在1.6THzτ1ps时拐点降到160GHz。如果你的工作频段在几个THzτ1ps的样品在太赫兹高端已经处于欧姆损耗主导的区域表面阻抗实部随频率上升明显而τ0.1ps时整个太赫兹窗口仍在惯性主导区域电导率随频率变化较平缓。这个趋势对判断器件是损耗受限还是惯量受限很重要。散射时间的取值很有讲究。高迁移率悬浮石墨烯的τ可以到几皮秒而跟衬底强耦合的CVD石墨烯通常只有0.1~0.3ps。做仿真时尽量根据自己样品的霍尔迁移率来估算τ μE_F/(e v_F²)其中v_F≈1e6m/s。迁移率10000cm²/(V·s)的样品E_F0.2eV时τ≈0.63ps迁移率5000cm²/(V·s)时τ≈0.31ps。不要闭着眼睛填τ。4.3 从复电导率到等效介电常数有些仿真场景不提供阻抗边界条件选项比如某些时域有限差分软件这时可以把石墨烯等效为厚度为d的介质层用等效相对介电常数来建模ε_eq(ω) 1 iσ(ω)/(ωε₀d)推导过程很简单把石墨烯看作厚度d的介质面电流J_s σE体电流J_3D σ_3D_E而σ_3D σ/d。另一方面介质中的电流J σ_3D E -iωε₀(ε_r-1)E所以上面那个公式就出现了。实现时取d1nm或者0.5nm取决于你想把多少额外厚度强加给模型。在太赫兹频段ε_eq的虚部会非常大因为σ大、频率低这会让时域仿真中的Courant稳定性条件变得苛刻。如果你用FDTD做石墨烯器件仿真强烈建议还是用阻抗边界条件而不是等效介质层否则网格需要做得很细计算资源会被吃掉一大块。等效介质法有一个坑厚度d是一个人为拟合参数换一个d等效介电常数就变了。不同文献里d取值从0.3nm到10nm都有这导致了文献间结果对不齐的乱象。如果你需要跟某篇文献对比仿真数据一定先看它用的d是多少。5. 常见问题与排查实录5.1 积分告警与不收敛最常遇到的报错是integral函数提示Maximum function count exceeded或者Minimum step size reached通常出现在频率接近带间阈值、温度较低、或者奇异点delta取太小的时候。我总结了一套排查顺序先看奇异点delta。如果delta取得过小比如E_h的1e-12倍左右两个积分区间几乎贴着奇异点自适应积分会在边界处疯狂细化导致函数调用次数爆炸。调到E_h的1e-8~1e-9相对稳妥。再看截断上限Emax。如果Emax取得不够大比如只取5kBT而E_F远大于k_BT那么G(ε)在截断点附近还没有完全饱和到1被积函数尾部截断误差会传进主值积分。建议Emax至少取max(20|E_F|, 20k_BT, 10Ep)。最后看频率采样。当你扫频时如果带间阈值恰好落在你的频率网格中间你会看到电导率实部出现异常的毛刺这不是物理问题而是数值分辨率不够。加密采样点或者把频率扫描设成对数间隔可以改善。5.2 低频段表面阻抗虚部跳变我在最初写代码时发现σ_total的实部在非常低的频率比如0.01THz以下会轻微震荡虚部则乱七八糟。排查后发现是带间积分在低频极限下出了问题。当ω→0时E_h→0奇异点ε0恰好落在积分左端点。此时被积函数在ε→0时本身有定义可以用洛必达法则算但integral函数不知道这个内部奇异点来自端点自适应过程容易被误导。解决方法是把左区间积分拆成[0, 某个极小值]和[极小值, E_h-delta]两段或者干脆在低频段比如ℏω0.01k_BT用解析的低频展开式代替数值积分。实际上低频段远小于带间阈值带间项贡献本身就很小实部趋近于零虚部有一个有限的负值。很多应用场景中在ω远小于2E_F/ℏ时直接忽略带间项也不会差太多。我通常设一个判断如果f 0.1×2E_F/h就跳过带间积分的精细计算只保留带内项。这样既加速又避免了奇异积分。5.3 量纲与单位自查清单这个领域最常见的翻车现场就是单位制混乱。我列一个自查清单每次跑代码前过一遍E_F和k_BT必须同为焦耳或同为eV。如果E_F用eVk_BT也用eV300K对应约0.0259eV那么ln项无量纲没问题。但交给MATLAB时e_charge²、hbar²、kB这些常数全是SI单位你混用eV就会在前面系数上出错。ℏ和h别搞混。Kubo公式里用的是约化普朗克常数ℏ不是普朗克常数h。两者差2π倍误差巨大。类似地频率转能量时用Ehf但如果是角频率ω则Eℏω。我的代码里用的是omega2πf所以光子能量Ephbar*omega没错。电导率单位是S不是S/m。如果你从某篇论文里看到σ1e-3单位标注S这是二维电导率可以直接代入Z_s1/σ。如果看到S/m那是三维体电导率必须搞清楚他们用的等效厚度是多少才能转换。表面阻抗的虚部符号依赖时间约定。CST、HFSS、FDTD方法、频域求解器各自的约定可能不同有的用e^{iωt}有的用e^{-iωt}对应的电导率虚部符号会翻转。在把MATLAB计算出的σ赋给仿真软件前务必确认软件的约定。通常CST的Drude模型选项会带你设置参数的imaginary part sign别默认瞎填。5.4 性能优化与代码加速上面的代码在Nf500时大概几秒钟跑完完全够用。如果你要做大规模参数扫描比如E_F扫100个点、τ扫50个点500×100×50250万个频率点多次积分就不太妙了。这时候有几个提速手段。第一个是并行循环。MATLAB的parfor可以直接把for loop改成并行前提是装有Parallel Computing Toolbox。带间项积分之间相互独立天然适合并行。第二个是减少重复计算。Gfun里cosh(E_F/(kB*T))是一个常数可以提到循环外面算一次。同理带内项的ln_term也是一个常量向量没必要在每次频率循环里重复调用。第三个是缓存积分结果或使用查表插值。如果不要求每一点的绝对精度可以先用稀疏频率网格比如50个点算完积分再用spline插值得到500个点的值误差通常能控制在0.5%以内。这在设计迭代中绰绰有余速度能提升一个数量级以上。第四个如果你用的是低温解析近似θ对数形式带间项根本不需要积分整个σ可以写成全解析表达式向量化直接跑速度最快。但记得只在T远小于E_F/k_B时使用否则结果跟完整Kubo对不齐。6. 写在最后的实操体会代码分享到这里最后说几点我自己的体会。第一石墨烯Kubo公式虽然看起来像一团数学符号但拆开之后物理非常清晰带内项就是载流子的摩擦输运带间项就是光生电子-空穴对的量子跃迁两者在频域上互补低频靠带内高频靠带间中间的过渡区域由费米能级和温度共同决定。理解了这个图像你在调试代码时至少不会走错方向。第二无论数值结果怎样一定要先过一个常识检查DC电导率是否跟E_Fτ成正比表面阻抗是否在合理的几十到几百Ω/sq范围带间吸收边是否落在2E_F的位置这三个检查能在五分钟内帮你抓住90%的错误。第三真实项目中石墨烯的射频/微波特性非常依赖衬底。SiO₂衬底上的石墨烯跟悬浮石墨烯的散射时间可以差一个数量级从而让表面阻抗谱完全不同。仿真时把τ当可调参数跟实验S参数做拟合反推样品实际的τ和E_F这是一个非常实用的标定思路。碰到测不准的样品参数就用这个办法来校准。如果用这套代码复现出问题建议先跑一遍T300K、E_F0.2eV、τ0.5ps的默认参数把得到的结果跟我上面描述的量级对照一下。确认基线正确后再去扫描其他参数这样就不会一头扎进数值迷雾里出不来。祝各位在石墨烯器件仿真的路上一次跑通不用像我当初那样在奇点积分里泡了整整两天。
