AN IMPROVED MULTIPOLE APPROXIMATION FOR SELF-GRAVITY AND ITS IMPORTANCE FOR CORE-COLLAPSE SUPERNOVA SIMULATIONS

Gravity is a key phenomenon in many astrophysical contexts, and, in particular, plays an essential role in the explosions of core-collapse supernovae (CCSNe). Accurate computation of self-gravity is therefore an important objective for astrophysical simulation codes. For self-gravitating Newtonian systems this requires solving Poisson's equation. Poisson's equation is an elliptic partial differential equation, which couples every part of the domain at each time step. The optimal solver strategy for an astrophysical Poisson problem in which gravity is coupled to a hydrodynamic flow depends on the typical mass configuration in the domain. Multigrid algorithms (e.g., Huang & Greengard 2000; Trottenberg et al. 2001; Ricker 2008) are popular for cosmological structure-formation and star formation simulations with Newtonian gravity (e.g., Yang et al. 2009; ZuHone et al. 2010; Latif et al. 2011; Federrath & Klessen 2012), since these algorithms work well with mass configurations spread out over a computational domain. In problems where a single, large condensation of mass arises, however, a multipole expansion using spherical harmonics is more appropriate. Solving Poisson's equation using spherical harmonic expansions is a common approach for computing the self-gravity of nearly-spherical mass distributions. Multipole approximations have been used in a number of astrophysical applications including N-body calculations (see Sellwood 1987) and grid-based hydrodynamics (Müller & Steinmetz 1995). The unstable collapse of the core of a massive star that precedes a CCSN is particularly sensitive to a highly dynamic gravitational potential. Many approaches have been adopted for treating self-gravity in CCSN simulations ranging from full general relativity (e.g., O'Connor & Ott 2010; Ott et al. 2007, 2013; Müller et al. 2010) to simplified one-dimensional (1D) "monopole" approximations (e.g., Hanke et al. 2012; Dolence et al. 2013; Couch 2013a).

The Newtonian potential of a spherically symmetric self-gravitating mass is trivial, of course, and is represented by the monopole term of the expansion. However, as departures from spherical symmetry accumulate, the mass distribution must be represented by an expansion of spherical harmonics beyond l = 0, the accuracy of which depends on the degree of non-sphericity of the mass distribution and the number of terms used in the expansion (Müller & Steinmetz 1995). Such multipole approaches for self-gravity have been used in a number of multidimensional CCSN simulations (e.g., Livne et al. 2004; Buras et al. 2006; Bruenn et al. 2013). Multipole approaches are suited for CCSNe because the gravitational potential is dominated by the monopole contribution, but the higher-order contributions due to significant non-spherical motions in the post-shock region can be important. Additionally, in non-spherical geometries wherein the central proto-neutron star is allowed to move, the physical kick imparted on the star by the requirement of momentum conservation—a model-constraining observable—is critically dependent on an accurate, momentum-conserving self-gravity computation (Wongwathanarat et al. 2010, 2012).

In this article, we investigate the multipole expansion approach to solving Poisson's equation numerically for the self-gravity of an approximately spherical mass distribution. We identify, and correct, two heretofore neglected sources of significant errors that arise in implementations of multipole self-gravity for non-spherical coordinates:

• 1.  
  The numerical approximation that effectively places grid cell mass at the central point of a computational cell, then computes the gravitational potential at that point, resulting in a convergence failure of the multipole expansion, so that larger choices of multipole cutoff value l_max actually make the potential computation less accurate.
• 2.  
  Sub-optimal choice of location for the expansion center, which results in angular power at high multipole l values in the gravitational field, requiring a high—and expensive—value of l_max.

We show here that source 1 of error can be eliminated by a collocation scheme that effectively staggers point mass placement and potential computation; and source 2 of error can be minimized by a careful, unique choice of expansion center, which we derive. We demonstrate that CCSN simulations are particularly sensitive to these details of the multipole approximation for gravity and show that our improvements result in dramatic improvements in important metrics such as momentum conservation and convergence with number of terms in the multipole expansion.

This paper is organized as follows. In Section 2.1 we briefly present the discretized multipole equations and exhibit the intrinsic error that can afflict solutions to these equations due to the singularity in the Green's function of Poisson's equation. We show that this error is eliminated by computing gravitational potentials on cell faces rather than at cell centers. In Section 2.2 we derive the optimal location for centering the expansion for general mass distributions, based on the minimization of an angular "spectral compactness" measure that characterizes the extent in l-space of the spherical-harmonic spectrum. In Section 3 we describe our implementation of fast, efficient multipole gravity solver in the FLASH simulation framework. We test our new solver, which includes the improvements we discuss, on static potentials in Section 4 and exhibit the effects of the errors described above, as well as the result of their correction. In Section 5 we test our new implementation with highly dynamical CCSN simulations in two dimensions and show that the results are highly sensitive to the centering of the multipole expansion and to the relative collocation of the mass and potential evaluation points. We discuss our conclusions in Section 6.

2.1. The Self-potential Error

The gravitational potential of an isolated distribution of mass with density ρ(x) is given by the well-known Green's function of the Poisson equation

$$\begin{equation} \Phi (\mathbf {x})=-G\int d^{3}\mathbf {x}^{\prime }\frac{\rho (\mathbf {x}^{\prime })}{\left|\mathbf {x}-\mathbf {x}^{\prime } \right|}. \end{equation} \tag{ 1 }$$

Direct numerical implementation of this formula in a simulation is inefficient, often necessitating approximate approaches. For mass distributions that can be described as spherical to lowest-order, multipole expansions of Equation (1) can be used to efficiently compute solutions. The multipole expansion version of the potential is given by the equally well-known formula

$$\begin{eqnarray} \Phi (\mathbf {x})&=&{-}G\sum _{l=0}^{\infty }\sum _{m=-l}^{l}\frac{4\pi }{2l+1} \nonumber \\ &&\times \int d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x} ^{\prime })\, Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n}^{\prime })^{*}g_{l}(r,r^{\prime }), \end{eqnarray} \tag{ 2 }$$

where r ≡ |x|, n ≡ x/|x|, and

$$\begin{equation} g_{l}(r,r^{\prime })\equiv \Theta (r-r^{\prime })\frac{r^{\prime l}}{r^{l+1}}+\Theta (r^{\prime }-r)\frac{r^{l}}{r^{\prime l +1}}. \end{equation} \tag{ 3 }$$

Here, Θ(x) is the usual Heaviside function.

In Eulerian hydrodynamic codes, a standard discretization strategy for this expansion (Müller & Steinmetz 1995) begins with a subdivision of the domain into N_R spherical shells bounded by radii R_t, t = 1, ..., N_R, chosen to suit the problem (and not necessarily uniformly spaced). A cell centered at the position x_q is ascribed a radius r_q that is the mean radius of the spherical shell containing x_q, where q is an index running over mesh cells. The discretized potential is then computed as

$$\begin{eqnarray} \Phi (\mathbf {x}_{q})&=&-G\sum _{l=0}^{l_{\rm max}}\sum _{m=-l}^{l}\frac{4\pi }{2l+1}Y_{lm}(\mathbf {n}_{q}) \sum _{q^{\prime }} \Delta _{q^{\prime }}^{3}\rho (\mathbf {x}_{q^{\prime }})Y_{lm}(\mathbf {n}_{q^{\prime }}^{\prime })^{*} \nonumber \\ &&\times \left\lbrace \tilde{\Theta } _{qq^{\prime }}\frac{r_{q^{\prime }}^{l}}{r_{q}^{l+1}}+ \tilde{\Theta }_{q^{\prime }q}\frac{r_{q}^{l}}{r_{q^{\prime }}^{l+1}}+ \tilde{\delta }_{qq^{\prime }}\frac{1}{r_{q}}\right\rbrace , \end{eqnarray} \tag{ 4 }$$

where l_max is some chosen cutoff value for the expansion, $\Delta _{q^{\prime }}^{3}$ is the volume of the cell centered at $\mathbf {x}_{q^{\prime }}$, and where

$$\tilde{\Theta }_{qq^{\prime }}\equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & r_{q}>r_{q^{\prime }}\\ 0 & r_{q}\le r_{q^{\prime }} \end{array}\right.;\quad \quad \tilde{\delta }_{qq^{\prime }}\equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & r_{q}=r_{q^{\prime }}\\ 0 & r_{q}\ne r_{q^{\prime }} \end{array}\right..$$

If one were directly implementing the potential using the expression of Equation (1), discretization in the presence of the singular Green's function |x − x'|⁻¹ might give rise to misgivings having to do with the delicate handling of gravitational self-interaction within a mesh cell. This issue of self-gravity appears superficially to magically cure itself in the passage to the discrete multipole expansion of Equation (4), wherein no short-distance singularities are explicitly visible. This miracle cure is illusory, unfortunately: the singularity still lurks in the expression, and manifests itself in the failure of the self-interaction terms in the expression to converge as l_max → ∞.

To see this, consider the self-interacting term q' = q in Equation (4):

$$\begin{eqnarray} \Phi _{\rm Self}(\mathbf {x}_{q}) &\equiv & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}\sum _{l=0}^{l_{\rm max}}\frac{4\pi }{2l+1} \sum _{m=-l}^{l}Y_{lm}(\mathbf {n}_{q})Y_{lm}(\mathbf {n}_{q})^{*}.\nonumber \\ \end{eqnarray} \tag{ 5 }$$

The addition theorem of spherical harmonics states that

$$\begin{equation} \sum _{m=-l}^{l}Y_{lm}(\mathbf {n}_{1})Y_{lm}(\mathbf {n}_{2})^{*}= \frac{2l+1}{4\pi }P_{l}(\mathbf {n}_{1}\cdot \mathbf {n}_{2}), \end{equation} \tag{ 6 }$$

where P_l(x) is a Legendre polynomial. We therefore have that

$$\begin{eqnarray} \Phi _{\rm Self}(\mathbf {x}_{q}) & \equiv & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}\sum _{l=0}^{l_{\rm max}}P_{l} (1)\nonumber \\ & = & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}\sum _{l=0}^{l_{\rm max}}1\nonumber \\ & = & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}\times (l_{\rm max}+1). \end{eqnarray} \tag{ 7 }$$

It follows that the discrete expression for Φ_Self(x_q) is not convergent with multipole order, and that the accuracy of the discrete scheme described above cannot be improved by increasing l_max. We note that Sellwood (1987) remarked upon related difficulties in the context of N-body simulations, but did not give the explicit form of this self-potential error nor expound on its origins in the discrete multipole expansion. In numerical simulations, this pathology manifests itself as a dramatic failure in accuracy of the potential calculation, which gets worse with increasing l_max. We also note that due to the factor r_q in the denominator of Equation (7) this error is worst near the origin of the multipole expansion. This error is also larger for computational zones containing large masses, $\Delta _{q}^{3}\rho (\mathbf {x}_{q})$. Both of this conditions are met in the extreme for CCSN simulations containing a proto-neutron star near the center of the domain.

A more deft handling of self-interaction is required if the scheme is to be rescued. We may begin by observing that the physical origin of the difficulty is that the scheme in effect treats all masses as points at the cell centers, then computes potentials at those same cell centers. If the points of potential computation were offset from the cell centers, the problem would go away. This is akin to the idea of a gravitational softening length. Mathematically, in the limit l_max → ∞, the self-gravity expression calculated at an offset point x_off near x_q is

$$\begin{eqnarray} \Phi _{\rm Self}(\mathbf {x}_{\rm off}) & = & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}\sum _{l=0}^{\infty }P_{l} (\mathbf {n}_{\rm off}\cdot \mathbf {n}_{q})\nonumber \\ & = & \frac{-G\Delta _{q}^{3}\rho (\mathbf {x}_{q})}{r_{q}}[2(1-\mathbf {n}_{\rm off}\cdot \mathbf {n}_{q}) ]^{-1/2}, \end{eqnarray} \tag{ 8 }$$

where we have used the generating function of the Legendre polynomials, $(1-2xt+t^{2})^{-1/2}=\sum _{l=0}^{\infty }t^{l}P_{l}(x)$ (Arfken & Weber 2005) with t = 1. This expression is obviously finite, so the expansion converges. Of course, we need the potential at cell centers to compute gravitational forces—momentum and energy fluxes—at cell faces. So we modify the basic scheme above by computing potentials at all cell faces, and ascribing to each cell center the average of the potentials on the faces bounding the cell. This should be a very accurate operation as the gravitational potential is generally a smooth function in space. As shown below, this scheme works well: it converges with multipole order, and provides excellent momentum conservation.

It is important to note that the self-potential error described above is a product of the discrete evaluation of Equation (2). In spherical coordinates, it is possible to compute Equation (2) analytically, assuming constant density within the zone (Müller & Steinmetz 1995). Such an approach is not subject to the self-potential error (A. Wongwathanarat 2013, private communication). Analytic evaluation of Equation (2) in general coordinate systems is more difficult, particularly in non-spherical curvilinear systems. Thus, in order to retain uniformity amongst different coordinate systems while avoiding the self-potential error, we choose to evaluate the potentials discretely at cell faces, as discussed above.

2.2. Optimal Centering of a Multipole Expansion

The issue of where a multipole expansion should be centered has received surprisingly little analytic attention, given its importance to accurate computation of the gravitational potential. A possible reason for this is that in many cases, a spherical coordinate system is adopted, obviating the ambiguity in the choice of expansion center. For other coordinate geometries, such as cylindrical and Cartesian, the optimal location of the expansion origin is not so obvious, and a careless choice can be costly to the accuracy of the gravity solve.

There exist intuitive arguments for different choices of expansion center. The center of the grid is the obvious choice in spherical coordinate meshes. The center of mass (CoM) is indicated, perhaps a little indirectly, on the basis of the importance that it plays as a diagnostic of linear momentum conservation, since motion of the CoM directly indicates a failure of momentum conservation. The CoM is also a good choice as centering the expansion there eliminates the l = 1 dipole term (e.g., Müller & Steinmetz 1995). McGlynn (1984) working in an N-body context, advocates an expansion center location a minimizing the sum ∑_n|x_n − a|^2k, with the parameter k chosen empirically to balance the relative weighting of inner and outer particles. McGlynn (1984) also points out that the truncated multipole expansion is not translationally invariant, a point that has significant consequence for the conservation of linear momentum in calculations relying on multipole gravity solvers. This feature of multipole expansions underscores the criticality of optimally centering the expansion so as to best maintain momentum conservation.

Sellwood (1987) stresses that the origin of the multipole expansion should be placed at the location of peak density, because failure to do so can result in errors in the gravitational force, and in energy non-conservation. The intuitive reason that the peak density makes sense as the expansion origin is that condensations at large radii subtend small angles at the origin, and, if massive, can show up as power in higher-l regions of the angular momentum spectrum than would be the case were they placed near the center. It is important that the angular power spectrum of the potential be concentrated to as low values of l as is practicable, because discrete multipole Poisson solvers truncate the expansion in spherical harmonics at some l_max. This cutoff should be as low as possible, for the sake of computational efficiency (the computational cost of the Poisson solve grows as $\mathcal {O}(l_{\rm max}{}^{2})$ in three dimensions), but higher than any substantial power in the spectrum.

In this section we give more rigorous arguments than have been offered to date for the choice of expansion center. We use angular spectral "compactness," as described informally above, as the criterion for making the choice. We show that the choice advocated by Sellwood (1987) is, for all intents and purposes, very close to optimal when there is a significant fraction of total mass in a condensed object.

2.2.1. Spectral Compactness Minimization

As adumbrated above, we need a way to characterize the global angular spectral distribution in the gravitational field, so as to have some way to discuss how well the spectrum is concentrated to low values of l.

The multipole expansion of the potential Φ(x), given in Equations (2) and (3), is not ideal for this purpose, since its spectral content varies in space. We may, however, average Φ(x) spatially, weighted by the density ρ(x), to obtain the binding energy,

$$\begin{eqnarray} \mathcal {E} & = & -\frac{1}{2}\int d^{3}\mathbf {x}\,\Phi (\mathbf {x})\rho (\mathbf {x})\nonumber \\ & = & \sum _{l=0}^{\infty }\mathcal {E}_{l}, \end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray} \mathcal {E}_{l} & \equiv & \frac{G}{2}\frac{4\pi }{2l+1} \sum _{m=-l}^{l}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\, \rho (\mathbf {x})\rho (\mathbf {x}^{\prime })\nonumber \\ &&\times Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n}^{\prime })^{*}g_{l}(r,r^{\prime }), \end{eqnarray} \tag{ 10 }$$

and g_l(r, r') is the function given in Equation (3).

We propose to use $f_l\equiv \mathcal {E}_{l}/\mathcal {E}$ as a global angular spectral density in what follows. In order for this to make sense, it is of course necessary to establish that $\mathcal {E}_{l}\ge 0$ for all l. We demonstrate that this is the case in Appendix A.

How can we measure the concentration to low l of the distribution f_l? A reasonable approach is to use a moment measure, such as the mean 〈l〉 ≡ ∑_llf_l, and examine its behavior as a function of expansion center location a. It is clear that as a moves very far away from the region where most of the mass resides, the mass distribution acquires very small angular scales, and the moment measure must increase without bound. The moment measure is also obviously bounded below by 0. We require that the choice of the expansion center location a should result in a value of that moment that is as small as possible.

From the point of view of practical computation, it turns out that the most convenient moment for this purpose is

$$\begin{eqnarray} \mu (\mathbf {a}) & \equiv & \langle l(l+1) \rangle (\mathbf {a})\nonumber \\ & = & \sum _{l=0}^{\infty }l(l+1)f_{l}(\mathbf {a}). \end{eqnarray} \tag{ 11 }$$

In order to find the ideal expansion origin, we seek to minimize this "spectral compactness parameter" with respect to expansion origin, a. In Appendix B we show that the location that minimizes μ(a) is approximately

$$\begin{eqnarray} \mathbf {a} & \approx & \frac{\int d^{3}\mathbf {x}\,\mathbf {x}\rho (\mathbf {x})^{2}}{\int d^{3}\mathbf {x}\, \rho (\mathbf {x})^{2}}\nonumber \\ & \equiv & \left\langle \mathbf {x}\right\rangle _{\rho ^{2}}. \end{eqnarray} \tag{ 12 }$$

It is clear that this "square-density-weighted mean location (SDML)" is more biased toward large condensations of mass than the ordinary CoM. It is instructive to consider a simple example to illustrate the behavior of $\left\langle \mathbf {x}\right\rangle _{\rho ^{2}}$. We imagine a cubic box of side L, centered at a location x_D and containing a uniform diffuse density ρ_D corresponding to a diffuse mass M_D = L³ρ_D. The box also contains a sphere of condensed mass of radius r ≪ L and uniform density ρ_C (and hence of mass M_C = 4πr³/3ρ_C) centered at a location x_C. It is straightforward to show that with this mass configuration, the square-density-weighted CoM is

$$\begin{equation} \left\langle \mathbf {x}\right\rangle _{\rho ^{2}}=\frac{M_{C}\left(\rho _{C}+2\rho _{D}\right)\mathbf {x}_{C}+M_{D}\rho _{D} \mathbf {x}_{D}}{M_{C}\left(\rho _{C}+2\rho _{D}\right)+M_{D}\rho _{D}}. \end{equation} \tag{ 13 }$$

If, for example, we assume the situation that prevails in CCSN simulations—that is, M_C ∼ M_D, ρ_C ≫ ρ_D, then this expression becomes

$$\begin{equation} \left\langle \mathbf {x}\right\rangle _{\rho ^{2}}=\mathbf {x}_{C}+\frac{\rho _{D}}{\rho _{C}}\frac{M_{D}}{M_{C}} \left(\mathbf {x}_{D}-\mathbf {x}_{C}\right)+\mathcal {O}\left(\left[\frac{\rho _{D}}{\rho _{C}}\right]^{2}\right). \end{equation} \tag{ 14 }$$

We can see that when the density contrast between ρ_D and ρ_c is of many orders of magnitude, the square-density-weighted CoM basically takes up residence at the center of the condensation. This is the reason that the peak-density prescription for the expansion center is so effective. By contrast, the usual CoM location is the mass-weighted average of x_D and x_C, which can be well-separated from x_C if M_D ∼ M_C. Any such separation can obviously lead to troublesome angular power at high values of l.

We use the FLASH hydrodynamic simulation framework (Dubey et al. 2009) to exhibit the effects of the self-potential correction and the expansion centering schemes described above. In this section, we outline the implementation of the multipole gravity solver in FLASH. A more complete technical description of the algorithm is supplied in the FLASH User's Guide.²

The discretized potential computation expressed in Equation (4) may be separated into two distinct computations: the computation of an array of moments, and the computation of the potential itself using the moments. For notational convenience, we introduce the solid harmonic functions

$$\begin{eqnarray} R_{l m}(\mathbf {x}) & = & \sqrt{{4\pi \over {2l+1}}}r^l Y_{l m}({\bf n}) \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} I_{l m}(\mathbf {x}) & = & \sqrt{{4\pi \over {2l+1}}}{Y_{l m}({\bf n})\over r^{l+1}}. \end{eqnarray} \tag{ 16 }$$

We will define multipole moments using a grid of concentric spheres of increasing radii r_μ, μ = 1, 2, .... These radii are chosen at runtime depending on the nature of the mass distribution, and are not necessarily uniformly spaced. The spacing between radii is always more than one grid cell width, so that the shells between successive spheres encompass multiple spherical layers of cells. Given this grid, we may define the "inner" and "outer" multipole moment functions

$$\begin{eqnarray} M^R_{l m}(r_\mu) & = & \sum _{\left|\mathbf {x}_{q^\prime }\right|\le r_\mu } R_{l m}(\mathbf {x}_{q^\prime })m(q^\prime) \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} M^I_{l m}(r_\mu) & = & \sum _{\left|\mathbf {x}_{q^\prime }\right|> r_\mu } I_{l m}(\mathbf {x}_{q^\prime })m(q^\prime), \end{eqnarray} \tag{ 18 }$$

where $m(q^\prime)=\Delta _{q^\prime }^3\rho (\mathbf {x}_{q^\prime })$ is the mass of the cell indexed by q'.

We further define μ⁺(r) as the index μ of the smallest of the r_μ exceeding r, and μ⁻(r) as the index μ of the largest of the r_μ not exceeding r, so that μ⁺(r) − μ⁻(r) = 1. We may then linearly interpolate the multipole moments:

$$\begin{eqnarray} \widetilde{M}^{R,I}_{l m}(r)&\equiv & \frac{r-r_{\mu ^-(r)}}{r_{\mu ^+(r)}-r_{\mu ^-(r)}}M^{R,I}_{l m}(r_{\mu ^+(r)})\nonumber \\ &&+\frac{r-r_{\mu ^+(r)}}{r_{\mu ^-(r)}-r_{\mu ^+(r)}}M^{R,I}_{l m}(r_{\mu ^-(r)}) \end{eqnarray} \tag{ 19 }$$

Using the interpolated moments, we write the discretized potential as

$$\begin{eqnarray} &&\Phi (\mathbf {x}_q) =\nonumber\\ &&\quad - G {\mathcal {R}e}\left[\sum _{l m}\widetilde{M}^R_{l m}(|\mathbf {x}_q|))I_{l m}^*(\mathbf {x}_q) + \sum _{l m}\widetilde{M}^{I*}_{l m} (|\mathbf {x}_q|))R_{l m}(\mathbf {x}_q)\right].\nonumber\\ && \end{eqnarray} \tag{ 20 }$$

The potential evaluation strategy is to first compute the multipole moments from Equations (17) and (18) using the chosen grid of concentric spheres of radii r_μ; then, at the second stage, use this array of moments to compute the potential using Equation (20).

The FLASH implementation of this strategy relies on explicitly real (sine and cosine) versions of these formulae, which are described in the FLASH User's Guide. The real solid harmonic functions that arise are computed by recurrence relations that follow from the Legendre function recurrence relations (Arfken & Weber 2005). The radial arguments of the solid harmonic functions are carefully scaled before the recursion relations are applied, to prevent over- and underflows in large, highly resolved domains.

The implementation allows for different choices of spacing functions for the sphere radii in different radial zones, so that, for example, the spacing could be linear in an inner zone and logarithmic in an outer zone. The range of possible choices is described in the FLASH User's Guide.

As discussed in Section 2.1, the potential evaluation described by Equation (20) is always carried out at cell faces. The cell-centered potential is then computed by averaging the potential of the faces bounding a cell. Again, for multipole gravity algorithms based in spherical coordinates that compute the cell-centered potentials analytically (Müller & Steinmetz 1995), rather than discretely, the self-potential error mitigated by our staggered computation approach should not be an issue.

The center of the multipole expansion is chosen by the FLASH solver to be the cell corner nearest the square-density-weighted mean position, Equation (12). This choice minimizes the spectral compactness, as described in Section 2.2.1, and also prevents any problematic potential evaluations at zero radius.

The analytic form of the gravitational potential of a stable, rotationally symmetric, hydrostatic, uniform-density spheroid is due to MacLaurin (see Chandrasekhar 1987, p. 77). Such "MacLaurin" spheroids are useful for the validation of self-gravity solvers as they provide an exact analytic solution against which to compare the approximate calculated potentials. Here we consider the accuracy of the multipole gravity solver for static MacLaurin spheroids. We compare the accuracy of the method using cell-centered potential solves to that of using face-centered solves.

The exact gravitational potential for a point within a MacLaurin spheroid of density ρ is

$$\begin{eqnarray} &&\Phi ({\bf x}) = \pi G \rho \left[ 2A_1 a_1^2 - A_1(x^2+y^2) + A_3\left(a_3^2-z^2\right)\right], \end{eqnarray} \tag{ 21 }$$

where a₁, a₂, and a₃ are the semi-major axes of the spheroid and a₁ = a₂ > a₃. Here

$$\begin{eqnarray} A_1 &=& \frac{\sqrt{1-e^2}}{e^3} \sin ^{-1}e - \frac{1-e^2}{e^2}\ , \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} A_3 &=& \frac{2}{e^2} - \frac{2\sqrt{1-e^2}}{e^3} \sin ^{-1}e\ , \end{eqnarray} \tag{ 23 }$$

where e is the ellipticity of a spheroid:

$$\begin{equation} e = \sqrt{1 - \left(\frac{a_3}{a_1} \right) ^2}. \end{equation} \tag{ 24 }$$

For a point outside the spheroid, potential is

$$\begin{eqnarray} \Phi ({\bf x})& =& \frac{2a_3}{e^2} \pi G \rho \left[a_1 e \tan ^{-1}h - \frac{1}{2} \left(\vphantom{\frac{1}{1}}(x^2+y^2)\right.\right. \nonumber\\ &&\left.\left. \times\!\left(\tan ^{-1}h - \frac{h}{1+h^2}\right)+ 2z^2 (h-\tan ^{-1}h)\right) \!\right], \qquad \end{eqnarray} \tag{ 25 }$$

where

$$\begin{equation} h = \frac{a_1 e}{\sqrt{a_3^2 + \lambda }}\ , \end{equation} \tag{ 26 }$$

and λ is the positive root of the equation

$$\begin{equation} \frac{x^2}{a_1^2 + \lambda } + \frac{y^2}{a_2^2 + \lambda } + \frac{z^2}{a_3^2 + \lambda } = 1. \end{equation} \tag{ 27 }$$

For the present tests we consider a spheroid of uniform density ρ = 1 g cm⁻³ embedded in a background of vanishing density, ρ ≈ 0. We use an eccentricity 0.9 in 2D cylindrical geometry and compare the L2-norm error of the cell-centered potential calculation with that of the face-centered potential calculation. Figure 1 shows the results. These tests span a very large range in l_max, from 0 to 384. We find that at every value of l_max the face-centered calculation yields a smaller L2-norm error, i.e., it is more accurate. And at high values of l_max, beyond about 24, the cell-centered calculation error increases with higher l_max. The character of this increase is very nearly linear, just as we would expect based on Equation (7). The face-centered calculation, on the other hand, results in an error that continues to decrease with l_max, i.e., the accuracy of the calculation converges with l_max.

Zoom In Zoom Out Reset image size

Further evidence that the self-potential error isolated and exhibited in Equation (7) is real and present in the cell-centered potential calculation is given by inspection of the normalized error in the potential. In Figure 2 we show pseudocolor plots of the normalized error in the potential for a MacLaurin spheroid with e = 0.9 for two different values of l_max, and compare cell-centered and face-centered potential calculations. The self-potential error of Equation (7) predicts that the largest errors occur near the center of the multipole expansion. In the case of the 2D cylindrical spheroid of Figure 2 this is R = 0, z = 0.5. We see that this is precisely the case. For the cell-centered calculation there is a large normalized error at the center of the spheroid that is absent in the face-centered calculation. Additionally we see that the magnitude of this error increases for larger l_max in the cell-centered case.

Zoom In Zoom Out Reset image size

Static potentials for which analytic solutions are known are useful in verifying the accuracy of the self-gravity solver but we also seek to test if our novel handling for the errors present in multipole approximations have a positive impact on dynamical simulations that hinge critically on self-gravity. For this we turn to CCSN simulations. Having established in Section 4 that face-centered potential calculations avoid the self-potential error, resulting in greater accuracy of the potential and convergence with increasing l_max, we focus only on the face-centered potential calculation approach for the CCSN simulations. We test the impact of different multipole expansion centering on the CCSN problem by running simulations with different values of l_max for three different expansion centers: the CoM, fixed at the coordinate origin, and the SDML.

In our finite-volume Eulerian approach, gravity is coupled to the hydrodynamic calculation via source terms on the right-hand sides of the momenta and energy equations. In FLASH, these source terms are included in the Riemann solver as corrections to the intermediate cell face states that are used in calculating time-centered face fluxes of conserved quantities. We have modified the coupling of gravity and hydro in FLASH in the following way. Previous versions of FLASH extrapolated the gravitational acceleration to the time step midpoint (n + 1/2) using the current (n) and previous (n − 1) time step accelerations. This approach is formally only first-order accurate in time. We have adopted instead the second-order accurate approach of interpolating the acceleration to the time step midpoint by first updating the density field via the continuity equation, then reevaluating the gravitational potential, then finishing the finite-volume update of momenta and energy with time-centered gravitational accelerations interpolated to n + 1/2 using the n and n + 1 state accelerations. This is the approach used in, e.g., CASTRO (Almgren et al. 2010). Since this approach still utilizes source terms, the scheme is not expected to conserve momenta and energy perfectly. Such conservation can be achieved by using the method of, e.g., Jiang et al. (2013).

For these simulations we use the approach of Couch (2013a, 2013b). We follow the evolution from the collapse phase through core bounce and into shock revival by neutrino heating. We assume simple local neutrino heating/cooling as introduced by Murphy & Burrows (2008) with an exponential cutoff of the neutrino source terms at high density. Deleptonization is accounted for using the density-dependent parameterization of Liebendörfer (2005), both pre- and post-bounce. The only modification we make to the method of Couch (2013a, 2013b) is to weight the density-dependent neutrino source term cutoff so that we achieve a critical luminosity for explosion closer to that of Murphy & Burrows (2008), as was also done in Hanke et al. (2012). All of our simulations are carried out in 2D cylindrical geometry with a maximum resolution of 0.5 km and we use the 15 M_☉ progenitor of Woosley & Weaver (1995). We use a fixed neutrino luminosity of 2.2 × 10⁵² erg s⁻¹.

In Figure 3 we graphically present the results of the CCSN simulations for several values of l_max and multiple expansion origins. We show as function of post-bounce time the z-coordinate of both the $\langle \mathbf {x} \rangle _{\rho ^2}$ and the CoM along with the total z-momentum and average shock radius. For these 2D axisymmetric calculations initialized from spherically symmetric initial conditions the CoM should remain fixed at the coordinate origin, which is simply a restatement of the conservation of total z-momentum. We find that for l_max>0 centering the multipole expansion on the SDML results in dramatically improved conservation of z-momentum. For other choices of expansion center the z-momentum non-conservation can be in excess of 400 M_☉ km s⁻¹. This spurious momentum is about the same as what is observed for typical neutron stars! The magnitude of the momentum non-conservation is indiscernible for the case of centering on $\langle \mathbf {x} \rangle _{\rho ^2}$, though conservation is not perfect as reflected by the slight drift in the CoM.

Zoom In Zoom Out Reset image size

These simulations result in non-symmetric explosions and so we expect that the proto-neutron star (PNS) will receive a kick. The $\langle \mathbf {x} \rangle _{\rho ^2}$ tracks very well the center of the PNS and so its motion can be regarded as that of the PNS. Much larger kicks are imparted to the PNS for the CoM and x = 0 cases, for which we measure large non-conservations of momenta. The PNS also begins its motion much earlier than the SDML case. The kick of the PNS is obviously affected by the momentum non-conservation. It is worth noting that for CoM centering, the conservation of momentum improves with increasing l_max, but even for l_max=16 the CoM still moves by about 3 km, or six numerical zones while the CoM barely moves by one zone for any l_max for SDML centering.

The SDML centering is obviously superior to other centering choices for l_max>0, but equally obvious is its utter failure for l_max= 0. In the case of monopole gravity centering the expansion on SDML allows the PNS to move too easily away from the CoM while not correctly accounting for the strong dipole term that would result and pull the PNS back. We also see that for centering at the coordinate origin and l_max= 0, the PNS is held fixed in place and does not receive a kick. The CoM still moves in this case, reflecting non-conservation of momentum. Centering on the coordinate origin also displays divergent behavior with increasing l_max: higher values result in greater non-conservation of momentum and greater spurious motion of the PNS.

The average shock radius histories for coordinate origin centering are also highly variable with respect to changes in l_max. The other expansion centering approaches yield highly consistent shock radius histories for all values of l_max, save for l_max=0 in the SDML case.

Our choice of the SDML for the multipole expansion centering is motivated by our minimization of the spectral compactness, μ, introduced in Section 2. This metric, defined as 〈l(l + 1)〉, measures the concentration of total gravitational potential energy at low multipole orders. Our analysis in Section 2 indicates that centering the expansion on the SDML should maximize the amount of total potential energy from low orders, i.e., yield the most spherical representation of the gravitational potential. To test this for the CCSN simulations we compute the normalized potential energy spectra, $f_l \equiv \mathcal {E}_l / \mathcal {E}$, for the three different expansion centering approaches at three different times, shown in Figure 4. All of the simulations in Figure 4 were run with l_max= 16 but we compute the spectra out to l = 64 and indicate l = 16 by the vertical dashed line. Prior to core bounce (t_pb = −2 ms) the spectra are highly concentrated at l = 0, as expected for the spherically-symmetry mass distribution, and the odd multipoles are much reduced due to the symmetry. By 200 ms post-bounce the spectra remain very similar except for the slightly reduced power in l = 1 and greater power in l = 0 for the SDML case. The shock radii are also very similar at this time (see Figure 3). The differences in the centering approaches are more obvious at 500 ms after runaway shock expansion has begun. The reduced power at multipoles greater than 0 is evident for the SDML case while centering on the coordinate origin results in a spreading of the spectrum to larger multipoles.

Zoom In Zoom Out Reset image size

We have identified and corrected two sources of error arising in general discretized multipole approximations to Poisson's equation. The first error results from assuming that all the mass in a computational zone resides at the cell center and then evaluating the potential at the same point. Inspection of the Green's function for the continuous Poisson equation makes obvious that this error has its origin in the divergent |x − x'|⁻¹ term. This term is explicitly absent from the discretized equations but the error it induces is still lurking in the method. We show that this error is proportional to the mass in a zone, divided by the distance of the zone center from the origin of the multipole expansion, multiplied by l_max+1. This error therefore grows rather than shrinking as the number of terms retained in the truncated expansion increases. We show that the "self-potential" error can be corrected by evaluating the gravitational potential at cell faces, where no mass has been located, rather than cell centers. The cell-centered potential is then found by averaging the potential at the cell-bounding faces. Using MacLaurin spheroids, for which exact analytic potential solutions are known, we show that this approach improves the accuracy of the potential calculation and leads to convergence of the solution with increasing l_max, i.e., the self-potential error is eliminated.

The second error we identify has to do with a poor selection of the multipole expansion origin. By suggesting a useful metric, the spectral compactness μ = 〈l(l + 1)〉, characterizing the symmetry of the potential we find that the optimal location for the origin that minimizes μ is the SDML, $\langle \mathbf {x} \rangle _{\rho ^2}$. For diffuse mass distributions, or distributions in which the total mass in the computational domain is dominated by a single condensation, this location is not too different from the CoM, the common choice for multipole expansion origin. For high-mass condensations embedded in high-mass diffuse flows, such as occur in CCSN simulations that include the proto-neutron star, the $\langle \mathbf {x} \rangle _{\rho ^2}$ is close to the peak density of the high-mass condensate. Using a series of CCSN simulations we demonstrate the superiority of locating the expansion center at the SDML: momentum conservation is dramatically improved resulting in significantly different kicks imparted to the PNS by the development of asymmetric explosions. CCSN simulations that include the PNS are especially susceptible the two errors we discuss because of the enormous mass density in few zones near the expansion origin.

Our computational approach is embedded in an Eulerian hydrodynamic framework. Nevertheless, the multipole approach to the solution of the Poisson equation is quite general, and its numerical implementation stands apart from the specific numerical hydrodynamic scheme employed here. It follows that the improvements we describe above to the discretized multipole approximation to Poisson's equation are generally applicable. In particular, the optimal choice of expansion center is relevant to all simulations that employ multipole approach for calculating self-gravity, and the face-centering of the potential calculation is relevant to all such approaches that are grid-based and do not evaluate potentials analytically, as can be done in spherical geometry (Müller & Steinmetz 1995).

It is important to note that the momentum non-conservation, and concomitant erroneous motion of the PNS, is due to the movement of the PNS away from the origin of the multipole expansion. In spherical geometry where the PNS is unable to move away from the origin, or in CCSN simulations that excise the PNS, we do not expect to see such bad momentum non-conservation. We, therefore, do not expect previous studies of PNS kicks that excise the PNS from the domain in spherical geometry (Scheck et al. 2004, 2006; Wongwathanarat et al. 2010, 2013) to suffer from the inaccuracies we here uncover and correct. Likewise our results have no bearing on PNS kick studies that do not utilize multipole gravity solvers (Nordhaus et al. 2010, 2012).

The multipole approach is appropriate for systems wherein the mass distribution is approximately spherical, so that a spherical harmonic expansion can be expected to reach high accuracy after a moderate number of terms. For such problems it has substantial benefits over other approaches for solving Poisson's equation, such as multigrid or tree methods, because it is comparatively inexpensive. For the time-dependent CCSN simulations described in Section 5 the multipole implementation we present in Section 3 requires less than 7% of the time to calculate the hydrodynamics. More exact multigrid and tree methods can dominate the computational expense of simulations utilizing them (cf. Ricker 2008). By incorporating the two essential reforms of the multipole algorithm we present the method can deliver on its promise of accurate calculation of self gravity while also retaining its efficient computation.

The authors thank Paul Ricker, Thomas Janka, and Annop Wongwathanarat for helpful conversations. We especially thank our referee, Ewald Müller, for a careful review that improved this article. S.M.C. is supported by NASA through Hubble Fellowship grant No. 51286.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. This work was supported in part by the National Science Foundation under grant AST-0909132. This work was supported in part at the University of Chicago by the U.S. Department of Energy (DOE) under contract B523820 to the NNSA ASC/Alliances Center for Astrophysical Thermonuclear Flashes. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research used computational resources at ALCF at ANL, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.

Our proposed angular decomposition of the global spectral of the gravitational field is proportional to $\mathcal {E}_{l}$, where

$$\begin{eqnarray} \mathcal {E}_{l} & \equiv & \!\frac{G}{2}\frac{4\pi }{2l+1}\!\sum _{m=-l}^{l}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\, \rho (\mathbf {x})\rho (\mathbf {x}^{\prime })\, Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n}^{\prime })^{*}g_{l}(r,r^{\prime }), \nonumber\\ \end{eqnarray} \tag{ A1 }$$

and g_l(r, r') is the function given in Equation (3).

In order for this choice of spectral decomposition of the field to give rise to a sensible distribution, it is necessary to establish that $\mathcal {E}_{l}\ge 0$ for all l. To do this, we write Equation (A1) as

$$\begin{equation} \mathcal {E}_{l}=\frac{G}{2}\frac{4\pi }{2l+1}\sum _{m=-l}^{l}\int _{0}^{\infty }\!r^{2}dr\int _{0}^{\infty }\!r^{\prime 2}dr^{\prime } g_{l}(r,r^{\prime })\,\eta _{lm}(r)\eta _{lm}(r^{\prime })^{*}, \end{equation} \tag{ A2 }$$

where we have defined the moments

$$\begin{equation} \eta _{lm}(r)\equiv \oint d^{2}\mathbf {n}\,\rho (r\mathbf {n})Y_{lm}(\mathbf {n}). \end{equation} \tag{ A3 }$$

From Equation (A2), we see that the requirement that $\mathcal {E}_{l}\ge 0$ is equivalent to requiring positive-semi-definiteness of the integral operator $\hat{g}_{l}$, whose action on a function η(r) is $[\hat{g}_{l}\circ \eta ](r)=\int _{0}^{\infty }r^{\prime 2}dr^{\prime }\, g_{l}(r,r^{\prime })\eta (r^{\prime })$. In other words, we must have $(\eta ,\hat{g}_{l}\circ \eta)=\int _{0}^{\infty }r^{2}dr\,\eta (r)[\hat{g}_{l}\circ \eta ](r)\ge 0$. But $(2l+1)^{-1}\hat{g}_{l}$ is the inverse of the radially-separated Laplacian differential operator $\mathcal {L}_{l}\equiv -{1}/{r^{2}}{\partial }/{\partial r}\left(r^{2}{\partial }/{\partial r}\right)+{l(l+1)}/{r^{2}}$, for which g_l(r, r') is the Green's function: $[\mathcal {L}_{l}\circ \hat{g}_{l}] (r,r^{\prime })={2l+1}/{r^2}\delta (r-r^{\prime })$. Furthermore, we may easily show that $\mathcal {L}_{l}$ is positive-definite, $(\eta ,\mathcal {L}_{l}\circ \eta )>0$, for the boundary conditions of interest here (finite at the origin, zero at infinity) by means of an integration by parts. Setting $\hat{g}_{l}\circ \eta \equiv (2l+1)\chi$, so that $\eta =\mathcal {L}_{l}\circ \chi$ we therefore have

$$\begin{equation} (\eta ,\hat{g}_{l}\circ \eta)=(2l+1)\times \left(\mathcal {L}_{l}\circ \chi ,\chi \right)>0. \end{equation} \tag{ A4 }$$

Since $\hat{g}_{l}$ is a positive-definite integral operator, it follows immediately from Equation (A2) that $\mathcal {E}_{l}>0$ for all l.

The normalized distribution over l $f_{l}\equiv \mathcal {E}_{l}/\mathcal {E}$ is therefore a sensible measure of the angular spectrum in a gravitating mass distribution. The total binding energy $\mathcal {E}$ is obviously independent of the expansion center position a. The individual terms $\mathcal {E}_{l}$ in the decomposition are certainly functions of a, however, so that the spectral distribution is also dependent on a. We will therefore write this dependence as $\mathcal {E}_{l}(\mathbf {a})$ explicitly below.

When calculating the spectrum f_l empirically from a mass distribution, as we do in Section 5, there is a subtle source of error to be guarded against, which is traceable to discretization noise. The effect comes about because, as remarked earlier, the mass of each cell, which represents a volume integral of some smooth, nearly constant mass density function over the cell, is represented in the numeric quadratures of the multipole algorithm as a Dirac-δ-function at the cell center. Obviously, an infinitely-narrow density peak is capable of contributing power to arbitrarily-high multipole orders l, whereas the cell's contribution to the angular spectrum due to the underlying, nearly constant density function should cut off rapidly above some angular scale. The error therefore manifests itself in the spectrum f_l as a noisy positive direct current (DC)-offset level at high l-values. In order to exhibit normalizable spectra, it is necessary to remove this error. This can be done by observing that a cell with index q, of size Δ_q, located at a distance r_q from the center of the expansion, subtends an angle θ_q ∼ Δ_q/r_q at the center. We should not expect such a cell to contribute anything but noise to multipoles of order l > 2π/θ_q. Discarding such terms from the multipole moment contribution of these cells, the DC offset is removed, and normalizable spectra such as the ones shown in Section 5 are recovered.

As asserted in Section 2.2.1, from the point of view of practical computation, it turns out that the most convenient moment for the purpose of quantifying angular spectrum compactness is

$$\begin{eqnarray} \mu (\mathbf {a}) & \equiv & \langle l(l+1) \rangle (\mathbf {a}) \nonumber \\ & = & \sum _{l=0}^{\infty }l(l+1)f_{l}(\mathbf {a}). \end{eqnarray} \tag{ B1 }$$

The reason this is convenient is because the term l(l + 1) arises naturally from the application of the Laplacian to the spherical harmonic expansion of the Green's function |x − x'|⁻¹:

$$\begin{eqnarray*} -4\pi \delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime }) &= & \nabla ^{2}\left|\mathbf {x}-\mathbf {x}^{\prime }\right| ^{-1}\\ & = & \nabla ^{2}\left\lbrace \sum _{l=0}^{\infty }\sum _{m=-l}^{l}\frac{4\pi }{2l+1}Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n} ^{\prime })^{*}g_{l}(r,r^{\prime })\right\rbrace \\ & = & \sum _{l=0}^{\infty }\sum _{m=-l}^{l}\frac{4\pi }{2l+1}\left[\frac{1}{r^{2}}\frac{\partial }{\partial r} \left(r^{2}\frac{\partial }{\partial r}\right)-\frac{l(l+1)}{r^{2}}\right]\nonumber\\ && \times\,Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n}^{\prime })^{*} g_{l}(r,r^{\prime })\\ & = & \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}\frac{\partial }{\partial r}\right)|\mathbf {x}-\mathbf {x} ^{\prime }|^{-1}\nonumber\\ &&-\frac{1}{r^{2}}\sum _{l=0}^{\infty }\sum _{m=-l}^{l}\frac{4\pi } {2l+1}\nonumber\\ &&\times l(l+1)Y_{lm}(\mathbf {n})Y_{lm} (\mathbf {n}^{\prime })^{*}g_{l}(r,r^{\prime }) \end{eqnarray*}$$

so that

$$\begin{eqnarray} &&\sum _{l=0}^{\infty }\sum _{m=-l}^{l}\frac{4\pi }{2l+1}l(l+1)Y_{lm}(\mathbf {n})Y_{lm}(\mathbf {n}^{\prime })^{*}g_{l} (r,r^{\prime }) \nonumber\\ &&\quad = \frac{\partial }{\partial r}\left(r^{2}\frac{\partial }{\partial r}\right)|\mathbf {x}-\mathbf {x} ^{\prime }|^{-1}+4\pi r^{2}\delta ^{3}(\mathbf {x}-\mathbf {x}^{\prime }).\quad \end{eqnarray} \tag{ B2 }$$

Combining Equations (A1) and (B2) with Equation (B1), and setting a = 0 temporarily, we obtain

$$\begin{eqnarray} \mu (\mathbf {0})&=&\frac{G}{2\mathcal {E}}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x} ^{\prime })\nonumber\\ &&\times\,\left\lbrace 4\pi r^{2}\delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })+\frac{\partial }{\partial r} \left(r^{2}\frac{\partial }{\partial r} \right)|\mathbf {x}-\mathbf {x}^{\prime }|^{-1}\!\right\rbrace .\qquad \end{eqnarray} \tag{ B3 }$$

At the cost of some algebra, we may evaluate the second term in Equation (B3). We obtain

$$\begin{eqnarray} \mu (\mathbf {0})&=&\frac{G}{2\mathcal {E}}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x} ^{\prime })\left\lbrace\vphantom{\frac{-(\mathbf {x}\cdot \mathbf {x} ^{\prime })^{2}-3r^{2}r^{\prime ^{2}}+2(r^{2}+r^{\prime 2})\mathbf {x}\cdot \mathbf {x}^{\prime }}{\left|\mathbf {x}-\mathbf {x} ^{\prime }\right|^{5}}} 4\pi r^{2}\delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })\right.\nonumber\\ &&\left. +\frac{-(\mathbf {x}\cdot \mathbf {x} ^{\prime })^{2}-3r^{2}r^{\prime ^{2}}+2(r^{2}+r^{\prime 2})\mathbf {x}\cdot \mathbf {x}^{\prime }}{\left|\mathbf {x}-\mathbf {x} ^{\prime }\right|^{5}}\right\rbrace . \end{eqnarray} \tag{ B4 }$$

After making the substitution x' = (x' − x) + x in the numerator of the second term, some further algebra yields

$$\begin{eqnarray} \mu (\mathbf {0}) & = & \frac{G}{2\mathcal {E}}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x}) \rho (\mathbf {x}^{\prime })\left \lbrace\vphantom{\sum _{i,k=1}^{3}x_{i} x_{k}} 4\pi r^{2}\delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })\right.\nonumber\\ &&+\left.\sum _{i,k=1}^{3}x_{i} x_{k}\frac{3(x_{i}^{\prime }-x_{i})(x_{k}^{\prime }-x_{k})-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{\left|\mathbf {x}-\mathbf {x}^{\prime }\right|^{5}}\right.\nonumber\\ &&+ \left.2\sum _{i=1}^{3}x_{i}\frac{(x_{i}^{\prime }-x_{i})}{\left|\mathbf {x}-\mathbf {x}^{\prime }\right|^{3}} \right \rbrace. \end{eqnarray} \tag{ B5 }$$

To obtain μ(a) from μ(0) all that is required is to make the replacements x → x − a, x' → x' − a inside the braces in Equation (B5). The result is

$$\begin{eqnarray} \mu (\mathbf {a}) & = & \frac{G}{2\mathcal {E}}\int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x}) \rho (\mathbf {x}^{\prime }) \left\lbrace\vphantom{\frac{3(x_{i}-x_{i}^{\prime }) (x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}-\mathbf {x}^{\prime } |^{5}}} 4\pi \delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })|\mathbf {x}-\mathbf {a}| ^{2}\right.\nonumber\\ && \left. +\!\sum _{i,k=1}^{3}(x_{i}\,{-}\,a_{i})(x_{k}\,{-}\,a_{k})\!\frac{3(x_{i}^{\prime }\,{-}\,x_{i})(x_{k}^{\prime }\,{-}\,x_{k}){-}|\mathbf {x}\,{-}\, \mathbf {x}^{\prime }|^{2}\delta _{ik}}{\left|\mathbf {x}\,{-}\,\mathbf {x}^{\prime }\right|^{5}}\right.\nonumber \\ && \left. +2\sum _{i=1}^{3}(x_{i}-a_{i})\frac{(x_{i}^{\prime }-x_{i})}{\left|\mathbf {x}-\mathbf {x}^{\prime } \right|^{3}}\right \rbrace. \end{eqnarray} \tag{ B6 }$$

The extremization with respect to a of this quadratic expression in a is straightforward, and leads to a 3 × 3 linear problem,

$$\begin{equation} \mathbf {Ma}=\mathbf {b}, \end{equation} \tag{ B7 }$$

with

$$\begin{eqnarray} \left[\mathbf {M}\right]_{ik} & \equiv & \int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x} ^{\prime }) \left\lbrace\vphantom{\frac{3(x_{i}-x_{i}^{\prime }) (x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}-\mathbf {x}^{\prime } |^{5}}} 4\pi \delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })\delta _{ik}\right.\nonumber\\ &&+ \left.\frac{3(x_{i}-x_{i}^{\prime }) (x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}-\mathbf {x}^{\prime } |^{5}}\right\rbrace \nonumber \\ & = & \int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x}^{\prime })\nonumber\\ &&\times \left\lbrace 4\pi \delta ^{3} (\mathbf {x}-\mathbf {x}^{\prime })\delta _{ik}+\frac{\partial ^{2}}{\partial x_{i}\partial x_{k}}\frac{1}{\left|\mathbf {x}-\mathbf {x}^{\prime }\right|}\right\rbrace \nonumber \\ & = & \int d^{3}\mathbf {x}\,\left\lbrace 4\pi \rho (\mathbf {x})^{2}\delta _{ik}-\rho (\mathbf {x})\, G^{-1}\frac{\partial ^{2}\Phi (\mathbf {x})}{\partial x_{i}\partial x_{k}}\!\right\rbrace , \end{eqnarray} \tag{ B8 }$$

and

$$\begin{eqnarray} \left[\mathbf {b}\right]_{i} & = & \int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x} ^{\prime }) \left\lbrace\vphantom{\frac{3(x_{i}-x_{i}^{\prime }) (x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}-\mathbf {x}^{\prime } |^{5}}} 4\pi \delta ^{3}(\mathbf {x}-\mathbf {x}^{\prime })x_{i}\right.\nonumber\\ &&+\sum _{k=1}^{3}x_{k}\frac{3(x_{i}- x_{i}^{\prime })(x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}- \mathbf {x}^{\prime }|^{5}}\nonumber\\ &&-\left.\frac{x_{i}-x_{i}^{\prime }}{|\mathbf {x}-\mathbf {x}^{\prime }|^{3}}\right\rbrace \nonumber \\ & = & \int d^{3}\mathbf {x}\, d^{3}\mathbf {x}^{\prime }\,\rho (\mathbf {x})\rho (\mathbf {x}^{\prime })\times \left\lbrace\vphantom{\frac{3(x_{i}-x_{i}^{\prime }) (x_{k}-x_{k}^{\prime })-|\mathbf {x}-\mathbf {x}^{\prime }|^{2}\delta _{ik}}{|\mathbf {x}-\mathbf {x}^{\prime } |^{5}}} 4\pi \delta ^{3}(\mathbf {x}-\mathbf {x}^{\prime })x_{i}\right.\nonumber\\ &&+\left.\sum _{k=1}^{3}x_{k}\frac{\partial ^{2}}{\partial x_{i}\partial x_{k}}\frac{1}{|\mathbf {x}-\mathbf {x}^{\prime }|}+\frac{\partial }{\partial x_{i}}\frac{1}{|\mathbf {x}- \mathbf {x}^{\prime }|}\right\rbrace \nonumber \\ & = & \int d^{3}\mathbf {x}\,\left\lbrace 4\pi \rho (\mathbf {x})^{2}x_{i}-\rho (\mathbf {x})\sum _{k=1}^{3}x_{k} G^{-1}\frac{\partial ^{2}\Phi (\mathbf {x})}{\partial x_{i}\partial x_{k}}\right.\nonumber\\ &&-\left.\rho (\mathbf {x})G^{-1}\frac{\partial \Phi (\mathbf {x})}{\partial x_{i}}\right\rbrace . \end{eqnarray} \tag{ B9 }$$

The last term in Equation (B9) yields, upon integration, G⁻¹ times the net self-force of the gravitating mass configuration, which is necessarily zero. We therefore have for b

$$\begin{equation} \left[\mathbf {b}\right]_{i}=\int d^{3}\mathbf {x}\,\left\lbrace 4\pi \rho (\mathbf {x})^{2}x_{i}-\rho (\mathbf {x})\sum _{k=1}^{3} x_{k}G^{-1}\frac{\partial ^{2}\Phi (\mathbf {x})}{\partial x_{i}\partial x_{k}}\right\rbrace . \end{equation} \tag{ B10 }$$

The integrands in Equations (B8) and (B10) feature the sum of a term proportional to the square of the density, and a term proportional to the tidal tensor ∂²Φ/∂x_i∂x_k. In the next section, we estimate the relative sizes of the two terms in each of the two integrals, and find that it is an acceptable approximation to drop the tidal terms in comparison with the square-density terms. Making this approximation, we obtain

$$\begin{eqnarray} \mathbf {M} & \approx & \mathbf {1}\times \int d^{3}\mathbf {x}\,4\pi \rho (\mathbf {x})^{2} \end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray} \mathbf {b} & \approx & \int d^{3}\mathbf {x}\,4\pi \rho (\mathbf {x})^{2}\,\mathbf {x}, \end{eqnarray} \tag{ B12 }$$

so that

$$\begin{eqnarray} \mathbf {a} & \approx & \frac{\int d^{3}\mathbf {x}\,\mathbf {x}\rho (\mathbf {x})^{2}}{\int d^{3}\mathbf {x}\, \rho (\mathbf {x})^{2}}\nonumber \\ & \equiv & \left\langle \mathbf {x}\right\rangle _{\rho ^{2}}. \end{eqnarray} \tag{ B13 }$$

That is to say, the optimal expansion center location is the average location weighted by the square of the density.

A useful result worth setting down is a formula for the spectral compactness, μ(0) that is convenient for numerical computation. Starting from Equation (B5), we may replace the dipole and quadrupole tensors with suitable derivatives of the Green's function, as we did in Equations (B8) and (B9). We find that

$$\begin{eqnarray} \mu (\mathbf {0})&=&\frac{G}{2\mathcal {E}}\int d^{3}\mathbf {x}\,\left\lbrace 4\pi \left|\mathbf {x}\right|^{2}\rho (\mathbf {x})^{2}+ \rho (\mathbf {x})G^{-1}\sum _{i,k=1}^{3}x_{i}x_{k}\frac{\partial ^{2}\Phi (\mathbf {x})}{\partial x_{i}\partial x_{k}}\right.\nonumber\\ &&{-} \left.2\rho (\mathbf {x})G^{-1}\sum _{i=1}^{3}x_{i}\frac{\partial \Phi (\mathbf {x})}{\partial x_{i}}\right\rbrace . \end{eqnarray} \tag{ B14 }$$

All the data required to compute this integral over the domain is available after the potential has been computed.

In this Appendix we estimate the relative sizes of the two terms in Equations (B8) and (B10).

By re-expressing the potential Φ(x) in terms of the density ρ(x), and making the change of variables x' → y = x' − x, Equations (B8) and (B10) may be written as

$$\begin{eqnarray} \left[\mathbf {M}\right]_{ik} & = & \int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\int d^{3}\mathbf {y}\rho (\mathbf {x}+\mathbf {y}) \nonumber\\ &&\times\left\lbrace 4\pi \delta ^{3}\left(\mathbf {y}\right)\delta _{ik}-\frac{y^{2}\delta _{ik}-3y_{i}y_{k}}{y^{5}}\right\rbrace \end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray} \left[\mathbf {b}\right]_{i} & = & \int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\int d^{3}\mathbf {y}\rho (\mathbf {x}+\mathbf {y}) \nonumber\\ &&\times\left\lbrace 4\pi \delta ^{3}\left(\mathbf {y}\right)x_{i}-\frac{\sum _{k=1}^{3}x_{k}\left[y^{2}\delta _{ik}-3y_{i}y_{k}\right]}{y^{5}}\right\rbrace. \qquad \end{eqnarray} \tag{ C2 }$$

We single out the rational terms in these integrals:

$$\begin{eqnarray} &&[\mathbf {M}^{R}]_{ik} \equiv -\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\int d^{3}\mathbf {y}\rho (\mathbf {x}+ \mathbf {y})\frac{y^{2}\delta _{ik}-3y_{i}y_{k}}{y^{5}}\qquad \end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray} \left[\mathbf {b}^{R}\right]_{i} & \equiv & -\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\int d^{3}\mathbf {y}\rho (\mathbf {x}+ \mathbf {y})\nonumber\\ &&\times\frac{\sum _{k=1}^{3}x_{k} [y^{2}\delta _{ik}-3y_{i}y_{k}]}{y^{5}}. \end{eqnarray} \tag{ C4 }$$

The numerators of the integrands contain the trace-free symmetric tensor y²δ_ik − 3y_iy_k. We recognize this as quadrupole tensor, and exploit its nature as a spherical tensor—a spherical harmonic in tensor guise—to reduce the order of the singularity in y.

We will require the following spherical integrals:

$$\begin{eqnarray} \oint d^{2}\mathbf {n}\, n_{i}n_{k} & = & \frac{4\pi }{3}\delta _{ik} \end{eqnarray} \tag{ C5 }$$

$$\begin{eqnarray} \oint d^{2}\mathbf {n}\, n_{i}n_{k}n_{l}n_{m} & = & \frac{4\pi }{15}\left(\delta _{ik}\delta _{ml}+\delta _{im}\delta _{kl}+ \delta _{il}\delta _{mk}\right). \end{eqnarray} \tag{ C6 }$$

These may be obtained by observing that the resulting tensors must be rotationally-invariant and totally symmetric under index interchange. Such tensors can only be constructed from the only tensor at hand—the identity tensor δ_ik—by the combinations indicated. The coefficients may then be calculated by setting i = k = l = m = 3 in the resulting expressions and performing the integrals in spherical coordinates. In addition, we observe that any similar integral featuring an odd number of components of n as factors in the integrand is necessarily zero, since it changes sign under the variable change n → −n. Note also that Equation (C5) implies that the spherical integral of the quadrupole tensor δ_ik − 3n_in_k is zero.

Since we are interested in the y → 0 behavior, we expand ρ(x + y) around x:

$$\begin{eqnarray} \rho (\mathbf {x}+\mathbf {y})&=&\rho (\mathbf {x})+\sum _{l=1}^{3}y_{l}\frac{\partial \rho (\mathbf {x})}{\partial x_{l}}\nonumber\\ &&+\frac{1}{2}\sum _{l=1}^{3}\sum _{m=1}^{3} [y_{l}y_{m}+\mathcal {O}(y^{3})]\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{l}\partial x_{m}}. \end{eqnarray} \tag{ C7 }$$

Inserting this expansion in Equation (C3), we obtain

$$\begin{eqnarray} [\mathbf {M}^{R}]_{ik} & \approx & \!-\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\frac{1}{2}\sum _{l=1}^{3}\sum _{m=1} ^{3}\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{l}\partial x_{m}}\int _{0}^{\infty }y^{2}\,dy\nonumber\\ &&\times\oint d^{2}\mathbf {n}\, [y_{l}y_{m}+\mathcal {O}(y^{3})]\frac{y^{2}\delta _{ik}-3y_{i}y_{k}}{y^{5}}\nonumber \\ & = & \!\!-\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\frac{1}{2}\sum _{l=1}^{3}\sum _{m=1}^{3}\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{l}\partial x_{m}}\int _{0}^{\infty }\![y+\mathcal {O}(y^{2})]\, dy\nonumber\\ &&\times\oint d^{2}\mathbf {n} \left(\delta _{ik}n_{l}n_{m}-3n_{i}n_{k}n_{l}n_{m}\right)\nonumber \\ & = & \!-\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\frac{1}{2}\sum _{l=1}^{3}\sum _{m=1}^{3}\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{l}\partial x_{m}}\int _{0}^{\infty }\![y+\mathcal {O}(y^{2})]\, dy\nonumber\\ &&\times\frac{4\pi }{15}\left(2\delta _{ik} \delta _{lm}-3\delta _{il}\delta _{km}-3\delta _{im}\delta _{kl}\right)\nonumber \\ & = & \!-\frac{4\pi }{15}\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\left(\delta _{ik} \nabla ^{2}\rho (\mathbf {x})-3\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{i}\partial x_{k}}\right)\nonumber\\ &&\times \int _{0}^{\infty } [y+\mathcal {O}(y^{2})] \,dy, \end{eqnarray} \tag{ C8 }$$

where in the first line we summarily dropped from the density expansion both the $\mathcal {O}(y^{0})$ term—because it results in a spherical integral of the quadrupole tensor, which is zero—and the $\mathcal {O}(y^{1})$ term—because it results in a spherical integral with an odd number of vector factors, which is also zero. In the inner integrand, we see that the dependence on y as y → 0 is a very benign $\mathcal {O}(y^{1})$.

We proceed similarly, inserting the expansion of Equation (C7) into Equation (C4). Again, only one term from the expansion survives, with the $\mathcal {O}(y^{2})$ term latching on to the quadrupole, as before. We obtain

$$\begin{eqnarray} [\mathbf {b}^{R}]_{i}&\approx &\! {-}\frac{4\pi }{15}\int d^{3}\mathbf {x}\,\rho (\mathbf {x})\,\left\lbrace \left(x_{i} \nabla ^{2}\rho (\mathbf {x})-3\sum _{k=1}^{3}x_{k}\frac{\partial ^{2}\rho (\mathbf {x})}{\partial x_{l}\partial x_{k}}\right) \right\rbrace\nonumber\\ &&\times\ \int _{0}^{\infty } [y+\mathcal {O}(y^{2})]\, dy. \end{eqnarray} \tag{ C9 }$$

We again find that the "singular" behavior of the integrand is in fact $\mathcal {O}(y^{1})$ as y → 0. This $\mathcal {O}(y^{1})$ behavior is no different from the short-distance behavior of the Poisson Green's function, which combines a y⁻¹ singularity with the d³y measure to produce an $\mathcal {O}(y^{1})$ dependence in the integrand.

We now use the expressions just derived to estimate the relative size of the rational terms and the δ-function terms in Equations (C1) and (C2). To do this, we assume a distribution of matter bounded to some region of size R. We estimate the term $\int dy[y+\mathcal {O}(y^{2})]\sim R^{2}/2$. We also assume the presence of a sharp peak in the density, so that the integral measure d³x ρ(x) places most of the action near the density peak. In this region, the linear term in the expansion in Equation (C7) is small compared to the $\mathcal {O}(y^{2})$ term, and may be neglected, and we may estimate 1/2R²|∂²ρ/∂x²| ∼ ρ. We use this estimate for the aggregate second-derivative terms in brackets in Equations (C8) and (C9). We also replace the factors x in Equation (C9) by a typical value R. By these means, we obtain for the size of the matrix elements of M^R

$$\begin{eqnarray} |[\mathbf {M}^{R}]_{ik}| & \sim & \frac{4\pi }{15}\int d^{3}\mathbf {x}\,\rho (\mathbf {x})^{2}. \end{eqnarray} \tag{ C10 }$$

Since the δ-function term in Equation (C1) is

$$\begin{equation} [\mathbf {M}^{\delta }]_{ik}=4\pi \delta _{ik}\int d^{3}\mathbf {x}\,\rho (\mathbf {x})^{2}, \end{equation} \tag{ C11 }$$

we obtain the ratio

$$\begin{equation} \left|\frac{[\mathbf {M}^{R}]_{ik}}{[\mathbf {M}^{\delta }]_{ik}}\right|\sim \frac{1}{15}, \end{equation} \tag{ C12 }$$

give or take a little slop. By the same means, we obtain

$$\begin{equation} |[\mathbf {b}^{R}]_{k}|\sim \frac{4\pi R}{15}\int d^{3}\mathbf {x}\,\rho (\mathbf {x})^{2}. \end{equation} \tag{ C13 }$$

The δ-function term in Equation (C2) is

$$\begin{eqnarray} [\mathbf {b}^{\delta }]_{k} & = & 4\pi \int d^{3}\mathbf {x}\,\rho (\mathbf {x})^{2}\, x_{k}\nonumber \\ & \sim & 4\pi R\int d^{3}\mathbf {x}\,\rho (\mathbf {x})^{2}, \end{eqnarray} \tag{ C14 }$$

so that the ratio of terms is again

$$\begin{equation} \left|\frac{[\mathbf {b}^{R}]_{k}}{[\mathbf {b}^{\delta }]_{k}}\right|\sim \frac{1}{15}. \end{equation} \tag{ C15 }$$

On the basis of these estimates it appears that the neglect of the rational terms in Equations (B8) and (B10) is a justifiable approximation.