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 “10digit 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] [slightly edited for a better fit]With some background in theoretical mechanics, the challenge amounts to evaluating the weakly singular sixfold integral
(1) 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 with^{1}^{1}1In a personal communication Fornberg kindly shared his recollections: “Playing with that was a fun about weeklong 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 typofree bookkeeping of terms. Actually using Mathematica (at least in how I tried that) was somewhat dangerous, as multivalued functions arose, and I recall some nonrelevant 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 ‘cleanup’. I was not at all sophisticated about that, mostly using the general technique of ‘messing around’ with basic log and arctan formulas.”
(2) 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,^{2}^{2}2In 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 selfenergy 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
(3) 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):
(4) 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;^{3}^{3}3See 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 ,
(5) 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 selfenergy of the unit cube,
(6) a result that has repeatedly been rediscovered, without reference though; see, e.g., (Hackbusch, p. 207) and (Ciftja3, Eq. (6)).
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,^{4}^{4}4
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)
(7) 
gives, by Fubini’s theorem,
(8) 
where specifies the following function (with , suppressed in the notation) for the parameters and variables with index :
(9) 
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 formwhere 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 form^{5}^{5}5Recall the notation that we introduced in Waldvogel’s formula (5).
where satisfies and is, up to a permutation, the componentwise 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 nonstandard manner for reasons of convenience:^{6}^{6}6The standard form of the error function is normalized as
(10) 
from which we directly obtain the asymptotics as , and , namely
(11) 
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 isProof
The recursions for and combine into the single one
(12a)  
(12b) 
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:
Proof
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
(13) 
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 form^{7}^{7}7Since , the second one can be dropped if .
(14) 
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.
Proof
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 nonzero 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:
(15) 
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 nonzero 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
(16a)  
(16b)  
(16c)  
(16d) 
where and the triple is, up to a permutation, the componentwise 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, termwise 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
(17) 
and the integral can, finally, be evaluated termwise.
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:
(18a)  
(18b)  
(18c)  