跳转至

线性因果耦合传播：marcher 与 GEM 专用 IF2

摘要

本文说明 CausalWave 的数学模型、离散方程、实现映射和误差来源。底层 marcher 面向式 (3) 所定义的一般线性耦合系统：一组场沿空间单向传播，另一组场沿时间演化；它在每个网格点联立求解梯形离散所得的小型复线性系统。当前公开接口与主要数值验证围绕一维弱场 GEM 展开。IF2 则是这个 GEM 模型的专用后端：用空间累积积分消去光场，精确处理区间内的局部对角相位，再对耦合项作二阶端点指数拟合。空间 Volterra 结构允许 IF2 每个时间步只解一个下双对角系统。这里的“精确”仅指指定离散式内的局部指数与线性代数求解，不意味着完整 PDE 的精确解。

梯形离散、指数求积与 Volterra 积分都是已有数值方法。CausalWave 提供这些方法的结构化实现、问题接口和数值验证，不把基础算法称为本项目原创。Gourlay（1970）和Hochbruck 与 Ostermann（2010）给出下文所用方法的直接背景；具体对应关系分别见 §2.4 与 §3.2。

关键词： 线性耦合传播；复振幅；Volterra 积分；梯形法；指数拟合；梯度回波记忆。

1. 模型、边界与适用范围

设 \(z\in[0,L]\)、\(t\geq0\)，探针场为 \(E\in\mathbb C\)，原子相干性为 \(\alpha\in\mathbb C\)：

\[ \partial_zE=iN\alpha,\qquad \partial_t\alpha=D(z,t)\alpha+igE,\qquad D(z,t)=\gamma+i\eta s(t)(z-L/2). \tag{1} \]

边界和初态为

\[ E(0,t)=u(t),\qquad \alpha(z,0)=0. \tag{2} \]

\(g,N,\eta,L\) 为实参数；它们必须与坐标、入射场采用一致的单位或无量纲化。代码不进行自动单位换算。\(\gamma\) 是式 (1) 的有符号系数：退相干写作 \(\gamma=-\Gamma\)、\(\Gamma\geq0\)；输入正值将产生增长，不能在报告中称作衰减。

常见 GEM 控制为写入段 \(s=+1\)、读出段 \(s=-1\)。这里允许任意有限的实值分段常数 \(s(t)\) 和复数分段线性 \(u(t)\)。这是弱场、单向传播的线性模型；不包含布居饱和、反向传播、扩散或横向衍射。低层 marcher 另有三能级 LambdaGEM，并不意味着本页 IF2 已推广到该体系。

GEM 的两能级背景可参阅Hétet 等（2008），Raman GEM 的模型与适用近似可参阅Hosseini 等（2012）。本页直接声明软件求解的方程；这些物理参考不是参数单位、衰减符号或实际器件性能已经核对的替代品。三能级 LambdaGEM 保留光学极化与基态相干性，弱场 EIT 的背景见Fleischhauer 与 Lukin（2000）；它不是包含布居饱和的完整密度矩阵模型。

层次 已提供能力 验证与适用边界
marcher.march_system 与 Protocol 式 (3) 的线性多场、已知时空系数及源项 有制造解、边界和小块线性代数检查；每个新物理模型仍需独立核对
GEM 与 LambdaGEM 两能级 GEM 和线性三能级协议；后者可设置 EIT 参数 包含两能级收敛及三能级解析透射检查；不等于任意协议已充分验收
causalwave.solve_gem 的 IF2 后端 式 (1) 的两复场、分段驱动、零初始相干性 没有通用多场 IF2、完整 Lindblad、双向或多维传播接口

1.1 原生驱动与真实事件

原生节点 \(\tau_k\) 定义输入：在 \([\tau_k,\tau_{k+1})\) 内，\(s(t)=s_k\)，\(u(t)\) 按两个端点线性插值。内部时间网格取原生节点、输出节点和零时刻的并集再逐区间加密；因此不会跨越一个未显式拆分的已知跳变。

到达翻转时刻的前一积分区间使用左侧控制；从该时刻出发的后一积分区间使用新控制。状态连续，导数可跳变。降低输出采样率不会改变驱动；但遗漏驱动事件会改变被求解的问题，增加求解器步数不能补回缺失的信息。

2. marcher：时空梯形块系统

2.1 一般形式

令输运场 \(U\in\mathbb C^P\)、局部演化场 \(V\in\mathbb C^Q\)：

\[ U_z=AV+s,\qquad V_t=BV+CU+r. \tag{3} \]

给定左边界 \(U(0,t)\) 与初态 \(V(z,0)\)。GEM 对应 \(P=Q=1\)、\(A=iN\)、\(B=D\)、\(C=ig\)、\(s=r=0\)。系数与源项可以依赖时空，但场的耦合必须是线性的局部零阶项。

在 \((z_i,t_j)\) 记 \(c_i=(z_i-z_{i-1})/2\)、\(d_j=(t_j-t_{j-1})/2\)。空间、时间方向分别用梯形公式，得到

\[ \begin{pmatrix} I&-c_i A_{ij}\\ -d_j C_{ij}&I-d_j B_{ij}^{-} \end{pmatrix} \begin{pmatrix}U_{ij}\\V_{ij}\end{pmatrix} = \begin{pmatrix}R^U_{ij}\\R^V_{ij}\end{pmatrix}, \tag{4} \]

其中

\[ R^U_{ij}=U_{i-1,j}+c_i\left(A_{i-1,j}V_{i-1,j}+s_{i-1,j}+s_{ij}\right), \tag{5} \]
\[ R^V_{ij}=V_{i,j-1}+d_j\left(B_{i,j-1}V_{i,j-1}+C_{i,j-1}U_{i,j-1}+r_{i,j-1}+r_{ij}\right). \tag{6} \]

\(B_{ij}^{-}\) 对声明为阶跃的分量取当前时间区间内部的左侧值；其他平滑分量取端点值。这与代码中的 B_step_mask 一致。时间端点的 \(C_{ij}U_{ij}\) 已进入式 (4) 左侧，不应在式 (6) 重复添加。

2.2 因果推进顺序

式 (4) 只依赖左侧点 \((i-1,j)\) 与下侧点 \((i,j-1)\)。按 \(i+j=1,2,\ldots\) 的顺序推进，同一对角线上各点可向量化求解。两个边缘必须分别处理：

  1. \(t=0\)：\(V\) 由初态固定，\(U\) 沿空间梯形积分。对于零相干性，应有整个初始截面 \(E(z,0)=u(0)\)。
  2. \(z=0\)：\(U\) 由边界固定，\(V\) 沿时间梯形积分。一般 \(\alpha(0,t)\) 不为零。
  3. 内部点：联立求解式 (4)，保留复振幅与相位。

当前实现对 \(2\times2\)、\(3\times3\) 块使用显式小矩阵公式，其他维度使用 NumPy 线性代数。它是对角线推进，不是沿某条特征线返回 PDE 解析解。

2.3 精度与成本

在分段光滑、网格充分解析且离散系统适定时，梯形空间与时间离散具有二阶一致性。必须检查网格收敛；不能把标量 Crank–Nicolson 的 A 稳定性直接推广成任意耦合系数、任意步长下的无条件准确或稳定。

固定 \(P,Q\) 时，工作量与内部格点数 \(n_z n_t\) 成正比；一般块求解另有约 \((P+Q)^3\) 的局部代价。当前实现保留细时空系数和状态，内存随 \(n_z n_t\) 增长。这是大加密倍数下需要注意的主要成本。

2.4 与既有梯形格式的关系

空间和时间方向各用一次端点梯形公式，是式 (4)–(6) 的基本构造。在标量、常数 \(A,C\)、\(B=0\) 且无源项的子类中，消去 \(V\) 可得 \(U_{zt}=ACU\)；同样消去离散未知量后，得到四角平均的 Goursat 梯形差分。这一子类与Gourlay（1970）§2讨论的经典 Goursat 梯形构造对应。

这是对本实现的代数归类，不表示该文直接给出了本项目的所有变系数、多场或事件处理代码。对角线推进利用的是网格依赖关系，小型块求解和向量化属于实现选择。本项目不把这些特征单独称作新数值原理，也不据此把一般格式直接等同于所有被称作 Keller box 的方法。

3. IF2：空间消元与二阶指数拟合

这里 IF2 是本项目对以下二阶指数拟合格式的简称，不是对所有文献中“IF2”命名方法的等同声明，也不是一个新的算法优先权主张。当前实现的空间消元和系数构造专用于式 (1)。

3.1 用 Volterra 积分消去光场

由式 (1) 的 Maxwell 方程与边界得

\[ E(z,t)=u(t)+iN\int_0^z\alpha(\zeta,t)\,d\zeta. \tag{7} \]

在空间网格上，用累积梯形算子 \(Q\)：

\[ (Qa)_j=\sum_{\ell=1}^{j}\frac{\Delta z_\ell}{2}(a_{\ell-1}+a_\ell),\qquad (Qa)_0=0. \tag{8} \]

于是空间半离散系统为

\[ \dot a=Da-\kappa Qa+igu(t)\boldsymbol1,\qquad \kappa=gN. \tag{9} \]

\(Q\) 虽可表示成下三角矩阵，代码并不构造稠密矩阵，而是做累积求和。

3.2 变常数公式与指数权重

在没有控制跳变的时间区间 \([t_n,t_{n+1}]\) 内，\(D\) 是固定对角算子，步长为 \(h\)。令 \(f=-\kappa Qa+igu\boldsymbol1\)，则

\[ a_{n+1}=e^{hD}a_n+\int_0^h e^{(h-\tau)D}f(t_n+\tau)\,d\tau. \tag{10} \]

对积分内的 \(f\) 作端点线性近似。定义逐元素计算的

\[ \varphi_1(x)=\frac{e^x-1}{x},\qquad \varphi_2(x)=\frac{e^x-1-x}{x^2},\qquad \varphi_1(0)=1,\quad\varphi_2(0)=\tfrac12, \]
\[ R=e^{hD},\qquad w_0=h[\varphi_1(hD)-\varphi_2(hD)],\qquad w_1=h\varphi_2(hD). \tag{11} \]

这些权重是已有指数梯形求积的端点权重。Hochbruck 与 Ostermann（2010），§2.2、Example 2.6对节点 \(0,1\) 给出 \(\varphi_1-\varphi_2\) 与 \(\varphi_2\)；取其线性算子为 \(-D\)，即得到式 (11)。在本问题中，右端耦合依赖未知的 \(a_{n+1}\)，因此还要解式 (12)–(13) 的隐式系统。文献中不同的积分因子、ETD 或 IIF 方法未必使用相同权重，不能仅凭简称移用它们的稳定性结论。

实现对小 \(|hD|\) 使用级数，避免相近数相减导致的消减误差。更新式为

\[ b=Ra_n-\kappa w_0Qa_n+ig(w_0u_n+w_1u_{n+1}),\qquad a_{n+1}=b-q p,\quad q=\kappa w_1,\quad p=Qa_{n+1}. \tag{12} \]

向量乘积 \(w_0Qa_n\)、\(qp\) 均为逐元素乘法。当 \(D=0\)，\(w_0=w_1=h/2\)，式 (12) 退化成空间半离散系统的梯形时间格式。

3.3 双对角系统

设 \(c_j=\Delta z_j/2\)。由 \(p_j-p_{j-1}=c_j(a_{j-1}+a_j)\)，将式 (12) 代入可得

\[ (1+c_jq_j)p_j-(1-c_jq_{j-1})p_{j-1}=c_j(b_{j-1}+b_j),\quad j\geq1;\qquad p_0=0. \tag{13} \]

这就是下双对角系统。实现用 SciPy solve_banded((1, 0), ...)，随后恢复 \(a_{n+1}=b-qp\)，重新累积积分并按式 (7) 恢复光场。重新计算积分使下一步状态与累积梯形算子保持一致；不是额外的训练或拟合。

式 (13) 的局部递推来自这里特定的空间累积梯形结构。它解释了当前实现为什么不需要存储稠密 \(Q\)；不意味着一般 Volterra 核、多场耦合或任意非线性系统也能直接使用同一个双对角求解器。

3.4 “指数处理”解决和没有解决的部分

  • 区间内局部对角相位 \(e^{hD}\) 被直接计算，减少普通梯形推进对高频局部相位的时间色散。
  • 当 \(N=0\) 且驱动按合同分段线性时，式 (10) 的区间积分可由这些权重精确表达，误差主要是浮点运算。
  • 有耦合时 \(Qa(t)\) 一般不是时间线性的，端点拟合仍有截断误差；IF2 不是完整耦合系统的精确指数。
  • 高频 \(\alpha(z,t)\) 仍需要足够的空间分辨率；精确处理时间相位不能恢复未解析的空间振荡。

在适当稳定与光滑条件下，时间全局误差为二阶，空间累积梯形误差也为二阶。跳变需要明确分段。是否达到 1% 必须由两个复场的独立收敛和参考比较决定。

3.5 缓存与复杂度

每步累积积分和双对角解均为 \(O(n_z)\)，总工作量为 \(O(n_z n_t)\)。内部状态为 \(O(n_z)\)，外加 \(O(n_z^{out}n_t^{out})\) 的请求输出。当前实现每次求解内部使用最多 64 MiB 的 LRU 指数权重缓存，以精确 (h, sw) 为键，不舍入步长或控制，也不跨不同参数/网格的请求复用。浮点步长看似相同但二进制不同会产生不同键；缓存命中率影响耗时，不改变离散公式。

4. 两个实现的比较边界

通过 solve_gem 使用时，两者求解相同的两能级模型，但误差常数、内存布局和相位处理不同。IF2 的当前实现利用了该 GEM 模型的 Volterra 结构；marcher 的一般块接口可表达多个局部演化场。接口能表达一个模型，与该模型已通过物理验证，是两个不同结论。不能只用代码行数、一个示例的速度或图形相似来判断哪个方法“更准确”。

公平比较应固定：物理参数、原生驱动、输出坐标、复场误差定义、参考精度、线程数、硬件和是否计入编译。先为每个后端找满足相同精度门槛的网格，再报告时间与峰值内存。当前仓库的验证页面把自包含解析证据与历史面板分开说明。

5. 实现与参考

  • 式 (4)–(6)：src/marcher/core.py 的 _solve_block 与 march_system。
  • 式 (11)：src/causalwave/_if2_reference.py 的 _exponential_weights。
  • 式 (12)–(13)：src/causalwave/_if2.py 的 solve_gem_cached。
  • 统一驱动/坐标合同：src/causalwave/api.py。

带状矩阵存储遵循 SciPy solve_banded 官方定义。XMDS2 是独立工具链，见Dennis 等（2013）。完整书目信息、物理背景和 CoOMBE、clerq、mbsolve 等相邻软件的范围见参考文献与相关软件。

本页的推导对应随仓库发布的离散实现；验证页记录实际通过的检查与数据边界。推导、实现测试和真实器件验证各自回答不同问题，不应互相替代。