The Challenge of Sixfold Integrals: The Closed-Form Evaluation of Newton Potentials between Two Cubes

The challenge of explicitly evaluating, in elementary closed form, the weakly singular sixfold integrals for potentials and forces between two cubes has been taken up at various places in the mathematics and physics literature. It created some strikingly specific results, with an aura of arbitrariness, and a single intricate general procedure due to Hackbusch. Those scattered instances were mostly addressing the problem heads on, by successive integration while keeping track of a thicket of primitives generated at intermediate stages. In this paper we present a substantially easier and shorter approach, based on a Laplace transform of the kernel. We clearly exhibit the structure of the results as obtained by an explicit algorithm, just computing with rational polynomials. The method extends, up to the evaluation of single integrals, to higher dimensions. Among other examples, we easily reproduce Fornberg's startling closed form solution of Trefethen's two-cubes problem and Waldvogel's symmetric formula for the Newton potential of a rectangular cuboid.


page 1

page 2

page 3

page 4


A closed-form formula for the Kullback-Leibler divergence between Cauchy distributions

We report a closed-form expression for the Kullback-Leibler divergence b...

An Elementary Proof of the Generalization of the Binet Formula for k-bonacci Numbers

We present an elementary proof of the generalization of the k-bonacci Bi...

Variance Analysis for Monte Carlo Integration: A Representation-Theoretic Perspective

In this report, we revisit the work of Pilleboue et al. [2015], providin...

Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping

The method of Helsing and co-workers evaluates Laplace and related layer...

On a seventh order convergent weakly L-stable Newton Cotes formula with application on Burger's equation

In this paper we derive 7^th order convergent integration formula in tim...

High Performance Evaluation of Helmholtz Potentials usingthe Multi-Level Fast Multipole Algorithm

Evaluation of pair potentials is critical in a number of areas of physic...

High Performance Evaluation of Helmholtz Potentials using the Multi-Level Fast Multipole Algorithm

Evaluation of pair potentials is critical in a number of areas of physic...

1 Introduction

Here are three openings for this paper, looking at some remarkable spotlights of closed form evaluations of Newton potentials and forces in reverse chronological order. Interestingly, the protagonists started their work each time from scratch, completely independent of each other.

  • In Oct. 2006, as yet another instance of his “10-digit problems”, Nick Trefethen cooked up the following challenge for the graduate students in Oxford’s Numerical Analysis “Problem Solving Squad” Tref2:

    Two homogeneous unit cubes of unit mass attract each other gravitationally according to Newton’s law with unit gravitational constant. Their centers are one unit apart, so the cubes are right up against each other, touching. What is the force, to ten digit accuracy?
    *[1mm]XXX [slightly edited for a better fit]

    With some background in theoretical mechanics, the challenge amounts to evaluating the weakly singular sixfold integral


    As a problem in numerical analysis, using a delightful subdivision idea the challenge was nailed by graduate student Alex Prideaux, see (Tref1, p. 124):

    While on sabbatical at Oxford in 2010, Bengt Fornberg learnt about the problem, thought it ought to be doable in closed form by elementary means and, after a week full of fun and intense work, surprised everybody with111In a personal communication Fornberg kindly shared his recollections: “Playing with that was a fun about week-long episode [ …] I do not any longer have any of my paper scribbles left from when working out the successive integrals, but I recall a lot of integrations by parts (but no fancier calculus ingredients than that). I used Mathematica in ‘calculator mode’, mostly to verify steps and to ensure typo-free bookkeeping of terms. Actually using Mathematica (at least in how I tried that) was somewhat dangerous, as multivalued functions arose, and I recall some non-relevant choices. I ended up about four days later with a pretty horrendous expression which you, in the attached [26 terms along the lines of and the like], can see my numerical verification of. From that point on, what remained was about another three days of ‘clean-up’. I was not at all sophisticated about that, mostly using the general technique of ‘messing around’ with basic log and arctan formulas.”


    Of course, to be really confident, it was checked against the numerical value.

    Asked by a reader whether Mathematica can do it, Michael Trott discussed the problem in the Oct. 2012 entry Trott of his blog. He explored, by a Laplace transform technique,222In the physics literature on this topic, the Laplace transform technique was previously used by Orion Ciftja in 2010–11 for the evaluation, in elementary finite terms, of the electrostatic self-energy of a homogeneous unit square Ciftja2 and a homogeneous unit cube Ciftja3. the Newton potential between two homogeneous unit cubes with centers at a distance of , then taking the derivative at ; in this way, basically, reproducing Fornberg’s solution.

    Understanding, and simplifying, the mathematical structure of that blog entry, and of Fornberg’s solution in the first place, motivated a preliminary version of this paper.

  • In 2001, with boundary element methods in mind, Wolfgang Hackbusch Hackbusch devised a direct method to calculate explicit expressions for gravitational or electrostatic potentials, their evaluation being a task often required in weakly singular integral equations. He considered the sixfold integral


    where and , are cuboids of the form

    The monomials in the numerator correspond, e.g., to polynomial approximations of inhomogeneous mass or charge distributions (cf. cubic), or to the possible choice of higher order finite elements in a Galerkin method.

    Hackbusch proceeded by successive integration, from the first integral to the sixth in the order , while carefully controlling the structure of the primitives as linear combinations of terms of a certain form (altogether there are 18 different such forms, see (HackbuschCode, Table (6.5))). The description of the terms, and the proofs of their recursion formulae, extend over 16 pages in his paper. To deal with concrete cases, he assigned the actual computations to a Pascal program HackbuschCode with about 3000 lines of code, including TeX output of the results. Here is a typical result, very much resembling the look and feel of Fornberg’s expression (2):


    As a verification, Hackbusch checked the examples in (Hackbusch, §3.14) against numerical values; whereas we enjoy the luxury to check against his code.

    Inspired by the exploration found in Trott’s blog entry Trott, we present a new general method for the direct evaluation of that is much easier to derive, to describe (about 5 pages), and to code (25 lines of basic Mathematica code;333See the Mathematica notebook coming with the sources at arXiv:2204.02793; it includes all the algorithmic calculations reported below and also some supplementary material. but pen and paper would suffice for all the examples given here).

  • In 1976, to simplify MacMillan’s (MacMillan, pp. 72–80) classical formula from 1930 which extends over pages, Jörg Waldvogel (MR442257, Eq. (15)) obtained that

    the potential of a homogeneous rectangular cuboid at a point , is given by a symmetric formula that has 48 terms if expanded: namely, short and crisp, writing and for sums cyclic in ,


    By nicely exploiting symmetry and some homogeneity (via Euler’s theorem), Waldvogel just had to evaluate a double integral. In addition, subsequently integrating over , he obtained the self-energy of the unit cube,


    a result that has repeatedly been rediscovered, without reference though; see, e.g., (Hackbusch, p. 207) and (Ciftja3, Eq. (6)).

    As a bonus, we show in Section 6 that the method presented here allows us to straightforwardly reproduce the lovely symmetric formula (5).

The main objective of this paper is to give these three spotlights a common frame that allows an understanding of their structure as well as a comparatively simple algorithm to obtain, among other things, all previous results and formulae.

The basic idea and structure of the results in 3D

The simple idea of this paper, quite classical in multivariate integration,444

As it deserves to be better known, Appendix A discusses two historic examples: one in multivariate calculus, attributed to Cauchy, and one in probability theory, due to Montroll.

is to turn a sixfold integral such as (1) into a single integral of a product of three double integrals by means of a Laplace transform (that is, by bringing yet another integral into the fold, followed by a change of order of integration).

Let us describe this for the evaluation of the Newton potential as given by (3). First, inserting the following Laplace transform (where the substitution simplifies things right at the beginning)


gives, by Fubini’s theorem,


where specifies the following function (with , suppressed in the notation) for the parameters and variables with index :


Next, as explained in Section 2, by introducing the error function

simply as the unique odd primitive (antiderivative) of the Gaussian

, we get an algorithmic description of the factors : they are linear combinations of terms of the form

where is a difference between one of the bounds of the outer integral and one of the bounds of the inner integral. Actually, this form allows right away the rather straightforward numerical evaluation of by just integrating the single integral in (8).

Finally, as explained in Section 3, the error function serves merely as an intermediary on the way to the elementary closed form of . We arrive there by a renormalization process, based on partial integration for removing cancelling singularities at of the integrand . In passing, this process also reduces the number of possible terms to eventually just three basic forms—for which, ultimately, there are entries found in a table of integrals such as the one by Prudnikov et al. MR950173 (for the sake of completeness, we give simple proofs of these entries in Appendix B).

As explained in Section 4, this algorithmic construction extends also to the evaluation of force components and yields, in summary:

In the 3D case, the potential and also the force components are linear combinations of terms of the form555Recall the notation that we introduced in Waldvogel’s formula (5).

where satisfies and is, up to a permutation, the component-wise difference between the coordinates of a vertex of and a vertex of ; the coefficients of the terms are rational polynomials in the integral bounds. Note that this theorem reveals a general structure that is already visible in Waldvogel’s formula (5)—an example that effectively corresponds to the case that is shrunk to a point while the potential is properly rescaled. By rewriting the inverse hyperbolic tangent in terms of logarithms, this structure is also present in the concrete cases (2), (4), and (6).

Extensions and limitations of the method

When extending the class of problems to higher dimensions, or other (weakly) singular integral kernels, we hit the limits of multiple integrals being evaluable in terms of elementary closed form expressions:

First, in dealing with force components in Section 4, we note the importance of factors such as in the numerator of (1) to make the method work; by means of an example we argue that general integrals with an kernel are highly unlikely to be evaluable in elementary finite terms at all.

Second, other dimensions than 3D are the topic of Section 5. Whereas integration in elementary finite terms extends easily to the 1D and 2D cases (actually even reducing the number of possible forms in the Main Theorem from three to just two), we argue that the 4D and higher dimensional cases are highly unlikely to be evaluable in elementary finite terms at all.

Therefore, in such cases we have to resort to the (straightforward) numerical evaluation of a single integral; with dimensions or higher posing no problem whatsoever.

2 Evaluation of the Factor

Since it serves only the role of an intermediary, we normalize the error function in a slightly non-standard manner for reasons of convenience:666The standard form of the error function is normalized as

which would yield a lot of distracting factors in our formulae; for see (MR2655347, §7.1). it is simply defined here as the unique odd primitive of the Gaussian, that is


from which we directly obtain the asymptotics as , and , namely


Writing primitives as indefinite integrals (with constants of integration determined implicitly by the given formulae), we get the following result:

Lemma 1

For  , define

with of degree and of degree recursively given by

and initial values taken as and . Then these functions are primitives of the

-moments of

and wrt , that is


The recursions for and combine into the single one


with initial value . Now a straightforward direct calculation confirms inductively, for  , that

Note that since , and are multiplied by with in their recursions, any value for them will do. ∎

A concrete example of the use of this lemma, with , can be found in Section 6, following Eq. (31). We proceed by considering the double indefinite integrals that are the basis for calculating the factor .

Corollary 1

For there holds

where are polynomials of degree and with

which are obtained by expanding and replacing according to the following rules:


Using Lemma 1 we calculate a primitive, written as a double integral,

By expanding in the first integral into powers of , another application of Lemma 1 shows that contributes

whereas expanding in the second integral shows that contributes

Combining the coefficients of and gives the table. The claim about the polynomial degrees follows from the degrees of and . ∎

Now, by substitution and rescaling we get the primitive


Using this primitive we calculate the definite integral defined in (9) as

Algorithmically, with the integral bounds handled as variables, the construction of this expression for by means of Corollary 1 uses simply some arithmetic of rational polynomials. Structurally, we get:

Lemma 2

The factor

is a linear combination of terms of the form777Since , the second one can be dropped if .


where and the exponent of , namely in the first case and in the second, is restricted to

The coefficient of a term belonging to with and is a homogeneous polynomial of degree in .

Note that Lemma 2 enumerates possible terms. In concrete cases, by partial symmetries, some of them may vanish, cancel, or combine.


By construction, is a linear combination of terms of the form

with and taken from the set of differences of the integral bounds; the coefficients are of the stated form. Since is an even function of , by linear independence of the even Gaussians and the odd error functions (listed by different values of the parameter ) over the field of rational functions, only terms that are themselves even functions of can have non-zero coefficients. Hence it suffices to take only those exponents into account which are even in the first case and odd in the second.∎

Remark 1

Though both forms of the terms in (14) behave as close to , all those singularities must cancel since obviously, as ,

Example 1

We consider two specific cases of the integral

Here, the exponents of are restricted to . The first specific case, with difference set , exhibits all of the 15 possible terms enumerated in Lemma 2:


The singularities of those 12 terms that behave as and close to must cancel since, as ,

In the second specific case, with the symmetric difference set , fewer terms than enumerated in Lemma 2 have non-zero coefficients (there are no terms with and ):

Once again, the singularities at cancel since

3 Evaluation of the Potential

By (8) the Newton potential (3) is given as

Lemma 2 shows that is a linear combination of terms of the form


where and the triple is, up to a permutation, the component-wise difference between the coordinates of a vertex of and a vertex of .

All of these terms are integrable at , but none of them (except the last one with ) is integrable at . So a direct, term-wise integration of is bound to fail. Yet, those singularities must cancel since, as ,

So, to remove the cancelling singularities, some kind of renormalization is called for. We write

and apply, following an idea explored by Trott Trott, repeated integration by parts as long as there are any singular terms left in the reminder integrals.

Specifically, if there is sufficient decay at , we have

Since the functions will be analytic at , the sum of all the boundary terms must cancel in the end: we therefore do not calculate them in the first place. This way, by repeated application of the replacement rule

until no term is left to be replaced, the integrand is eventually transformed into a renormalized one, denoted by , such that


and the integral can, finally, be evaluated term-wise.

Because the derivatives of and wrt are and , these replacement rules do not introduce any new forms of terms; for  , there are effectively just four rules: