Summation by parts methods for the spherical harmonic decomposition of the wave equation in arbitrary dimensions

We investigate numerical methods for wave equations in $n+2$ spacetime dimensions, written in spherical coordinates, decomposed in spherical harmonics on $S^n$, and finite-differenced in the remaining coordinates $r$ and $t$. Such an approach is useful when the full physical problem has spherical symmetry, for perturbation theory about a spherical background, or in the presence of boundaries with spherical topology. The key numerical difficulty arises from lower-order $1/r$ terms at the origin $r=0$. As a toy model for this, we consider the flat space linear wave equation in the form $\dot\pi=\psi'+p\psi/r$, $\dot\psi=\pi'$, where $p=2l+n$, and $l$ is the leading spherical harmonic index. We propose a class of summation by parts (SBP) finite differencing methods that conserve a discrete energy up to boundary terms, thus guaranteeing stability and convergence in the energy norm. We explicitly construct SBP schemes that are second and fourth-order accurate at interior points and the symmetry boundary $r=0$, and first and second-order accurate at the outer boundary $r=R$.