企业官网建设流程全解析

搞电磁仿真的朋友应该都遇到过这个尴尬：在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_F（E_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(1+e^{-E_F/k_BT})

于是带内项趋近于：

σ_intra ≈ e²E_F/(πℏ²(ω - i/τ)) × (-i) ≈ e²E_Fτ/(πℏ²(1-iωτ))

这就是标准的Drude模型，DC电导率σ_DC = e²E_Fτ/πℏ²。随便代入几个数：E_F=0.2eV、τ=0.5ps，σ_DC≈11.8mS，对应表面阻抗约85Ω/sq，跟实验上测得的高质量石墨烯直流方块电阻量级完全一致。这个验证很重要，说明公式和代码没写错。

在室温300K下，E_F=0.2eV时，k_BT≈0.0259eV，E_F/(2k_BT)≈3.86，ln(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时G=0，当ε > E_F时G=1。这时带间项可以解析积分出来：

σ_inter = e²/(4ℏ) × [ θ(ℏω-2E_F) + (i/π) ln |(ℏω+2E_F)/(ℏω-2E_F)| ]

θ是单位阶跃函数。这个公式是大家最常用的"低温解析近似"，非常方便手算验证代码。它清楚说明了一个物理事实：当ℏω < 2E_F时，带间吸收被Pauli阻塞抑制，实部很小；当ℏω > 2E_F时，实部突变为e²/4ℏ≈60.8μS，这是石墨烯在红外波段的"统一透明率"对应的电导率量子值。

但做室温仿真时，T=300K下热展宽会把这条突变边抹平成平滑过渡，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里用体电导率建石墨烯模型，要把二维电导率除以等效厚度。比如σ_2D=0.1mS，等效厚度取1nm，则体电导率应为σ_3D=0.1mS/1nm=1e5 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_F=0.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)，a=E_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_F=0.2eV, T=300K)'); 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验证。在f=0.1THz附近（近似静态极限），看σ_total的实部是否接近理论值。E_F=0.2eV、τ=0.5ps、T=300K时，σ_intra在低频段实部应该在几十mS量级，表面阻抗实部在几十到上百Ω/sq量级。

二是带间边验证。E_F=0.2eV时，带间跃迁阈值在ℏω=0.4eV，对应频率f=2E_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，其他参数固定（T=300K，τ=0.5ps），扫描频率0.1~10THz。

结果非常直观：E_F越高，带内电导率实部越大，表面阻抗越小。在1THz处，E_F=0.1eV时表面电阻可能接近几百Ω/sq，E_F=0.5eV时可以降到十几Ω/sq。这跟DC方块电阻的E_F成正比关系是一致的。

更关键的是虚部。低频段石墨烯表面阻抗虚部为负，呈电容性还是感性？这个需要按你用的时间约定来判断。在我的符号约定（e^{-iωt}）下，表面阻抗虚部为负对应容抗，表现为电场能量存储。随着频率升高，虚部的绝对值先增大后减小。这说明石墨烯在太赫兹频段的行为像一个损失较大的RC分布式电路，而不是简单的电阻层。做天线匹配时，这个虚部必须纳入考量，否则你设计的阻抗匹配点会整体偏移。

我还发现一个规律：E_F越高，带间阈值频率越高。阈值频率f_inter = 2E_F/h，E_F=0.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_F=0.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，所以上面那个公式就出现了。

实现时，取d=1nm（或者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用eV，k_BT也用eV（300K对应约0.0259eV），那么ln项无量纲，没问题。但交给MATLAB时，e_charge²、hbar²、kB这些常数全是SI单位，你混用eV就会在前面系数上出错。

ℏ和h别搞混。Kubo公式里用的是约化普朗克常数ℏ，不是普朗克常数h。两者差2π倍，误差巨大。类似地，频率转能量时用E=hf，但如果是角频率ω，则E=ℏω。我的代码里用的是omega=2πf，所以光子能量Ep=hbar*omega，没错。

电导率单位是S不是S/m。如果你从某篇论文里看到σ=1e-3（单位标注S），这是二维电导率，可以直接代入Z_s=1/σ。如果看到S/m，那是三维体电导率，必须搞清楚他们用的等效厚度是多少才能转换。

表面阻抗的虚部符号依赖时间约定。CST、HFSS、FDTD方法、频域求解器各自的约定可能不同，有的用e^{+iωt}，有的用e^{-iωt}，对应的电导率虚部符号会翻转。在把MATLAB计算出的σ赋给仿真软件前，务必确认软件的约定。通常CST的Drude模型选项会带你设置参数的"imaginary part sign"，别默认瞎填。

5.4 性能优化与代码加速

上面的代码在Nf=500时大概几秒钟跑完，完全够用。如果你要做大规模参数扫描（比如E_F扫100个点、τ扫50个点），500×100×50=250万个频率点多次积分就不太妙了。这时候有几个提速手段。

第一个是并行循环。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，这是一个非常实用的标定思路。碰到测不准的样品参数，就用这个办法来校准。

如果用这套代码复现出问题，建议先跑一遍T=300K、E_F=0.2eV、τ=0.5ps的默认参数，把得到的结果跟我上面描述的量级对照一下。确认基线正确后，再去扫描其他参数，这样就不会一头扎进数值迷雾里出不来。祝各位在石墨烯器件仿真的路上，一次跑通，不用像我当初那样在奇点积分里泡了整整两天。