Fast multipole expansions and error

For a cell $C$ let $h(C)$ be its remaining tree height, and for a pair $(C,D)$ use $h(C)+h(D)$. Every nonterminal call split by [F1] selects a nonleaf cell (both coordinates when $C=D$); each child call strictly lowers the height of every split coordinate and leaves any unsplit coordinate fixed. The measure is a nonnegative integer, so every branch terminates. The recursion is finitely branching (at most sixteen child pairs at a simultaneous equal-cell split), hence its finite-depth recursion tree is finite.

If the source cell is split, its child label sets are disjoint and cover the source labels, so the child Cartesian products are a disjoint cover of $I_C\times I_D$. The same statement holds when the target cell is split. When $C=D$ is internal, splitting both ordered coordinates gives all ordered child products, also disjoint. Thus each recursive replacement preserves exactly the parent pair set.

A terminal far cell pair records its entire Cartesian product $I_C\times I_D$. A terminal near pair of leaves records exactly that product with the diagonal labels removed. By step 2.1, induction from the root pair shows these terminal records cover the off-diagonal part of the root product without overlap.

A far cell pair has disjoint circumscribed disks because $R_{CD}>a_C+a_D$; hence its label sets cannot intersect and its product contains no diagonal pair. Removing the diagonal only at terminal leaf pairs therefore leaves exactly all ordered pairs $(j,i)$ with $j\mathrm{\ne }i$, each once.

Let ${\mathcal{P}}_{\mathrm{\ell }}$ be the homogeneous degree-$\mathrm{\ell }$ polynomials in three variables; it has $\left(\genfrac{}{}{0px}{}{\mathrm{\ell }+2}{2}\right)$ monomials. For a homogeneous polynomial $g$ of degree $d$, direct differentiation and Euler's identity give $\mathrm{\Delta }(\mathrm{\mid }x{\mathrm{\mid }}^{2s}g)=2s(2d+2s+1)\mathrm{\mid }x{\mathrm{\mid }}^{2s-2}g+\mathrm{\mid }x{\mathrm{\mid }}^{2s}\mathrm{\Delta }g$. If $\mathrm{\ell }\ge 2$, put $n=\mathrm{\ell }-2$, $J=\lfloor n\mathrm{/}2\rfloor$, $A_j=2(j+1)(2n-2j+3)$, $c_0=1\mathrm{/}{A}_0$, and $c_j=-c_{j-1}\mathrm{/}{A}_j$ for $1\le j\le J$. For $f\in {\mathcal{P}}_n$, the degree-$\mathrm{\ell }$ polynomial $q={\sum }_{j=0}^J{c}_j\mathrm{\mid }x{\mathrm{\mid }}^{2(j+1)}{\mathrm{\Delta }}^{j}f$ satisfies $\mathrm{\Delta }q=f$: the coefficients of $\mathrm{\mid }x{\mathrm{\mid }}^{2j}{\mathrm{\Delta }}^{j}f$ telescope by $c_0{A}_0=1$ and $c_j{A}_j+c_{j-1}=0$, and ${\mathrm{\Delta }}^{J+1}f=0$. Hence $\mathrm{\Delta }:{\mathcal{P}}_{\mathrm{\ell }}\to {\mathcal{P}}_{\mathrm{\ell }-2}$ is onto; for $\mathrm{\ell }=0,1$ the target is zero.

Restriction of a homogeneous polynomial to ${\mathbb{S}}^2$ is injective, because homogeneity recovers its value at every nonzero $x$ from its value at $x\mathrm{/}\mathrm{\mid }x\mathrm{\mid }$. Rank-nullity therefore gives $\dim {\mathcal{H}}_{\mathrm{\ell }}=\left(\genfrac{}{}{0px}{}{\mathrm{\ell }+2}{2}\right)-\left(\genfrac{}{}{0px}{}{\mathrm{\ell }}{2}\right)=2\mathrm{\ell }+1$ for $\mathrm{\ell }\ge 2$, and direct counting gives $1$ and $3$ for $\mathrm{\ell }=0,1$. Gram-Schmidt on any finite spanning list gives an orthonormal basis; no choice axiom is used.

For an orthonormal basis define $K_{\mathrm{\ell }}(\omega ,{\omega }^{\mathrm{\prime }})={\sum }_m{Y}_{\mathrm{\ell }m}(\omega )\overline{Y_{\mathrm{\ell }m}({\omega }^{\mathrm{\prime }})}$. A unitary basis change leaves this sum unchanged, and rotations carry one orthonormal basis of ${\mathcal{H}}_{\mathrm{\ell }}$ to another, so $K_{\mathrm{\ell }}$ is rotation-invariant. Thus $K_{\mathrm{\ell }}(\omega ,\omega )$ is constant in $\omega$; integrating it over the sphere gives $2\mathrm{\ell }+1$, so this constant is $(2\mathrm{\ell }+1)\mathrm{/}(4\pi )$.

Fix ${\omega }^{\mathrm{\prime }}$. After rotating it to the third coordinate axis, use $w=x_1+i{x}_2$ and $\bar{w}=x_1-i{x}_2$. A rotation sends these to $e^{i\theta }w$ and $e^{-i\theta }\bar{w}$, so an invariant monomial $w^a{\bar{w}}^b{x}_3^c$ must have $a=b$. The invariant homogeneous polynomials are therefore sums of $(x_1^2+x_2^2{)}^j{x}_3^{\mathrm{\ell }-2j}$, whose restrictions to the sphere are polynomials $f(t)$ of degree at most $\mathrm{\ell }$, with $t=x_3$. The polar-coordinate Laplacian is ${\mathrm{\partial }}_{rr}+2r^{-1}{\mathrm{\partial }}_r+r^{-2}{\mathrm{\Delta }}_{{\mathbb{S}}^2}$; applying it to $r^{\mathrm{\ell }}f(\cos \theta )$ gives ${\mathrm{\Delta }}_{{\mathbb{S}}^2}f=-\mathrm{\ell }(\mathrm{\ell }+1)f$. For a zonal $f$, the angular term is $(\sin \theta {)}^{-1}{\mathrm{\partial }}_{\theta }(\sin \theta {\mathrm{\partial }}_{\theta }f)=(1-t^2)f^{\mathrm{\prime }\mathrm{\prime }}-2t{f}^{\mathrm{\prime }}$ by $t=\cos \theta$, so $(1-t^2)f^{\mathrm{\prime }\mathrm{\prime }}-2t{f}^{\mathrm{\prime }}+\mathrm{\ell }(\mathrm{\ell }+1)f=0$. Writing $f(t)={\sum }_{k=0}^{\mathrm{\ell }}{a}_k{t}^k$ yields $(k+2)(k+1)a_{k+2}+(\mathrm{\ell }(\mathrm{\ell }+1)-k(k+1))a_k=0$. For $\mathrm{\ell }\ge 1$, the equation at $k=\mathrm{\ell }-1$, using $a_{\mathrm{\ell }+1}=0$, gives $2\mathrm{\ell }{a}_{\mathrm{\ell }-1}=0$; descending through $k=\mathrm{\ell }-3,\mathrm{\ell }-5,\dots$ forces every coefficient of parity opposite to $\mathrm{\ell }$ to vanish. For $k=\mathrm{\ell }-2,\mathrm{\ell }-4,\dots$, the coefficient $\mathrm{\ell }(\mathrm{\ell }+1)-k(k+1)$ is positive, so each lower coefficient of parity $\mathrm{\ell }$ is determined by the next one. Thus the zonal solution space has dimension at most one (for $\mathrm{\ell }=0$ it is the constants). The function $K_{\mathrm{\ell }}(\cdot ,{\omega }^{\mathrm{\prime }})$ is a zonal harmonic, and $K_{\mathrm{\ell }}({\omega }^{\mathrm{\prime }},{\omega }^{\mathrm{\prime }})=(2\mathrm{\ell }+1)\mathrm{/}(4\pi )>0$, so a nonzero such solution exists. Define $P_{\mathrm{\ell }}(t)=K_{\mathrm{\ell }}(\omega ,{\omega }^{\mathrm{\prime }})\mathrm{/}{K}_{\mathrm{\ell }}({\omega }^{\mathrm{\prime }},{\omega }^{\mathrm{\prime }})$ when $t=\omega \cdot {\omega }^{\mathrm{\prime }}$. Rotation invariance makes this depend only on $t$, and $P_{\mathrm{\ell }}(1)=1$. The diagonal value in step 3.1 now gives the displayed addition theorem.

Cauchy-Schwarz applied to the finite basis sum gives $\mathrm{\mid }{K}_{\mathrm{\ell }}(\omega ,{\omega }^{\mathrm{\prime }})\mathrm{\mid }\le \sqrt{K_{\mathrm{\ell }}(\omega ,\omega )K_{\mathrm{\ell }}({\omega }^{\mathrm{\prime }},{\omega }^{\mathrm{\prime }})}=(2\mathrm{\ell }+1)\mathrm{/}(4\pi )$. Dividing by the diagonal constant proves $\mathrm{\mid }{P}_{\mathrm{\ell }}(t)\mathrm{\mid }\le 1$.

Fix ${\omega }^{\mathrm{\prime }}\in {\mathbb{S}}^2$ and let $u(x)=\mathrm{\mid }{\omega }^{\mathrm{\prime }}-x{\mathrm{\mid }}^{-1}$ for $\mathrm{\mid }x\mathrm{\mid }<1$. The radial formula $\mathrm{\Delta }\mathrm{\mid }x-{\omega }^{\mathrm{\prime }}{\mathrm{\mid }}^{-1}=f^{\mathrm{\prime }\mathrm{\prime }}(r)+2f^{\mathrm{\prime }}(r)\mathrm{/}r=0$ away from ${\omega }^{\mathrm{\prime }}$ shows $u$ is harmonic in the unit ball. Its homogeneous Taylor term of degree $\mathrm{\ell }$ is therefore harmonic, since the homogeneous terms in $\mathrm{\Delta }u=0$ of each degree vanish separately; it is also invariant under rotations fixing ${\omega }^{\mathrm{\prime }}$. By step 4.1 it has the form $b_{\mathrm{\ell }}{r}^{\mathrm{\ell }}{P}_{\mathrm{\ell }}(\omega \cdot {\omega }^{\mathrm{\prime }})$. On the ray $x=r{\omega }^{\mathrm{\prime }}$, $u(r{\omega }^{\mathrm{\prime }})=(1-r{)}^{-1}={\sum }_{\mathrm{\ell }\ge 0}{r}^{\mathrm{\ell }}$, so coefficient comparison gives $b_{\mathrm{\ell }}=1$. Finally $u(r\omega )=(1-2r(\omega \cdot {\omega }^{\mathrm{\prime }})+r^2{)}^{-1\mathrm{/}2}$, whose Taylor series in $r$ has radius at least $1$: its possible denominator zeros are $t\pm i\sqrt{1-t^2}$, both of modulus one. This proves the generating identity for $0\le r<1$.

For $\mathrm{\mid }w\mathrm{\mid }<1$, the geometric series $1\mathrm{/}(1+t)={\sum }_{r=0}^{\mathrm{\infty }}(-t{)}^r$ converges uniformly on each closed disk $\mathrm{\mid }t\mathrm{\mid }\le r_0<1$. Integrating it from $0$ to $w$ gives the analytic branch identity $\log (1+w)={\sum }_{m=1}^{\mathrm{\infty }}(-1{)}^{m+1}{w}^m\mathrm{/}m$.

For a nonempty child cell $E$ of $C$, write $z_j-c_C=(z_j-c_E)+(c_E-c_C)$ and apply the finite binomial theorem to its $n$th power. Summing over the sources in each nonempty child and then over the nonempty children gives $M_n(C)={\sum }_E{\sum }_{r=0}^n\left(\genfrac{}{}{0px}{}{n}{r}\right)(c_E-c_C{)}^{n-r}{M}_r(E)$ for every $n\le p$; empty children contribute zero. This is exactly the M2M pass.

If a parent local polynomial is ${\sum }_{k=0}^p{L}_k(z-c_D{)}^k$, put $z-c_D=(z-c_E)+(c_E-c_D)$ and expand each power. The coefficient of $(z-c_E{)}^j$ is ${\sum }_{k=j}^p\left(\genfrac{}{}{0px}{}{k}{j}\right)L_k(c_E-c_D{)}^{k-j}$; hence L2L exactly re-centers the polynomial without adding truncation error.

For a far pair $(C,D)$, put $z-c_C=D_0+x$ with $D_0=c_D-c_C$ and $\mathrm{\mid }x\mathrm{\mid }\le a_D$. Then $\mathrm{\mid }{D}_0+x\mathrm{\mid }\ge R-a_D>a_C$. Substituting $w=-(z_j-c_C)\mathrm{/}(z-c_C)$ into the branch-normalized series from step 1.1 and summing the finite source set gives $M_0(C)L_{CD}(z-c_C)-{\sum }_{n\ge 1}{M}_n(C)\mathrm{/}[n(z-c_C{)}^n]$. The branch is $L_{CD}$ because both $z-c_C$ and $z-z_j$ lie in its disk and the difference of their logarithms is the series that vanishes at $z_j=c_C$. The convergence is absolute because $\mathrm{\mid }{z}_j-c_C\mathrm{\mid }\mathrm{/}\mathrm{\mid }z-c_C\mathrm{\mid }<1$.

For a far pair put $D=c_D-c_C$, $x=z-c_D$ and $y=z_j-c_C$. By [F1], $\mathrm{\mid }x-y\mathrm{\mid }\le \mathrm{\mid }x\mathrm{\mid }+\mathrm{\mid }y\mathrm{\mid }<R=\mathrm{\mid }D\mathrm{\mid }$. The branch-normalized logarithm therefore has expansion $L_{CD}(D+x-y)=L_{CD}(D)+{\sum }_{m\ge 1}(-1{)}^{m+1}(x-y{)}^m\mathrm{/}(m{D}^m)$.

Expand $(x-y{)}^m$ by the binomial theorem and set $m=n+k$. For $n\ge 1,k\ge 0$, the coefficient of $x^k{y}^n$ is $(-1{)}^{k+1}\left(\genfrac{}{}{0px}{}{n+k-1}{k}\right)\mathrm{/}(n{D}^{n+k})$, since $\left(\genfrac{}{}{0px}{}{n+k}{k}\right)\mathrm{/}(n+k)=\left(\genfrac{}{}{0px}{}{n+k-1}{k}\right)\mathrm{/}n$. The terms with $n=0$ are $M_0{L}_{CD}(D)$ for $k=0$ and $(-1{)}^{k+1}{M}_0\mathrm{/}(k{D}^k)$ for $k\ge 1$; the terms with $k=0,n\ge 1$ are $-M_n\mathrm{/}(n{D}^n)$. Summing $y^n$ over sources replaces it by $M_n$, and retaining $0\le n,k\le p$ yields precisely the M2L formulas in [F2].

Substitution of a target point into the re-centered degree-$p$ polynomial is exactly L2P, and termwise differentiation gives its derivative formula. These are finite polynomial identities and require no additional convergence assumption.

At each nonempty leaf, P2M accumulates $p+1$ moments per resident source, so the total is at most $N(p+1)$ operations up to a constant. The nonempty cells form a rooted connected subtree, so enumerating them by inspecting four children per nonempty internal cell costs $O(T)$ and yields exactly $T-1$ parent-child edges whose two cells are nonempty. On each such edge, M2M or L2L evaluates triangular sums with ${\sum }_{n=0}^p(n+1)=(p+1)(p+2)\mathrm{/}2$ terms, hence at most $O((p+1{)}^2)$ operations. Empty-child moments are zero and their edges are skipped. Initializing the local vectors costs $O(T(p+1))$, and both this and the $O(T)$ enumeration cost are absorbed in the displayed bound.

For each of the $K$ far pairs, M2L forms $p+1$ coefficients, each using at most $p+1$ terms, and evaluates the fixed logarithm once at the pair center; this costs $O(K(p+1{)}^2)$ because $p\ge 1$. L2P evaluates one degree-$p$ polynomial and its derivative at each of $N$ targets in $O(N(p+1))$ operations.

The recursive list decomposition uses $V$ cell-pair inspections. Each of the $J$ direct near particle-pair evaluations uses a constant number of arithmetic operations and one fixed-logarithm evaluation. By [F3] no required pair is counted twice or omitted; the work to dispatch the finite lists is included in these counts.

Adding the counts from steps 1.1, 1.2, and 1.3 gives the displayed operation bound after increasing one absolute constant $C_0$. The assumptions $T,V\le CN$, $K\le CT$, and $J\le CN$ then give $O(N(p+1{)}^2)$; since $p\ge 1$, $(p+1{)}^2\le 4p^2$, giving $O(N{p}^2)$. No bound on the cost of constructing the input quadtree or assigning its points was used.

Fix a far interaction $(C,D)$ contributing to target $i$, and write $D_0=c_D-c_C$, $x=z_i-c_D$, $y_j=z_j-c_C$, and $R=\mathrm{\mid }{D}_0\mathrm{\mid }$. Then $\mathrm{\mid }x\mathrm{\mid }\le a_D$, $\mathrm{\mid }{y}_j\mathrm{\mid }\le a_C$, and $\rho =(a_C+a_D)\mathrm{/}R<1$ by [F1]. The series in [F3] for each source is the double binomial expansion of $L_{CD}(D_0+x-y_j)$. The algorithm retains all terms with source index $n\le p$ and target degree $k\le p$; every omitted term has total degree $m=n+k\ge p+1$.

For total degree $m$, the sum of absolute values of the potential terms for one source is at most $\mathrm{\mid }{q}_j\mathrm{\mid }{\sum }_{k=0}^m\left(\genfrac{}{}{0px}{}{m}{k}\right)\mathrm{\mid }x{\mathrm{\mid }}^k\mathrm{\mid }{y}_j{\mathrm{\mid }}^{m-k}\mathrm{/}(m{R}^m)=\mathrm{\mid }{q}_j\mathrm{\mid }(\mathrm{\mid }x\mathrm{\mid }+\mathrm{\mid }{y}_j\mathrm{\mid }{)}^m\mathrm{/}(m{R}^m)\le \mathrm{\mid }{q}_j\mathrm{\mid }{\rho }^m\mathrm{/}m$. Summing sources and omitted degrees gives an interaction error at most $Q_C{\sum }_{m=p+1}^{\mathrm{\infty }}{\rho }^m\mathrm{/}m\le Q_C{\rho }^{p+1}\mathrm{/}((p+1)(1-\rho ))$, since $1\mathrm{/}m\le 1\mathrm{/}(p+1)$ for $m\ge p+1$.

For each total degree $m\ge 1$, the sum of absolute differentiated monomials for one source, before multiplying by $q_j$, is ${\sum }_{k=1}^m\frac{k}{m}\left(\genfrac{}{}{0px}{}{m}{k}\right)\frac{\mathrm{\mid }x{\mathrm{\mid }}^{k-1}\mathrm{\mid }{y}_j{\mathrm{\mid }}^{m-k}}{R^m}=\frac{(\mathrm{\mid }x\mathrm{\mid }+\mathrm{\mid }{y}_j\mathrm{\mid }{)}^{m-1}}{R^m}\le \frac{{\rho }^{m-1}}{R}\mathrm{.}$ Here a zeroth power is $1$, including at zero. Every subset of these monomials, in particular the omitted subset from the rectangular truncation $n,k\le p$, has absolute sum bounded by the same nonnegative full sum. The series and its differentiated monomials converge absolutely and uniformly on $\mathrm{\mid }x\mathrm{\mid }\le a_D$ because $\rho <1$, so termwise differentiation is valid. By step 1.1, omitted terms have $m\ge p+1$; multiplying by $\mathrm{\mid }{q}_j\mathrm{\mid }$ and summing sources and degrees bounds the interaction derivative error by $Q_C{R}^{-1}{\sum }_{m=p+1}^{\mathrm{\infty }}{\rho }^{m-1}=Q_C{R}^{-1}{\rho }^p\mathrm{/}(1-\rho )$. The derivative of any additive branch constant is zero.

By [F4], each off-diagonal particle pair is either one exact direct near interaction or belongs to exactly one far cell product. M2M and L2L introduce no additional truncation error by [F2]–[F3], and L2P evaluates the stored polynomial exactly. Therefore the total potential and derivative errors for target $i$ are bounded by the sums of the finite per-interaction bounds from steps 2.1 and 2.2; all near terms contribute zero error.

For a source $y_j$, put $v_j=y_j-c$, $r_j=\mathrm{\mid }{v}_j\mathrm{\mid }\le a$, and ${\omega }_j=v_j\mathrm{/}{r}_j$ when $r_j>0$. For a target put $r=\mathrm{\mid }x-c\mathrm{\mid }>a$ and $\omega =(x-c)\mathrm{/}r$. Applying the generating identity and then the addition identity in [F2] with ratio $r_j\mathrm{/}r<1$ gives $\frac{1}{\mathrm{\mid }x-y_j\mathrm{\mid }}={\sum }_{\mathrm{\ell }=0}^{\mathrm{\infty }}\frac{r_j^{\mathrm{\ell }}}{r^{\mathrm{\ell }+1}}{P}_{\mathrm{\ell }}(\omega \cdot {\omega }_j)=4\pi {\sum }_{\mathrm{\ell }=0}^{\mathrm{\infty }}\frac{1}{2\mathrm{\ell }+1}\frac{r_j^{\mathrm{\ell }}}{r^{\mathrm{\ell }+1}}{\sum }_{m=-\mathrm{\ell }}^{\mathrm{\ell }}{Y}_{\mathrm{\ell }m}(\omega )\overline{Y_{\mathrm{\ell }m}({\omega }_j)}$. For $r_j=0$, the same formula follows by the limit: the degree-zero term is $1\mathrm{/}r$ and all positive-degree regular solid harmonics vanish at the center.

The dimension result in [F2] gives ${\sum }_{\mathrm{\ell }=0}^p(2\mathrm{\ell }+1)=(p+1{)}^2$ coefficients. At degree $\mathrm{\ell }$, each homogeneous polynomial has at most $\left(\genfrac{}{}{0px}{}{\mathrm{\ell }+2}{2}\right)$ monomials and there are $2\mathrm{\ell }+1$ basis polynomials. Evaluating all basis values at one point by the prescribed monomial sums therefore uses at most a constant times ${\sum }_{\mathrm{\ell }=0}^p(2\mathrm{\ell }+1)\left(\genfrac{}{}{0px}{}{\mathrm{\ell }+2}{2}\right)=O((p+1{)}^4)$ operations.

Summing the finite source set with weights $q_j$ and grouping the finite sums over $m$ produces exactly the moments and the truncated evaluator in [F3]. The remainder magnitude is at most ${\sum }_j\mathrm{\mid }{q}_j\mathrm{\mid }{\sum }_{\mathrm{\ell }=p+1}^{\mathrm{\infty }}\frac{r_j^{\mathrm{\ell }}}{r^{\mathrm{\ell }+1}}\mathrm{\mid }{P}_{\mathrm{\ell }}(\omega \cdot {\omega }_j)\mathrm{\mid }\le \frac{Q}{r}{\sum }_{\mathrm{\ell }=p+1}^{\mathrm{\infty }}(a\mathrm{/}r{)}^{\mathrm{\ell }}=\frac{Q}{r}\frac{(a\mathrm{/}r{)}^{p+1}}{1-a\mathrm{/}r}$, using [F2] and the geometric-series formula.

If $\mathrm{\mid }x-d\mathrm{\mid }\le b$, the triangle inequality gives $r=\mathrm{\mid }x-c\mathrm{\mid }\ge R-b>a$. The function $a^{p+1}\mathrm{/}(r^{p+1}(r-a))$, equal to $r^{-1}(a\mathrm{/}r{)}^{p+1}\mathrm{/}(1-a\mathrm{/}r)$, is nonincreasing for $r>a$ because its denominator is nondecreasing; substituting $r=R-b$ in step 2.1 gives the stated uniform target-region bound.

Repeating this basis evaluation at $S$ source points and $T_0$ target points costs $O((S+T_0)(p+1{)}^4)$. The remaining moment and target accumulations use $(p+1{)}^2$ coefficients per point, hence $O((S+T_0)(p+1{)}^2)$ additional operations; if basis values are supplied, this latter count is the full coefficient-stage cost.