麦克斯韦全方程绘制TE/TM模式N-厚度曲线异常：仅显示m=4

核心问题分析 #

1. 模式方程完全错误：原代码中根号内表达式写反——有效折射率$n_{eff}$必然满足$\max(n_s, n_c) < n_{eff} < n_f$（芯层折射率最高，模式被限制在芯层内），原代码中$\sqrt{n_{eff}^2 - n_f^2}$会得到虚数，直接导致数值求解失败，仅m=4偶然收敛到错误值。
2. 初始猜测值不合理：所有模式都用$n_f$作为初始猜测，但不同模式的$n_{eff}$分布范围不同（m越大，$n_{eff}$越接近衬底折射率$n_s$），单一初始值无法覆盖所有模式的收敛区间。
3. 忽略包层折射率$n_c$：原代码未考虑包层的倏逝波衰减，非对称波导的模式方程需要同时结合衬底和包层的影响。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root_scalar

# 给定参数
wavelength = 1550  # 波长（nm）
nc = 1.10          # 包层折射率
nf = 2.55          # 芯层折射率
ns = 1.94          # 衬底折射率
T_range = np.linspace(0.01, 5, 500)  # 波导厚度范围（μm）
mode_numbers = np.arange(0, 5)       # 模式m=0到4
k0 = 2 * np.pi / wavelength          # 自由空间波数（nm⁻¹）

# 初始化结果数组，默认填充无效值
TE_modes = np.full((len(mode_numbers), len(T_range)), np.nan)
TM_modes = np.full((len(mode_numbers), len(T_range)), np.nan)

# TE模式色散方程：返回值为0时满足模式条件
def te_dispersion(neff, m, T):
    k0T = k0 * T * 1000  # 转换厚度单位：μm→nm，与波长单位统一
    # 芯层横向波数（实数，模式在芯层内振荡）
    h = k0 * np.sqrt(nf**2 - neff**2)
    # 衬底、包层倏逝波衰减常数（实数，模式在衬底/包层内指数衰减）
    q_s = k0 * np.sqrt(neff**2 - ns**2)
    q_c = k0 * np.sqrt(neff**2 - nc**2)
    
    # TE模式相位匹配条件：芯层内相位变化 + 界面相位跃变 = mπ
    left = h * k0T
    right = m * np.pi + np.arctan(q_s / h) + np.arctan(q_c / h)
    return left - right

# TM模式色散方程：考虑Goos-Hänchen位移的修正因子
def tm_dispersion(neff, m, T):
    k0T = k0 * T * 1000
    h = k0 * np.sqrt(nf**2 - neff**2)
    q_s = k0 * np.sqrt(neff**2 - ns**2)
    q_c = k0 * np.sqrt(neff**2 - nc**2)
    
    # TM模式的界面相位跃变需要加入折射率平方比的修正
    factor_s = (nf**2 / ns**2) * (q_s / h)
    factor_c = (nf**2 / nc**2) * (q_c / h)
    
    left = h * k0T
    right = m * np.pi + np.arctan(factor_s) + np.arctan(factor_c)
    return left - right

# 遍历所有模式和厚度求解
for m_idx, m in enumerate(mode_numbers):
    # 为每个模式设置合理的求解区间：m越大，neff越接近ns
    neff_low = ns + 0.01
    neff_high = nf - 0.01
    # 调整初始猜测中心，m越大越靠近衬底折射率
    guess_center = neff_high - (m / (max(mode_numbers)+1)) * (neff_high - neff_low)
    
    for t_idx, T in enumerate(T_range):
        # 求解TE模式
        try:
            te_sol = root_scalar(te_dispersion, args=(m, T), 
                                method='brentq', bracket=[neff_low, neff_high])
            if te_sol.converged:
                TE_modes[m_idx, t_idx] = te_sol.root
        except ValueError:
            # 厚度不足时该模式不支持，保留无效值
            pass
        
        # 求解TM模式
        try:
            tm_sol = root_scalar(tm_dispersion, args=(m, T), 
                                method='brentq', bracket=[neff_low, neff_high])
            if tm_sol.converged:
                TM_modes[m_idx, t_idx] = tm_sol.root
        except ValueError:
            pass

# 绘制TE模式曲线
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
for m_idx, m in enumerate(mode_numbers):
    valid_idx = ~np.isnan(TE_modes[m_idx])
    plt.plot(T_range[valid_idx], TE_modes[m_idx][valid_idx], 
             linestyle='-', label=f'TE模式 m={m}')
plt.xlabel('波导厚度 (μm)')
plt.ylabel('有效折射率 N')
plt.title('TE模式有效折射率随波导厚度变化')
plt.legend()
plt.grid(True)
plt.ylim(ns-0.1, nf+0.1)

# 绘制TM模式曲线
plt.subplot(1, 2, 2)
for m_idx, m in enumerate(mode_numbers):
    valid_idx = ~np.isnan(TM_modes[m_idx])
    plt.plot(T_range[valid_idx], TM_modes[m_idx][valid_idx], 
             linestyle='--', label=f'TM模式 m={m}')
plt.xlabel('波导厚度 (μm)')
plt.ylabel('有效折射率 N')
plt.title('TM模式有效折射率随波导厚度变化')
plt.legend()
plt.grid(True)
plt.ylim(ns-0.1, nf+0.1)

plt.tight_layout()
plt.show()

关键修正点说明 #

1. 修正模式方程：
   • 正确区分芯层横向波数（实数）和衬底/包层的倏逝波衰减常数（实数），确保所有根号内表达式为正。
   • TM模式加入折射率平方比的修正因子，符合非对称波导的Goos-Hänchen位移效应。
2. 优化数值求解：
   • 使用root_scalar的brentq区间二分法，比root更适合单变量单调函数的求解，且指定了严格的求解区间（$n_s < n_{eff} < n_f$）。
   • 为每个模式设置差异化的初始猜测中心，提高收敛成功率。
3. 单位统一：将波导厚度从μm转换为nm，与波长单位统一，避免计算误差。
4. 无效值处理：当波导厚度不足以支持某模式时（薄厚度无法容纳高次模），保留nan并仅绘制有效数据点。

内容的提问来源于stack exchange，提问作者Abel Stewart