如何解决质子-氘离子大散射角碰撞的弛豫过快问题？

我基于Nanbu算法编写了Fortran代码模拟质子-氘离子碰撞，核心碰撞逻辑在calculations模块的collision子routine中实现带电粒子速度交换。目前已经修正了时间步长与库仑对数logCou的计算，但模拟结果存在两个明显不符合物理预期的问题：

• 粒子速度弛豫速度远快于理论预测值
• 计算得到的散射角χ普遍过大，违背了离子-离子碰撞散射角应较小的物理规律

核心代码片段（collision子routine） #

subroutine collision(particles, logCou, dt)
    use particle_mod
    implicit none
    type(Particle), intent(inout) :: particles(:)
    real(8), intent(in) :: logCou, dt
    real(8) :: v_rel, b_max, b, chi, mu, q1q2, pre_factor
    integer :: i, j

    q1q2 = particles(1)%q * particles(2)%q  ! 质子与氘核的电荷乘积
    mu = (particles(1)%m * particles(2)%m) / (particles(1)%m + particles(2)%m)  ! 约化质量

    do i = 1, size(particles)-1
        do j = i+1, size(particles)
            ! 计算相对速度
            v_rel = sqrt( (particles(i)%vx - particles(j)%vx)**2 + &
                          (particles(i)%vy - particles(j)%vy)**2 + &
                          (particles(i)%vz - particles(j)%vz)**2 )
            
            ! 计算最大碰撞参数b_max
            b_max = (q1q2) / (4.0d0 * acos(-1.0d0) * epsilon0 * mu * v_rel**2) * logCou
            
            ! 随机采样碰撞参数
            call random_number(b)
            b = b * b_max

            ! 计算散射角χ
            chi = 2.0d0 * atan( (q1q2) / (4.0d0 * acos(-1.0d0) * epsilon0 * b * mu * v_rel**2) )

            ! 速度变换与更新（Nanbu算法核心步骤）
            ! ... 此处为速度分量旋转与交换的代码
            ! 示例：将相对速度旋转χ角后更新粒子速度
            ! ...
        end do
    end do
end subroutine collision

结果示例 #

图中可见散射角峰值集中在远大于预期的区间，同时粒子速度分布在极短时间内就趋近于平衡态。

已排查的点 #

1. 确认时间步长dt的取值满足碰撞频率的约束（dt << 1/ν_coll）
2. 库仑对数logCou采用了等离子体物理中的标准公式计算（log(λ)，λ为德拜长度与b_max的比值）
3. 粒子的质量、电荷参数赋值正确（质子m=1.67e-27kg，q=1.6e-19C；氘核m=2*质子质量，q=1.6e-19C）

希望能得到关于代码实现错误或Nanbu算法应用误区的排查方向，谢谢。

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