ECC3 程序说明

1. 概述与程序思想

ECC3 是经验耦合道（ECC）模型的升级版：进一步计算蒸发剩余（ER）截面 $\sigma_{\mathrm{ER}}$。在 ECC2 俘获模型（BW91* 势 + 不对称高斯势垒分布）的基础上，ECC3 加入简化统计衰变存活几率 $W_{\mathrm{sur}}$ 与复合核形成几率 $P_{\mathrm{CN}}$，构成"俘获 × 形成 × 存活"三段式框架。

程序思想一句话：蒸发剩余截面 =（俘获截面，ECC2 解析势垒分布）×（复合核形成几率 $P_{\mathrm{CN}}$，准裂变经验抑制）×（中子链存活几率 $W_{\mathrm{sur}}$，简化统计衰变）——三因子各司其职，分别回答"撞上了没有""粘住了没有""活下来了没有"三个物理问题。

2. 程序框架

2.1 总体计算式

先从整体图像出发：蒸发剩余截面由俘获截面、复合核形成几率与中子链存活几率三个因子相乘得到，

其中俘获截面 $\sigma_{\mathrm{cap}}$ 由 ECC2 解析势垒分布给出（§4.1），$P_{\mathrm{CN}}$ 为复合核形成几率（§4.5），$W_{\mathrm{sur}}$ 为存活几率（§4.4）。由于存活几率依赖角动量 $J=L$，需按分波写出：

$k^2=2\mu E_{\mathrm{c.m.}}/\hbar^2$，$E^*=E_{\mathrm{c.m.}}+Q_0$。程序在 sigma_ER.DAT 中同时输出 $\sigma_{\mathrm{ER}}$（第 2 列）、$W_{\mathrm{sur}}$（第 3 列）与 $P_{\mathrm{CN}}$（第 4 列），供拟合脚本使用；各中子蒸发道截面输出到 sigma_XN.DAT（§4.7、§6.3）。

2.2 主程序执行流程

3. 子程序一览

4. 物理模型细节

4.1 ECC2 俘获截面

BW91* 势 + 不对称高斯势垒分布 + Hill–Wheeler 透射：

势垒宽度 $W = 0.63(1+W_{\mathrm{def}}) + 0.015\,V_B\sqrt{\omega_1^2+\omega_2^2+\omega_0^2}$（壳隙项 $\Delta_B$ 禁用）。对 $A\le16$ 的轻弹核使用简化式 $W=0.63+0.015V_B$。

4.2 原子核基态壳修正

对 WS4 质量表做 6 参数液滴拟合（$k_s=2.2$，对称能系数 $a_{\mathrm{sym}}$ 取值位于合理范围），拟合常数由程序内置：

壳修正 $\delta_{\mathrm{sh}}=B_{\mathrm{WS4}}-B_{\mathrm{LD}}$。该方案对每个核素的壳修正都确定、连续、不依赖位垒表是否覆盖，且与后续裂变位垒完全解耦。

4.3 能级密度、分离能与裂变位垒

平滑能级密度参数 $\tilde{a}$ 的典型值为 A/12 ~ A/10 （上式中可裂变参数 $x$ 的液滴系数 $a_s,\ a_c,\ \kappa$ 取自 WS4 质量模型）。此处 $x$ 的表面能含 $(1-\kappa I^2)$ 不对称修正，用于能级密度参数 $K=6(1+x)$。

裂变位垒：优先取 WS4 表 Bf_WS4.txt（两列 $B_f$ 与 $B_f+\Delta B$，取较小值）；表中缺失时用宏观裂变位垒代替（基于$0.38(0.75-x)E_S$ || $0.88(1-x)^3E_S$ 计算）。有效位垒随激发能衰减并扣除转动能：

4.4 简化统计衰变：核心思想与存活几率递推

理论框架：ECC3 的统计衰变部分建立在 Hauser–Feshbach 统计模型（复合核退激）框架之上，但采用 Weisskopf–Ewing 型简化——各道宽度由能级密度相空间积分给出，不引入光学模型透射系数与自旋-宇称耦合。竞争道包括 中子 n、质子 p、α 粒子、γ 射线与裂变 五个出射道。

• 数值积分而非蒙卡：各道宽度（$I_n,I_p,I_\alpha,I_\gamma,I_f$）均由 50 点梯形积分计算（中子道在 NEUTRON_INTEGS 中一次算全部，p/α/γ/裂变道为 CPINT/GINT/FINT）；存活几率先在激发能–角动量网格（131×101）上递推并存储，运行时按 $(E^*,J)$ 直接查表；
• 采用近似：p/α 用库仑位垒阶跃 $\max(1-V_C/\varepsilon,0)$ 近似透射，γ 用统计 E1（∝ $C_\gamma\varepsilon^3$），裂变用 Bohr–Wheeler 型鞍点积分, 核子道的自旋多重度统一取 1；
• 纯中子链：仅沿中子链（$A\to A-1$）递推存活几率；p/α/γ 道只参与宽度竞争、不改变核种，其子核存活取 1；
• 角动量处理：转动能 $E_{\mathrm{rot}}(J)$ 从激发能中扣除，中子发射后 $J\to J-1$。

沿纯中子链（$A_0\to A_0-1\to\cdots$）从末端反推：

对于$Z \ge 82$ 的裂变核，鞍点能级密度参数带增强因子，其取值介于下限 $r_{\min}$ 与上限 $r_{\max}$ 之间，并按经典可裂变参数 $x_0$ 以 $(x_0/0.75)^3$ 幂律增强（这里 $x_0= \dfrac{a_c}{2 a_s} Z^2/A $，不含表面对称能项）：

其中下限 $r_{\min}=r_f$（$x_0\le 0.75$ 时增强因子取 1，$r_{\mathrm{eff}}=r_f$），上限 $r_{\max}=r_f\cdot r_{f,\max}$（增强因子被 RFMAX 截断，防止热熔合体系裂变竞争过度）。该规则用于复合核及其各代子核的裂变宽度。注意此处的 $x_0$ 与 §4.3 中用于能级密度参数 $K=6(1+x)$ 的可裂变参数 $x$【表面能含 $(1-\kappa I^2)$ 项】不同。

4.4.1 各道宽度积分的具体表达式

五道积分均由 50 点梯形公式在 $\varepsilon\in[0,u_0]$ 上数值求积（中子道 NEUTRON_INTEGS，带电/γ/裂变道 CPINT/GINT/FINT）。记复合核 $(A,Z)$ 的激发能 $E^*$、转动能 $E_{\mathrm{rot}}(J)$，扣除转动后的热激发能 $\tilde{E}=E^*-E_{\mathrm{rot}}(J)$；对能背移中 $\Delta_c=\Delta_{\mathrm{pair}}(A,Z)$ 为复合核对能（PAIRING_TABLE），$\Delta_d=\Delta_{\mathrm{pair}}(A-1,Z)$ 为发射中子后子核对能。各道表达式如下（共同因子在比例中约去）。

① 中子蒸发宽度

$S_n$ 为中子分离能（SEPN），$\rho_d$ 为子核 $(A-1,Z)$ 的能级密度（按子核质量数 $A-1$ 与壳修正计算），$W_{A-1}$ 为子核存活几率（网格查询，角动量取 $J-1$，能量取发射后剩余激发能）。

② 质子蒸发宽度

$S_p$ 为质子分离能（SE_P），$V_C^{(p)}$ 为质子-子核库仑位垒（Weisskopf–Ewing 阶跃近似）。

③ α 蒸发宽度

$S_\alpha$ 为 α 分离能（SE_A），$V_C^{(\alpha)}$ 为 α-子核库仑位垒，因子 $4$ 为 α 质量数（$g=1$，$\mu=4$）。

④ γ 发射宽度（统计 E1）

⑤ 裂变宽度（Bohr–Wheeler 型）

${B_f}^{\mathrm{eff}}$ 由 §4.3 给出（含壳修正阻尼与转动能扣除）；鞍点能级密度参数取 $a_{\mathrm{saddle}}=\dfrac{A}{K}\,r_{f,\mathrm{eff}}$（壳修正置零，use_shell=.FALSE.），$r_{f,\mathrm{eff}}$ 按本节上文的增强公式取值。

能级密度（五道共用）

各子核能级密度参数按各自质量数 $A'$、壳修正 $\delta_{\mathrm{sh}}(Z',A'-Z')$ 计算；存活几率 $W_A = (I_n+I_p+I_\alpha+I_\gamma)/(I_n+I_p+I_\alpha+I_\gamma+I_f)$。

4.5 $P_{\mathrm{CN}}$ 复合核形成几率

总体表达式（v3.0.30 统一形式）——所有体系的 $P_{\mathrm{CN}}$ 写成"渐近形成几率 $P_{0}(Z_1Z_2)$ × 能量因子 $W(E^*)$"两部分之积：

• $P_{0}$ 为归一化 Fermi（幅度 $A_0=1$），冷、热熔合同一形式，只用不同的 $\zeta,\tau$（$Z\!Z_0/DZ$ 与 $Z\!Z_{0c}/DZ_c$）；随 $Z_1Z_2$ 增大的准裂变抑制即由此体现（参考【Zagrebaev, Greiner, PRC 78 (2008) 034610】）。
• 热熔合能量因子 $W$ 为 Fermi 窗 × 深垒因子 $S$：$W_1(E^*)$ 乘以 $S=1/[1+\exp((B_s-E_{\mathrm{cm}})/C_{\mathrm{urv}})]$。$P_{\mathrm{CN}}$ 在高能端趋于 $P_{0}$，故此时 $P_{0}$ 是渐近值。$B_{\mathrm{fus}}$：$Z_1+Z_2\ge104$ 取 $V_B-B_{\mathrm{cap}}+\overline{E}_{\mathrm{def}}$，否则取 $-Q_0$。
• 冷熔合能量因子 $W$ 为 高斯，峰位 $E_0=g B_f$、峰高 $P_0\cdot A_{0c}$、宽度 $w=B_{\mathrm{cap}}-1$（分母中为 $2w$，1 为零点能）；两端下降。此时 $P_0$ 是峰值因子（不是渐近值）。奇 $Z_0$ 再乘 $F_{\mathrm{odd}}$。
• 深垒因子 $S(B_s)$ 的取舍：入射能量低于势袋底部时（$B_s=V_B-B_{\mathrm{cap}}$，$V_B$ 为最可几垒高），势垒厚度显著增加，导致熔合概率锐减。热熔合保留 （对总 ER 的 RMSD 影响 <0.02，但会改变超重体系分道峰值能量），冷熔合去掉 （高斯 $W$ 已含低 $E^*$ 下降）。
• 冷熔合判据（ICOLD）：$Z\ge104$ 且靶为 Pb/Bi（$\max(Z_1,Z_2)=82/83$）。低激发能存活增强（IRFEDC，§4.6 规则 3）对所有 $(B_m+Q_0)x<23$ 体系生效。

俘获势袋判据（v3.0.30）：当 $B_{\mathrm{cap}}\le1$ MeV 时，俘获势袋基本消失（深度已低于零点振动能 $\sim$1 MeV），熔合概率趋于 0，程序直接终止并提示，不再输出截面。例如 $^{82}$Se+$^{209}$Bi（$B_{\mathrm{cap}}\approx0.67$ MeV）、$^{86}$Kr+$^{208}$Pb（0.20）均被终止。

4.6 ${B_f}^{\min}$ 处理与鞍点增强（针对极端复合核：极度缺中子/极低势垒/极低激发能/极大可裂变参数；v3.0.33 增规则 5）

针对极端复合核 ER 截面的系统性偏离，ECC3 内置以下五条判据以改善模型结果：

规则 3，对于低激发熔合，考虑到低激发能下复合核存活概率增强、防止反推 $P_{\mathrm{CN}}>1$。

规则 4（v3.0.29 新增）针对 x>1 的极端不稳定超重核 ：这里的 $x$ 为含表面对称能项的液滴可裂变参数（见 §4.4）。与规则 2、3 的 $\pm0.05$ 加性修正不同，规则 4 是乘法增强，作用顺序为“基础增强 $\to$ 规则 2/3 的 $\pm0.05$ $\to$ 乘 $x^3$”。它把最重体系的鞍点能级密度进一步增强（存活下降），修正了此前对 $Z=118$（$^{48}$Ca+$^{249}$Cf）高估约 0.9 个数量级的问题，而 $Z\le115$ 的体系基本不受影响。

规则 5（v3.0.33 新增）针对极缺中子锕系复合核（$89<Z<104$、$S_{2p}<6.9$ MeV）：以 $S_{2p}$ 为标度把 $r_{\mathrm{eff}}$ 从标准值平滑推向 $r_{\max}=1.831$（$S_{2p}\to0$ 的滴线核取 $r_{\max}$），并令其所有蒸发子核沿用同一 $r_{\mathrm{eff}}$。该规则显著改善 $^{36}$Ar+$^{181}$Ta、$^{40}$Ca+$^{175}$Lu、$^{48}$Ca+$^{204}$Pb 等极端缺中子体系此前 2–5 个数量级的高估。

4.7 中子蒸发分道 $\sigma(x{\rm n})$ 截面

程序在总蒸发剩余截面之外，同时输出 1n–9n 各中子蒸发道的部分截面 $\sigma_{x\mathrm{n}}$，默认采用 3Wsur 的实现概率 prob_realiz(x)。各道之和在每个入射能量上严格等于总截面 $\sigma_{\mathrm{ER}}$，结果写入 sigma_XN.DAT（§6.3）。

4.7.1 默认分道算法（3Wsur prob_realiz）

子程序 WSURSHAPE 在能量扫描结束后、写出 sigma_XN.DAT 之前调用。对每个 $E_{\mathrm{c.m.}}$：

1. 取总截面 $\sigma_{\mathrm{ER}}$（含俘获、$P_{\mathrm{CN}}$ 与裂变存活）及 $E^*=Q_0+E_{\mathrm{c.m.}}$；中子分离能链 $S_{n,k}=S_n(Z_0,A_0-k)$ 由 SEPN 预先计算。
2. 复合核温度 $T=(1+\sqrt{1+4aE^*})/(2a)$，能级密度参数取 3Wsur 默认形式（体积+表面+曲率，$r_0=1.16$）。
3. 3Wsur 多重数分布： $$P_{1}=\exp\!\left[-\frac{(E^*-S_{n,0}-2T)^2}{2\cdot2.5^2}\right]$$ $$P_{x}=I(D_N,2x-3)-I(D_M,2x-1)\quad(x>1),$$ $$D_N=\frac{E^*-\sum_{k<x}S_{n,k}}{T},\qquad D_M=\frac{E^*-\sum_{k\le x}S_{n,k}}{T},$$ 其中 $I(z,m)=P(m+1,z)$ 为正则化下不完全伽马函数（GAMMP）。
4. 归一化并生成分道截面： $$\sigma_{x}=\sigma_{\mathrm{ER}}\,\frac{P_x}{\sum_{y=1}^{9}P_y},\qquad \sum_{x\ge1}\sigma_x=\sigma_{\mathrm{ER}} .$$ ±0.05 不确定度带用 $\sigma_{\mathrm{ER,hi}}$、$\sigma_{\mathrm{ER,lo}}$ 按同一公式生成。

该算法是确定性的，不使用实验数据、不含可调参数；各道为峰位随 $x$ 右移、两翼单调的光滑窄峰，无双峰或高能回折。

4.7.2 输出与守恒

• sigma_XN.DAT 每行 Ecm, σ_ER, σ(1n), …, σ(7n)（mb），共 9 列；
• 每个能量点严格满足 $\sum_{x\ge1}\sigma(x{\rm n})=\sigma_{\mathrm{ER}}$；
• 在线程序 ECC3.0.33.html 采用同一默认算法：由 WN-only 级联得到 $\sigma_{\mathrm{ER}}$（及 ±0.05 带），再调用 wsurShape9 生成各道。

5. 参数取值

5.1 可调参数（当前取值）

当前 ecc3_params.dat 取值为（ECC3 v3.0.30 起未变）：

5.2 固定参数

5.3 系统计算结果（60 个反应体系，ECC3 v3.0.33）

数据集：doc/ER-data.json 中全部 60 个反应、412 个实验点。ECC3.0.33（源码基于 3.0.32）的原生扫描窗为 $V_B+60$ MeV，其中 391 点参与统计；为与旧版可比，同时给出 $V_B+30$ 窗口（361 点）结果。参数取 §5.1 的值。

与 ECC3.0.30 的可比对比（同为 $V_B+30$ 窗口、361 点）：全部 $0.2759\to$ 0.2724；$Z_0<104$ $0.2410\to$ 0.2334；冷熔合基本不变（$\approx0.327$）。v3.0.33 的规则 5 显著改善了极端缺中子体系。

对照：ECC3.0.29 为 hot 0.2674 / cold 0.8849 / all 0.4265；ECC3.0.27 为 hot 0.2973 / all 0.4429。

6. 数据文件与输入输出

6.1 数据文件

6.2 输入

标准输入一行：A1, Z1, A2, Z2；参数文件 ecc3_params.dat 启动时读取。

6.3 输出

7. 编译与运行

# Linux / gfortran
gfortran -O2 -o ECC3.exe ECC3.0.33.for
./ECC3.exe

# Linux / Intel（结果与 gfortran 完全一致；必须带 -fp-model=precise）
ifx -O2 -fp-model=precise -o ECC3.exe ECC3.0.33.for

# Windows / Intel（源码文件按实际版本号命名，如 ECC3.0.33.for）
ifx  /O2 /fp:precise ECC3.0.33.for
ifort /O2 /fp:precise ECC3.0.33.for

# 输入反应体系，例如： 48 20 248 96
# 生成 sigma_cap.DAT、sigma_ER.DAT 与 sigma_XN.DAT

8. ECC3 与 HIVAP 的主要差异

一句话总结：ECC3 是 HIVAP 的"降维 + 补物理"实现——用 ECC2 解析势垒分布代替光学模型、用 WS4 全表 + 液滴壳修正代替人工核数据、用 6 参数代替 ~20 参数；同时用 $P_{\mathrm{CN}}$ 经验因子 补上超重核区准裂变抑制、用 ${B_f}^{\min}$ 补上缺中子复合核位垒欠收缩问题、用 $x>1$ 门控补上最重核（$Z\ge116$）鞍点能级密度的额外增强。两者物理内核一致，ECC3 换来了全局普适性与可拟合性，代价是丢失逐体系微观细节。

8.1 ECC3.0.33 与 3Wsur统计模型 分道截面计算的异同

8.2 ECC3.0.33 与 HIVAP 分道截面处理方法的不同

两者都把 ER 过程分为俘获 → 复合核形成 → 统计衰变存活三步，总截面 $\sigma_{\mathrm{ER}}=\sigma_{\mathrm{capture}}\,P_{\mathrm{CN}}\,W_{\mathrm{sur}}$ 相同；但在"如何把总存活分解到各 xn 分道"这一层，实现根本不同：HIVAP 在 $(Z,N)$ 二维布居网格上做逐级五道竞争的主方程演化，每个残核网格点就是一个 xn/pxn/αxn 道；ECC3.0.33 只沿一维中子链反推总存活 $\sigma_{\mathrm{ER}}$，再用 Wsur 解析多重数（prob_realiz）给定各 xn 道道形。

1. HIVAP：二维 $(Z,N)$ 布居 + 逐级 $\gamma/n/p/\alpha/$裂变五道竞争
   • 复合核 $(E,J)$ 布居 POPN 归一到融合截面（hivap.f:637-638）。
   • 每步按 $\gamma/n/p/\alpha/$裂变分流，分支比 $=\Gamma_i/\sum_j\Gamma_j$。
   • 分道截面 = 该残核网格点的布居；经 LINKFUS 乘 Wong 俘获截面后写出 sigxpn_fus.dat（列：Ecm, σ, $Z_p$, $N_n$）。即 xn 道 = $Z_{\mathrm{CN}}$、$N_{\mathrm{CN}}{-}x$ 网格点；pxn/αxn 是其他 $Z$ 网格点。
   • 全 $(E,J)$ 分辨：扣 yrast/转动能 $U=E-YJ$、移除 yrast 以下、发射后移能、布居截断 CUTPOP（eva.f:229、cutoff.f）；能级密度 Reisdorf 体积+表面+曲率（seqin.f:98-99，$r_0=1.12$）。
   • 闭合：$\sum_{\text{道}}\sigma=\sigma_{\mathrm{ER}}$（例 16O+197Au, Ecm=100: $\sum=111.37$ mb）。
2. ECC3.0.33：一维中子链 + Wsur 解析多重数
   • 只沿中子链（$Z$ 固定 $=Z_{\mathrm{CN}}$）反推总存活 $W_N$，得到 $\sigma_{\mathrm{ER}}=\sigma_{\mathrm{capture}}P_{\mathrm{CN}}W_N$；各 xn 道由 3Wsur 实现概率按 $\sigma_x=\sigma_{\mathrm{ER}}P_x/\sum_y P_y$ 生成（§4.7.1）。
   • 裂变只作整体存活 $W_N$；$p/\alpha/\gamma$ 只作宽度竞争进入 $W_N$，不产生 pxn/αxn 道，也不改变核种。
   • 当前默认：$\sigma_x=\sigma_{\mathrm{ER}}\,P_{xn}/\sum_y P_y$，其中 $P_{xn}$ 为 Wsur 解析实现概率（$x=1$ 高斯、$x>1$ Pearson 不完全伽马），道形来自"多重数分布假设"，总 $\sigma_{\mathrm{ER}}$ 才来自级联（WSURSHAPE）。
   • 能级密度 $a=A/[6(1+x)]$；分道道形用 Wsur 的 $a$（$r_0=1.16$）。

本质区别与后果：ECC3 的存活递推是"HIVAP 型"的，但 HIVAP 用全维度布居演化得到分道，ECC3 用一维中子链反推总存活，再用 Wsur 解析多重数定道形；因此 ECC3 的分道不含逐道的 $p/\alpha/\gamma/$裂变竞争——这是与 HIVAP 最本质的不同。HIVAP 更严格（Hauser–Feshbach 型主方程，含全部出射道与逐道裂变）；ECC3 更快、参数少，代价是在裂变强烈、$p/\alpha$ 出射显著的超重缺中子区会有系统差异。对"竞争如何进入分道"的对照：HIVAP 为 $P_{xn}\prod_i\Gamma_n/(\Gamma_n+\Gamma_f)$（且含 $p/\alpha/\gamma$）；Wsur 同式（§8.1）；ECC3.0.33 为 $\sigma_{\mathrm{ER}}P_{xn}/\sum P$（无逐级连乘）。

9. 进一步完善建议

1. $P_{\mathrm{CN}}$ 物理化：当前 $P_{0}$ 与能量项是经验 Fermi 函数，可尝试基于TDHF、ImQMD等微观动力学模型计算$P_{\mathrm{CN}}$，减少自由参数。
2. 误差加权拟合：当前为无权重 RMSD(log10)，可加入实验误差倒数加权（超重核不对称误差大）。
3. 分道截面的进一步物理化：当前 xn 分道由 3Wsur prob_realiz 生成（§4.7），已与 HIVAP 的 4n/5n 峰位一致；可进一步在分道解析中引入逐道裂变竞争（Hauser–Feshbach 分道种群演化）并加入 γ 与 pxn 分支，以获得更严格的中重核同位素产额。
4. 与 HIVAP 交叉验证：对中等体系（$P_{0}\approx1$ 区）同时跑 ECC3 与 HIVAP，确认无 $P_{\mathrm{CN}}$ 区域两者一致。

10. 相关链接