1 Introduction

There is a vast and multi-disciplinary literature on nonlocal energies, as they are at the crossover of different mathematical fields, and of different applications. Nonlocal energies arise as the macroscopic limit of long-range discrete interactions in the many-particle limit. The expression of the interaction kernel and its properties depend on the particle system of interest: it can model attraction, repulsion or a combination of both; it can be bounded or singular; it can be radial or anisotropic.

The mathematical literature on nonlocal energies has been mainly focused on the case of radial potentials, which model interactions depending on the mutual distance between the particles only (see, e.g., [1, 3, 4, 6,7,8,9, 14, 29]). In many applications, however, radial potentials are not realistic, and interactions may depend not only on the inter-particle distance, but also on the angle between their position vector and a given, preferred direction (see, e.g., [2, 5, 21, 28]). This is for instance the case for many biological systems, e.g., crowds, flocks of birds, schools of fish.

In materials science, some defects in metals, like dislocations of edge type, interact via an anisotropic potential, and this is the particle system that motivates our work. The anisotropy of the interactions is due to the motion of each dislocation being restricted to a given direction (the Burgers’ vector, \({\mathbf {b}}\), of the dislocation), which is reminiscent of the metal’s microscopic lattice structure. Under the simplifying assumption that dislocations are all parallel to each other, and have \({\mathbf {b}}=\mathbf {e_1}\), they can be modelled as point defects in two dimensions, and their interaction potential is

$$\begin{aligned} W_{\text {edge}}(x) = -\log |x| + \frac{x_1^2}{|x|^2}, \quad x\ne 0, \, x=(x_1, x_2)\in {\mathbb {R}}^2, \end{aligned}$$

see, e.g., [18]. The anisotropy of the interactions results into anisotropic low-energy dislocation structures (LEDS), dislocation walls in particular. The minimality of vertical one-dimensional structures (walls) was a long-standing conjecture in the engineering literature, and was recently proved in [25] for the nonlocal dislocation energy

$$\begin{aligned} I_{\text {edge}}(\mu ) = \iint _{{\mathbb {R}}^2\times {\mathbb {R}}^2} W_{\text {edge}}(x-y) \,d\mu (x) \,d\mu (y) + \int _{{\mathbb {R}}^2} |x|^2 \,d\mu (x) \end{aligned}$$
(1.1)

defined on probability measures \(\mu \in {\mathcal {P}}({\mathbb {R}}^2)\) representing the density of defects. For the derivation of an interaction energy related to (1.1) (but in a bounded domain) from a semi-discrete strain energy we refer to [24]. We also mention the recent work [22] where, starting from the nonlinear version of the strain-energy model considered in [24], ‘low-angle grain boundaries’ (like vertical walls) are shown to have the optimal energy scaling in accordance with the celebrated Read-Shockley formula.

In this paper we study an n-dimensional generalisation of (1.1). More precisely, we consider the family of nonlocal energies

$$\begin{aligned} I_{\alpha }(\mu ) = \iint _{{\mathbb {R}}^n\times {\mathbb {R}}^n} W_{\alpha }(x-y) \,d\mu (x) \,d\mu (y) + \int _{{\mathbb {R}}^n} |x|^2 \,d\mu (x) \end{aligned}$$
(1.2)

defined on \(\mu \in {\mathcal {P}}({\mathbb {R}}^n)\), where the interaction potential \(W_{\alpha }\) is given by

$$\begin{aligned} W_{\alpha }(x) = W_0(x) + \alpha \frac{x_1^2}{|x|^n}, \qquad W_0(x)= {\left\{ \begin{array}{ll} \displaystyle -\log |x| \quad &{} \text {if } n=2,\\ \displaystyle \frac{1}{|x|^{n-2}} &{}\text {if } n \ge 3, \end{array}\right. } \end{aligned}$$
(1.3)

for \(x\ne 0\), \(W_\alpha (0)=+\infty \), with \(x=(x_1,\dots , x_n)\in {\mathbb {R}}^n\), and \(\alpha > -1\). Note that \(W_{1}=W_{\text {edge}}\) for \(n=2\). The two-dimensional case was considered in [10, 25] for every \(\alpha \in {\mathbb {R}}\).

The main result of this paper is Theorem 3.1, where we show that, if \(n\ge 3\) and \(\alpha \in (-1,n-2]\), then the functional \(I_\alpha \) has a unique minimiser \(\mu _\alpha \) in \({\mathcal {P}}({\mathbb {R}}^n)\) and \(\mu _\alpha \) is of the form

$$\begin{aligned} \mu _\alpha := \frac{1}{| \Omega _\alpha |} \, \chi _{ \Omega _\alpha }, \quad \Omega _\alpha = \left\{ x=(x_1,\dots ,x_n)\in {\mathbb {R}}^n: \ \frac{x_1^2}{a(\alpha )^2} + \frac{1}{b(\alpha )^2}\sum _{i=2}^nx_i^2 < 1\right\} , \end{aligned}$$

for some \(a(\alpha ), b(\alpha )>0\). In other words, the minimiser of \(I_\alpha \) is the (normalised) characteristic function of a n-dimensional spheroid. Moreover, the spheroid is prolate for \(\alpha \in (-1,0)\) (that is, \(a(\alpha )>b(\alpha )\)) and oblate for \(\alpha \in (0,n-2]\) (that is, \(a(\alpha )<b(\alpha )\)).

Note that unlike in the case \(\alpha =0\), where the radial symmetry of the kernel and of the confinement strongly suggests that the ball is the natural candidate for the minimisation (see [11, 16, 26]), the situation for \(\alpha \ne 0\) is less clear. Indeed, while heuristically one can expect the anisotropy to cause an elongation (for \(\alpha <0\)) or a contraction (for \(\alpha >0\)) of the ball in the \(x_1\)-direction, it is surprising that the minimiser is still a characteristic function, and that its support is a spheroid.

1.1 Our approach and discussion

In most of the cases treated in the mathematical literature on nonlocal systems, the interaction kernel is assumed to be radial, and one of the goals is to show that the corresponding minimiser is radially symmetric (or to show that the minimiser is unique, which trivially implies its radial symmetry), see e.g. [1, 4, 6,7,8,9, 14, 29]. Radial symmetry is paramount in the identification of the minimiser in the classical case of purely Coulomb interactions, corresponding in our setting to \(\alpha =0\) (see [16, 26] for \(n=2\), and [11] for \(n\ge 3\)). Explicitly characterising the minimiser, or even understanding its shape and general properties, is therefore much more challenging in the case of anisotropic interactions. To the best of our knowledge, this has been previously done only in [10, 25], in dimension \(n=2\).

Our result generalises to any dimension \(n\ge 3\) the work [10] and is another paradigmatic example of the role of the anisotropy of the kernel on the shape of minimisers. Our approach here is however completely different from the approach in [10], which was based on complex-analysis techniques, clearly no longer available in higher dimensions.

A key step in [10] was to compute exactly the (gradient of the) potential \(W_\alpha *\mu _{a,b}\), with \(\mu _{a,b}\) being the (normalised) characteristic function of an ellipse of semi-axes a and b, and to impose the Euler–Lagrange conditions associated to \(I_\alpha \). Also for \(n\ge 3\) we prove the minimality of spheroids via the Euler–Lagrange conditions [see (3.1), (3.2)], in the range \(\alpha \in (-1,n-2]\) for which \(I_\alpha \) is both well-defined and strictly convex (hence, the Euler–Lagrange conditions are necessary and sufficient for minimality). To do so, we need to compute the potential \(W_\alpha *\mu _{a,b}\), with \(\mu _{a,b}\) being the (normalised) characteristic function of a spheroid \(\Omega (a,b)\) of semi-axis a in the \(x_1\)-direction and b in all the other directions. The computation of \(W_\alpha *\mu _{a,b}\) for \(n\ge 3\) is substantially different from the two-dimensional case in [10]. For \(n=2\) it was crucial to rewrite the potential in complex variables and recognise that \(\nabla W_0*\mu _{a,b}\) is the Cauchy transform of the ellipse, which had been computed for instance in [19] for rotating vortex patches in fluid dynamics. Then the gradient of the anisotropic part of the potential in complex coordinates was computed by noting that it could be written as a suitable complex derivative of the fundamental solution of the operator \(\partial ^2\), where \(\partial = \partial /\partial z\).

Such complex-analysis techniques are clearly not available in the higher-dimensional case, which we tackle here by means of the following strategy. We write \(\Phi _\alpha :=W_\alpha *\mu _{a,b}\), and \(\alpha \Psi := \Phi _\alpha -\Phi _0\). First of all, the expression of the Coulomb potential \(\Phi _0\) of a spheroid in \({\mathbb {R}}^n\) is well-known (see, e.g., [13, 20]). The challenge is to express the anisotropic potential \(\Psi \) in terms of the known potential \(\Phi _0\). We do it differently in \(\Omega (a,b)\) and outside \(\Omega (a,b)\). In \(\Omega (a,b)\) we show that \(\Psi \) can be obtained by differentiating \(\Phi _0\) with respect to the aspect ratio \(a^2/b^2\) of the spheroid. In \({\mathbb {R}}^n{\setminus } \Omega (a,b)\), instead, \(\Psi \) is obtained by differentiating \(\Phi _0\) with respect to a parameter spanning a family of spheroids confocal with \(\Omega (a,b)\), using the fact that the expression of \(\Phi _0\) is invariant on confocal spheroids (see (3.33)).

With the expression of \(\Phi _\alpha \) at hand, we then impose the Euler–Lagrange conditions (3.1)–(3.2). We find that the first condition is satisfied, for \(\alpha \in (-1,n-2]\), by at least a pair \((a(\alpha ),b(\alpha ))\) of semi-axes, with \(a(\alpha )>b(\alpha )>0\) for \(\alpha \in (-1,0)\) and \(0<a(\alpha )<b(\alpha )\) for \(\alpha \in (0,n-2]\). Hence there is at least one stationary, non-degenerate spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for the energy \(I_\alpha \). We then show that, for any spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for which the stationarity condition (3.1) is satisfied, also the unilateral condition (3.2) is satisfied. Since (3.1)–(3.2) are necessary and sufficient conditions for minimality, this implies that any spheroid \(\Omega (a(\alpha ),b(\alpha ))\) satisfying (3.1) is in fact a minimiser for \(I_\alpha \). The strict convexity of the energy then gives uniqueness of the minimiser, and in particular implies that there is only one spheroid satisfying (3.1). This approach can be carried out also in the two-dimensional case (see [27]). We note that in two dimensions there are several methods that work, see for instance [23] where a maximum-principle argument is used.

1.1.1 Dimensionality of minimisers for \(n=2\) and \(n\ge 3\)

For the energy (1.1), it was shown in [25] that the unique minimiser is one dimensional, and is given by the semi-circle law on the vertical axis,

where denotes the restriction of the one-dimensional Hausdorff measure to the interval \((-\sqrt{2},\sqrt{2}).\) The semi-circle law also arises as the unique minimiser of the one-dimensional logarithmic energy with quadratic confinement (see [31]), and represents in that case the optimal positions of the eigenvalues of a Hermitian random matrix with Gaussian entries.

We recall that the minimiser of the Coulomb-gas energy \(I_0\) for \(n=2\) is the two-dimensional measure \(\mu _0 = \frac{1}{\pi }\chi _{B_1(0)}\), the so-called circle-law, also well-known in the context of random matrices. In fact, the minimiser of \(I_0\) is the normalised characteristic function of a ball in any dimension. The change of dimension of the minimiser of the energy \(I_\alpha \), for \(n=2\), between \(\alpha =0\) and \(\alpha =1\) was investigated by the authors in [10]. In [10] it was shown in particular that the minimiser of \(I_\alpha \), for \(n=2\) and \(\alpha \in (-1,1)\), is the two-dimensional measure

$$\begin{aligned} \mu _\alpha = \frac{1}{\pi }\frac{1}{\sqrt{1-\alpha ^2}} \chi _{\Omega (\sqrt{1-\alpha },\sqrt{1+\alpha })}, \end{aligned}$$

with

$$\begin{aligned} \Omega (\sqrt{1-\alpha },\sqrt{1+\alpha })=\left\{ x=(x_1,x_2)\in {\mathbb {R}}^2:\ \frac{x_1^2}{1-\alpha }+\frac{x_2^2}{1+\alpha }<1\right\} . \end{aligned}$$

Hence the minimiser of \(I_\alpha \) has full dimension for \(\alpha \in (-1,1)\), and is one-dimensional for both \(\alpha \le -1\) and \(\alpha \ge 1\), being respectively the semi-circle law on the horizontal or the vertical axis.

A question left open in [10] was to understand why there is a change of dimension of the minimiser \(\mu _\alpha \) at \(\alpha =\pm 1\), for \(n=2\). The relation between the dimensionality of the minimiser of a nonlocal energy and the singularity of the interaction kernel is a fascinating and subtle problem. The available results in the literature are usually of the form of a lower bound for the dimension of the measure (see, e.g., [1]), which is helpful if the goal is to prove that the dimension is full, but less so to prove that there is a loss of dimension.

For the energy \(I_\alpha \) in dimension \(n=2\), since the Fourier transform of \(W_\alpha \) changes sign exactly at the values \(\alpha =\pm 1\), it was natural to conjecture that the change of dimension could be due to the change of sign of \(\widehat{W}_\alpha \). Similarly, since the Laplacian of \(W_\alpha \) is

$$\begin{aligned} \Delta W_\alpha (x) = -(1-\alpha )\partial _{x_1}^2 \log |x| - (1+\alpha )\partial _{x_2}^2 \log |x|, \end{aligned}$$

it was reasonable to expect that the singular behaviour exhibited by \(\mu _\alpha \) at \(\alpha =\pm 1\) was a consequence of the ‘degeneracy’ of \(\Delta W_\alpha \) at those values.

The analysis done in this paper demonstrates that the situation is more delicate. While it is still plausible to expect that a positive Fourier transform (or a non-degenerate Laplacian) results into a fully-dimensional minimiser, the contrary is not true, at least for \(n\ge 3\). Indeed while for \(n\ge 3\) the Fourier transform \(\widehat{W}_\alpha \) of the interaction kernel changes sign at \(\alpha =n-2\) (see (2.5)), and similarly

$$\begin{aligned} \Delta W_\alpha (x) =\bigg (1-\frac{\alpha }{n-2}\bigg ) \partial _{x_1}^2\bigg (\frac{1}{|x|^{n-2}}\bigg ) + \bigg (1+\frac{\alpha }{n-2}\bigg ) \sum _{i=2}^n\partial _{x_i}^2\bigg (\frac{1}{|x|^{n-2}}\bigg ), \end{aligned}$$

we prove that the minimiser of \(I_{\alpha }\) is the characteristic function of a non-degenerate spheroid also for the limit value \(\alpha =n-2\).

1.1.2 The shape of minimisers for \(\alpha >0\) and \(\alpha <0\)

In the two-dimensional case \(n=2\), changing sign to \(\alpha \) corresponds to swapping \(x_1\) and \(x_2\) (up to a constant in the energy), due to the zero-homogeneity of the energy. Hence it is sufficient to characterise the minimisers of \(I_\alpha \) for \(\alpha >0\), which is what we did in [10].

This is no longer true for \(n\ge 3\), since in this case there is only one privileged coordinate. Intuitively, configurations elongated on \(x_1\) are penalised for \(\alpha >0\), and preferred for \(\alpha <0\) (see, e.g., (2.1)), hence the minimisers for \(\alpha >0\) and \(\alpha <0\) cannot be congruent up to a rotation, which is the case in dimension two. More precisely, for \(n\ge 3\) we may write

$$\begin{aligned} W_{\alpha }(x) = (1+\alpha ) W_0(x) -\alpha \frac{1}{|x|^n} \sum _{i=2}^n x_i^2. \end{aligned}$$

Thus, changing sign to \(\alpha \) corresponds not only to a change in the anisotropy, but also to a rescaling of the Coulomb kernel. This also suggests that a different behaviour of the energy should be expected at \(\alpha =-1\) in dimension \(n=2\) and \(n\ge 3\), as discussed next.

1.1.3 The limit case \(\alpha =-1\)

In dimension \(n\ge 3\), for \(\alpha =-1\) the anisotropy cancels completely the \(x_1\)-component of the Coulomb potential, since from (1.3)

$$\begin{aligned} W_{-1}(x) = \frac{1}{|x|^{n-2}} - \frac{x_1^2}{|x|^n} = \frac{x_2^2+\ldots +x_n^2}{|x|^n}, \quad x\ne 0. \end{aligned}$$

As a consequence, there is a discrepancy with the situation for \(\alpha \in (-1,n-2]\): While \(I_\alpha \) is lower semicontinuous for any \(\alpha \in (-1,n-2]\), the functional \(I_{-1}\) with the kernel \(W_{-1}\) above (and \(W_{-1}(0)=+\infty \)) is not lower semicontinuous (see Remark 2.2). In particular, \(I_{-1}\) does not describe the asymptotic behaviour of the functionals \(I_\alpha \), as \(\alpha \rightarrow -1^+\) (see Remark 4.2). We resolve this issue in Sect. 4, where we characterise the \(\Gamma \)-limit \(J_*\) of \(I_\alpha \), as \(\alpha \rightarrow -1^+\), in Fourier space, on probability measures with compact support. This partial representation allows us to show strict convexity of \(J_*\) on a class of measures that is the relevant one for minimisation, and hence to deduce uniqueness of the minimiser for \(J_*\). Moreover, we show that the minimiser is, also in this limit case, the (normalised) characteristic function of an n-dimensional (prolate) spheroid.

1.2 Open questions and future work

There are several questions that we will address in future work. We believe that ellipses, or spheroids, arise as minimisers of more general anisotropic energies. A first step would be to consider interaction kernels of the form \(W_0 + \alpha W_{\text {aniso}}\), with

$$\begin{aligned} W_{\text {aniso}}(x) = \frac{x_1^{2\kappa }}{|x|^{n-2+2\kappa }}, \quad \kappa \in {\mathbb {N}}. \end{aligned}$$

It is plausible to expect that in the two-dimensional case \(n=2\) and in a suitable range of \(\alpha \) the minimisers are ellipses. It would be interesting to understand whether, as for \(\kappa =1\), they shrink to a segment for some special value \(\alpha ^*\) within, or at the boundary of, the interval of strict convexity of the corresponding energy (or of positivity of the Fourier transform of \(W_0 + \alpha W_{\text {aniso}}\)). This analysis would help to shed more light on what causes the loss of dimension of the minimising measure.

1.3 Outline

The plan of the paper is as follows. In Sect. 2 we show that the energy \(I_\alpha \) admits a unique global minimiser \(\mu _\alpha \), for \(\alpha \in (-1,n-2]\). Section 3 is devoted to the main result of the paper, Theorem 3.1, in which \(\mu _\alpha \) is identified as the (normalised) characteristic function of a spheroid, for \(\alpha \in (-1,n-2]\). Finally, in Sect. 4, we discuss the limit case \(\alpha \rightarrow -1^+\).

2 Existence and uniqueness of the minimiser of \(I_{\alpha }\) for \(\alpha \in (-1,n-2]\), \(n\ge 3\)

In this section we prove that for every \(\alpha \in (-1,n-2]\) the nonlocal energy \(I_{\alpha }\) defined in (1.2), for \(n\ge 3\), has a unique minimiser \(\mu _\alpha \in {\mathcal {P}}({\mathbb {R}}^n)\), and that the minimiser has a compact support.

Proposition 2.1

Let \(n\ge 3\), and let \(\alpha \in (-1,n-2]\). Then the energy \(I_{\alpha }\) is well defined on \({\mathcal {P}}({\mathbb {R}}^n)\), is strictly convex on the class of measures with compact support and finite interaction energy, and has a unique minimiser in \({\mathcal {P}}({\mathbb {R}}^n)\). Moreover, the minimiser has compact support and finite energy.

Proof

The case \(n =2\) has been proved in [25, Section 2] and [10, Proposition 2.1]. For \(n\ge 3\) the proof follows by a similar argument. The main novelty is Step 3.2, where we present a new, more transparent way of proving condition (2.3) below.

Step 1: Well definiteness of \(I_\alpha \). Since \(\alpha > -1\), if we write \(W_\alpha \), for \(x\ne 0\), as

$$\begin{aligned} W_{\alpha }(x) = \frac{1}{|x|^n}\bigg ((1+\alpha ) x_1^2 + \sum _{i=2}^n x_i^2\bigg ), \end{aligned}$$
(2.1)

we can immediately see that the energy is well-defined and non-negative on \({\mathcal {P}}({\mathbb {R}}^n)\).

Step 2: Existence of a compactly supported minimiser. First of all, it is easy to see that \(I_\alpha (\mu _B)<+\infty \), where \(\mu _B= \frac{1}{|B_1(0)|} \chi _{B_1(0)}\). This implies that \(\inf _{{\mathcal {P}}({\mathbb {R}}^n)} I_\alpha < +\infty \). Moreover, we have that

$$\begin{aligned} W_{\alpha }(x-y) + \frac{1}{2}(|x|^2+|y|^2) \ge \frac{1}{2}(|x|^2+|y|^2). \end{aligned}$$
(2.2)

This lower bound provides tightness and hence compactness with respect to narrow convergence for minimising sequences, that, together with the lower semicontinuity of \(I_\alpha \), guarantees the existence of a minimiser. As in [25, Section 2.2] (see also [3]), one can show that any minimiser of \(I_{\alpha }\) has compact support, again by (2.2).

Step 3: Strict convexity of \(I_\alpha \) and uniqueness of the minimiser. We prove that

$$\begin{aligned} \int _{{\mathbb {R}}^n} W_{\alpha }*(\nu _1-\nu _2) \,d(\nu _1-\nu _2) > 0 \end{aligned}$$
(2.3)

for every \(\nu _1, \nu _2 \in {\mathcal {P}}({\mathbb {R}}^n)\), \(\nu _1\ne \nu _2\), with compact support and finite interaction energy, namely such that \(\int _{{\mathbb {R}}^n} (W_{\alpha }*\nu _i) \, d\nu _i < + \infty \) for \(i=1,2\). Condition (2.3) implies strict convexity of \(I_{\alpha }\) on the set of probability measures with compact support and finite interaction energy and, consequently, uniqueness of the minimiser.

To prove (2.3), we follow again the same strategy as in [25, Section 2.3]. The idea consists in rewriting the interaction energy of \(\nu :=\nu _1-\nu _2\) in Fourier space, as

$$\begin{aligned} \int _{{\mathbb {R}}^n} W_{\alpha }*\nu \,d \nu = \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |{\widehat{\nu }}(\xi )|^2\, d\xi , \end{aligned}$$
(2.4)

and proving that \(\widehat{W}_\alpha \) is a positive distribution.

Step 3.1: Computation of \(\widehat{W}_\alpha \). Note that \(W_{\alpha }\) is a tempered distribution, namely \(W_{\alpha }\in {{\mathcal {S}}}'\), where \({\mathcal {S}}\) denotes the Schwartz space; hence also \({\widehat{W}}_{\alpha }\in {{\mathcal {S}}}'\). We recall that \({\widehat{W}}_{\alpha }\) is defined by the formula

$$\begin{aligned} \langle {\widehat{W}}_{\alpha }, \varphi \rangle := \langle W_{\alpha }, {\widehat{\varphi }}\,\rangle \qquad \text { for every }\, \varphi \in {{\mathcal {S}}} \end{aligned}$$

where, for \(\xi \in {\mathbb {R}}^n\),

$$\begin{aligned} {\widehat{\varphi }}(\xi ):=\int _{{\mathbb {R}}^n}\varphi (x)e^{-2\pi i\xi \cdot x}\, dx. \end{aligned}$$

To compute \(\widehat{W}_\alpha \) it is convenient to rewrite \(W_\alpha \) as

$$\begin{aligned} W_\alpha (x) = \Big (1+\frac{\alpha }{n}\Big ) \frac{1}{|x|^{n-2}} + \frac{\alpha }{n}\,\frac{1}{|x|^n} \bigg ((n-1) x_1^2 -\sum _{i=2}^n x_i^2\bigg ), \end{aligned}$$

namely as the sum of the n-dimensional Coulomb potential (up to a multiplicative constant) and the ratio between a homogeneous harmonic polynomial of degree two and a power of |x|. By [30, eq. (32), p. 73] and by [15, Exercise 1, p. 154] we have that the Fourier transform \({\widehat{W}}_{\alpha }\) of \(W_{\alpha }\) is given by

$$\begin{aligned} \langle {\widehat{W}}_{\alpha }, \varphi \rangle&= \Big (1+\frac{\alpha }{n}\Big )\frac{(n-2)\pi ^{\frac{n}{2}-2}}{2\Gamma (\frac{n}{2})} \int _{{\mathbb {R}}^n}\frac{1}{|\xi |^2}\varphi (\xi ) d\xi - \frac{\alpha }{n} \,\frac{\pi ^{\frac{n}{2}-2}}{\Gamma (\frac{n}{2})} \int _{{\mathbb {R}}^n} \frac{(n-1) \xi _1^2 - \sum _{i=2}^n\xi _i^2}{|\xi |^4} \, \varphi (\xi ) d\xi \nonumber \\&= \frac{\pi ^{\frac{n}{2}-2}}{2\Gamma (\frac{n}{2})} \int _{{\mathbb {R}}^n} \frac{(n-2-\alpha ) \xi _1^2 + (n-2+\alpha ) \sum _{i=2}^n \xi _i^2}{|\xi |^4} \, \varphi (\xi ) d\xi \end{aligned}$$
(2.5)

for every \(\varphi \in {{\mathcal {S}}}\). Thus, (2.5) implies that \({\widehat{W}}_{\alpha }\) is a positive function in \(L^1_{\text {loc}}({\mathbb {R}}^n)\) for every \(\alpha \in (-1, n-2]\), hence in particular a positive tempered distribution.

Step 3.2: Proof of (2.4). We start by proving that (2.4) holds when \(\nu \) is a non-negative finite Borel measure with compact support, where we understand that the two sides of the formula are either both finite and coincide, or both equal to \(+\infty \).

We proceed by regularisation. Let \(\varphi \in C^\infty _c(B_1(0))\) be non-negative, radial, and with \(\int _{{\mathbb {R}}^n} \varphi (x)\,dx =1.\) For \(\varepsilon >0\) we define

$$\begin{aligned} \varphi _\varepsilon (x):= \frac{1}{\varepsilon ^n} \varphi \left( \frac{x}{\varepsilon }\right) \qquad \text { and} \qquad \nu _\varepsilon : = \nu *\varphi _\varepsilon . \end{aligned}$$
(2.6)

We claim that

$$\begin{aligned} \int _{{\mathbb {R}}^n} (W_\alpha *\nu _\varepsilon )(x) \nu _\varepsilon (x) \,dx = \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |\widehat{\nu _\varepsilon }(\xi )|^2 \, d\xi . \end{aligned}$$
(2.7)

To show this, let us set \(f:=W_\alpha *\nu _\varepsilon \) and \(g:=\nu _\varepsilon \), and note that \(g \in C^{\infty }_c({\mathbb {R}}^n)\) and \(f \in C^{\infty }({\mathbb {R}}^n)\). Moreover, since \(\widehat{\nu _\varepsilon }\in {{\mathcal {S}}}\) and \(\widehat{W}_\alpha \in L^1_{\text {loc}}({\mathbb {R}}^n)\) behaves as \(1/|\xi |^2\) at infinity by (2.5), we have that \(\widehat{f}=\widehat{W}_\alpha \,\widehat{\nu _\varepsilon }\in L^1({\mathbb {R}}^n)\). Let \(\psi \in C^{\infty }_c({\mathbb {R}}^n)\) be such that \(\psi = 1\) on \(B_1(0)\) and let \(R>0\) be such that the support of g is contained in \(B_R(0)\). If \(\tau >0\) is such that \(\tau R < 1\), then, by Parseval formula,

$$\begin{aligned} \int _{{\mathbb {R}}^n} f(x)g(x)\,dx= & {} \int _{{\mathbb {R}}^n} \psi (\tau x)f(x)g(x) \,dx \nonumber \\= & {} \int _{{\mathbb {R}}^n} (\widehat{\psi (\tau \, \cdot )} *\widehat{f}\, )(\xi ) \,\overline{\widehat{g}(\xi )}\, d\xi = \int _{{\mathbb {R}}^n} ({\widehat{\psi }}_\tau *\widehat{f}\,)(\xi )\, \overline{\widehat{g}(\xi )}\, d\xi , \end{aligned}$$
(2.8)

where \({\widehat{\psi }}_\tau (x):= \tau ^{-n} {\widehat{\psi }}(x/\tau )\). Using that \({{\widehat{\psi }}}\in {\mathcal {S}}\subset L^1({\mathbb {R}}^n)\), \(\int _{{\mathbb {R}}^n}{{\widehat{\psi }}}(\xi )\,d\xi =\psi (0)=1\), and that \(\widehat{f}\in L^1({\mathbb {R}}^n)\), it is easy to see that \({\widehat{\psi }}_\tau *\widehat{f}\) converges to \({\widehat{f}}\) in \(L^1({\mathbb {R}}^n)\), as \(\tau \rightarrow 0\). In fact,

$$\begin{aligned} \Vert {\widehat{\psi }}_\tau *\widehat{f} -\widehat{f} \Vert _{L^1}\le \int _{{\mathbb {R}}^n}|{{\widehat{\psi }}}(\eta )|\int _{{\mathbb {R}}^n} |\widehat{f}(\xi -\tau \eta )-\widehat{f}(\xi )|\,d\xi d\eta \end{aligned}$$

and one can conclude by the continuity of translations in the \(L^1\)-norm and by the Dominated Convergence Theorem. Since \(\widehat{g}\in L^{\infty }({\mathbb {R}}^n)\), we deduce that

$$\begin{aligned} \lim _{\tau \rightarrow 0} \int _{{\mathbb {R}}^n} ({\widehat{\psi }}_\tau *\widehat{f}\,)(\xi )\, \overline{\widehat{g}(\xi )}\, d\xi = \int _{{\mathbb {R}}^n} \widehat{f}(\xi )\, \overline{\widehat{g}(\xi )}\, d\xi , \end{aligned}$$

which, together with (2.8), proves (2.7).

We now let \(\varepsilon \rightarrow 0\) in (2.7). For the right-hand side we observe that for every \(\xi \in {\mathbb {R}}^n\)

$$\begin{aligned} \widehat{\varphi _\varepsilon }(\xi ) = {\widehat{\varphi }}(\varepsilon \xi ) \rightarrow {\widehat{\varphi }}(0) =1, \end{aligned}$$

as \(\varepsilon \rightarrow 0\), and that \(\Vert \widehat{\varphi _\varepsilon }\Vert _{L^\infty }\le \Vert \varphi \Vert _{L^1} =1\) for every \(\varepsilon >0\). Therefore, either by the Dominated Convergence Theorem or the Fatou Lemma, we have

$$\begin{aligned} \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |\widehat{\nu _\varepsilon }(\xi )|^2 \, d\xi = \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |{\widehat{\nu }}(\xi )|^2 |\widehat{\varphi _\varepsilon }(\xi ) |^2 \, d\xi \ \rightarrow \ \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |{\widehat{\nu }}(\xi )|^2 \, d\xi , \end{aligned}$$
(2.9)

as \(\varepsilon \rightarrow 0\), even if the right-hand side is infinite.

To deal with the left-hand side of (2.7) we note that, since \(\alpha >-1\), there exists a positive constant \(C=C(\alpha )\) such that

$$\begin{aligned} \frac{1}{C}\, W_0(x) \le W_\alpha (x) \le C\, W_0(x) \qquad \text { for every }\, x\in {\mathbb {R}}^n. \end{aligned}$$
(2.10)

Hence,

$$\begin{aligned} (W_\alpha *\varphi _\varepsilon )(z) \le C(W_0*\varphi _\varepsilon )(z) \qquad \text { for every }\, z\in {\mathbb {R}}^n. \end{aligned}$$

Since \(W_0\) is superharmonic and \(\varphi \) is radial with integral 1, the mean value property on spheres yields

$$\begin{aligned} (W_0*\varphi _\varepsilon )(z) \le W_0(z) \qquad \text { for every }\, z\in {\mathbb {R}}^n. \end{aligned}$$

Thus, combining the three previous inequalities,

$$\begin{aligned} (W_\alpha *\varphi _\varepsilon )(z) \le C\, W_0(z) \le C^2\, W_\alpha (z) \qquad \text { for every }\, z\in {\mathbb {R}}^n. \end{aligned}$$
(2.11)

Note that for every \(z \in {\mathbb {R}}^n\)

$$\begin{aligned} (W_\alpha *\varphi _\varepsilon )(z) \rightarrow W_\alpha (z), \end{aligned}$$
(2.12)

as \(\varepsilon \rightarrow 0\), since \(W_\alpha \) is continuous as a function with values into \([0,+\infty ]\). Owing to the convergence (2.12) and the domination (2.11), we can either apply the Dominated Convergence Theorem or the Fatou Lemma to deduce that

$$\begin{aligned} \iint _{{\mathbb {R}}^n\times {\mathbb {R}}^n}(W_\alpha *\varphi _\varepsilon )(x-y) \, d\nu (x)\,d\nu (y) \ \rightarrow \ \iint _{{\mathbb {R}}^n\times {\mathbb {R}}^n} W_\alpha (x-y) \, d\nu (x)\, d\nu (y),\qquad \end{aligned}$$
(2.13)

as \(\varepsilon \rightarrow 0\), even if the right-hand side is infinite.

We now go back to the left-hand side of (2.7) and observe that

$$\begin{aligned} \int _{{\mathbb {R}}^n} (W_\alpha *\nu _\varepsilon )(x)\nu _\varepsilon (x) \, dx = \iint _{{\mathbb {R}}^n\times {\mathbb {R}}^n}(W_\alpha *\varphi _\varepsilon *\varphi _\varepsilon )(x-y) \, d\nu (x)\,d\nu (y). \end{aligned}$$

Note that \((\varphi _\varepsilon *\varphi _\varepsilon )(x) = \varepsilon ^{-n}(\varphi *\varphi )(x/\varepsilon )\) and that \(\varphi *\varphi \) inherits the properties of \(\varphi \): it is radial, belongs to \(C^\infty _c({\mathbb {R}}^n)\), and \(\int _{{\mathbb {R}}^n} (\varphi *\varphi )(x)\,dx =1\). Therefore, (2.13) holds with \(\varphi _\varepsilon \) replaced by \(\varphi _\varepsilon *\varphi _\varepsilon \). This concludes the proof of (2.4) for a non-negative measure \(\nu \).

We now prove (2.4) for a signed and neutral measure \(\nu :=\nu _1-\nu _2\), where \(\nu _1, \nu _2 \in {\mathcal {P}}({\mathbb {R}}^n)\) have compact support and finite interaction energy. First of all, by using (2.4) for \(\nu _1+\nu _2\) we have that

$$\begin{aligned} \int _{{\mathbb {R}}^n} (W_\alpha *(\nu _1+\nu _2))(x)\,d(\nu _1+\nu _2)(x)= & {} \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi ) |{\widehat{\nu }}_1(\xi )+{\widehat{\nu }}_2(\xi )|^2\,d\xi \\\le & {} 2 \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi )(|{\widehat{\nu }}_1(\xi )|^2 +|{\widehat{\nu }}_2(\xi )|^2 )\, d\xi < +\infty . \end{aligned}$$

By expanding both sides of the identity above and using (2.4) for \(\nu _1\) and \(\nu _2\) we get

$$\begin{aligned} \int _{{\mathbb {R}}^n} (W_\alpha *\nu _1)(x)\,d\nu _2(x)= \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi )\,\mathrm{Re}\big ({\widehat{\nu }}_1(\xi )\overline{{\widehat{\nu }}_2(\xi )}\,\big ) \,d\xi , \end{aligned}$$

which, by using again (2.4) for \(\nu _1\) and \(\nu _2\), gives

$$\begin{aligned} \int _{{\mathbb {R}}^n} (W_\alpha *(\nu _1-\nu _2))(x)\,d(\nu _1-\nu _2)(x) = \int _{{\mathbb {R}}^n} \widehat{W}_\alpha (\xi )|{\widehat{\nu }}_1(\xi )-{\widehat{\nu }}_2(\xi )|^2\,d\xi . \end{aligned}$$

Since the right-hand side is strictly positive for \(\nu _1\ne \nu _2\), Eq. (2.3) is proved. \(\square \)

Remark 2.2

In the case \(n\ge 3\), unlike in the two-dimensional case, the energy is not well-defined for \(\alpha <-1\). Indeed, writing \(W_\alpha \) as in (2.1), we can see that the two terms in the right-hand side of (2.1) are both unbounded for |x| close to zero, and have opposite sign.

Even for \(\alpha =-1\) the situation is subtle. The functional \(I_{-1}\) with kernel \(W_{-1}\) defined as in (1.3) is not lower semicontinuous with respect to narrow convergence: Indeed, the probability measures converge narrowly to the Dirac delta at 0, but

$$\begin{aligned} 0=\lim _{k\rightarrow +\infty } I_{-1}(\mu _k) < I_{-1}(\delta _0) = +\infty . \end{aligned}$$
(2.14)

So, in particular, \(I_{-1}\) is not the \(\Gamma \)-limit of \(I_\alpha \) for \(\alpha \rightarrow -1^{+}\). Moreover, (2.14) implies that the relaxed functional \(\overline{I_{-1}}\) of \(I_{-1}\) is equal to 0 at \(\delta _0\), which is therefore a minimiser of \(\overline{I_{-1}}\). Note that the kernel \(W_{-1}\) is also not lower semicontinuous, since \(W_{-1}(0)=+\infty \) and

$$\begin{aligned} \liminf _{x\rightarrow 0} W_{-1}(x)=0. \end{aligned}$$

If we define a new kernel \({\widetilde{W}}_{-1}\) to be as \(W_{-1}\) for \(x\ne 0\), and \({\widetilde{W}}_{-1}(0):=0\), then \({\widetilde{W}}_{-1}\) is lower semicontinuous, and the corresponding functional \(\widetilde{I}_{-1}\) has a unique minimiser, which is simply the Dirac delta at 0.

The case \(\alpha =-1\) will be discussed in detail in Sect. 4.

3 Minimality of spheroids

It is a standard computation in potential theory to show that any minimiser \(\mu \) of \(I_\alpha \) must satisfy the following Euler–Lagrange conditions: There exists \(C\in {\mathbb {R}}\) such that

$$\begin{aligned}&(W_{\alpha }*\mu )(x) + \frac{|x|^2}{2} = C \quad \text {for}\, \mu \text {-a.e. }\, x\in {{\,\mathrm{supp}\,}}\mu , \end{aligned}$$
(3.1)
$$\begin{aligned}&(W_{\alpha }*\mu )(x) + \frac{|x|^2}{2} \ge C \quad \text {for q.e.}\,x\in {\mathbb {R}}^n, \end{aligned}$$
(3.2)

where quasi everywhere (q.e.) means up to sets of zero capacity (see [26, Chapter I, Theorem 1.3] or [25, Section 3.1]). The Euler–Lagrange conditions (3.1) and (3.2) are in fact equivalent to minimality for \(\alpha \in (-1,n-2]\) due to Proposition 2.1. We refer to [25, Section 3.1] for details.

Let \(a,b>0\), and let \(\Omega (a,b)\subset {\mathbb {R}}^n\) denote the ellipsoid with semi-axis a in the \(x_1\) direction and the other semi-axes of the same length b, namely

$$\begin{aligned} \Omega (a,b):= \left\{ x=(x_1,\dots ,x_n)\in {\mathbb {R}}^n: \ \frac{x_1^2}{a^2} + \frac{1}{b^2}\sum _{i=2}^nx_i^2 < 1\right\} . \end{aligned}$$

This special ellipsoid is called oblate spheroid if \(a<b\) and prolate spheroid if \(a>b\).

The main result is the following.

Theorem 3.1

Let \(n\ge 3\) and \(\alpha \in (-1,n-2]\). There exist \(a(\alpha ), b(\alpha )>0\) such that the measure

$$\begin{aligned} \mu _\alpha := \frac{1}{| \Omega _\alpha |} \, \chi _{ \Omega _\alpha }, \qquad \Omega _\alpha := \Omega (a(\alpha ), b(\alpha )), \end{aligned}$$

is the unique minimiser of the functional \(I_\alpha \) in \(\mathcal P({\mathbb {R}}^n)\) and satisfies the Euler–Lagrange conditions

$$\begin{aligned} (W_{\alpha }*\mu _\alpha )(x) + \frac{|x|^2}{2}&= C_\alpha \quad \text {for every}\, x\in \Omega _\alpha , \end{aligned}$$
(3.3)
$$\begin{aligned} (W_{\alpha }*\mu _\alpha )(x) + \frac{|x|^2}{2}&\ge C_\alpha \quad \text {for every}\,x\in {\mathbb {R}}^n, \end{aligned}$$
(3.4)

with \(C_\alpha = I_\alpha (\mu _\alpha )- \frac{1}{2}\int _{{\mathbb {R}}^n} |x|^2 \,d\mu _\alpha (x)\). Moreover, the spheroid \(\Omega _\alpha \) is prolate for \(\alpha \in (-1,0)\) and oblate for \(\alpha \in (0,n-2]\).

Remark 3.2

For \(\alpha =0\), \(n\ge 3\), it is well-known that the unique minimiser of the Coulomb energy \(I_0\) is the normalised characteristic function of the ball centred at 0 with radius \((n-2)^{1/n}\) (see, e.g., [11, Corollary 1.3]). In other words, \(a(0)=b(0)=(n-2)^{1/n}\). In the proof below we will focus only on the case \(\alpha \ne 0\).

We split the proof of Theorem 3.1 into Sect. 3.1, where we prove the stationarity condition (3.3), and Sect. 3.2, where we prove (3.4). The heart of the proof consists in the exact evaluation of the convolution

$$\begin{aligned} \Phi _\alpha := W_\alpha *\frac{1}{|\Omega (a,b)|}\chi _{\Omega (a,b)} \end{aligned}$$

both in \(\Omega (a,b)\) (in Sect. 3.1) and in \({\mathbb {R}}^n{\setminus } \Omega (a,b)\) (in Sect. 3.2), for \(a,b>0\). Note that \(\Phi _\alpha \in C^1({\mathbb {R}}^n)\), as can be seen by e.g. adapting the proof of [17, Lemma 4.1]).

We write

$$\begin{aligned} \Phi _\alpha (x)= \Phi _0(x) + \alpha \Psi (x), \quad \text {with } \, \Psi (x):= \frac{1}{|\Omega (a,b)|}\int _{\Omega (a,b)}\frac{(x_1-y_1)^2}{|x-y|^n} \,dy. \end{aligned}$$
(3.5)

3.1 The condition (3.3) on spheroids

Let \(\alpha \in (-1,n-2]\). We claim that there exist \(a(\alpha ), b(\alpha )>0\) (with \(a(\alpha )< b(\alpha )\) for \(\alpha \in (0,n-2]\) and \(a(\alpha )> b(\alpha )\) for \(\alpha \in (-1,0)\)) such that

$$\begin{aligned} (W_{\alpha }*\mu _\alpha )(x) + \frac{|x|^2}{2} = C_\alpha \quad \text {for every}\, x\in \Omega _\alpha , \qquad \mu _\alpha := \frac{1}{| \Omega _\alpha |} \, \chi _{ \Omega _\alpha }, \end{aligned}$$
(3.6)

where we recall that \(\Omega _\alpha = \Omega (a(\alpha ), b(\alpha ))\).

In the two-dimensional case studied in [10], we computed the semi-axes of \(\Omega _\alpha \) in terms of \(\alpha \), and deduced the explicit values \(a=\sqrt{1-\alpha }\) and \(b=\sqrt{1+\alpha }\). For \(n\ge 3\) we do not have explicit expressions for the semi-axes in terms of \(\alpha \).

3.1.1 The potential inside a spheroid

In this section, for \(\alpha \in (-1,n-2]\) and \(a,b>0\), we evaluate \(\Phi _\alpha (x)\) with \(x\in \Omega (a,b)\). We start by recalling the case \(\alpha =0\) of the Coulomb potential, namely

$$\begin{aligned} \Phi _0(x)= \frac{1}{|\Omega (a,b)|}\int _{\Omega (a,b)}\frac{1}{|x-y|^{n-2}} \,dy =:\mathop {\,{--}\int }\nolimits _{\!\!\Omega (a,b)}\frac{1}{|x-y|^{n-2}} \,dy. \end{aligned}$$

From here onwards we use a barred integral to denote the mean over the domain of integration. For \(x\in \Omega (a,b)\) we have that

$$\begin{aligned} \Phi _0(x)&= \frac{n(n-2)}{4}\int _0^\infty \left( 1-\frac{x_1^2}{a^2+s} - \frac{r^2}{b^2+s}\right) \frac{ds}{\sqrt{a^2+s}(b^2+s)^{\frac{n-1}{2}}}\nonumber \\&= - \frac{n(n-2)}{4b^n} \left( x_1^2 \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}} + r^2 \int _0^\infty \frac{d\sigma }{\sqrt{t+\sigma }(1+\sigma )^{\frac{n+1}{2}}}\right) \nonumber \\&+ C(a^2,b^2), \end{aligned}$$
(3.7)

where \(r^2=\sum _{i=2}^n x_i^2\), \(C(a^2,b^2)\) is a constant that depends smoothly on \(a^2\) and \(b^2\), and in the last step we set \(\sigma :=s/b^2\) and denoted with t the aspect ratio \(t:=a^2/b^2\), \(t>0\) (see, e.g., [13]). In particular, \(\Phi _0\) in \(\Omega (a,b)\) is a second-degree polynomial with no linear terms.

We now obtain the anisotropic term \(\Psi \) of \(\Phi _\alpha \) on \(\Omega (a,b)\) [see (3.5)] by differentiating \(\Phi _0\) with respect to the aspect ratio t, in the spirit of [10, Section 4]. First of all note that, by the definition of \(\Phi _0\) and by a change of variables,

$$\begin{aligned} \Phi _0(b\sqrt{t} u_1,bu') = \frac{1}{b^{n-2}}\mathop {\,{--}\int }\nolimits _{\!\!B_1(0)} \frac{dv}{\big (t(u_1-v_1)^2+|u'-v'|^2\big )^{\frac{n-2}{2}}}, \end{aligned}$$
(3.8)

where \(u':=(u_2,\dots ,u_n)\), and \(u=(u_1,u')\in B_1(0)\). By differentiating (3.8) with respect to the aspect ratio t we obtain

$$\begin{aligned} \frac{\partial }{\partial t} \Big (\Phi _0(b\sqrt{t} u_1,bu')\Big ) = -\frac{1}{b^{n-2} }\frac{n-2}{2} \mathop {\,{--}\int }\nolimits _{\!\!B_1(0)} \frac{ (u_1-v_1)^2}{\big (t(u_1-v_1)^2+|u'-v'|^2\big )^{n/2}}\,dv , \end{aligned}$$
(3.9)

and since

$$\begin{aligned} \Psi (b\sqrt{t} u_1,bu') = \frac{1}{b^{n-2}} \mathop {\,{--}\int }\nolimits _{\!\!B_1(0)} \frac{t (u_1-v_1)^2}{\big (t(u_1-v_1)^2+|u'-v'|^2\big )^{n/2}}\, dv, \end{aligned}$$

it follows by (3.9) that

$$\begin{aligned} \Psi (b\sqrt{t} u_1,bu') =-\frac{2t}{n-2} \frac{\partial }{\partial t} \Big (\Phi _0(b\sqrt{t} u_1,bu')\Big ). \end{aligned}$$
(3.10)

On the other hand, by the explicit expression (3.7) we have that for every \(u=(u_1,u')\in B_1(0)\)

$$\begin{aligned} \Phi _0(b\sqrt{t} u_1,bu')&= - u_1^2 \, \frac{n(n-2) t}{4 b^{n-2}} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}}\\&- |u'|^2 \, \frac{n(n-2)}{4 b^{n-2}} \int _0^\infty \frac{d\sigma }{\sqrt{t+\sigma }(1+\sigma )^{\frac{n+1}{2}}} + C(tb^2,b^2), \end{aligned}$$

and hence

$$\begin{aligned}&\frac{\partial }{\partial t} \Big (\Phi _0(b\sqrt{t} u_1,bu') \Big )\nonumber \\&\quad = - u_1^2 \, \frac{n(n-2)}{4 b^{n-2}} \left( \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}} - \frac{3 t}{2} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{5/2}(1+\sigma )^{\frac{n-1}{2}}}\right) \nonumber \\&\qquad + |u'|^2\, \frac{n(n-2)}{8 b^{n-2}} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}} +{\tilde{C}}(tb^2,b^2), \end{aligned}$$
(3.11)

where \({\tilde{C}}(tb^2,b^2)\) is another constant. So, by (3.10) and (3.11) we obtain the expression of \(\Psi \) for \(x\in \Omega (a,b)\), namely

$$\begin{aligned} \Psi (x)&= x_1^2 \, \frac{n}{2b^n} \left( \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}} - \frac{3t}{2} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{5/2}(1+\sigma )^{\frac{n-1}{2}}}\right) \nonumber \\&\quad - r^2 \, \frac{nt}{4 b^n} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}} -\frac{2t}{n-2} {\tilde{C}}(tb^2,b^2). \end{aligned}$$

In conclusion, for \(x\in \Omega (a,b)\),

$$\begin{aligned}&\Phi _\alpha (x)\nonumber \\&\quad = x_1^2 \frac{n}{4 b^n} \bigg ((2\alpha -(n-2))\int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}} - 3\alpha t\int _0^\infty \frac{d\sigma }{(t+\sigma )^{5/2}(1+\sigma )^{\frac{n-1}{2}}}\bigg )\nonumber \\&\qquad + r^2 \frac{n}{4 b^n} \bigg (-(n-2)\int _0^\infty \frac{d\sigma }{\sqrt{t+\sigma } (1+\sigma )^{\frac{n+1}{2}}} - \alpha t\int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}}\bigg ) + C, \end{aligned}$$
(3.12)

where C denotes a constant and \(t>0\).

3.1.2 The condition (3.3) on spheroids

We now use the expression (3.12) of the potential on spheroids to verify that there is a spheroid for which the first Euler–Lagrange condition (3.3) is satisfied.

We start by establishing some relations among the integrals appearing in the expressions of the coefficients of \(x_1^2\) and \(r^2\). Note that, by defining

$$\begin{aligned} H(t):=\int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n-1}{2}}} \end{aligned}$$
(3.13)

for \(t>0\), we have that

$$\begin{aligned} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{5/2}(1+\sigma )^{\frac{n-1}{2}}}&=-\frac{2}{3} H'(t), \end{aligned}$$
(3.14)
$$\begin{aligned} \int _0^\infty \frac{d\sigma }{\sqrt{t+\sigma }(1+\sigma )^{\frac{n+1}{2}}}&=\frac{2}{n-1}\frac{1}{\sqrt{t}} - \frac{1}{n-1} H(t), \end{aligned}$$
(3.15)
$$\begin{aligned} \int _0^\infty \frac{d\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}}&=\frac{2}{n-1}\frac{1}{t^{3/2}} + \frac{2}{n-1}H'(t). \end{aligned}$$
(3.16)

Note also that

$$\begin{aligned} -n H(t)+2(1-t)H'(t) + \frac{2}{t^{3/2}}=0. \end{aligned}$$
(3.17)

To prove it, we rewrite (3.16) as

$$\begin{aligned}&\frac{2}{n-1}H'(t) = \frac{1}{t} \int _0^\infty \frac{t+\sigma -\sigma }{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}}\,d\sigma - \frac{2}{n-1}\frac{1}{t^{3/2}} \\&\quad = \frac{1}{t} \left( \frac{2}{n-1}\frac{1}{\sqrt{t}} - \frac{1}{n-1} H(t) \right) - \frac{1}{t} \int _0^\infty \frac{1+ \sigma -1}{(t+\sigma )^{3/2}(1+\sigma )^{\frac{n+1}{2}}}\,d\sigma - \frac{2}{n-1}\frac{1}{t^{3/2}} \\&\quad = - \frac{1}{n-1} \frac{1}{t} H(t) - \frac{1}{t} H(t) +\frac{1}{t} \left( \frac{2}{n-1}\frac{1}{t^{3/2}} + \frac{2}{n-1}H'(t) \right) , \end{aligned}$$

hence obtaining the relation (3.17).

Integrating (3.17) we deduce that

$$\begin{aligned} H(t)=\frac{1}{|1-t|^{\frac{n}{2}}}\int _t^1 \frac{|1-s|^{\frac{n}{2}}}{s^{3/2}(1-s)}\, ds \qquad \text { for}\, t\ne 1. \end{aligned}$$

By integration by parts we obtain

$$\begin{aligned} H(t)= \frac{2}{\sqrt{t}(1-t)}- \frac{n-2}{|1-t|^{n/2}}\int _t^1 \frac{|1-s|^{\frac{n}{2}-2}}{\sqrt{s}}\, ds \qquad \text { for}\, t\ne 1, \end{aligned}$$
(3.18)

and by differentiation

$$\begin{aligned} H'(t)= \frac{(n+1)t-1}{t^{3/2}(1-t)^2}- \frac{n(n-2)}{2}\frac{1-t}{|1-t|^{\frac{n}{2}+2}}\int _t^1 \frac{|1-s|^{\frac{n}{2}-2}}{\sqrt{s}}\, ds \qquad \text { for}\, t\ne 1.\qquad \end{aligned}$$
(3.19)

Condition (3.3) on spheroids is equivalent to the following two equations, obtained by equating the coefficients of \(x_1^2\):

$$\begin{aligned} -\frac{1}{2} = \frac{n}{4b^n}\left( (2\alpha -(n-2))H(t)+2\alpha t H'(t)\right) , \end{aligned}$$
(3.20)

and of \(r^2\):

$$\begin{aligned} -\frac{1}{2} = \frac{n}{4(n-1)b^n}\left( -2(n-2+\alpha )\frac{1}{\sqrt{t}} +(n-2) H(t)-2\alpha t H'(t)\right) . \end{aligned}$$
(3.21)

For what follows it is more convenient to reduce to two alternative equivalent conditions: the first one is obtained by adding (3.20) to \((n-1)\)-times (3.21):

$$\begin{aligned} b^n = \frac{n-2+\alpha }{\sqrt{t}} - \alpha H(t), \end{aligned}$$
(3.22)

and the second condition is obtained by subtracting (3.20) from (3.21):

$$\begin{aligned} \alpha \left( (n-1)H(t)+nt H'(t)+\frac{1}{\sqrt{t}}\right) + (n-2)\left( -\frac{n}{2} H(t)+\frac{1}{\sqrt{t}}\right) =0. \end{aligned}$$
(3.23)

We claim that, for every \(\alpha \in (-1,n-2]\) there exists a spheroid satisfying (3.22) and (3.23), and hence satisfying the stationarity condition (3.3). More precisely, we will show that there exists a pair \((a(\alpha ), b(\alpha ))\) (or equivalently \((t(\alpha ), b(\alpha ))\)) satisfying (3.22) and (3.23).

3.1.3 The equation (3.23)

We denote with \(F(t,\alpha )\), for \(\alpha \in {\mathbb {R}}\) and \(t>0\) (we recall that \(t=a^2/b^2\)), the left-hand side of (3.23). Then (3.23) is of the form \(F(t,\alpha )=0\). We write

$$\begin{aligned} F(t,\alpha ) = \frac{1}{\sqrt{t}}\big (A(t)\alpha +B(t)\big ), \quad (t,\alpha )\in (0,+\infty )\times {\mathbb {R}}, \end{aligned}$$

where

$$\begin{aligned} A(t):= (n-1)\sqrt{t} H(t)+nt^{3/2} H'(t)+1, \qquad B(t):= -\frac{n(n-2)}{2} \sqrt{t} H(t) + (n-2). \end{aligned}$$
(3.24)

We claim that for every \(\alpha \in (-1,n-2]\) there exists \(t=t(\alpha )>0\) such that \(F(t(\alpha ),\alpha )=0\). We start by analysing the behaviour of A and B for t close to zero. By (3.18) we have that

$$\begin{aligned} \sqrt{t} H(t) = \frac{2}{1-t}- \frac{(n-2)\sqrt{t}}{(1-t)^{n/2}}\int _t^1 \frac{(1-s)^{\frac{n}{2}-2}}{\sqrt{s}}\, ds \qquad \text { for}\, 0<t<1, \end{aligned}$$
(3.25)

and since

$$\begin{aligned} \int _0^1 \frac{(1-s)^{\frac{n}{2}-2}}{\sqrt{s}}\, ds<+\infty , \end{aligned}$$
(3.26)

we immediately deduce that \(\lim _{t\rightarrow 0^+} \sqrt{t} H(t) = 2\). Moreover, by (3.19),

$$\begin{aligned} t^{3/2} H'(t) = \frac{(n+1)t-1}{(1-t)^2}- \frac{n(n-2)t^{3/2}}{2(1-t)^{\frac{n}{2}+1}}\int _t^1 \frac{(1-s)^{\frac{n}{2}-2}}{\sqrt{s}}\, ds \qquad \text { for}\, 0<t<1, \qquad \end{aligned}$$
(3.27)

thus, by (3.26) we deduce that \(\lim _{t\rightarrow 0^+} t^{3/2} H'(t) = -1\). This implies that

$$\begin{aligned} \lim _{t\rightarrow 0^+} A(t) = \lim _{t\rightarrow 0^+} \left( (n-1)\sqrt{t} H(t)+ n t^{3/2} H'(t) + 1\right) = 2(n-1) - n +1 =n-1, \end{aligned}$$

and

$$\begin{aligned} \lim _{t\rightarrow 0^+} B(t) = (n-2)\lim _{t\rightarrow 0^+} \left( -\frac{n}{2} \sqrt{t} H(t) + 1\right) = -(n-1)(n-2). \end{aligned}$$

Hence, if \(\alpha \ne n-2\), \(\lim _{t \rightarrow 0^+} F(t,\alpha ) = -\infty \). If \(\alpha =n-2\), by (3.24) we have

$$\begin{aligned} \sqrt{t}F(t,n-2)= \frac{(n-2)^2}{2} (\sqrt{t} H(t)-2) + n(n-2)(t^{3/2}H'(t)+1). \end{aligned}$$

Using (3.25) and (3.27), we obtain

$$\begin{aligned} F(t,n-2)= & {} \frac{(n-2)\sqrt{t}}{(1-t)^2}(2t+n^2-2)\\&\quad - \frac{(n-2)^2}{2(1-t)^{\frac{n}{2}+1}}\left( (n-2)(1-t)+n^2t\right) \int _t^1 \frac{(1-s)^{\frac{n}{2}-2}}{\sqrt{s}}\, ds \end{aligned}$$

for \(0<t<1\). By (3.26) we deduce that

$$\begin{aligned} \lim _{t \rightarrow 0^+} F(t, n-2)= -\frac{(n-2)^3}{2}\int _0^1 \frac{(1-s)^{\frac{n}{2}-2}}{\sqrt{s}}\, ds< 0. \end{aligned}$$

By a direct computation from (3.13) and (3.14) we have that \(H(1) = \frac{2}{n}\) and \(H'(1)=-\frac{3}{(n+2)}\), therefore \(F(1,\alpha ) = \frac{4(n-1)}{n(n+2)}\alpha \). Finally, one can check directly that

$$\begin{aligned} \lim _{t\rightarrow +\infty }A(t)=1 \qquad \text { and} \qquad \lim _{t\rightarrow +\infty }B(t)=n-2. \end{aligned}$$
(3.28)

Now fix \({\bar{\alpha }}\in (0,n-2]\): then, since \(\lim _{t \rightarrow 0^+} F(t,{\bar{\alpha }})< 0\) and \(F(1,{\bar{\alpha }})>0\), and \(F(\cdot ,{\bar{\alpha }})\) is continuous on \((0,\infty )\), it follows that there exists at least one \(t({\bar{\alpha }})\in (0,1)\) such that \(F(t({\bar{\alpha }}),{\bar{\alpha }})=0\). Similarly, fixing \({\bar{\alpha }}\in (-1,0)\), since \(F(1,{\bar{\alpha }})<0\) and \(\lim _{t\rightarrow +\infty } F(t,{\bar{\alpha }}) = 0^+\), it follows that there exists at least one \(t({\bar{\alpha }})>1\) such that \(F(t({\bar{\alpha }}),{\bar{\alpha }})=0\).

In conclusion, for every \(\alpha \in (-1,n-2]\) there exists at least one \(t(\alpha )> 0\) such that \(F(t(\alpha ),\alpha ) = 0\); in other words, for every \(\alpha \in (-1,n-2]\) there exists a solution \(t(\alpha )>0\) of (3.23).

3.1.4 The equation (3.22)

Now we solve (3.22) for the \(t(\alpha )\) found above by solving (3.23) for spheroids and we compute the corresponding b. Then the spheroid will be the one with semi-axes \(b(\alpha )\) and \(a(\alpha )\), with \(t(\alpha )=a(\alpha )^2/b(\alpha )^2\). From (3.22) and by (3.15)

$$\begin{aligned} 0< \sqrt{t} H(t) =2 - (n-1)\sqrt{t} \int _0^\infty \frac{d \sigma }{\sqrt{t+\sigma }(1+\sigma )^{\frac{n+1}{2}}} < 2 \qquad \text { for}\, t>0, \end{aligned}$$

we have that, for \(\alpha \in (0,n-2]\)

$$\begin{aligned} b^n = \frac{1}{\sqrt{t}}\big (n-2+\alpha - \alpha \sqrt{t} H(t)\big ) > \frac{1}{\sqrt{t}}\big (n-2 - \alpha \big )\ge 0, \end{aligned}$$
(3.29)

and for \(\alpha \in (-1,0)\)

$$\begin{aligned} b^n = \frac{1}{\sqrt{t}}\big (n-2+\alpha - \alpha \sqrt{t} H(t)\big )> \frac{1}{\sqrt{t}}\big (n-2+\alpha ) >0. \end{aligned}$$
(3.30)

3.2 The condition (3.4) outside spheroids

In this section we show that for \(\alpha \in (-1,n-2]\) and for any spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for which the stationarity condition (3.3) is satisfied, also the unilateral condition (3.4) is satisfied. This implies that any spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for which the stationarity condition (3.3) is satisfied is in fact a minimiser for the functional \(I_\alpha \), and by Proposition 2.1 it is the unique minimiser (which in particular implies that there is only one spheroid satisfying (3.3)).

We do it in several steps. We start by evaluating the Coulomb potential, which corresponds to \(\alpha =0\), namely

$$\begin{aligned} \Phi _0(x)= \mathop {\,{--}\int }\nolimits _{\!\!\Omega (a,b)}\frac{1}{|x-y|^{n-2}} \,dy, \end{aligned}$$

for \(x\in {\mathbb {R}}^n{\setminus }\Omega (a,b)\) and \(a,b>0\). For \(x\in {\mathbb {R}}^n{\setminus }\Omega (a,b)\) and \(a,b>0\),

$$\begin{aligned} \Phi _0(x) = \frac{n(n-2)}{4}\int _{\lambda (x)}^\infty \left( 1-\frac{x_1^2}{a^2+s} - \frac{r^2}{b^2+s}\right) \frac{ds}{\sqrt{a^2+s}(b^2+s)^{\frac{n-1}{2}}}, \end{aligned}$$
(3.31)

where \(r^2=\sum _{i=2}^n x_i^2\) and \(\lambda (x)\) is the largest root of the equation

$$\begin{aligned} \frac{x_1^2}{a^2+\lambda } +\frac{r^2}{b^2+\lambda }=1, \end{aligned}$$

(see, e.g., [13]). By straightforward computations one can see that

$$\begin{aligned} \lambda (x)= {\left\{ \begin{array}{ll} \displaystyle \frac{1}{4}\left( \sqrt{x_1^2+(r+c)^2} + \sqrt{x_1^2+(r-c)^2}\right) ^2 - b^2 \quad &{}\text {if } a<b,\\ \displaystyle \frac{1}{4}\left( \sqrt{(x_1+c)^2+r^2} + \sqrt{(x_1-c)^2+r^2}\right) ^2 - a^2 \quad &{}\text {if } a>b, \end{array}\right. } \end{aligned}$$

where

$$\begin{aligned} c^2:={\left\{ \begin{array}{ll} b^2-a^2 \quad &{}\text {if } a<b,\\ a^2-b^2 \quad &{}\text {if } a>b. \end{array}\right. } \end{aligned}$$

3.2.1 The anisotropic potential outside a spheroid

In this section we prove that the anisotropic term \(\Psi \) of \(\Phi _\alpha \), both for oblate and prolate spheroids, is related to the Coulomb potential \(\Phi _0\) by the relation

$$\begin{aligned}&\Psi (x) = \frac{a^2}{a^2-b^2} \, \Phi _0(x)+ \frac{1}{(n-2)(a^2-b^2)} \nabla \Phi _0(x)\cdot \left( b^2 x_1, a^2 x'\right) , \end{aligned}$$
(3.32)

\(\text {for } x=(x_1,x')\in {\mathbb {R}}^n{\setminus } \Omega (a,b)\).We obtain (3.32) by an ingenious differentiation of \(\Phi _0\). While in Sect. 3.1\(\Psi \) on \(\Omega (a,b)\) was obtained by differentiating \(\Phi _0\) with respect to the aspect ratio \(a^2/b^2\) of the spheroid, the geometric quantity that is relevant in this case is the parameter spanning a family of spheroids confocal with \(\Omega (a,b)\), and surrounding it from the outside.

We prove (3.32) in the case of oblate spheroids, but the case of prolate spheroids is completely analogous. For oblate spheroids we set \(a^2=t\) and \(b^2=t+c^2\), where c is fixed, and set

$$\begin{aligned} \Phi ^t_0(x):=\mathop {\,{--}\int }\nolimits _{\!\!\Omega _t}\frac{1}{|x-y|^{n-2}} \,dy, \quad \Psi ^t(x):=\mathop {\,{--}\int }\nolimits _{\!\!\Omega _t}\frac{(x_1-y_1)^2}{|x-y|^n} \,dy, \quad \Omega _t:=\Omega (\sqrt{t}, \sqrt{t+c^2}). \end{aligned}$$

From (3.31) one can easily rewrite the Coulomb potential on a spheroid as

$$\begin{aligned} \Phi _0(x) = \frac{n(n-2)}{4} \int _{\ell (x)}^\infty \Big (1-\frac{x_1^2}{\sigma -c^2} - \frac{r^2}{ \sigma }\Big ) \frac{d\sigma }{\sqrt{\sigma -c^2}\sigma ^{\frac{n-1}{2}}}, \end{aligned}$$

where \(\ell (x) = \lambda (x) + b^2\) and \(\sigma :=s+b^2\), and \(\ell (x)\) depends only on \(c^2\). Hence \(\Phi _0^t\) depends on its semi-axes a and b only via c, and not on a and b separately, namely

$$\begin{aligned} 0=\frac{\partial }{\partial t} \Phi ^t_0(x) = \frac{\partial }{\partial t}\mathop {\,{--}\int }\nolimits _{\!\!B_1(0)}\!\! \bigg (t\bigg (\frac{x_1}{\sqrt{t}}-v_1\bigg )^2\!\!\!+(t+c^2)\sum _{i\ne 1}\bigg (\frac{x_i}{\sqrt{t+c^2}}-v_i\bigg )^2\bigg )^{-\frac{n-2}{2}}\!\!dv. \end{aligned}$$

By expanding the derivative above and rewriting

$$\begin{aligned}&\bigg (\frac{x_1}{\sqrt{t}}-v_1\bigg )^2 + \sum _{i\ne 1}\bigg (\frac{x_i}{\sqrt{t+c^2}}-v_i\bigg )^2 \\&\quad = \frac{1}{t+c^2}\bigg (t\bigg (\frac{x_1}{\sqrt{t}}-v_1\bigg )^2\!\!\!+(t+c^2)\sum _{i\ne 1}\bigg (\frac{x_i}{\sqrt{t+c^2}}-v_i\bigg )^2\bigg ) + \frac{c^2}{t(t+c^2)} t \bigg (\frac{x_1}{\sqrt{t}}-v_1\bigg )^2, \end{aligned}$$

we obtain

$$\begin{aligned} 0=\frac{1}{t+c^2} \, \Phi _0^t(x) + \frac{c^2}{t(t+c^2)}\,\Psi ^t(x) + \frac{1}{n-2}\nabla _x \Phi _0^t(x)\cdot \left( \frac{x_1}{t}, \frac{x'}{t+c^2}\right) . \end{aligned}$$
(3.33)

The expression (3.33) gives the (unknown) expression of the anisotropic term \(\Psi ^t\) in terms of the (known) Coulomb potential \(\Phi _0^t\) and its spatial gradient. Substituting t and c in terms of a and b in (3.33) and rearranging the terms we then have (3.32).

With the expression of the Coulomb potential \(\Phi _0\) [see (3.31)] and a closed formula for the anisotropic potential \(\Psi \) [in (3.32)] in the outer region \({\mathbb {R}}^n{\setminus } \Omega (a,b)\) at hand, we now prove (3.4). More precisely, we prove that for every \(\alpha \in (-1,n-2]\),

$$\begin{aligned} \Phi _\alpha (x) + \frac{|x|^2}{2} \ge C_\alpha \quad \text {for every } x\in {\mathbb {R}}^n{\setminus } \Omega (a(\alpha ),b(\alpha )), \end{aligned}$$
(3.34)

where \(\Omega (a(\alpha ),b(\alpha ))\) is a stationary point, satisfying (3.3). To prove (3.34) we first rewrite the potentials \(\Phi _0\) and \(\Psi \) (and hence \(\Phi _\alpha \)) outside a spheroid in a more convenient way, in terms of a set of coordinates—oblate or prolate spheroidal—alternative to the Euclidean ones, and more suitable for the geometry of the problem.

We deal with the oblate and prolate case separately.

3.2.2 The condition (3.34) outside an oblate spheroid

To prove (3.34) for an oblate spheroid, we rewrite the potentials \(\Phi _0\) in (3.31) and \(\Psi \) in (3.32) outside a spheroid in terms of the oblate spheroidal coordinates. By the symmetry of \(\Phi _\alpha \) and of the confinement, it is sufficient to reduce to computations in the \(x_1 x_2\)-plane. In terms of oblate spheroidal coordinates we have

$$\begin{aligned} {\left\{ \begin{array}{ll} x_1= cz\rho \\ x_2= c\sqrt{(1+z^2)(1-\rho ^2)} \end{array}\right. } \quad z\ge 0, \, \rho \in [-1,1], \end{aligned}$$

where we recall that \(c^2= b^2-a^2\). Note that the outer region \({\mathbb {R}}^n{\setminus } \Omega (a,b)\) corresponds to \(z\ge \frac{a}{c}\).

For \(z\ge \frac{a}{c}\) and \(\rho \in [-1,1]\), the expression of the Coulomb potential (3.31) in oblate spheroidal coordinates reads as

$$\begin{aligned} \Phi _0(z,\rho ) = \frac{n(n-2)}{4} \rho ^2\int _{c^2z^2}^\infty \frac{c^2(\sigma -c^2z^2)}{\sigma ^{3/2}(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma + \frac{n(n-2)}{4}\int _{c^2z^2}^\infty \frac{\sigma -c^2z^2}{\sqrt{\sigma }(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma , \end{aligned}$$

where we used that \(r^2=x_2^2\), \(\lambda (x)=c^2z^2-a^2\) and the change of variables \(\sigma = a^2+s\). We recall that the gradient of the oblate spheroidal coordinates with respect to Cartesian coordinates is given by the following formulas:

$$\begin{aligned} \nabla \rho (x)= & {} \dfrac{1}{c(z^2+\rho ^2)}\big (z(1-\rho ^2), -\rho \sqrt{(1+z^2)(1-\rho ^2)}\big ), \\ \nabla z (x)= & {} \dfrac{1}{c(z^2+\rho ^2)}\big (\rho (1+z^2), z\sqrt{(1+z^2)(1-\rho ^2)}\big ). \end{aligned}$$

Since

$$\begin{aligned} \partial _z \Phi _0(z,\rho )&= -\frac{n(n-2)}{2} c^2 z\int _{c^2z^2}^\infty \frac{c^2\rho ^2+\sigma }{\sigma ^{3/2}(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma ,\\ \partial _\rho \Phi _0(z,\rho )&= -\frac{n(n-2)}{2} c^2 \rho \int _{c^2z^2}^\infty \frac{c^2 z^2-\sigma }{\sigma ^{3/2}(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma , \end{aligned}$$

we deduce that

$$\begin{aligned} \nabla \Phi _0(x) = -\frac{c \, n (n-2)}{2}\left( \int _{c^2z^2}^\infty \frac{z\rho \, d\sigma }{\sigma ^{3/2}(\sigma +c^2)^{\frac{n-1}{2}}}, \int _{c^2z^2}^\infty \frac{\sqrt{(1+z^2)(1-\rho ^2)} \, d\sigma }{\sqrt{\sigma }(\sigma +c^2)^{\frac{n+1}{2}}}\right) , \end{aligned}$$

and

$$\begin{aligned}&- \frac{1}{c^2} \nabla \Phi _0(x)\cdot (b^2x_1,a^2x_2) \\&\quad = \frac{n(n-2)}{4} \left( \int _{c^2z^2}^\infty \frac{2b^2z^2\rho ^2}{\sigma ^{3/2}(\sigma +c^2)^{\frac{n-1}{2}}}\, d\sigma + \int _{c^2z^2}^\infty \frac{2a^2(1+z^2)(1-\rho ^2)}{\sqrt{\sigma }(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma \right) . \end{aligned}$$

Hence the anisotropic potential is

$$\begin{aligned} \Psi (x)&= \frac{n}{4} \rho ^2 \int _{c^2z^2}^\infty \frac{-na^2(\sigma -c^2z^2)+2c^2z^2(\sigma +c^2)}{\sigma ^{3/2}(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma \\&\quad + \frac{n}{4} \int _{c^2z^2}^\infty \frac{-na^2(\sigma -c^2z^2)+2a^2(\sigma +c^2)}{c^2\sqrt{\sigma }(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma . \end{aligned}$$

Finally, the confinement term \(|x|^2/2\), in terms of the spheroidal coordinates, is

$$\begin{aligned} \frac{|x|^2}{2} = \frac{c^2}{2}(1-\rho ^2+z^2). \end{aligned}$$

Note that \(\Phi _0\), \(\Psi \) and the confinement are all quadratic functions in the variable \(\rho \). More precisely, for \(z\ge \frac{a}{c}\) and \(\rho \in [-1,1]\),

$$\begin{aligned} \Phi _\alpha (x)+\frac{|x|^2}{2} = A_\alpha (z)+B_\alpha (z)\rho ^2, \end{aligned}$$
(3.35)

where

$$\begin{aligned} A_\alpha (z)&:= \frac{n}{4} \int _{c^2z^2}^\infty \frac{\big ( (n-2)c^2 - n\alpha a^2\big )(\sigma -c^2z^2) + 2\alpha a^2 (\sigma +c^2)}{c^2\sqrt{\sigma }(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma + \frac{c^2}{2}(1+z^2),\\ B_\alpha (z)&:= \frac{n}{4} \int _{c^2z^2}^\infty \frac{\big ( (n-2)c^2 - n\alpha a^2\big )(\sigma -c^2z^2) + 2\alpha c^2z^2 (\sigma +c^2)}{\sigma ^{3/2}(\sigma +c^2)^{\frac{n+1}{2}}}\, d\sigma - \frac{c^2}{2}. \end{aligned}$$

We now prove (3.34) for any oblate spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for which the stationarity condition (3.3) is satisfied. To avoid burdening the reader with heavy notation, we will write a and b instead of \(a(\alpha )\) and \(b(\alpha )\).

By (3.35), proving (3.34) is now equivalent to show that

$$\begin{aligned} A_\alpha (z)+B_\alpha (z)\ge C_\alpha \quad \text {and } \quad A_\alpha (z)\ge C_\alpha \quad \text {for } z\ge \frac{a}{c}, \end{aligned}$$
(3.36)

namely to check the inequality for \(\rho =0\) and \(\rho =1\). Moreover, note that the functions \(A_\alpha \) and \(B_\alpha \) are well-defined and smooth at every \(z>0\). Since, as observed before, the potential \(\Phi _\alpha + |\cdot |^2/2\) belongs to \(C^1({\mathbb {R}}^2)\) and satisfies (3.3), Eq. (3.35) implies that

$$\begin{aligned} A_\alpha \Big (\frac{a}{c}\Big )=A_\alpha \Big (\frac{a}{c}\Big )+ B_\alpha \Big (\frac{a}{c}\Big )=C_\alpha \qquad \text { and } \qquad A_\alpha '\Big (\frac{a}{c}\Big )=A_\alpha '\Big (\frac{a}{c}\Big )+ B'_\alpha \Big (\frac{a}{c}\Big )=0.\nonumber \\ \end{aligned}$$
(3.37)

We show that (3.36) is satisfied by proving that

$$\begin{aligned} \Big (\frac{1}{z} \big (A'_\alpha (z)+B'_\alpha (z)\big )\Big )'\ge 0 \quad \text {and } \quad \Big (\frac{1}{z} A'_\alpha (z)\Big )'\ge 0 \quad \text {for } z\ge \frac{a}{c}. \end{aligned}$$
(3.38)

Clearly (3.38), together with the second condition in (3.37), implies

$$\begin{aligned} A'_\alpha (z)+ B'_\alpha (z)\ge 0 \quad \text {and } \quad A'_\alpha (z)\ge 0 \quad \text {for } z\ge \frac{a}{c}, \end{aligned}$$

which, in turn, gives (3.36) owing to the first condition in (3.37). While (3.38) may look more complicated, it actually gives rise to simpler computations. Indeed we have

$$\begin{aligned} \Big (\frac{1}{z} A'_\alpha (z)\Big )' = \frac{n}{c^{n-2} z^2(1+z^2)^{\frac{n+1}{2}}}\Big ((n-2)z^2+\alpha \frac{a^2}{c^2}\Big ) \ge 0 \quad \text {for every } z\ge \frac{a}{c}, \end{aligned}$$

since \(n\ge 3\) and \(\alpha \ge 0\). Moreover,

$$\begin{aligned} \Big (\frac{1}{z} \big (A'_\alpha (z)+B'_\alpha (z)\big )\Big )' = \frac{n}{c^{n-2} z^2(z^2+1)^{\frac{n+1}{2}}} \Big ( (n-2)(1+\alpha )z^2 +n-2-\alpha -\alpha \frac{a^2}{c^2}(n-1)\Big ). \end{aligned}$$

Hence to prove the claim (3.38) it is sufficient to show that

$$\begin{aligned} (n-2)(1+\alpha )z^2 +n-2-\alpha -\alpha \frac{a^2}{c^2}(n-1)\ge 0 \qquad \text {for}\, z\ge \frac{a}{c}. \end{aligned}$$

This condition is satisfied since

$$\begin{aligned} (n-2)(1+\alpha ) \frac{a^2}{c^2} +n-2-\alpha -\alpha \frac{a^2}{c^2}(n-1) = (n-2-\alpha )\left( 1+\frac{a^2}{c^2}\right) \ge 0. \end{aligned}$$

3.2.3 The condition (3.34) outside a prolate spheroid

To prove (3.34) for a prolate spheroid, we rewrite the potentials \(\Phi _0\) in (3.31) and \(\Psi \) in (3.32) outside a spheroid in terms of the prolate spheroidal coordinates. By the symmetry of \(\Phi _\alpha \) and of the confinement, it is sufficient to reduce to computations in the \(x_1 x_2\)-plane. In terms of prolate spheroidal coordinates we have

$$\begin{aligned} {\left\{ \begin{array}{ll} x_1= cz\rho \\ x_2= c\sqrt{(z^2-1)(1-\rho ^2)} \end{array}\right. } \quad z\ge 1, \, \rho \in [-1,1], \end{aligned}$$

where we recall that now \(c^2= a^2-b^2\). The outer region \({\mathbb {R}}^n{\setminus } \Omega (a,b)\) corresponds also in this case to \(z\ge \frac{a}{c}\).

For \(z\ge \frac{a}{c}\) and \(\rho \in [-1,1]\), the expression of the Coulomb potential (3.31) in prolate spheroidal coordinates reads as

$$\begin{aligned} \Phi _0(z,\rho ) = \frac{n(n-2)}{4} \rho ^2\int _{c^2z^2}^\infty \frac{c^2(c^2z^2-\sigma )}{\sigma ^{3/2}(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma + \frac{n(n-2)}{4}\int _{c^2z^2}^\infty \frac{\sigma -c^2z^2}{\sqrt{\sigma }(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma , \end{aligned}$$

where we used that \(r^2=x_2^2\), \(\lambda (x)=c^2z^2-a^2\) and a change of variables. We recall that the gradient of the prolate spheroidal coordinates with respect to Cartesian coordinates is given by the following formulas:

$$\begin{aligned} \nabla \rho (x)= & {} \dfrac{1}{c(z^2-\rho ^2)}\big (z(1-\rho ^2), -\rho \sqrt{(z^2-1)(1-\rho ^2)}\big ), \\ \nabla z(x)= & {} \dfrac{1}{c(z^2-\rho ^2)}\big (\rho (z^2-1), z\sqrt{(z^2-1)(1-\rho ^2)}\big ). \end{aligned}$$

Since

$$\begin{aligned} \partial _z \Phi _0(z,\rho )&= \frac{n(n-2)}{2} c^2 z\int _{c^2z^2}^\infty \frac{c^2\rho ^2-\sigma }{\sigma ^{3/2}(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma ,\\ \partial _\rho \Phi _0(z,\rho )&= \frac{n(n-2)}{2} c^2 \rho \int _{c^2z^2}^\infty \frac{c^2 z^2-\sigma }{\sigma ^{3/2}(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma , \end{aligned}$$

we deduce that

$$\begin{aligned} \nabla \Phi _0(x) = -\frac{c \, n (n-2)}{2}\left( \int _{c^2z^2}^\infty \frac{z\rho \, d\sigma }{\sigma ^{3/2}(\sigma -c^2)^{\frac{n-1}{2}}}, \int _{c^2z^2}^\infty \frac{\sqrt{(z^2-1)(1-\rho ^2)} \, d\sigma }{\sqrt{\sigma }(\sigma -c^2)^{\frac{n+1}{2}}}\right) , \end{aligned}$$

and

$$\begin{aligned}&\frac{1}{c^2} \nabla \Phi _0(x)\cdot (b^2x_1,a^2x_2) \\&\quad = -\frac{n(n-2)}{4} \left( \int _{c^2z^2}^\infty \frac{2b^2z^2\rho ^2}{\sigma ^{3/2}(\sigma -c^2)^{\frac{n-1}{2}}}\, d\sigma + \int _{c^2z^2}^\infty \frac{2a^2(z^2-1)(1-\rho ^2)}{\sqrt{\sigma }(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma \right) . \end{aligned}$$

Hence the anisotropic potential is

$$\begin{aligned} \Psi (x)&= \frac{n}{4} \rho ^2 \int _{c^2z^2}^\infty \frac{-na^2(\sigma -c^2z^2)+2c^2z^2(\sigma -c^2)}{\sigma ^{3/2}(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma \\&\quad + \frac{n}{4} \int _{c^2z^2}^\infty \frac{na^2(\sigma -c^2z^2)-2a^2(\sigma -c^2)}{c^2\sqrt{\sigma }(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma . \end{aligned}$$

Finally, the confinement term \(|x|^2/2\), in terms of the spheroidal coordinates, is

$$\begin{aligned} \frac{|x|^2}{2} = \frac{c^2}{2}(-1+\rho ^2+z^2). \end{aligned}$$

Note that, as before, \(\Phi _0\), \(\Psi \) and the confinement are all quadratic functions in the variable \(\rho \). More precisely, for \(z\ge \frac{a}{c}\) and \(\rho \in [-1,1]\),

$$\begin{aligned} \Phi _\alpha (x)+\frac{|x|^2}{2} = A_\alpha (z)+B_\alpha (z)\rho ^2, \end{aligned}$$
(3.39)

where

$$\begin{aligned} A_\alpha (z)&:= \frac{n}{4} \int _{c^2z^2}^\infty \frac{\big ( (n-2)c^2 + n\alpha a^2\big )(\sigma -c^2z^2) - 2\alpha a^2 (\sigma -c^2)}{c^2\sqrt{\sigma }(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma + \frac{c^2}{2}(z^2-1),\\ B_\alpha (z)&:= \frac{n}{4} \int _{c^2z^2}^\infty \frac{\big (-(n-2)c^2 - n\alpha a^2\big )(\sigma -c^2z^2) + 2\alpha c^2z^2 (\sigma -c^2)}{\sigma ^{3/2}(\sigma -c^2)^{\frac{n+1}{2}}}\, d\sigma + \frac{c^2}{2}. \end{aligned}$$

We now prove (3.34) for any prolate spheroid \(\Omega (a(\alpha ),b(\alpha ))\) for which the stationarity condition (3.3) is satisfied. Again, to avoid burdening the notation, we will write a and b instead of \(a(\alpha )\) and \(b(\alpha )\). Note that, by (3.39), proving (3.34) is equivalent to show that

$$\begin{aligned} A_\alpha (z)+B_\alpha (z)\ge C_\alpha \quad \text {and } \quad A_\alpha (z)\ge C_\alpha \quad \text {for } z\ge \frac{a}{c}, \end{aligned}$$

namely to check the inequality for \(\rho =0\) and \(\rho =1\). Arguing as in the case of an oblate spheroid, it is in fact sufficient to prove that

$$\begin{aligned} \Big (\frac{1}{z} \big (A'_\alpha (z)+B'_\alpha (z)\big )\Big )'\ge 0 \quad \text {and } \quad \Big (\frac{1}{z} A'_\alpha (z)\Big )'\ge 0 \quad \text {for } z\ge \frac{a}{c}. \end{aligned}$$
(3.40)

We have

$$\begin{aligned} \Big (\frac{1}{z} A'_\alpha (z)\Big )' = \frac{n}{c^{n-2} z^2(z^2-1)^{\frac{n+1}{2}}}\Big ((n-2)z^2+\alpha \frac{a^2}{c^2}\Big ) \ge 0 \quad \text {for every } z\ge \frac{a}{c}, \end{aligned}$$

since \(n\ge 3\) and \(\alpha >-1\). Moreover,

$$\begin{aligned} \Big (\frac{1}{z} \big (A'_\alpha (z)+B'_\alpha (z)\big )\Big )' = \frac{n}{c^{n-2} z^2(z^2-1)^{\frac{n+1}{2}}} \Big ( (n-2)(1+\alpha )z^2 -(n-2)+\alpha -\alpha \frac{a^2}{c^2}(n-1)\Big ). \end{aligned}$$

Hence to prove the claim (3.40) it is sufficient to show that

$$\begin{aligned} (n-2)(1+\alpha )z^2 -(n-2)+\alpha -\alpha \frac{a^2}{c^2}(n-1)\ge 0 \qquad \text { for } z\ge \frac{a}{c}. \end{aligned}$$

This condition is satisfied because \(\alpha >-1\) and

$$\begin{aligned} (n-2)(1+\alpha ) \frac{a^2}{c^2} -n+2+\alpha -\alpha \frac{a^2}{c^2}(n-1) = (n-2-\alpha )\left( \frac{a^2}{c^2}-1\right) \ge 0. \end{aligned}$$

This concludes the proof of Theorem 3.1.

Remark 3.3

For every \(\alpha \in (-1,n-2]\), let \(t(\alpha )>0\) be the solution of the equation \(F(t,\alpha )=0\) found in Sect. 3.1.3. Note that this solution is unique. Indeed, by (3.29) and (3.30), for a given \(\alpha \in (-1,n-2]\), any solution \(t(\alpha )>0\) of \(F(t,\alpha )=0\) identifies a spheroid satisfying the stationarity condition (3.3). Moreover, in Sect. 3.2 we show that any stationary spheroid satisfies also condition (3.4), so it is a minimiser of \(I_\alpha \), and hence is unique by strict convexity. Moreover, the function \(t: \alpha \in (-1,n-2] \mapsto t(\alpha )\in (0,+\infty )\) is continuous and strictly decreasing. To prove it, let \(\alpha _0\in (-1,n-2]\) and let \(\alpha \rightarrow \alpha _0\). Let \((\alpha _k)\) denote a subsequence converging monotonically to \(\alpha _0\), as \(k\rightarrow +\infty \); note that the subsequence \((I_{\alpha _k})\) is also monotone. Since \(I_{\alpha _k}\) is lower semicontinuous for every k, we have that \(I_{\alpha _k}\) \(\Gamma \)-converges to \(I_{\alpha _0}\), as \(k\rightarrow +\infty \), with respect to narrow convergence (see, e.g., [12, Propositions 5.4 and 5.7]). By the Urysohn property of \(\Gamma \)-convergence (see for instance [12, Proposition 8.3]) the whole sequence \((I_\alpha )\) \(\Gamma \)-converges to \(I_{\alpha _0}\), as \(\alpha \rightarrow \alpha _0\), since the space of probability measures endowed with the narrow convergence is metrisable. Moreover, the equilibrium measures \(\mu _\alpha \) satisfy

$$\begin{aligned} \int _{{\mathbb {R}}^n}|x|^2\, d\mu _\alpha (x)\le I_\alpha (\mu _\alpha )\le I_\alpha (\mu _{n-2})\le I_{n-2}(\mu _{n-2}) \qquad \text { for every}\, \alpha \in (-1,n-2], \end{aligned}$$

where we used the minimality of \(\mu _\alpha \). In other words, the equilibrium measures are tight. By the Fundamental Theorem of \(\Gamma \)-convergence (see, e.g., [12, Corollary 7.17]) and the uniqueness of the minimiser of \(I_{\alpha _0}\) we deduce that \(\mu _\alpha \) converges narrowly to \(\mu _{\alpha _0}\), as \(\alpha \rightarrow \alpha _0\). The characterisation of \(\mu _\alpha \) given in Theorem 3.1 allows us to conclude that \(t(\alpha )\rightarrow t(\alpha _0)\), as \(\alpha \rightarrow \alpha _0\), that is, the function \(\alpha \mapsto t(\alpha )\) is continuous. It is also injective by the following argument: If \(\alpha _1,\alpha _2\in (-1,n-2]\), \(\alpha _1\ne \alpha _2\), are such that \(t(\alpha _1)=t(\alpha _2)=t_0>0\), then \(F(t_0,\alpha _1)=0=F(t_0,\alpha _2)\), which implies \(A(t_0)=0\), with A defined in (3.24). Since \(F(t_0,\alpha _1)=0\), we deduce that \(B(t_0)=0\), but this would imply that \(F(t_0,\alpha )=0\) for every \(\alpha \in (-1,n-2]\), that is, the function \(\alpha \mapsto t(\alpha )\) is constant. This is not possible, since \(t(\alpha )>1\) for \(\alpha <0\) and \(t(\alpha )<1\) for \(\alpha >0\). This last property, together with continuity and injectivity, implies that \(\alpha \mapsto t(\alpha )\) is strictly decreasing.

4 The limit case, as \(\alpha \rightarrow -1\)

In this section we discuss the behaviour of the nonlocal energies \(I_\alpha \), as \(\alpha \rightarrow -1\).

Theorem 4.1

As \(\alpha \rightarrow -1^+\), the functionals \(I_\alpha \) \(\Gamma \)-converge, with respect to narrow convergence, to a functional \(J_*:\mathcal P({\mathbb {R}}^n)\rightarrow [0,+\infty ]\), whose unique minimiser is a measure \(\mu _*\) that is the normalised characteristic function of a prolate spheroid \(\Omega _*\). Moreover, the following representation holds:

$$\begin{aligned} J_*(\mu )= \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\mu }}}(\xi )|^2\, d\xi + \int _{{\mathbb {R}}^n} |x|^2 \,d\mu (x) \end{aligned}$$
(4.1)

for every \(\mu \in {\mathcal {P}}({\mathbb {R}}^n)\) with compact support, where \({\widehat{W}}_*\in L^1_{\text {loc}}({\mathbb {R}}^n)\) is the function given by

$$\begin{aligned} {\widehat{W}}_*(\xi ):=\frac{\pi ^{\frac{n}{2}-2}}{2\Gamma (\frac{n}{2})} \frac{(n-1) \xi _1^2 + (n-3) \sum _{i=2}^n \xi _i^2}{|\xi |^4} \end{aligned}$$

for a.e. \(\xi \in {\mathbb {R}}^n\).

Remark 4.2

The functional \(J_*\) does not coincide with the functional \(I_{-1}\) defined in Remark 2.2, since \(J_*\) is lower semicontinuous, whereas \(I_{-1}\) is not. Moreover \(J_*\) does not coincide with the functionals \(\overline{I_{-1}}\) or \(\widetilde{I}_{-1}\) either, since the Dirac delta at 0 is a minimiser for both \(\overline{I_{-1}}\) and \({\widetilde{I}}_{-1}\).

Proof of Theorem 4.1

The proof is subdivided into several steps.

Step 1: Convergence of the equilibrium measures. We prove that, as \(\alpha \rightarrow -1^+\), the equilibrium measures \(\mu _\alpha \) converge to a measure \(\mu _*\), that is the normalised characteristic function of a prolate spheroid \(\Omega _*\). First of all, we claim that

$$\begin{aligned} \lim _{\alpha \rightarrow -1^+} t(\alpha )=t_*\in {\mathbb {R}}\end{aligned}$$
(4.2)

for some \(t_*>1\). Since \(\alpha \mapsto t(\alpha )\) is strictly decreasing by Remark 3.3 and \(t(\alpha )>1\) for \(\alpha <0\), the limit as \(\alpha \rightarrow -1^+\) exists and is strictly greater than 1. Assume by contradiction that \(t_*=+\infty \). By (3.28) we can pass to the limit in the equation

$$\begin{aligned} 0=\sqrt{t(\alpha )} F(t(\alpha ), \alpha )= A(t(\alpha ))\alpha + B(t(\alpha )) \end{aligned}$$
(4.3)

and deduce that \(n-3=0\), which gives a contradiction for \(n>3\). If \(n=3\), Eq. (4.3), together with the assumption that \(t_*=+\infty \), implies that

$$\begin{aligned} -A(t(\alpha ))+ B(t(\alpha ))<0 \end{aligned}$$
(4.4)

for \(\alpha +1>0\) small enough (note that \(A(t(\alpha ))>0\) for \(\alpha +1>0\) small enough by (3.28)). By (3.18) and (3.19), for \(t>1\) and \(n=3\) we obtain

$$\begin{aligned} -A(t)+ B(t)= -\frac{5t+4}{(t-1)^2}+\frac{\sqrt{t}}{2(t-1)^{5/2}}(2t+7)\int _1^t \frac{1}{\sqrt{s(s-1)}}\, ds. \end{aligned}$$

Since the right-hand side is positive as \(t\rightarrow +\infty \), this contradicts (4.4). Claim (4.2) is thus proved for every \(n\ge 3\). Passing to the limit in (3.30) we also obtain that

$$\begin{aligned} \lim _{\alpha \rightarrow -1^+} b(\alpha ) = b_*:=(t_*)^{-\frac{1}{2n}}\big (n-3+\sqrt{t_*}H(t_*)\big )^{\frac{1}{n}}>0. \end{aligned}$$

This implies that the equilibrium measures \(\mu _\alpha \) converge narrowly, as \(\alpha \rightarrow -1^+\), to the normalised characteristic function of the prolate spheroid \(\Omega _*=\Omega (\sqrt{t_*}b_*,b_*)\).

Step 2: \(\Gamma \)-convergence. Let \((\alpha _k)\) be a sequence such that \(\alpha _k\rightarrow -1^+\), as \(k\rightarrow +\infty \). Note that we can extract a decreasing subsequence \((\alpha _{k_j})\searrow -1^+\) along which also the sequence of functionals \((I_{\alpha _{k_j}})\) is decreasing, and hence \(\Gamma \)-convergent as \(j\rightarrow +\infty \) to the functional

$$\begin{aligned} J_*(\mu ):={\overline{J}}(\mu ), \qquad J(\mu ):=\inf _{\alpha \in (0,-1)} I_\alpha (\mu ) \end{aligned}$$

for every \(\mu \in {\mathcal {P}}({\mathbb {R}}^n)\), where \({\overline{J}}\) denotes the lower semicontinuous envelope of J with respect to narrow convergence. By the Urysohn property of \(\Gamma \)-convergence (see for instance [12, Proposition 8.3]) the whole sequence \((I_\alpha )\) \(\Gamma \)-converges to \(J_*\), as \(\alpha \rightarrow -1^+\). By the Fundamental Theorem of \(\Gamma \)-convergence we deduce that \(\mu _*\) is a minimiser of \(J_*\) and

$$\begin{aligned} \lim _{\alpha \rightarrow -1^+} I_{\alpha }(\mu _{\alpha })= J_*(\mu _*). \end{aligned}$$
(4.5)

Step 3: Representation formula for \(J_*\). Using formula (2.4), which holds for every \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\) with compact support and every \(\alpha \in (-1,n-2]\), we can now prove the representation formula (4.1).

We first observe that by (2.5)

$$\begin{aligned} {\widehat{W}}_\alpha (\xi )\rightarrow {\widehat{W}}_*(\xi ) \qquad \text { for a.e.}\, \xi \in {\mathbb {R}}^n, \end{aligned}$$
(4.6)

as \(\alpha \rightarrow -1^+\), and there exists a constant C, independent of \(\alpha \), such that

$$\begin{aligned} 0\le {\widehat{W}}_\alpha \le C {\widehat{W}}_0 \end{aligned}$$
(4.7)

for every \(\alpha \in (-1,n-2]\). Let now \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\) be a measure with compact support. Let \((\alpha _k)\subset (-1,n-2]\) be any sequence converging to \(-1\), and let \((\nu _k)\) be a sequence in \({\mathcal {P}}({\mathbb {R}}^n)\) converging narrowly to \(\nu \) and such that

$$\begin{aligned} \liminf _{k\rightarrow \infty } I_{\alpha _k}(\nu _k)<+\infty . \end{aligned}$$

Up to subsequences, we can assume that \(\sup _k I_{\alpha _k}(\nu _k)<+\infty \). By (2.2) there exists a compact set \(K\subset {\mathbb {R}}^n\), containing the support of \(\nu \), such that \(\nu _k(K)>0\) and

$$\begin{aligned} W_{\alpha _k}(x-y) + \frac{1}{2}(|x|^2+|y|^2)\ge I_{\alpha _k}(\nu _k) +1 \qquad \text { for}\, (x,y)\not \in K\times K \end{aligned}$$

for every k. If we define

(4.8)

for every k, then \(\mu _k\in {\mathcal {P}}({\mathbb {R}}^n)\) has compact support and

$$\begin{aligned} I_{\alpha _k}(\nu _k)= & {} \big (\nu _k(K)\big )^2 I_{\alpha _k}(\mu _k) + \iint _{(K\times K)^c}\Big ( W_{\alpha _k}(x-y) + \frac{1}{2}(|x|^2+|y|^2)\Big )\, d\nu _k(x)\,d\nu _k(y) \\\ge & {} \big (\nu _k(K)\big )^2 I_{\alpha _k}(\mu _k) + \big (1- \big (\nu _k(K)\big )^2\big )(I_{\alpha _k}(\nu _k)+1). \end{aligned}$$

This inequality implies that

$$\begin{aligned} I_{\alpha _k}(\mu _k)\le I_{\alpha _k}(\nu _k) -\frac{1- \big (\nu _k(K)\big )^2}{\big (\nu _k(K)\big )^2}\le I_{\alpha _k}(\nu _k) \end{aligned}$$
(4.9)

for every k. Since \((\mu _k)\) converges narrowly to \(\nu \), as \(k\rightarrow \infty \), we have that \(({{\widehat{\mu }}}_k)\) pointwise converges to \({{\widehat{\nu }}}\) and by the Fatou Lemma

$$\begin{aligned} \liminf _{k\rightarrow \infty } \int _{{\mathbb {R}}^n} \widehat{W}_{\alpha _k}(\xi ) |{{\widehat{\mu }}}_k(\xi )|^2\, d\xi \ge \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi . \end{aligned}$$

By (4.9), (2.4), and the continuity of the confinement term with respect to narrow convergence (on measures with compact support), we obtain

$$\begin{aligned} \liminf _{k\rightarrow \infty } I_{\alpha _k}(\nu _k) \ge \liminf _{k\rightarrow \infty } I_{\alpha _k}(\mu _k) \ge \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi + \int _{{\mathbb {R}}^n} |x|^2 \,d\nu (x). \end{aligned}$$

Since by definition of \(\Gamma \)-convergence

$$\begin{aligned} J_*(\nu )=\min \Big \{ \liminf _{k\rightarrow \infty } I_{\alpha _k}(\nu _k): \ (\nu _k)\rightharpoonup \nu \text { narrowly}, (\alpha _k)\rightarrow -1^+\Big \}, \end{aligned}$$
(4.10)

we deduce that

$$\begin{aligned} J_*(\nu )\ge \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi + \int _{{\mathbb {R}}^n} |x|^2 \,d\nu (x) \end{aligned}$$
(4.11)

for every \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\) with compact support.

On the other hand, to prove the opposite inequality in (4.11), let us first consider \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\cap C^\infty _c({\mathbb {R}}^n)\). By (4.7) we have that

$$\begin{aligned} 0\le \widehat{W}_{\alpha _k}(\xi ) |{{\widehat{\nu }}}(\xi )|^2\le C\widehat{W}_0(\xi ) |{{\widehat{\nu }}}(\xi )|^2, \end{aligned}$$

which gives a domination in \(L^1({\mathbb {R}}^n)\), since \({{\widehat{\nu }}}\in {\mathcal {S}}\) and \(\widehat{W}_0\in L^1_{\text {loc}}({\mathbb {R}}^n)\) behaves as \(1/|\xi |^2\) at infinity. Therefore, by (4.6) and the Dominated Convergence Theorem

$$\begin{aligned} \lim _{k\rightarrow \infty } \int _{{\mathbb {R}}^n} \widehat{W}_{\alpha _k}(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi = \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi \end{aligned}$$

for every \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\cap C^\infty _c({\mathbb {R}}^n)\). By (4.10) and (2.4) this implies that

$$\begin{aligned} J_*(\nu )\le \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi + \int _{{\mathbb {R}}^n} |x|^2 \,d\nu (x) \end{aligned}$$

for every \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\cap C^\infty _c({\mathbb {R}}^n)\). Let now \(\nu \in {\mathcal {P}}({\mathbb {R}}^n)\) be a measure with compact support and let \(\nu _\varepsilon \) be defined as in (2.6). Using the inequality above and the lower semicontinuity of \(J_*\) with respect to narrow convergence, we obtain that

$$\begin{aligned} J_*(\nu )\le \liminf _{\varepsilon \rightarrow 0} J_*(\nu _\varepsilon ) \le \liminf _{\varepsilon \rightarrow 0} \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |\widehat{\nu _\varepsilon }(\xi )|^2\, d\xi + \int _{{\mathbb {R}}^n} |x|^2 \,d\nu (x). \end{aligned}$$
(4.12)

Arguing exactly as in (2.9), we have that

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0} \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |\widehat{\nu _\varepsilon }(\xi )|^2\, d\xi = \int _{{\mathbb {R}}^n} \widehat{W}_*(\xi ) |{{\widehat{\nu }}}(\xi )|^2\, d\xi , \end{aligned}$$
(4.13)

even if the right-hand side is infinite. Combining (4.11)–(4.13), we finally obtain (4.1).

Step 4: Uniqueness of the minimiser of \(J_*\). We now use the representation (4.1) to show that the minimiser of \(J_*\) is in fact unique.

First of all, we prove that any minimiser of \(J_*\) must have compact support. Indeed, assume by contradiction that \(\nu _*\) is a minimiser of \(J_*\) not compactly supported. Let \((\alpha _k)\subset (-1,n-2]\) be a sequence converging to \(-1\) and let \((\nu _k)\subset {\mathcal {P}}({\mathbb {R}}^n)\) be a recovery sequence for \(\nu _*\), that is, such that \((\nu _k)\) converges to \(\nu _*\) narrowly and

$$\begin{aligned} \lim _{k\rightarrow \infty } I_{\alpha _k}(\nu _k)=J_*(\nu _*). \end{aligned}$$

In particular, \(\sup _k I_{\alpha _k}(\nu _k)<+\infty \). We argue in a similar way as in (4.9). By (2.2) there exists a compact set \(K\subset {\mathbb {R}}^n\) such that \(0<\nu _*(K)<1\), \(\nu _k(K)>0\) and

$$\begin{aligned} W_{\alpha _k}(x-y) + \frac{1}{2}(|x|^2+|y|^2)\ge I_{\alpha _k}(\nu _k) +1 \qquad \text { for} \,(x,y)\not \in K\times K \end{aligned}$$

for every k. If we define \(\mu _k\) as in (4.8), then

$$\begin{aligned} I_{\alpha _k}(\mu _k)\le I_{\alpha _k}(\nu _k) -\frac{1- \big (\nu _k(K)\big )^2}{\big (\nu _k(K)\big )^2} \end{aligned}$$

for every k. Since by narrow convergence \(\nu _k(K)\rightarrow \nu _*(K)\), as \(k\rightarrow \infty \), we obtain

$$\begin{aligned} \liminf _{k\rightarrow \infty } I_{\alpha _k}(\mu _k)\le J_*(\nu _*) -\frac{1- \big (\nu (K)\big )^2}{\big (\nu (K)\big )^2} < J_*(\nu _*). \end{aligned}$$
(4.14)

On the other hand, by the minimality of \(\mu _{\alpha _k}\) and by (4.5)

$$\begin{aligned} \liminf _{k\rightarrow \infty } I_{\alpha _k}(\mu _k) \ge \liminf _{k\rightarrow \infty } I_{\alpha _k}(\mu _{\alpha _k})= J_*(\mu _*). \end{aligned}$$
(4.15)

Since both \(\mu _*\) and \(\nu _*\) are minimisers of \(J_*\), (4.14) and (4.15) give a contradiction. On measures with compact support the representation (4.1) holds and the right-hand side of (4.1) is strictly convex as a function of \(\mu \). We, thus, conclude that \(\mu _*\) is the only minimiser of \(J_*\). \(\square \)