Computation of Power Law Equilibrium Measures on Balls of Arbitrary Dimension

We present a numerical approach for computing attractive-repulsive power law equilibrium measures in arbitrary dimension. We prove new recurrence relationships for radial Jacobi polynomials on $d$-dimensional ball domains, providing a substantial generalization of the work based on recurrence relationships of Riesz potentials on arbitrary dimensional balls. Among the attractive features of the numerical method are good efficiency due to recursively generated banded and approximately banded Riesz potential operators and computational complexity independent of the dimension $d$, in stark contrast to the widely used particle swarm simulation approaches for these problems which scale catastrophically with the dimension. We present several numerical experiments to showcase the accuracy and applicability of the method and discuss how our method compares with alternative numerical approaches and conjectured analytical solutions which exist for certain special cases. Finally, we discuss how our method can be used to explore the analytically poorly understood gap formation boundary to spherical shell support.

1. Introduction. In this paper we present a novel banded and approximately banded spectral method for the computation of equilibrium measures with power law kernels on ball domains of arbitrary dimension, with computational cost independent of the dimension. Equilibrium measure problems naturally arise as the continuous limit of particle swarm systems and as such find various applications in the physics and biology, ranging from modelling cellular movement and interactions of charged particles to the flocking behaviour of birds [76,19,50,14,32,12,11,5]. Consider the problem of finding stationary states of a discrete system of N particles with a pairwise interaction potential Kp|x i´xj |q: where the f term correponds to self-propulsion and friction, see [15]. In the continuous limit N Ñ 8 and in the absence of friction the evolution of the system is governed by the aggregation equation with the constant E{2 given by the minimal total energy, which constitutes a more tractable problem. From a numerical perspective this minimization problem remains extremely difficult, however, as naively it would require us to minimize over the space of all positive densities of unknown radius. For the important special case problem of power law kernels, that is Kp|x´y|q " |x´y| α α´| x´y| β β , (1.2) there are some existence and uniqueness results [35,10,18,17] and for a few specific parameter combinations even analytic solutions [15,16]. The radial symmetry of these kernels suggests that the equilibrium states are themselves radially symmetric. However, this is not true and the loss of radial symmetry for local and global minimizers is very subtle and not theoretically shown even in the case of the power-law kernels (1.2), see [5]. Only very recently [17], we have been able to show that the explicit radially symmetric solutions found in particular ranges in [15] are in fact the global minimizers of the total interaction energy with K given by (1.2) with parameters inside the uniqueness range in [35]. The authors in [17] also give a generic family of kernels leading to non radially symmetric global minimizers.
In the simplest cases, the support of ρ is just a single interval in one-dimension [15,27], a disk in two dimensions and a d-dimensional ball in higher dimensions [15]. For different, in particular higher, powers α and β different support domains for the equilibrium measure may be observed when performing particle simulations, such as two interval systems in one-dimension [15,27], annuli in two dimensions and hyperspherical shells in higher dimensions [15]. While this is known from numerical particle swarm simulations, there are as of now very few analytic results on this gap formation phenomenon and to our knowledge the only numerical results dealing with the continuous gap formation is the one-dimensional ultraspherical spectral method described in [27]. The present paper sets the foundations for studying these phenomena in higher dimensions and can be seen as an arbitrary dimensional generalization of [27], with many of the obtained results mirroring the one-dimensional case. As such, the method is also related to [42], in which the last author introduced a Chebyshev spectral method for the computation of potential theory equilibrium measures with kernels K 0 " logp|x´y|q. The path to justifying similar methods in higher dimensions is fundamentally different because the one-dimensional method was able to leverage the structure of Riemann-Liouville fractional integrals. In constrast, their higher dimensional analogues called Riesz potentials are significantly more restrictive.
A lot of the material required to develop this method is dispersed across different subdisciplines and fields. Equilibrium measure problems, themselves often discussed in biology or physics publications due to interest in their applications [13,37], have close ties to mathematical potential theory [65]. The power law potentials used in this paper are in particular related to Riesz potentials and fractional Laplacians on balls [15,21]. Much of the groundwork for the modern theory of these concepts was laid by Rubin and Samko [66,64] and while some of their papers and books are available in English much of their earlier important work appears to remain untranslated from their Russian originals, see e.g. [31,62]. Furthermore, the theory of Zernike polynomials, a disk generalization of Jacobi polynomials, has been developed in tandem with applied fields such as optical abberrations research [9,36,72,73], which has led to a number of different notations and conventions for radially shifted Jacobi polynomials for d-dimensional balls which feature prominently in our numerical method. As a result of all of the above, we provide a detailed overview of the relevant results in the early sections of the paper, while providing references to both standard monographs and important recent results.
The sections of this paper are organized as follows: Section 1 serves as an introduction, giving an overview of the required concepts for power law equilibrium measures, sparse spectral methods using Jacobi polynomials in one and higher dimensions as well as the theory of fractional Laplacians and Riesz potentials. After collecting auxilliary lemmas and theorems section 2 presents new results and proofs for recurrence relationships and explicit formulae for power law integrals of Jacobi polynomials and Gaussian hypergeometric functions on d-dimensional balls. Section 3 contains a discussion of how these results may be applied to compute equilibrium measures, discusses the banded and approximately banded structures found for appropriate basis choices and finishes with a discussion of regularization and stability for this approach. Section 4 contains extensive numerical validation by means of comparison to special case analytic solutions as well as particle simulations and explores conjectures and other analytically unknown properties such as the uniqueness of solutions and the boundary of the gap forming behaviour for high parameter regimes. We close with a review of the method's properties and future research pathways.

Power law equilibrium measures in higher dimensions.
Definition 1.1 (Equilibrium measure). Given a kernel K : R d Ñ R, x, y P R d we define the equilibrium mesaure to be the measure dρ " ρpxqdx such that the expression When M " 1 the equilibrium measure is a probability measure.
Note that the energy as defined in (1.3) is a Lyapunov functional of the evolution equation stated in (1.1). For the power law kernel in (1.2) one can derive the Euler-Lagrange conditions as outlined above, more details and references on this approach may be found in [4,10,15]. In case the global minimizers are unique [35,17], then the equilibrium measure may instead be found by solving the following governing equation 1 α ż supppρq |x´y| α ρpyqdy´1 β ż supppρq |x´y| β ρpyqdy " E.
In any case, they are candidates for being global minimizers of the interaction energy (1.3). By the radial symmetry of the problem, assuming that the measure ρpyq is supported on the ball B R " tx P R d , |x| ď Ru, the d dimensional governing equation takes the form [15]: as before subject to a mass condition Our spectral method will be developed to find the minimizer among radial densities supported on a ball. General conditions on the interaction potential leading to the radial symmetry of local/global minimizers of (1.3) have been recently shown in [17].
1.2. Jacobi polynomials and sparse spectral methods. We say polynomials p n pxq are orthogonal on a given domain Ω Ď R with respect to some weight function wpxq if they satisfy ż Ω wpxqp n pxqp m pxqdx " c n,m δ n,m , where c n,m are constants and δ n,m is the Kronecker delta. A complete orthogonal polynomial basis may be used to expand sufficiently smooth functions via with constant coefficients f n and This expansion can be used to numerically approximate functions by only summing over the first N number of terms. There has been much interest in the mathematical properties of such polynomials, leading to a multitude of high quality formula collections [41] and monographs [20], to which we refer for details beyond the scope of this paper. Various complete orthogonal polynomial bases, most notably among them the ultraspherical, Jacobi and Chebyshev polynomials, are used to numerically approximate functions and solve differential and integral equations using sparse, in particular banded, operators. A recent in-depth survey may be found in [44]. These methods converge fast with competitive computational cost due to sparsity. Open source software implementations are available to perform such computations in the form of Ap-proxFun.jl [43] and FastTransforms.jl [67,68], which also form the framework within which the method introduced in this paper has been implemented for the numerical experiments in section 4.
In the remainder of this section we will describe properties of one of the most ubiquitously used one-dimensional bases for these methods, the Jacobi polynomials, which also often feature as components in various higher dimensional bases, including ball domains of arbitrary dimension. The Jacobi polynomials are conventionally denoted P pa,bq n pxq and are orthogonal on the interval r´1, 1s with respect to the weight function with a, b ą´1. The special case of a " b corresponds to the ultraspherical or Gegenbauer polynomials with slightly different normalization, while the case a " b " 0 corresponds to the Legendre polynomials [41,Table 18.3.1]. We will find the following explicit representation of the n-th Jacobi polynomial to be useful for our purposes [41, 18.5.7]: pn`1qpn`a`b`1qp2n`a`bq . This recurrence may be used to define a tridiagonal multiplication operator X acting on a function expanded in the Jacobi polynomials Together with element-wise addition and subtraction, as well as the raising and lowering operators defined below, these properties make up the basic toolbox for function approximation and arithmetic used in sparse spectral methods. The range of recurrences and symmetry properties of the Jacobi polynomials is extensive and we will not attempt to reproduce it fully here. A large collection of these may be found in [41, 18.3-18.18] and a number of additional ladder operators are listed in [45]. Instead we selectively state some results which we will be using in the development of our method, beginning with the symmetry relation [41, 18.6.1]: (1.7) P pa,bq n p´xq " p´1q n P pb,aq n pxq.
The following two relationships between parameters and degree allow one to take care of additional weight terms [41, 18.9.6]: pn`a 2`b 2`1 q P pa,bq n`1 pxq`p n`a`1q pn`a 2`b 2`1 q P pa,bq n pxq . (1.9) Basis conversions between different parameter Jacobi polynomials may also be accomplished with recurrence relations such as [41, 18.9.5]: p2n`a`b`1q P pa`1,bq n pxq´p n`bq p2n`a`b`1q P pa`1,bq n´1 pxq , (1.10) pxq . (1.11) Boundary evaluations for the Jacobi polynomials can be made using [41, pxq, (1.14) Sometimes applications may instead want to make use of shifted Jacobi polynomials, which via straightforward variable substitutions move the domain of orthogonality to e.g. r0, 1s instead. As we will see in the next section, a certain type of shifted Jacobi polynomial, namely P pa,bq n p2r 2´1 q, r P p0, 1q, will feature prominently in the development of our method for equilibrium measures. We thus also reproduce variants of the above-listed properties for the shifted radial basis P pa,bq n p2r 2´1 q in Appendix 6.2 as a reference.
1.3. Sparse spectral methods on unit balls of arbitrary dimension. Substantial progress has been made over the past years in expanding the tools described in the previous section to various higher dimensional domains such as triangles [45,46], disks [78], wedges [47], disk slices and trapeziums [70], quadratic curves [49] as well as surfaces of revolution [48] using multivariate orthogonal polynomials [20]. The most ubiquitous set of orthogonal polynomials on the disk are the so-called Zernike polynomials [83,82], a set of bivariate orthogonal polynomials typically given in terms of polar coordinates r and θ. The radial part of the Zernike polynomials can be written in terms of shifted Jacobi polynomials and as a result they inherit many of their useful properties. This turns out to be the case for higher dimensional ball domains as well. Due to the Zernike polynomials' numerous applications [79] especially in the study of optical abberations [9,58,36,72,80,73] and atmospheric wavefronts [40,59], one finds that their treatment in the scientific literature is dispersed across different fields in multiple disciplines, with normalization and ordering of the polynomials often inconsistent between different publications, see e.g. [26,40,80]. The spectral method we propose in this paper uses a generalized set of polynomials on the d dimensional unit ball. In two dimensions, this basis is usually referred to as the generalized Zernike polynomials [20,81,77,2,3], or sometimes simply Jacobi polynomials on the unit disk D " tx P R 2 , |x| ď 1u. They are orthogonal on D with respect to the weight function The radial part [78] of the generalized Zernike polynomials reads Q pk,mq n p|x|q " |x| m P pk,mq n p2|x| 2´1 q, where P pk,mq n p2|x| 2´1 q are shifted one-dimensional Jacobi polynomials as introduced in the previous section. Typically, an additional normalization constant is also present. The angle component of the generalized Zernike polynomials corresponds to Fourier modes, see e.g. [78,20]. For radially symmetric functions only the m " 0 modes of the basis are relevant, meaning that an appropriate basis of orthogonal polynomials on the d-dimensional unit ball is simply given by the polynomials P pk, d´2 2 q n p2|x| 2´1 q, which are orthogonal with respect to p1´|x| 2 q k . For k " 0 the above definitions are equivalent to the standard Zernike polynomials, making 'generalized Zernike polynomials' an appropriate description. We note that the ordinary Zernike polynomials thus have a similar relationship to the generalized Zernike polynomials as the Legendre polynomials do to the ultraspherical or Jacobi polynomials, cf. [6,81,78]. As a direct consequence of the generalized Zernike polynomials' defining relationship with shifted Jacobi polynomials, they satisfy a recurrence relationship and various other properties in the radial parameter, see Appendix 6.2.
1.4. The fractional Laplacian of radial functions. The Laplace operator and its non-local generalization in fractional calculus are important tools for the development of our method, so we give a brief overview of the relevant results. First, a standard definition: Definition 1.2 (Fractional Laplace operator). We define the negative fractional Laplace operator p´∆q γ 2 for γ P p0, 2q via the following singular integral where B " Bp0, q denotes a ball of radius around the origin. Equivalently [33] with range of validity γ P p0, dq we can write the fractional Laplacian as the inverse of the Riesz potential, thus denoted p´∆q´γ 2 : Kwaśnicki recently published a mostly self-contained and comprehensive survey of many equivalent definitions of the fractional Laplacian in [33]. The Riesz potential will be introduced and discussed in more detail in section 1.5. Without stating the construction explicitly, we note that the definition of the fractional Laplace operator as a singular integral can be extended for higher parameters γ P p0, kq where k is an even integer using certain centered differences of f , see [71,21].
As we will primarily be interested in functions f pxq " f prq with r " |x| with supppf q " B R Ă R d , we now shift focus from the whole space to ball domains and rotationally symmetric functions. For the ordinary Laplace operator ∆ acting on a radially symmetric function one has A closely related formula for the negative fractional Laplacian is known but due to nonlocality involves several non-trivial integrals, see [23,Lemma 7.1]. Dyda, Kuznetsov and Kwaśnicki recently made substantial progress in the theory of fractional Laplacians in [21], in particular also for ball domains. They proved the following result valid for γ ą´d and x P R d : If f pxq :" V pxqG mn pq p |x| 2 ; a; bq, then where G is the so-called Meijer G-function [41, 16.17.1] defined as a significant further generalization of the generalized hypergeometric functions p F q and V pxq is a solid harmonic polynomial of order l, i.e. a solution to Laplace's equation which is polynomial in |x|. In particular, if we set V pxq " 1 then l " 0. This very general formula can be used to derive several more specialized results. Of interest to us is the following statement which holds true for the unit ball x P B 1 for certain parameter ranges for γ [21, Theorem 3]: with V pxq as above and where P pa,bq n pxq are the Jacobi polynomials. While the authors in [21] state the range of validity of this derived formula as γ ą 0, their proof of it can be modified to also work when γ P p´2, 0q. As we will be needing this additional range we present a proof of this fact on the basis of the general results in [21] in Appendix 6.1. Note that (1.18) implies that the Jacobi polynomials with the stated parameter choices form a complete orthogonal basis of eigenfunctions for the fractional Laplacian p´∆q γ 2 weighted by p1´|x| 2 q γ 2 . 1.5. Riesz and power law potentials on balls. Riesz potentials pervade the wider mathematical context of equilibrium measures, with many results in Riesz potential theory going back to the original study of these constructions by the eponymous M. Riesz [56,57]. Some standard surveys which include modern and historical aspects are [55,34,66,1]. In the context of mathematical potential theory, a potential U ρ K pxq is generally an integral operator over the support of a measure ρ, with a given kernel Kpx, yq and where supppρq Ď R d . When dρpyq " ρpyqdy we will instead write it as acting on the density, i.e.
In many applications what one is primarily interested in are so-called convolution kernels Kpx, yq " Kpx´yq or even more specifically kernels which depend only on the Euclidean distance Kpx, yq " Kp|x´y|q. The Riesz potential is of this kind and includes the well-known logarithmic and classical Newtonian kernels as limiting and special cases respectively, cf. [34].
Definition 1.4 (Riesz potential). The Riesz potential of a function ρpyq with y P supppρq Ď R d is defined as the integral We have already seen in section 1.4 that the Riesz potential is the inverse of the negative fractional Laplacian given that γ P p0, dq [1,33]. The normalization constant in the Riesz kernel is chosen such that the following semi-group property holds for Riesz potentials given the condition pγ`βq, γ, β P p0, dq, cf. [34]: (1.19) or alternatively in the notation of the fractional Laplace operator: Most of the standard results for Riesz potentials, such as those found in [34,66], are only valid for parameters γ P p0, dq. As with the fractional Laplacian we have started with definitions and properties valid for the whole space R d . In the context of equilibrium measures we are interested in Riesz pontentials on balls of radius R, i.e. B R " tx P R d , |x| ď Ru or equivalently with measures where supppρq " B R . Such potentials have several unique properties compared to those on the whole space and were extensively covered in the seminal works of Rubin and Samko, see e.g. [64,61,66,60]. Of particular note is [61], in which Rubin showed that for γ P p0, dq and radially symmetric functions the Riesz potential on a ball of arbitrary dimension can be reduced to two nested one-dimensional fractional Riemann-Liouville integrals and [63,64] in which he introduced two new types of fractional integrals on balls, c.f. [22,66], the composition of which yields the Riesz potential. The power law kernels we are interested in for the purposes of equilibrium measures are Riesz kernels on balls but in the equilibrium measure literature, c.f. [15,27], the following notation with different normalization is sometimes preferred and will also be used in this paper henceforth: Definition 1.5 (Power law kernel). We define the power law kernel K α p|x´y|q, with x, y P R d by Definition 1.6 (Power law potential). The power law potential of a function ρpyq with y P supppρq Ď R d is defined as the integral Remark 1.7. The power law kernel as defined in Definition 1.5 corresponds to the Riesz potential in Definition 1.3 when α " γ´d with different normalization. Note that in contrast to the conventional Riesz potential the normalization in Definition 1.5 is chosen such that a kernel K α p|x´y|q will always result in an attractive potential while a negative sign´K α p|x´y|q will result in a repulsive potential, even when the sign of α changes.
Simple variable substitutions and multiplications with the appropriate constants allow for the conversion of all valid statements for Riesz potentials when γ´d " α P p´d, 0q. As many interesting equilibrium measure phenomena occur exlusively in high parameter regimes, compare the high parameter results in [15,27], we want to loosen the restriction of α P p´d, 0q or equivalently γ P p0, dq but whether this can be done depends fundamentally on the function the potential acts on. We will see important special cases in which this can be done in the next section.
2. Results for power law potential of radial functions. In [27] the authors derived recurrence relations for left and right-handed Riemann-Liouville fractional integrals in weighted ultraspherical polynomial bases and used them to prove a recurrence relationship for power law kernel integrals in the one-dimensional case. A natural higher dimensional analogue of Riemann-Liouville fractional integrals is found in the Riesz potential. Attempting a straightforward higher dimensional version of the approach in [27] fails because the Riemann-Liouville integrals satisfy properties, e.g. the semi-group property in (1.19), which for Riesz potentials are only valid in very limited parameter ranges. In this section we thus take a different approach, which nevertheless leads to a true generalization of the results in [27] to arbitrary dimension d. Throughout this section, we write Bpx, yq " ΓpxqΓpyq Γpx`yq to denote the Beta function and B R " tx P R d , |x| ď Ru for a d-dimensional ball of radius R centered at the origin. A straightforward rescaling allows us to set the radius of the ball domain to R " 1 without loss of generality, which significantly simplifies our notation when dealing with Jacobi polynomials.
2.1. Auxiliary results. The purpose of this section is to build and gather the mathematical tools needed to prove our main theorems. The following finite sum expression for hypergeometric p F q functions can be found in the extensive works of Prudnikov, Brychkov and Marichev [54, 5.3.6.3]: [54, p.798]. The proof of Theorem 2.16 will be using the case p " q " m " 1, which simplifies to: Next, we discuss some important lemmas which follow from results in potential theory.
The following general power law integral with x " 0 will be useful in determining certain integration constants and previously been used in the derivation of certain analytical power law equilibrium measure results, cf. [15, Appendix A]: Proof. The d " 1 variant of this formula may be obtained after splitting the expression into two integrals from 0 to R and´R to 0 respectively, then using the general formula [53, 3.3.2.4]: The higher dimensional variant then follows via a reduction formula to the one- We prove a generalization of Lemma 2.1 for Jacobi polynomials on the unit ball.
Lemma 2.2. On the d-dimensional unit ball B 1 the power law potential of the weighted Jacobi polynomial p1´|y| 2 q m P pa,bq n p2|y| 2´1 q at x " 0 evaluates: Proof. Via the explicit representation of Jacobi polynomials in (6.2) we can expand this integral into a form from which the result follows via Lemma 2.1.
The hypergeometric function variant then follows by definition of the beta and 3 F 2 hypergeometric function, see e.g. the properties in [41, 5.12, 16.2].
The following lemma is fundamental for the method introduced in this paper.
Lemma 2.3. On B R the power law potential of the function pR 2´| y| 2 q´α`d 2 , with power α P p´d, 2´dq, evaluates to an explicitly known constant: Proof. A proof of this result based on point inversion (Kelvin transforms) can be found in [34,Appendix].
The next lemma is a direct generalization of Lemma 2.3. While the general case proof of this generalization requires some additional work and we thus only provide a reference to its proof, it is nevertheless worthwhile for our later results to sketch how incremental generalizations of Lemma 2.3 can be obtained as was discussed in [15]. Keeping α P p´d, 2´dq, we note the action of the ordinary Laplace operator ∆ on a power law integral with increased power α`2: Inverting the Laplace operator here is straightforward and leads to where the integration constant was fixed by evaluation at x " 0 using Lemma 2.1 above. Finally replacing α`2 with α to retain the same kernel structure leads us to the following result which is now valid for α P p2´d, 4´dq: Repeating this process we can successively find solutions for higher powers. One finds that these are all special cases of the following general lemma: Lemma 2.4. On B R the power law potential of the function pR 2´| y| 2 q ´α`d 2 , with power α P p´d, 2`2 ´dq, evaluates as follows: If furthermore P N 0 , this hypergeometric function reduces to a polynomial. Explicitly: Proof. This generalization of Lemma 2.3 was derived in [28,8,7] on the basis of the above-mentioned connection to the theory of fractional Laplace operators. The Jacobi polynomial variant with order dependence in the parameter follows from their explicit hypergeometric function expression in (1.5) combined with (1.7).
The fact that Lemma 2.4 does not require P Z means that we can separate the power of the weight from that of the kernel, a fact we make explicit in the following corollary: Corollary 2.5. On B R the power law potential of the function pR 2´| y| 2 q m´α`d 2 with power α P p´d, 2`2m´dq, m P N 0 and β ą´d, evaluates as follows: Proof. Starting with ş B R |x´y| β pR 2´| y| 2 q ´β`d 2 dy and setting " m`β´α 2 , the result follows directly from Lemma 2.4.

2.2.
Laplacians of shifted Jacobi polynomials. This section collects further useful lemmas related to recurrences to compute the ordinary Laplace operator of various radial hypergeometric functions, including shifted Jacobi polynomials.
Lemma 2.6. The Laplace operator ∆ acting on the shifted Jacobi polynomials P pa,bq n p2|y| 2´1 q on the d-dimensional unit ball, y P B R " tx P R d , |x| ď Ru, evaluates: Proof. This is easily obtained via the representation for the Laplace operator on radial functions in (1.16) in combination with the derivative relations for Jacobi polynomials in (1.14-1. 15).
Note that the terms on the right-hand side of Lemma 2.6 can be transformed into a basis with consistent parameters by using shift operators such as those listed in (1.10-1.11). Unfortunately, when doing so the expressions do not generally simplify further and end up involving several Jacobi polynomials of different polynomial orders. That being said, we can obtain significant further simplification for certain choices of a and b. Specifically, this is the case for the radially symmetric Jacobi polynonimals with basis parameters as seen in section 1.3.
Lemma 2.7. The Laplace operator ∆ acting on the shifted Jacobi polynomials P pa, d´2 2 q n p2|y| 2´1 q on the d-dimensional unit ball B 1 " tx P R d , |x| ď 1u evaluates: Proof. This can be obtained from Lemma 2.6 via (1.10-1.11).
For the case a " 0 and d " 2, i.e. for the case of ordinary Zernike polynomials on the unit disk, similar results to Lemma 2.7 for the Laplacian and inverse Laplacian were discussed in [30] while also including their angular components.
2.3. Power law potentials of Jacobi polynomials on the unit ball. With all of the above results, we are now ready to investigate the action of the ball Riesz potential operator on radially symmetric functions expanded in orthogonal polynomials on arbitrary dimensional balls, where the radial part is expanded in terms of one-dimensional Jacobi polynomials P pa,bq n p2r 2´1 q with r 2 " |y| 2 as discussed in sections 1.2 and 1.3. First we prove the following result in low parameter ranges: Theorem 2.8. On the d-dimensional unit ball B 1 " tx P R d , |x| ď 1u the power law potential of the n-th radial Jacobi polynomial P p´α`d 2 , d´2 2 q n p2|y| 2´1 q weighted with p1´|y| 2 q´α`d 2 evaluates explicitly as follows in the parameter range α P p´d, 2´dq: Proof. There are a number of ways to derive this result from what we have so far discussed, the most straightforward of which is to use results based on those of Dyda, Kuznetsov and Kwaśnicki [21, Theorem 3] which we reproduced in (1.18) and extended in Appendix 6.1. When γ P p0, dq the Riesz potential is the inverse of the fractional Laplace operator, motivating the use of (1.18) to find a related result for the Riesz potential. With V pxq " 1 and thus l " 0 the expression reads Setting γ "´pα`dq then leads to The inherited validity of this expression is α P p´d, 2´dq, see Appendix 6.1. The above result holds for Riesz potentials, which as we saw in section 1.5 has different normalization constants to the equivalent power law potentials we are working with here. Taking a glance at Definitions 1.2 and 1.4, this difference can be accounted for by simply multiplying both sides in the above equation by Carrying out the appropriate cancellations we are then left with the first stated result. To obtain the second variation from the previous result, note that Note that this is the same factor that already appeared in Lemma 2.3.
Corollary 2.9. On the d-dimensional unit ball B 1 " tx P R d , |x| ď 1u the power law potential of p1´|y| 2 q´α`d 2 P p´α`d 2 , d´2 2 q n`1 p2|y| 2´1 q can be computed using a three term recurrence relationship: Proof. This immediately follows from expanding the Jacobi polynomial on the right-hand side of the formula given in Theorem 2.8 using the classical three-term recurrence for Jacobi polynomials in (1.6).
Remark 2.10. We can easily convince ourselves that the diagonality result in Corollary 2.9 cannot in general be true for higher parameters α ą 2´d, i.e. that diagonal operators only exist for α P p´d, 2´dq. First, note that in order for solutions to exist when α ą 2´d we need the weight parameter to satisfy ą 0 in Lemma 2.4. Obviously one necessary condition for the operator to be diagonal is that the 0-th order polynomial is mapped to a constant but the right-hand side of Lemma 2.4 tells us that this only occurs when " 0, as otherwise the hypergeometric function (or equivalently the Jacobi polynomial) is not constant with respect to |x|. Corollary 2.9 covers precisely the case " 0.
While the operators cannot be diagonal in higher parameter regimes, we can nevertheless obtain banded operators with explicit elements which is similarly efficient for computing purposes. The following theorem shows how to obtain higher bandwidth operators for higher parameter ranges and is the fundamental result for our sparse spectral method: Theorem 2.11. On the d-dimensional unit ball B 1 " tx P R d , |x| ď 1u the power law potential of the n-th radial Jacobi polynomial P p1´α`d 2 , d´2 2 q n p2|y| 2´1 q weighted with p1´|y| 2 q 1´α`d 2 evaluates as follows in the parameter range α P p2´d, 4´dq: where the constants are given by Proof. First we note how the d-dimensional ordinary Laplace operator ∆ x in x acts on power law integrals with increased power α`2, where α P p´d, 2´dq as before: The final step makes use of Theorem 2.8. Using Lemma 2.7 to invert the Laplacian in this form is ill-advised, as it will reduce the first parameter of the Jacobi polynomial on the right-hand side by 2 with no guarantee that the first Jacobi parameter will remain greater than´1. To avoid this happening, we first raise the basis parameters using the shift operators in (1.10-1.11). This first yields 2 q n´1`2 |x| 2´1˘, and in a second application With this conversion and because of the radial symmetry of the integrals we can invert the Laplacian obtaining ż Finally, the constant term c d may be fixed by evaluating the integral expression at x " 0, which can be done using Lemma 2.2. Using the evaluation operations in (1.12-1.13) then leaves us with the following equation: which can be solved for c d for all n. The result is as nice as one could hope for: Now replacing α`2 with α to retain the same form for the kernel as before yields the desired result, with inherited range of validity α P p2´d, 4´dq.
Remark 2.12. The method used to prove Theorem 2.11 can be used in a straightforward way to find banded operators for all P N 0 for the power law potential where the bandwidth increases with , starting from the diagonal case in Theorem 2.8 with " 0. The inherited range of validity for the successively obtained expressions is p2 ´d, 2`2 ´dq and the respective operators each have exactly 2 `1 bands in our basis. We omit these further explicit derivations, as the involved expressions become lengthy and do not require any new ideas.
2.4. Decoupling the weight and kernel power. When dealing with both an attractive and a repulsive power law integral at the same time, we need additional tools to treat both simultaneously in one consistent basis. We begin by proving a result for power law potentials of weighted monomials on balls: Lemma 2.13. On B R the power law potential of |y| 2k p1´|y| 2 q ´α`d 2 with k ą 0 and l P N 0 satisfies the following recurrence relation: Proof. Expanding one of the powers of the weight in the right-most term with the highest weight leads directly to Remark 2.14. Lemma 2.13 tells us that on B R the Riesz potential of finite sums of terms of the form |y| 2k pR 2´| y| 2 q ´α`d 2 for varying k can be reduced to a sum depending exclusively on solutions of the k " 0 case in Lemma 2.3 with increasing weight parameters. For example, the case k " 1 may be evaluated as follows: Note that even in the general case, we can evaluate all the terms on the right-hand side given α P p´d, 2`2 ´dq and that if is chosen such that the k " 0 term is a polynomial, i.e. P N, then all higher order terms are also polynomials.
The above results could in principle be used to design a dense spectral method for equilibrium measures in a basis of weighted monomials but such methods would be computationally expensive as well as significantly less robust. The above is easily modified to instead give a generic Jacobi polynomial recurrence: Lemma 2.15. On the unit ball B 1 , the Jacobi polynomials P p ´α`d 2 , d´2 2 q n p2|y| 2´1 q with weight p1´|y| 2 q ´α`d 2 satisfy the following two-term recurrence relationship: Proof. This is a consequence of the well known weight reduction recurrence relationship of the Jacobi polynomials in (1.8-1.9). In particular, with x Ñ`2|y| 2´1t he recurrence in (1.9) takes the following form: 2pn`a 2`b 2`1 q P pa,bq n p2|y| 2´1 q, which with some straightforward rearranging turns into P pa,bq n`1 p2|y| 2´1 q " pn`a`1q pn`1q P pa,bq n p2|y| 2´1 q 2pn`a 2`b 2`1 q pn`1q p1´|y| 2 qP pa`1,bq n p2|y| 2´1 q.
The result then follows from linearity after plugging this into the weighted power law integral with appropriate parameters pa, bq.
As before, advancing to step n`1 using this recurrence uses knowledge of the previous step n in the same basis as well as one with higher weight parameter, meaning that the computational cost of this recurrence does not scale linearly with n. Normally such a general recurrence relationship which holds true irrespective of the structure of the kernel would be of little use, as the solution for high orders requires the luxury of extensive knowledge about initial results with higher weight terms, in particular for the 0-th order terms. In the present case of the power law potential, however, exactly this knowledge is given by Lemma 2.4. In conjunction with Theorem 2.11 we could thus in principle use this to compute equilibrium measures even in the attractiverepulsive case where two potentials with different powers are present, while remaining in a consistent polynomial basis. In practice however this recurrence would also be a bad choice as it scales poorly with n and suffers from numerical instability as the polynomial orders increase. To circumvent this problem, we first take a detour and look for a direct generic solution to the power law integral of Jacobi polynomials: Theorem 2.16. On the d-dimensional unit ball B 1 the power law potential, with power α P p´d, 2`2m´dq, m P N 0 and β ą´d, of the n-th weighted radial Jacobi polynomial reduces to a Gaussian hypergeometric function as follows: Proof. Using the explicit representation found on the second line of (6.2), the left hand side power law integral can be rewritten Γpm`n´α`d 2`1 q n!Γpm`n´α 2 q n ÿ k"0 p´1q kˆn k˙Γ pm`n`k´α 2 q Γpm`k´α`d 2`1 q ż B1 |x´y| β p1´|y| 2 q k`m´α`d 2 dy.
Applying Corollary 2.5, we can resolve the integral despite the power of the weight and kernel being independent to obtain: n the above we expanded the Beta function into its Gamma function representation and performed the then obvious cancellation of terms. With some careful algebraic manipulations this can then be cast into the following form In general, any hypergeometric p F q function with an equal upper and lower parameter reduces to a p´1 F q´1 function, a fact that is straightforwardly seen from their series representation. Explicitly for the above appearing 3 F 2 we have: where we remind ourselves that p F q functions are symmetric when exchanging upper with other upper parameters and likewise for the lower parameters. Applying this to (2.2) and simplifying the constants concludes the proof of the theorem.
Remark 2.17. The specific form of Theorem 2.16 was chosen to suit our needs but it is easy to see that a simple substitution yields the generic weight parameter form without any change to the proof.
The importance of Theorem 2.16 for our method is that it can be used in conjunction with highly efficient and accurate implementations for the computation of the Gaussian hypergeometric functions, as implemented in HypergeometricFunctions.jl [69], to compute the operator entries even when α ‰ β. The specific hypergeometric function algorithm we relied on for the numerical experiments in this paper is the one discussed in [51]. This can be done either by directly expanding the n-th variant of the formula or preferably by using the following recurrence relationship where the basis and kernel powers are now decoupled: Corollary 2.18. On the unit ball B 1 , the power law integral of the Jacobi polynomials P pm´α`d 2 , d´2 2 q n p2|y| 2´1 q with weight p1´|y| 2 q m´α`d 2 , α P p´d, 2`2m´dq and β ą´d satisfies the following three term recurrence relationship: Proof. This is obtained by plugging the appropriate parameters obtained in Theorem 2.16 into the following recurrence relationship for Gaussian hypergeometric functions: The fact that there exist recurrence relationships for Gaussian hypergeometric functions of form 2 F 1ˆa` 1 n, b` 2 n c` 3 n ; z˙with i P t´1, 0, 1u which can be derived from contiguous relations [41, 15.5.11-15.5.18] is widely known, cf. [25,24,52] and the references therein, and Mathematica code for the symbolic generation of such recurrences is available [29].
Remark 2.19. The special case where the kernel and weight terms are related via α " β results in a polynomial right hand side, meaning by Theorem 2.16 the operator must have finitely many non-zero subdiagonal bands and by the three term recurrence in Corollary 2.18 it then also has finitely many superdiagonal non-zero bands. The bandwidth of the operator in the polynomial case for any specific parameter α " β can be determined by the order of the right hand side polynomial when n " 0. The minimal bandwidth that is achievable using Jacobi polynomials is determined by the range of validity β P p´d, 2`2m`β´α´dq, i.e. for α " β and α P p´d, 2´dq we can get diagonal operators, for α P p2´d, 4´dq we can at best get tridiagonal operators and so on, consistent with the results of the previous sections.
2.5. The mass condition in arbitrary dimension. As mentioned in section 1.1, equilibrium measure problems of the form discussed in this paper only become well-posed with the addition of a mass condition, i.e. by demanding that the equilibrium measure must satisfy ż B R ρpyqdy " M on its domain of support B R , the d-dimensional ball of radius R. In the computational setting this will require us to be able to evaluate the unit ball integral of sufficiently well-behaved functions expanded in a weighted basis of radial Jacobi polynomials.
Lemma 2.20. Let ρpyq " ρp|y| 2 q, y P B 1 be a function such that ş B1 ρpyqdy " M 1 , with M 1 constant. Furthermore, we assume that its expansion ρpyq " 8 ÿ n"0 ρ n p1´|y| 2 q a P pa, d´2 2 q n p2|y| 2´1 q in weighted radial Jacobi polynomials on the d-dimensional unit ball B 1 exists and satisfies ş B1 ř 8 n"0 |ρ n p1´|y| 2 q a P pa, d´2 2 q n p2|y| 2´1 q|dy ă 8. Then the value of M 1 is determined entirely by the 0-th coefficient, that is: Proof. The domain and radial symmetry of this problem suggests the use of hyperspherical coordinates: where σpS d´1 q " 2π Γpa`1qΓ`d 22 where the last equality relies on the classical orthogonality condition of the Jacobi polynomials in Equation (1.4) where we remind ourselves that P 0 ptq " 1. Combining this with the prior expression yields the stated result.

Description of the numerical method.
In the previous sections we derived recurrence relationships and explicit representations of the power law potential of Jacobi polynomials on d dimensional balls. We now detail how these results can be used to produce banded and approximately banded operators to efficiently compute the power law potential of a given rotationally symmetric function f prq and then discuss how these tools can be used to solve equilibrium measure problems.
3.1. Operator sparsity structure for classical Jacobi polynomials. We have seen in section 2.3 that we can achieve banded power law integral operators acting on coefficient vectors of weighted radial Jacobi polynomials by making an appropriate basis choice. In particular we have seen that for a given power α P p2 ´d, 2`2 ´dq with P N 0 in d dimensions, the choice of basis  gives banded operators with exactly 2 `1 bands. We present some illustrative examples of this in terms of matrix spy plots in Figure 1. Using Theorem 2.8, Theorem 2.11 and the straightforward generalizations one can generate these operators with very high efficiency. An important and interesting special case occurs when α is an even integer, as can e.g. be seen from the results in Theorem 2.16 and Corollary 2.18, namely the power law operator only has finitely many non-zero values, the exact number depending only on the value of α, see the example matrix spy plots in Figure 2. This is consistent with what was shown in one-dimension for ultraspherical polynomials in [27]. In the attractive-repulsive case this special case result means that we can obtain a simultaneously banded operator for any β for which solutions exist. When considering the general attractive-repulsive equilibrium measure problem, one has to consider both power law integrals in a consistent basis in order to apply spectral methods. To do this, we choose the basis such that the attractive α-operator is banded as above, then generate the attractive β-operator in the same basis by using Theorem 2.16 and Corollary 2.18. Again consistent with what was previously observed in one-dimension for ultraspherical polynomials [27], these operators are approximately banded and decay off the main band, see the example matrix spy plots in Figure 3. Improving on the understanding of the one-dimensional case, the general expression obtained in Theorem 2.16 explains the approximate bandedness as each column of the operator represents the coefficients of an expansion of the polynomial order dependent Gaussian hypergeometric function in a Jacobi polynomial basis. As we have obtained a recurrence relationship for this general case it is possible in practice to compute a (a) α "´π, β "´2.4, d " 5 (b) α " 1.9, β " 0.5, d " 2 (c) α " 3.76, β "´2 3 , d " 3 banded approximation for a chosen bandwidth without the wasteful step of having to compute the dense operator first. That being said, suitable convergence is generally obtained already for sizes were dense matrix computations are still feasible.
3.2. Application to equilibrium measure problems. With the established operator structure and recurrences to efficiently generate them, the application of these results to higher dimensional equilibrium measure problems is mostly analogous to that of the one-dimensional ultraspherical method in [27]. As such, we will slightly abbreviate this discussion by omitting the question of solving equilibrium measure problems with only one power and an external potential, which works analogously to what is discussed in [27, section 3.1], and instead focus on the significantly more challenging problem of attractive-repulsive systems with vanishing external potential.
We remind ourselves of the problem we intend to solve: Find the positive density ρpȳq which minimizes the scalar energy Note both ρ and the radius R of its support B R are unknown. To be able to use our Jacobi polynomial results we thus first normalize this expression to the unit ball: wherex " Rx andȳ " Ry and |y|, |x| ď 1. With its support normalized to the unit ball B 1 , the density ρpRyq is now assumed to be given in its expansion in weighted Jacobi polynomials ρpRyq " 8 ÿ n"0 ρ n p1´|y| 2 q ´α`d 2 P p ´α`d 2 , d´2 2 q n p2|y| 2´1 q.

Plugging this into the governing equation we find
Choosing instead to represent this as the action of operators on coefficient vectors as described in section 1.2, we obtain the linear system where the upper index in U a b denotes the power of the kernel and the bottom index the power of the operator which is made banded through the choice of basis. We have seen in the previous section that U α α is banded and U β α is approximately banded. This system can be uniquely solved for the coefficient vector ρ if R is known. Since R is the radius of support of the equilibrium measure, it is in general not known. However, this approach nevertheless allows us to compute a unique corresponding measure for any R, meaning that we have reduced a complicated minimization problem over the space of positie measures to a minimization in the single scalar value R. There is one further complication to address: obtaining the unique ρ corresponding to a given R appears to require knowledge of Epρq, which is certainly not known. This hurdle is overcome by the mass condition. Rearranging the above tô allows us to solve for a coefficient vector which corresponds to ρpRyq E . We can get rid of the constant E by normalizing our measure according to the given mass constraint. This can be done with very high efficiency using the approach described in section 2.5, ultimately allowing us to compute both a unique ρpyq and its corresponding energy E for a given radius R.
Remark 3.1. U α α and U β α are independent of R. Thus, when minimizing the energy by varying R, the operator does not need to be re-computed for each R. Instead, we store U α α and U β α after the first generation and then simply computé Remark 3.2. The dimension d has no direct impact on the computational cost, meaning that due to the rotational symmetry this is a rare case where calculations on the d-hypersphere are just as efficient as those on the interval. This is in stark constrast to the widely used methods of discrete particle swarm simulations to approach these problems which scale catastrophically with d, making even two or three dimensional computations almost impossible in interesting parameter ranges, not to mention higher dimensional problems.
3.3. Notes on convergence and stability. In this section we address the concern of numerical instability due to the fact that Fredholm integral equations of first kind posed on Banach spaces, such as the ones appearing in these equilibrium measure problems, are Hilbert-Schmidt and compact and thus cannot simply be inverted. In [27] we used a Tikhonov regularization [75,74] approach to overcome this problem and this also works for the d-dimensional generalization in this paper. We will only briefly sketch the idea behind Tikhonov regularization and leave the rest to the dedicated literature, see e.g. [39,38]. The idea behind this approach is to instead solve a well-posed second kind Fredholm integral equation which is in a well-defined sense 'adjacent' to the actual problem we intend to find solutions for. Denoting the full operator we wish to invert simply by F, the most straightforward Tikhonov regularization is to instead solve the problem where s is small. As mentioned in [27], we note that the error of the n-th order s-regularized approximation ρ s,n incurred from such a modification can be estimated via |ρ s,n´ρ | ď |ρ s,n´ρs |`|ρ s´ρ |, splitting into an error originated only from the Thikonov projection and an error caused by the spectral expansion of the regularized solution. We discuss the convergence properties of our implementation in section 4.3, where we also showcase the substantially improved stability gained from this regularization.

Numerical validation, convergence and experiments.
4.1. Uniqueness of the obtained equilibrium measures. Before comparing our method with analytic results and other numerical approaches, we explore whether our method produces a unique result for each set of parameters. It is currently not analytically known whether power law equilibrium measures are unique for general parameter ranges, so this investigation also provides numerical evidence that this is indeed the case. The most straightforward way to explore this question is to plot the energy of the measures obtained as a function of the radius input, which we do in Figure 4 for two generic combinations of dimensions and parameters. Figure 4 also shows the corresponding computed equilibrium measures. The minimal energy positive measure is unique in each tested case, including all of the other problems considered in the later sections. The consistently observed lower energy states for higher radius values than the local minimum are not positive measures and thus not admissible for the problem. Similar structures were observed in [27] for the one-dimensional ultraspherical method. In practice the problem of finding the minimal positive measure is nice enough for a simple constrained optimization to succeed when started on a lower radius than the solution. Automatically probing the energy on a grid before setting up the constraints of the optimization problem can lead to significant performance improvements.

4.2.
Comparison with analytic solutions in special cases. In [15], the authors described a number of special case solutions in arbitrary dimension to the power law equilibrium measure problem where one of the powers is an even integer. In this section we use these solutions to verify the accuracy of our method and to provide a practical demonstration of its convergence. We first consider the power law equilibrium measure problem defined by α " 2. The solution to this problem was found in [15] to be given by In Figure 6 we plot the numerically computed energy with varying R, showing that the analytic solutions are exactly the local minima obtained via our method. Additionally, in Figure 5 we plot the computed equilibrium measure for a two dimensional example on the disk and include an absolute error heatmap. Now for α " 4, the analytic solution as found in [15] may be written as with the constants ff , . As before, we plot the numerically computed energy with varying R along with the analytic radius in Figure 7. The reason these special cases have known analytic solutions is that they are particularly well-behaved -in fact, using our method both of these solutions may be computed to high precision using arbitrary precision floating point calculations with only one or two polynomial orders of approximation, meaning that these solutions can be computed almost instantaneously. Note how the energy as a function of R in Figure 6 and 7 show clear local minima at the analytic radius, further highlighting how these are special well-behaved cases compared to the general case. A glance at the form of the solutions tells us most of the work is already done by making the (a) pα, β, d, M q " p4, 0.5, 2, 1q (b) pα, β, d, M q " p4,´1.1, 3, 1q (c) pα, β, d, M q " p4,´3.9, 6, 2q   Figure 8 for completeness for an example where α " 4. The problem of visualizing convergence meaningfully is thus left to the next section. The fact that an analytic radius is known for these special cases does however give us a tool to visualize the error incurred by the precision chosen for the optimization. In Figure 9 we thus plot the maximum absolute error incurred for the obtained measure when deviating from the true radius solutions. As was observed in [27] for the one-dimensional case, errors in the computation of the radius propagate linearly to the error in the measure. This can be used as a guiding principle for the convergence conditions of the optimization method to obtain a result with a particular desired accuracy. 4.3. Numerical convergence properties and regularization. As seen in the previous section, convergence tests with known analytic solutions result in almost immediate convergence, making it difficult to sensibly visualize the properties of our method. We thus supplement the previous discussion with two example problems with no known analytic solution chosen to not fall into the well-behaved special cases. We furthermore explore the effect of Tikhonov regularization on the stability of the method. A natural way to think about convergence for solutions computed in terms of polynomial expansions is to measure whether and how fast the coefficients of the computed solution decay as the order of approximation increases. For particularly well-behaved functions the decay of such coefficients is exponential. In cases where only a single power law integral with an external potential are present we obtain exponential convergence in n as the right hand side hypergeometric functions are just polynomials when choosing the appropriate basis. Likewise, in the special attractive-repulsive case in which one of the powers is an even integer, the right-hand side of both power law integrals can simultaneously be brought into polynomial form, again leading to exponential convergence.
The convergence behaviour of the coefficients computed for ρ may be changed by choosing a different basis under the integrals, in particular via Remark 2.17. The natural choices are those which lead to one of the two operators being banded such that the bandwidth is minimized. In Figure 10 The regularized method as well as the early behaviour of the non-regularized method shows linear convergence in the coefficients for the generic attractive-repulsive. For large n, the instability of the coefficients in the non-regularized method is significant enough to be apparent in the resulting measures themselves. As seen in the example in Figure 10   Finally, it should be noted that our method via Theorem 2.16 and 2.18 does not inherently require the basis on the right hand side to be a particular set of Jacobi polynomials. As long as the operator mapping is accurately and consistently taken care of using the above results, the measure ρ will be obtained in a weighted Jacobi polynomial basis with any choice of right-hand side basis functions. As the basis in which the measure ρ is obtained does not change when doing this, only the computational efficiency in generating the operators may be affected by changing the right-hand side basis, not the n-convergence rate of the coefficients of ρ itself.
4.4. Comparison with particle swarm simulations. As mentioned in the introduction, the current go-to methods for numerically approximating equilibrium measure solutions are particle swarm simulations, see e.g. [14] where this approach and the observed phenomenon of gap formation for certain parameter ranges are explored. In the interesting parameter ranges, these simulations require some thousands to ten thousands of particles in one-dimension to meaningfully converge. Already in two dimensions, particle swarm simulations are close to the limit of what is computationally feasible. In contrast, our method converges independent of the dimension, which allows investigating problems far out of reach of the conventional methods. In this section, we compare the results of our method with a standard version of these particle swarm simulations in two dimensions. We also present some comparisons with lower density particle simulations in three dimensions, which have to be interpreted with care. In all of these cases, the results of our method are consistent with what is obtained via particle simulations. In Figure 11 we compare the particle solution to our spectral method for an example intentionally chosen to be a special case in which an analytic solution is also known via (4.1-4.2), allowing a three-way comparison via Figure 5. We present a further generic case one-dimensional and two dimensional example in Figure 12, without known analytic solutions.
(a) pα, β, dq " p2,´0.44, 2q (b) 2D histogram based on particles in (a) Fig. 11: (a) shows the results of a particle simulation for the indicated equilibrium measure problem with 1000 initially randomly distributed particles. (b) shows a twodimensional histogram approximating the density based on the particle simulation.
Note that the problem parameters here are intentionally the same as in Figure 5 and that an analytic solution is known in this special case.
(a) pα, β, dq " p1.3, 1.1, 2q (b) 2D histogram based on (a) (c) computed measure Fig. 12: (a) shows the results of a particle simulation for the indicated equilibrium measure problem with 1000 initially randomly distributed particles. (b) shows a twodimensional histogram approximating the density based on the particle simulation. (c) shows the computed equilibrium measure using the method we introduced in this paper.

4.5.
Tracing the gap formation boundary. It has been observed before [4,15,27] that for certain values of α and β, the assumption of the support of the equilibrium measure being a ball breaks down. For example, when α " 4, and β ą 2`2d´d 2 d`1 the analytic solution when assuming ball shaped support shows negative values at the origin [15], making the solution inadmissible. This gap formation phenomenon has been confirmed more generally by particle simulations in low dimensions [4] and recently via an ultraspherical spectral method in one-dimension in [27]. In Figure 13 we show an example of this gap formation phenomenon in a two dimensional particle simulation. In [27] the authors developed a two interval approach method, which was able to obtain admissible solutions beyond the gap formation boundary. Likewise, it is expected that in higher dimensions the correct support assumption for the measure is an annulus in two dimensions and hyperspherical shells in higher dimensions. While in one-dimension this problem can be solved with a two interval approach, this  [4,15,27].
is significantly more complicated in higher dimensions as hyperspherical shells are not simply comprised of two d-dimensional balls.
As is, the method proposed in this paper can thus not solve equilibrium measure problems past the gap formation boundary, although we intend to work on an extension of this approach for that case in the future. The present method can, however, help with understanding the shape of said gap formation boundary similarly to what was done in [27] in one-dimension, as the obtained measures will begin to show negative values at the origin. While the boundary point is known for a very limited number of special cases, such as α " 4 mentioned above, the general form of this boundary is presently not understood and is very difficult to impossible to explore with particle simulations even in just one or two dimensions due to slow convergence in the regions of low density around the origin immediately preceeding gap formation. In Figure 14 we show an exploration of the boundary using our method which can be performed in any dimension.
5. Discussion. In this paper, we have introduced a new approach to numerically solving power law equilibrium measures in arbitrary dimension for the case of ballshaped support. Our method is based on both existing as well as to our knowledge new recurrence results for radial Jacobi polynomials and Gaussian hypergeometric functions under the action of Riesz potentials. Aspects of these results, such as Theorem 2.16 may eventually be useful in analytic or computer-based proof aproaches to power law equilibrium measure problems, cf. [15]. In constrast to the conventional approach -particle simulations -our method's complexity is independent of the dimension d, making convergent arbitrary precision arithmetic computations possible in any dimension and for all parameter ranges where the assumption of ball-shaped support holds. Our method also allows the exploration of where this assumption breaks down, by tracing the gap formation boundary for parameters for which negative values around the origin appear. Future research will generalize this approach to solve equilibrium measure problems past the gap formation boundary on annuli and hyperspherical shells. 6. Appendices.
6.1. A -Extension of a Dyda-Kuznetsov-Kwaśnicki formula for fractional Laplacians. As mentioned in section 1.4 Dyda, Kuznetsov and Kwaśnicki proved that certain weighted Jacobi polynomials form a complete basis for the eigenfunctions of the fractional Laplacian on the d-dimensional unit ball [21]. We stated their result in (1.18). The domain of validity they prove for their formula is γ ą 0 but as we will show in this section, the formula also holds in the range γ P p´2, 0q. To prove this, we follow effectively the same proof Dyda, Kuznetsov and Kwaśnicki used, but at a critical point use a different theorem of theirs. The theorem we prove in this section is thus effectively an extension of Theorem 3 in [21]. As in [21], in the following V l,m pxq denotes a linear basis of the finite dimensional space with dimension M d,l " d`2l´2 d`l´2ˆ2`l´2 lṡ panned by solid harmonic polynomials with degree l ě 0. See section 1.3 as well as [21,20] for more details.  for all x such that |x| ă 1 and where δ :" d`2l.
Remark 6.2. Just as noted for [21,Theorem 3] for the fractional Laplacian, this result implies that the stated Jacobi polynomials form a complete orthogonal system of eigenfunctions for the weighted Riesz operator in the stated ranges.
6.2. B -Derived properties of radially shifted Jacobi polynomials. For ease of reference, this section concisely lists some of the basic properties for the radial Jacobi polynomials P pa,bq n p2|x| 2´1 q, where |x| P p0, 1q. They follow directly from the respective properties of the ordinary Jacobi polynomials, cf. [41, 18.9].
Classical recurrence relationship.