大家好，我现在在实现库仑相互作用的快速多极子方法（FMM），参考的是Leslie Greengard的《A Short Course on the Fast Multipole Method》。我的目标是精确计算点电荷系统的电势，目前Particle-to-Multipole（P2M）和Multipole-to-Particle（M2P）模块都验证正常，但卡在Multipole-to-Multipole（M2M）平移步骤上了，想请大家帮忙找找问题。

目前工作正常的模块（P2M与M2P） #

这部分的实现已经通过直接求和验证完全正确：

• 电势计算基于公式：Phi(P) = sum_n=0^inf sum_m=-n^n (M_n^m / r^(n+1)) * Y_n^m(theta, phi)
• 球谐函数使用scipy.special.sph_harm(m, n, phi, theta)，返回正交归一化形式的$Y_n^m(\theta, \phi)$
• P2M计算多极矩$M_{P2M}$的方式：sum q_i * r_i^n * conj(Y_n^m(theta_i, phi_i))
• 电势评估时，针对正交归一化球谐函数，已正确应用(4 * pi / (2n + 1))修正因子

验证结果如下，子单元多极展开计算的电势与直接求和完全吻合：

Potential from child multipole: 3.275435467795e+00
Potential from direct sum: 3.275435467795e+00

遇到的问题（M2M平移） #

按照论文定理5.3实现M2M平移后，平移得到的父单元多极矩计算出的电势与直接结果偏差极大：

Potential from parent multipole: 1.123673185858e+00
Potential from direct sum: 3.275435467795e+00

同时，我直接用父单元包含的所有粒子通过P2M计算的多极矩，和子单元M2M平移得到的父单元多极矩数值完全不匹配，说明M2M的实现肯定存在错误。我反复检查了公式对应关系，但始终找不到问题点。

以下是我的代码实现，麻烦大家帮忙排查：

from scipy.special import sph_harm, factorial
import numpy as np

def A_nm(n, m):
    """Compute A_n^m coefficient as defined in Greengard's paper."""
    return (-1)**n / np.sqrt(factorial(n - abs(m)) * factorial(n + abs(m)))

def M2M_theorem53(child_multipole, child_center, parent_center, p):
    """
    Perform M2M translation based on Theorem 5.3 (Greengard, 1997).
    Translates multipole expansion from child center to parent center (origin).
    Parameters:
    child_multipole: np.array of shape (p+1, 2p+1), with child moments O_n^m at [n, m+n]
    child_center: Point object
    parent_center: Point object
    p: order of expansion
    Returns:
    parent_multipole: np.array of shape (p+1, 2p+1), with translated M_j^k
    """
    parent_multipole = np.zeros((p + 1, 2 * p + 1), dtype=complex)
    # Vector from parent to child center
    dx = child_center.x - parent_center.x
    dy = child_center.y - parent_center.y
    dz = child_center.z - parent_center.z
    rho = np.sqrt(dx**2 + dy**2 + dz**2)
    if rho < 1e-14:
        return np.copy(child_multipole)  # No translation needed
    alpha = np.arctan2(dy, dx)
    beta = np.arccos(dz / rho)

    for j in range(p + 1):
        for k in range(-j, j + 1):
            acc = 0.0j
            A_jk = A_nm(j, k)
            for n in range(j + 1):
                for m in range(-n, n + 1):
                    j_minus_n = j - n
                    k_minus_m = k - m
                    if abs(k_minus_m) > j_minus_n:
                        continue  # Spherical harmonic is 0 outside its domain
                    A_nm_val = A_nm(n, m)
                    A_jn_km = A_nm(j_minus_n, k_minus_m)
                    # Complex exponential factor: i^{|k| - |m| - |k - m|}
                    phase = 1j ** (abs(k) - abs(m) - abs(k - m))
                    # Y_n^{-m}(alpha, beta)
                    Y = sph_harm(-m, n, alpha, beta)
                    term = (child_multipole[j_minus_n, k_minus_m + j_minus_n] * phase * A_nm_val * A_jn_km * (rho ** n) * Y / A_jk)
                    acc += term
            parent_multipole[j, k + j] = acc
    return parent_multipole

import numpy as np
from scipy.special import sph_harm

# --- Point and Particle classes ---
class Point():
    """The class for a point."""
    def __init__(self, coords=[], domain=1.0):
        if coords:
            assert len(coords) == 3, "the size of coords should be 3."
            self.x = coords[0]
            self.y = coords[1]
            self.z = coords[2]
        else:
            self.x = domain * np.random.random()
            self.y = domain * np.random.random()
            self.z = domain * np.random.random()
    def distance(self, other):
        return np.sqrt((self.x-other.x)**2 + (self.y-other.y)**2 + (self.z-other.z)**2)

class Particle(Point):
    """The derived class for a particle."""
    def __init__(self, coords=[], domain=1.0, m=1.0):
        Point.__init__(self, coords, domain)
        self.m = m  # Mass of the particle (charge for potential)
        self.phi = 0.  # Gravitational potential, initialized for direct sum

def P2M(sources, center, p):
    # 注：该函数完整实现已验证正确，此处保留用户提供的定义开头
    pass

我已经核对了论文中M2M平移的所有公式细节，包括球谐函数的共轭、系数因子、相位项等，但就是找不到错误点。有没有熟悉FMM或者球谐函数平移的朋友能帮我指出问题所在呀？

内容来源于stack exchange