Discrete minimizers of the interaction energy in collective behavior: a brief numerical and analytic review

Abstract

We consider minimizers of the $N$-particle interaction potential energy and briefly review numerical methods used to calculate them. We consider simple pair potentials which are repulsive at short distances and attractive at long distances, focusing on examples which are sums of two powers. The range of powers we look at includes the well-known case of the Lennard-Jones potential, but we are also interested in less singular potentials which are relevant in collective behavior models. We report on results using the software GMIN developed by Wales and collaborators for problems in chemistry. For all cases, this algorithm gives good candidates for the minimizers for relatively low values of the particle number $N$. This is well-known for potentials similar to Lennard-Jones, but not for the range which is of interest in collective behavior. Standard minimization procedures have been used in the literature in this range, but they are likely to yield stationary states which are not minimizers. We illustrate numerically some properties of the minimizers in 2D, such as lattice structure, Wulff shapes, and the continuous large-$N$ limit for locally integrable (that is, less singular) potentials.

1 Introduction

Finding configurations with minimum potential energy for a set of particles is a common mathematical problem found across several fields. In its simplest version we consider a potential $V\colon\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ and a number $N\geq 2$ of particles. The problem then consists in finding the configurations $X=(x_{1},\dots,x_{N})\in\mathbb{R}^{Nd}$ such that the energy

$$E_{N}(X)=\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}V(x_{i}-x_{j})$$ (1.1)

is as small as possible. These configurations are called minimizers or ground states. In a broader sense, the problem also involves understanding the structure of these minimizers: their shape, the distribution of their points, and their asymptotic properties for large $N$. The energy (1.1) is clearly invariant by a translation of the points $(x_{1},\dots,x_{N})$, and consequently any translation of a minimizer is also a minimizer. One usually removes this invariance by looking for minimizers with some fixed property, most often by fixing the mean position $\sum_{i=1}^{N}x_{i}$ to be $0$. If $V$ is radially symmetric, then $E_{N}$ is also invariant by rotations of the points $(x_{1},\dots,x_{N}$, and by their inversion $(-x_{1},\dots,-x_{N})$ (and hence by any rigid motion of these points). In all cases, minimizers are the same if $V$ is substituted by $V+C$ for some constant $C$, so some normalization on $V$ is usually assumed (often, either $\lim_{|x|\to+\infty}V(x)=0$ or $V(0)=0$, but not necessarily).

In this paper we give a short review of numerical methods to calculate these minimizers, and report on some results in dimension $d=2$ using the publicly available GMIN software developed by David J. Wales (see Wales). We carry out numerical simulations in 2D for potentials of the form

$$V(x)=\frac{|x|^{a}}{a}-\frac{|x|^{b}}{b},$$ (1.2)

where $a>b$. Due to our choice of coefficients, the term with the $a$ power always represents an attractive potential, and the term with the $b$ power is repulsive. We will also give a short review of the available rigorous results to date.

In the previous statement of the problem, the particles are all identical: they interact with others in an identical way through the pair potential $V$. More complicated versions of the problem may allow for particles of different types, interacting in different ways with one another. One of the better known contexts where this problem is relevant is chemical physics, where one often wants to find the configuration of minimum possible energy for a given molecule. This is important for example in the study of protein folding, where the natural structure of a protein may be related to its minimum-energy configuration. In the context of molecular structure, the Lennard-Jones interaction potential

$$V(x)=\frac{A}{|x|^{12}}-\frac{B}{|x|^{6}},$$ (1.3)

with $A,B>0$ given constants, is often taken as a test case. The study of its minimizers (in dimension $d=3$) has received a lot of attention, and is a notable example of a hard global optimization problem.

In addition to these problems in physics, this paper is motivated by several collective behavior models where these minimizers are relevant. Among them, perhaps the simplest is the so-called aggregation equation,

$$\partial_{t}u=\operatorname{div}(u(\nabla u*V)),$$ (1.4)

which is a partial differential equation for $u=u(t,x)$, which represents the density of a certain population which interacts through the potential $V$. Equilibria of (1.4) (if any) must satisfy $\nabla u*V=0$ on the support of $u$, which is a property shared by the minimizers of the continuous interaction energy

$$E(u):=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}u(x)u(y)V(x-y)\,\mathrm{d}x\,\mathrm{d}y$$ (1.5)

in the set of probability measures (or positive measures with a fixed total mass). Hence probability minimizers of the energy (1.5) are particular examples of stationary solutions to (1.4). To avoid repetition, whenever we speak of minimizers of $E$ it is understood that we always refer to minimizers in the set of probability measures. Finding minimizers of $E$ is a “continuous version” of the problem of minimizing $E_{N}$. The link between these problems is interesting, especially if we consider that numerical studies of the continuous problem must inevitably carry out some sort of discretization of the equation, which often leads to some version of the discrete minimization problem of $E_{N}$. Hence a proper understanding of the large-$N$ behavior of minimizers of $E_{N}$ is also useful to justify numerical methods for calculating minimizers of the continuous energy $E$.

There are also many individual-based models whose “organized states” share similarities with minimizers of $E_{N}$. As an example we mention a model from D’Orsogna et al. 2006: We consider $N$ individuals with positions and velocities given by $(x_{i},v_{i})$ for $i=1,\dots,N$, and the system of ordinary differential equations

$$\frac{\mathrm{d}}{\mathrm{d}t}x_{i} =v_{i},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}v_{i} =(\alpha-\beta|v_{i}|^{2})v_{i}-\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}\nabla V(x_{i}-x_{j}).$$

This models a set of individuals which interact through a potential $V$, but also have a preferred movement speed $|v|=(\alpha/\beta)^{1/2}$. It is easy to see that if $X=(y_{1},\dots,y_{N})$ is a critical point of the energy $E_{N}$, then $(x_{i}(t),v_{i}(t)):=(y_{i}+tv,v)$ is a solution to the previous system of equations. In particular, minimizers of $E_{N}$ yield coherent states or flocking states for this model. Minimizers (or critical points) of $E_{N}$ appear similarly in other models, and this has motivated interest in types of potentials different from the ones usually encountered in physics problems.

Existence of minimizers of $E_{N}$

For a very general family of potentials it is easy to show that minimizers must exist. If one assumes that $V$ is lower semicontinuous and radially strictly increasing whenever $|x|>R$ for some given $R>0$, then one may carry out an argument like that in Cañizo and Patacchini 2018 to show their existence: it is not hard to argue that there cannot be “gaps” of size larger than $R$ between the particles of a minimizer, which gives a simple estimate on their diameter. A different argument is given in Carrillo et al. 2014 for power-law potentials, and in Blanc and Lewin 2015 for potentials satisfying $\lim_{|x|\to+\infty}V(x)=0$.

2 Numerical methods

In general the energy $E_{N}$ from (1.1), as a function $E_{N}\colon\mathbb{R}^{Nd}\to\mathbb{R}\cup\{+\infty\}$, is not convex, and the global optimization problem of finding its minima cannot use the comparatively simpler methods of convex optimization. As mentioned before, there’s the problem of invariance by translations, easily avoided by fixing the mean position of particles. The problem of invariance by rotations is harder to avoid, and hence algorithms usually find one of any possible rotations of the minimizer. It can then be hard to decide whether two runs of the same algorithm have given a rotation of the same minimizer, or a fundamentally different minimizer. However, the deeper problem is that there are usually a huge number of stationary states for the energy $E_{N}$; that is, points $(x_{1},\dots,x_{N})$ at which

$$\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}\nabla V(x_{i}-x_{j})=0\qquad\text{for all $i=1,\dots,N$.}$$ (2.1)

Many of these are local minima, and runs of standard local minimization algorithms will most often fall into one of these states, and not get to a minimizer. Hence simple descent algorithms such as variants of conjugate gradient descent, or Newton-type methods like LBFGS (Liu and Nocedal 1989), will almost surely not find minimizers even for small $N$. Global minimization algorithms must somehow avoid these “traps” by more clever search techniques. We will give a short overview of numerical methods used for this kind of problems, and then describe the algorithm we use, implemented in the program GMIN by D. J. Wales (Wales).

Most numerical methods found in the literature have been tested with the Lennard-Jones potential (1.3), and most often in 3D, which is a challenging problem of physical interest and has been the subject of many publications in chemistry journals. A database with all currently available minima of $E_{N}$ for several potentials in 3D (and molecular minimization problems of other kinds) is available online (Wales et al.), among them Lennard-Jones and Morse potentials. This database contains references of the first publication for each minimizer, and the reader may also check the unpublished paper Zhang 2010 for a good list of historical references.

Figure 1: Minimizer for the Lennard-Jones potential with 91 particles.

Early methods such as those in Hoare 1979; Hoare and Pal 1972; Hoare and Pal 1971b; Hoare and Pal 1971a are biased methods, in the sense that they work by “guessing” that minimizers have some kind of regular crystal structure. Initial guesses for a minimizer are obtained by placing particles at the sites of an appropriate lattice, and then running local minimization algorithms to optimize the positions. Then one carries out some kind of random perturbation of these positions in order to find “nearby” local minima which may have a lower energy. In addition to guessing initial candidates by using the sites of a lattice, one may also start by finding minimizers for small $N$, and add particles to build candidates for larger $N$. These methods were quite successful in finding many candidates for minimizers which are still standing today, but there are several “hard” values of $N$ for which better (lower energy) configurations were found later (Northby 1987; Wales and Doye 1997). An important step of many of these methods, used for the “random perturbation” we mentioned, is simulated annealing (Kirkpatrick et al. 1983), a global optimization method which has an appealing physical motivation for chemistry problems. Roughly speaking, in simulated annealing a random perturbation $X^{\prime}$ of a given configuration $X$ is considered, and then:

1. It is accepted if $X^{\prime}$ has lower energy than $X$, or

2. It is accepted with a probability proportional to $\exp(-\frac{E_{N}(X^{\prime})-E_{N}(X)}{T})$ if $X^{\prime}$ has higher energy than $X$.

Here $T$ is a parameter which has an analogy to temperature, and governs how probable it is to jump “far” from the current state. In the literature on this, performing this algorithm any number of times for a fixed $T$ is sometimes referred to as a Monte Carlo simulation, and each step as a Monte Carlo step. Simulated annealing methods then work by running many steps of this random perturbation, and gradually reducing $T$ until departures from the current state become very unlikely. Methods focused on this technique have had some success (Wille 1987b; Wille 1987a), and are also part of more recent methods.

Another family of methods is that of evolutionary algorithms, which have also been quite successful in finding minima for Lennard-Jones potentials, even for the values of $N$ considered harder (Deaven et al. 1996; Niesse and Mayne 1996). Genetic algorithms are a particular kind of evolutionary algorithms, though the naming is not completely precise. These methods are based on an analogy with natural selection, and are also an example of biased methods which make informed guesses on the structure of minimizers. We refer to the more recent paper by Dittner and Hartke 2017 for a list of references on the use of this kind of algorithms for Lennard-Jones energy minimization.

The paper Ben Haj Yedder et al. 2003 contains a good presentation of related minimization problems, and describes a simple algorithm for minimizing $E_{N}$ for the Lennard-Jones potential, which is tested in 2D. It is also a biased method in the sense that particles are initially placed at the sites of a triangular lattice, and later suitable symmetrization methods are used.

As discussed, these biased methods have been tried mainly for potentials of Lennard-Jones type, motivated by problems in chemistry and crystallography. Since our aim is to use a good numerical algorithm for a wider range of potentials whose minimizers may behave very differently, an adaptation of these methods does not seem straightforward, especially for the cases which have a continuous limit such as those discussed in Section 4, where particles do not arrange as large-scale crystals. Instead, we use a technique known as basin-hopping, first reported for the Lennard-Jones potential in Wales and Doye 1997, and before used by Li and Scheraga 1987 for a different problem in chemistry. This technique has been in particular quite successful for Lennard-Jones type potentials. It is an unbiased algorithm in the sense that it does not assume any particular structure of the minimizers, and it seems a more appropriate algorithm for the range of potentials we are interested in. It works as follows: first, the energy $E_{N}(X)$ is substituted by an energy $\widetilde{E}_{N}(X)$, obtained by performing a suitable descent algorithm starting from the configuration $X$ (the Polak-Ribiere variant of the conjugate gradient algorithm is mentioned in Wales and Doye 1997, although it currently uses a customized LBFGS method). That is: if we denote by $Y=M(X)$ the result of running a descent algorithm from the starting configuration $X$, then

$$\widetilde{E}_{N}(X):=E_{N}(M(X)).$$

Hence the new energy landscape $\widetilde{E}_{N}$ is a piecewise constant function, which is constant in a region where our chosen descent algorithm stops at the same local minimum $Y$, and whose value in this region is $E_{N}(Y)$. This new landscape can be studied theoretically, taking $M(X)$ to be the point to which the descent algorithm converges, and it can be approximated numerically in practice by considering the point at which the numerical descent stops. Now several steps are followed:

1. A Monte Carlo exploration of the landscape $\widetilde{E}_{N}$ is carried out with a fixed temperature, chosen so that the acceptance rate of new points is about $1/2$ (as explained in our brief earlier discussion of simulated annealing).

2. During this exploration, we also observe the individual energy (or just energy) of each particle $i$, defined by

   $$E_{N,i}=\sum_{j\neq i}V(x_{i}-x_{j}).$$

   We consider the particle $i_{max}$ which has a maximum energy among all $N$ particles, and the one $i_{min}$ which has minimum energy. If this maximum energy is larger than a certain multiple of the lowest energy, then an angular move is performed for the particle $i_{max}$: this particle is removed from the cluster, and placed at a random position close to the boundary of the cluster.

In addition to this, a further Monte Carlo exploration is carried out with the following special starting positions:

1. The best configuration found by the algorithm with $N-1$ particles, adding one particle at the best position found by a fixed number of angular moves.

2. The best configuration found by the algorithm with $N+1$ particles, removing the particle with highest energy.

We have applied this algorithm for different potentials, otherwise with no modifications. Some of the results are explored in the next sections, together with a short review of the expected behavior of minimizers. We would like to emphasize that while obtaining local minimizers is straightforward in almost all cases, obtaining good candidates for a minimizer is a hard numerical problem. As an example, the hexagonal-shaped minimizers for very singular potentials such as Lennard-Jones (see Fig 1) require some kind of global search algorithm as described above. Or consider for example a minimizer for the locally integrable potential (1.2) with $a=2$ and $b=0.5$ (Figure 2). Asymptotically as $N\to+\infty$, this density is known to converge to the density of a minimizer of the continuous energy (1.5) (see Section 4), but the way this happens shows several interesting phenomena. One sees that locally particles still seem to be arranged in triangular patterns; but since the overall shape must be round (since the support of the minimizer is strongly suspected to be a ball), these local triangular arrangements must undergo some kind of dislocation patterns in order to fit the global shape. In Figure 2 one can clearly appreciate dislocation curves along which triangular patterns with different orientations meet. These patterns are quite sensitive to the minimization algorithm used, and we do not know whether the plot in Figure 2 actually represents a minimizer. But an understanding of these patterns is certainly an interesting problem, and even its numerical study requires powerful algorithms which can handle a large number of particles and yield good approximations to the minimizer. For particle numbers over a few thousand it is likely that this is still out of reach, even numerically. In Figures 2 and 3 we show minimizer candidates for 10,000 particles, which have required large computational power and several weeks of computer time.

Figure 2: Our best numerical candidate for a minimizer of the potential (1.2) with $a=2$, $b=0.5$ and $10^{4}$ particles. (We thank David J. Wales for providing this result)

3 Crystallization phenomena

The wider problem of understanding the reason why matter is often arranged in crystalline patterns is known as the problem of crystallization. The reviews by Ventevogel 1978a; Radin 1987 are a good introduction to this, giving a motivation and related mathematical problems. It happens that minimizers of $E_{N}$ have a close relationship to the observed arrangement of atoms or molecules at low temperature, and hence an aspect of the crystallization problem which is particularly simple to state is: prove rigorously that (or whether) minimizers of $E_{N}$ approach a periodic arrangement as $N\to+\infty$. For some potentials $V$ which have a singularity stronger than $|x|^{-d}$ at $x=0$, minimizers of $E_{N}$ numerically seem to do exactly this, but proofs are surprisingly difficult even in dimension $1$! The review by Blanc and Lewin 2015 is a nice source of references on this problem, and we would like to give a short summary with our current understanding.

Let us first state the problem a bit more precisely, essentially following Blanc and Lewin 2015. In the context of crystallization, we usually consider potentials $V=V(x)$ which are positive at $x=0$ (or have a singularity where the potential diverges to $+\infty$ as $x\to 0$), have a strictly negative minimum value, and satisfy

$$\lim_{|x|\to+\infty}V(x)=0.$$ (3.1)

Further conditions are needed for minimizers to show any kind of crystallization behavior; see more on this below. The canonical example of the potentials which exhibit this kind of phenomena is the Lennard-Jones potential from equation (1.3).

For a particle distribution $X=(x_{1},\dots,x_{N})\in\mathbb{R}^{Nd}$ we consider the associated empirical measure given by

$$\mu_{X}:=\sum_{i=1}^{N}\delta_{x_{i}}.$$ (3.2)

Notice that we are not dividing by $N$, so this measure has total mass equal to $N$. In many cases we expect this measure to approach some periodic arrangement as $N\to+\infty$, and in order to define this precisely we recall that a lattice in $\mathbb{R}^{d}$ is a subgroup of $\mathbb{R}^{d}$ of the form

$$L=\{k_{1}v_{1}+k_{2}v_{2}+\dots+k_{d}v_{d}\mid k_{1},\dots,k_{d}\in\mathbb{Z}\},$$ (3.3)

where $(v_{1},\dots,v_{d})$ form a vector space basis of $\mathbb{R}^{d}$. Lattices are a simple form of periodic structures, and are a basic concept in crystallography, where they are often referred to as Bravais lattices. The reason for this special naming is that in crystallography the name lattice is often broadened to include periodic structures which are not lattices in the sense of the above definition; notable examples are the honeycomb arrangement in 2D, or the hexagonal close packed arrangement in 3D. Here we will always use the word lattice in the mathematical sense stated in (3.3). Consistently with (3.2), for a lattice $L$ we denote by $\mu_{L}$ the measure

$$\mu_{L}:=\sum_{p\in L}\delta_{p}.$$ (3.4)

Notice that $\mu_{L}$ has infinite mass (is an unbounded measure), but it is finite in compact sets. Here is a very strong formulation of the problem of crystallization:

Q1 (lattice crystallization)

Does there exist a lattice $L\subseteq\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$\mu_{X_{N}}\overset{*}{\rightharpoonup}\mu_{L}\qquad\text{as $N\to+\infty$?}$$

The convergence here is the weak-$*$ convergence of measures (convergence when tested against continuous, compactly supported functions). Or even stronger: does this take place for all sequences of minimizers, up to rigid movements of these minimizers?

Crystalline structures may well be periodic structures which are not lattices. A very natural formalization of a periodic structure is a measure $\mu$ on $\mathbb{R}^{d}$ such that $\mu(\cdot+p)=\mu(\cdot)$ for all $p$ in some lattice $L$ (that is, a measure which is invariant by translations of vector $p$, for all $p$ in a given lattice). If this happens, we say that $\mu$ is invariant by the lattice $L$. Following this, here is a weaker form of the above problem:

Q2 (periodic crystallization)

Does there exist a lattice $L\subseteq\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$\mu_{X_{N}}\overset{*}{\rightharpoonup}\mu\qquad\text{as $N\to+\infty$,}$$

where $\mu$ is a measure invariant by the lattice $L$? Or even stronger: does this take place for all sequences of minimizers, up to rigid movements of these minimizers?

If these questions are too hard, one may consider the large-$N$ behavior of some macroscopic quantities like the energy per particle. Given a lattice $L$ in $\mathbb{R}^{d}$ and a potential $V$ we define its energy per particle by

$$e(L):=\sum_{p\in L\setminus\{0\}}V(p),$$ (3.5)

when this sum converges absolutely. Also, given a configuration $X=(x_{1},\dots,x_{N})\in\mathbb{R}^{Nd}$ we define its energy per particle by

$$e_{N}(X):=\frac{1}{N}E_{N}(X).$$

Here is a related, usually weaker, statement of the problems Q1 and Q2:

Q3 (energy crystallization)

Does there exist a lattice $L\subseteq\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$e_{N}(X_{N})\to e(L)\qquad\text{as $N\to+\infty$?}$$

Or even stronger: does this take place for all sequences of minimizers?

In a related problem, one can try to give good estimates for the convergence of $e_{N}(X_{N})$ to $e(L)$, for example giving bounds on its rate of convergence. The fact that $e_{N}(X_{N})$ has some limit is known for a wide range of potentials.

In the context of crystallization phenomena one always assumes (3.1). As a consequence, a sequence of configurations with $N$ particles such that the interparticle distances all diverge to $+\infty$ has energy as close to $0$ as we wish. Hence for a minimizer $X_{N}$ it must happen that $E_{N}(X_{N})\leq 0$, so if Q3 holds one must have $e(L)\leq 0$. In any case, if Q3 holds for the full sequence of minimizers (and not up to a subsequence), in particular it implies that for some $C>0$

$$E_{N}(X_{N})\geq-CN,\quad\text{for all}\;N\geq 2.$$ (3.6)

Potentials $V$ satisfying (3.6) are called H-stable, or sometimes just stable in the literature; see for example Ruelle 1969. Often (3.6), though seemingly weaker, actually implies that $e_{N}(X_{N})$ does have a limit as $N\to+\infty$. Since the terminology comes from statistical mechanics and thermodynamics, interest is placed on potentials for which thermodynamical quantities (such as the energy-per-particle) have a finite limit, which justifies the word stable. Minimizers for potentials which are not stable may still have different scaling limits which are also interesting, though their study has different motivations and is more recent; see Section 4 below.

Finally, there is the related question of the formation of a macroscopic object. The diameter of a minimizer is defined as the maximum distance between two of its points. If we rescale minimizers in order to fix their diameter, does the resulting measure converge to something? If we look at the regime where one expects that particles approach a periodic distribution, we would expect the diameter of a minimizer to be roughly proportional to $N^{1/d}$, and hence we pose the following question:

MO (formation of a macroscopic object)

Does there exist a nontrivial measure $\mu$ on $\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$\frac{1}{N}\mu_{X_{N}/N^{1/d}}\overset{*}{\rightharpoonup}\mu\qquad\text{as $N\to+\infty$?}$$ (3.7)

Or even stronger: does this take place for all sequences of minimizers, up to rigid movements of these minimizers?

The empirical measure $\mu_{X_{N}/N^{1/d}}$ corresponds to the scaled minimizer

$$\frac{1}{N^{1/d}}X_{N}:=\left(\frac{x_{1}}{N^{1/d}},\dots,\frac{x_{N}}{N^{1/d}}\right),$$

and we divide the measure by $N$ to ensure a total mass equal to $1$. This question concerns behavior at a different scale from questions Q1–Q3. An answer to MO seems to not give much information on any of Q1–Q3, and reciprocally.

We have focused on the above questions in order to narrow the discussion to a representative case, since the field is too broad for a short review. However, there are many related aspects of these problems. Some of them involve infinite distributions of particles. Among these we mention the following:

1. (Lattice minimization) Among all lattices in $\mathbb{R}^{d}$, which ones minimize the energy per particle (3.5)?

2. (Infinite particle distributions) Given a countably infinite set of points in $\mathbb{R}^{d}$, it is possible to define their mean energy per particle in a reasonable way, so that it is consistent with the energy per particle of a lattice given in (3.3). Among all infinite particle distributions, which ones minimize the energy per particle?

We do not elaborate on different mathematical interpretations of the problem, namely: minimization problems with fixed density; modified models with a background energy; statistical mechanical or quantum mechanical models. We refer the reader to the reviews mentioned above for more on these matters (Radin 1987; Blanc and Lewin 2015). For minimization on compact submanifolds of $\mathbb{R}^{d}$ we refer to Hardin and Saff 2005; Borodachov et al. 2019.

Minimizers for a specific case of potential $V$ have links to sphere packing problems. If we consider the following sticky hard sphere potential

$$V(x):=\begin{cases}+\infty\qquad&\text{if $|x|<1$,}\\ -1\qquad&\text{if $|x|=1$,}\\ 0\qquad&\text{if $|x|>1$,}\end{cases}$$ (3.8)

then minimizers cannot have any two particles at a distance closer than $1$ (otherwise their energy would be $+\infty$), and there is a reward for spheres which are touching. Hence a minimizer of this potential is just a configuration of $N$ hard spheres in $d$-dimensional space with the maximum total number of contact points among them. In dimension $2$ this problem was solved by Heitmann and Radin 1980; see also Luca and Friesecke 2018 for a more recent proof. But the problem is far from understood in dimensions $3$ and higher, with open questions even for low values of $N$ (above $11$). A popular review of the problem can be found in Hayes 2012, and see Meng et al. 2010; Hoy et al. 2012 for some recent results. Attempts have been made to relate the structure of these minimizers to those for the Lennard-Jones potential (Trombach et al. 2018), but a precise link remains unclear. It is also tempting to make a connection of the sticky hard sphere problem with that of sphere packing, which is the problem of finding a distribution of identical hard spheres in $\mathbb{R}^{d}$ with maximum possible density. This link has often been remarked, but as far as we know there are no rigorous results backing this up. There have been some surprising recent results in sphere packing by Viazovska 2017; Cohn et al. 2017, who give specific lattices in dimensions $8$ and $24$ which can be proved to be the densest possible sphere distributions in these dimensions. The proof relies on the existence of specially symmetric lattices in these dimensions, which are natural candidates for this densest packing.

Now that we have made an effort to state rigorously the results we are interested in, we can try to summarize the available rigorous results:

One-dimensional results

In dimension $1$ there were several early results by Ventevogel 1978b; Ventevogel and Nijboer 1979b; Ventevogel and Nijboer 1979a. These results concern mainly infinite particle distributions (as mentioned in points 1,2 above) but specifically Ventevogel 1978b contains a proof of the strong version of Q3 for a wide class of power-law potentials of the form (1.2) with $b<1$ and a certain restriction on $a$: the energy per particle of any sequence of minimizers converges to the energy per particle of a given lattice; this lattice is the unique one with minimum energy per particle. This is a rather complete answer to Q3 for some power-law potentials. However, as far as we know there are no satisfying general conditions on $V$ which ensure a similar answer.

In contrast, problem Q1 is not so well understood even in dimension $1$. The best available result seems to be the one by Gardner and Radin 1979, which does show the strong version of Q1 for some potentials including the Lennard-Jones one. The conditions for this to hold are more restrictive, and it seems a hard problem to generalize them to a wide class of potentials.

In dimension $1$ the minimization problem for the sticky hard sphere potential (3.8) is trivial, and minimizers are distributions of $N$ equally spaced particles at interparticle distance $1$.

Two-dimensional results

For the sticky hard sphere potential (3.8) in dimension $2$ results are already interesting. A first proof that minimizers must be placed at the sites of a triangular lattice was given by Heitmann and Radin 1980, giving an answer to Q1 (and hence Q2) and Q3. For a specific piecewise linear potential not too far from (3.8), a positive answer to Q3 was also given by Radin 1981. Then MO was answered by Au Yeung et al. 2011 for a family of short-range potentials with $V$ compactly supported and close enough to the sticky sphere potential (3.8): the convergence (3.7) takes place up to a subsequence, and the measure $\mu$ is $\frac{2}{\sqrt{3}}\mathbbm{1}_{E}$, where $\mathbbm{1}_{E}$ is the characteristic of a finite-perimeter set $E$ of area ${\sqrt{3}}/{2}$. In Theil 2006, Q3 is answered positively for a certain family of potentials which allow similar growth conditions as the Lennard-Jones one. In addition, it is proved that with periodic boundary conditions, the triangular lattice is a minimizer. Recently in Bétermin 2023 a computer-assisted method is used to show that the triangular lattice minimizes energy among lattices for the Lennard-Jones potential in 2D (hence giving an answer to the related lattice minimization problem we mentioned earlier).

Adding a three-body potential (and hence for an energy not strictly of the form (1.1)), and still with periodic boundary conditions, E and Li 2008 followed the ideas in Theil 2006 to show that minimizers must have particles placed at the sites of a hexagonal periodic arrangement (the vertices of a honeycomb, which are not a lattice in the sense we are using in this paper). Hence in this case, a version of Q2 holds, but not Q1. Problem Q3 has been solved by Bétermin et al. 2021 for some specific piecewise-constant potentials which give rise to a square lattice.

Three-dimensional results

From the point of view of rigorous results, the situation in three dimensions is almost completely open. As far as we know there are no results regarding any of Q1, Q2, Q3 or MO, even for the sticky hard sphere potential (3.8). The problem of optimal packing, however, was solved by (Hales 2005) with a non-conventional proof which requires several computer checks. One of the difficulties of the 3D problem is that there are many non-equivalent sphere distributions in $\mathbb{R}^{3}$ which maximize density, including non-periodic structures.

Higher dimensional results

In higher dimensions there have been important breakthroughs on optimal sphere packing, for which the optimal lattices have been recently identified in dimension $8$ (Viazovska 2017) and $24$ (Cohn et al. 2017). As a consequence of these ideas, in Cohn et al. 2022 it was proved that these optimal lattices are also universally optimal: they minimize the interaction energy among infinite configurations of points with fixed density $1$, for any potentials $V$ which are completely monotone functions of the squared distance. It is conjectured that the same holds in dimension $2$ for the usual triangular lattice, but this is not proved as far as we know. In Petrache and Serfaty 2020 it is shown that this universal optimality also implies a corresponding result for Riesz and Coulomb interaction energies, which require an appropriate renormalized definition for infinite particle distributions. Hence this latter result is known in dimensions $8$ and $24$, and is a conjecture in dimension $2$.

All of these results concern infinite particle configurations. As far as we know, no results are known concerning any of Q1, Q2, Q3 or MO in dimensions $d\geq 3$.

4 Continuous limits

Consider a sequence $(X_{N})_{N\geq 2}$ of minimizers of $E_{N}$, and a locally integrable potential $V$ (in the case of potentials of the form (1.2) this holds if and only if $b>-d$). Then it may happen that

$$\nu_{X_{N}}:=\frac{1}{N}\mu_{X_{N}}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}^{N}}\overset{*}{\rightharpoonup}\nu\qquad\text{as $N\to+\infty$},$$

for some measure $\nu$, where $X_{N}=(x_{1}^{N},\dots,x_{N}^{N})$. This type of limit is of mean-field type due to the factor $1/N$, and notice that there is no rescaling as there was in (3.7). In this case the discrete interaction energy $E_{N}$ is an approximation of the continuous one $E$,

$$\frac{1}{N^{2}}E_{N}(X)=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}V(x_{i}-x_{j})=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}V(x-y)\,\mathrm{d}\nu_{X_{N}}(x)\,\mathrm{d}\nu_{X_{N}}(y)\\ \approx\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}V(x-y)\,\mathrm{d}\nu(x)\,\mathrm{d}\nu(y)=E(\nu),$$ (4.1)

when this can be rigorously justified. Then we expect the measure $\nu$ to be a minimizer of the continuous interaction energy (1.5). Proving this is another interesting question which we will state explicitly:

MF (Mean-field limit)

Does there exist a nontrivial measure $\nu$ in $\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$\frac{1}{N}\mu_{X_{N}}\overset{*}{\rightharpoonup}\nu\qquad\text{as $N\to+\infty$?}$$

Does this take place for all sequences of minimizers, up to rigid movements of these minimizers? In this case, is $\nu$ a minimizer of the continuous interaction energy (1.5)?

In analogy to Q3, one may also look for good asymptotics of the energy-per-particle $E_{N}(X_{N})/N$ as $N\to+\infty$. For locally integrable potentials, if MF holds we expect (4.1) to hold and hence we also expect that

$$E_{N}(X_{N})\sim N^{2}E(\nu).$$

It is interesting to study this approximation independently, and improve it to more precise approximations. This problem was studied by Petrache and Serfaty 2017 when the potential $V$ is repulsive and there is an additional external confining potential. For the particular case $-d<b<a=2$ in (1.2), and when continuous minimizers are regular enough, this gives quite precise asymptotics for $E_{N}(X_{N})$.

The only answer we know to question MF was given in Carrillo et al. 2014 for power-law potentials (1.2) with $b\geq 1$, and in Cañizo and Patacchini 2018 for general potentials which include the power-laws (1.2) with $b>2-d$. In this case, the answer to all questions in MF is positive for strictly unstable potentials $V$, that is, potentials $V$ for which there exists a probability distribution $\rho$ on $\mathbb{R}^{d}$ such that

$$E(\rho)<\lim_{|x|\to+\infty}V(x)$$

(assuming this limit exists or is $+\infty$). For a wide range of potentials, this concept is actually shown to imply that the classical H-stability condition (3.6) does not hold (Simione et al. 2015; Cañizo and Patacchini 2018). We believe MF should also be true in the range $-d<b\leq 2-d$, but as far as we know there is no available proof. For the strategy in Cañizo and Patacchini 2018 two ingredients are needed:

1. One needs to show that the sequence of discrete energies $E_{N}$ Gamma-converges in a suitable sense to the continuum energy $E$.

2. One needs to show that the diameter of a minimizer $X_{N}$ stays bounded as $N\to+\infty$. As a consequence, the empirical measures $\mu_{X_{N}}/N$ converge up to a subsequence, in the weak-$*$ sense, to a certain probability measure $\mu$. By using the previous Gamma-convergence result (and possibly some additional uniform estimates on other properties of the minimizers) one may then deduce that $\mu$ must be a minimizer of $E$.

In relation to these questions, a description of the minimizers of $E$ is also an interesting problem. As remarked in the introduction, these minimizers are of interest for collective behavior models, and are also quite interesting from the point of view of calculus of variations. A general result on existence of these minimizers was given in Cañizo et al. 2015; Simione et al. 2015, but there remain many questions on their uniqueness, symmetry and regularity properties. There has been a lot of recent activity regarding these, and we will try to give here a short account of the literature on this.

An early discussion on the problem of minimizing $E$ can be found in Bavaud 1991. In Raoul 2012, for very regular potentials it is proved that stationary states (and in particular minimizers, if any) must be sums of delta functions. A link between the repulsive singularity of the potential at $x=0$ and the dimension of the support of minimizers was given by Balagué et al. 2013a, essentially by noticing that for any minimizer $\rho$ the convolution $\Delta V*\rho$ must be bounded. A consequence of this is that the support of a minimizer for $V$ of the form (1.2) with $2-d\leq b<2$ must have dimension at least $2-b$. This is not optimal since numerical simulations seem to indicate that, at least for “usual” potentials, the support always has integer dimension, and some effort has been devoted to understanding this effect. See however Carrillo and Shu 2021 for an example of fractal behavior for a specifically tailored potential $V$.

For “strongly repulsive” potentials with $-d<b<2-d$ (but still locally integrable), regularity of minimizers can be obtained via obstacle problems (Carrillo et al. 2016). When $a=2$ and $-d<b<2-d$, these minimizers correspond to the “fractional Barenblatt solutions” of Caffarelli and Vázquez 2011b; Caffarelli and Vázquez 2011a; candidates for minimizers can be found for $-d<b<a$ with $b\leq 2$, when $a$ or $b$ are even integers (Carrillo and Huang 2017). In dimension $d\geq 2$, when $a=2$ and $2-d<b<4-d$, minimizers are radially symmetric, unique, and were obtained in Carrillo and Shu 2021.

For $b=2$ and $a\in(2,4)$, the unique minimizer (up to translations) is the uniform measure on an $(n-1)$-dimensional sphere of appropriate radius (Davies et al. 2021). This has recently been extended by Frank and Matzke 2023 to the case $d\geq 2$, $2\leq a\leq 4$, $b_{*}(a)\leq b\leq 2$ with $b<a$, for a certain decreasing function $b_{*}=b_{*}(a)$ which is always between $-d+3$ and $-d+4$.

This is in contrast with minimizers for mildly repulsive potentials (potentials which are $\mathcal{C}^{2}$ at $0$), for which minimizers must be supported on a finite number of points under quite general conditions (Carrillo et al. 2017). This result is true even for local minimizers in the topology of the $\infty$-Wasserstein transport distance. This applies to power-law potentials of the form (1.2) whenever $2<b<a$. For $a\geq 4$ and $b\geq 2$ (excluding the case $(a,b)=(4,2)$) the minimizer is unique up to rigid movements, and it must be a sum of equal delta functions supported at the vertices of a regular $d$-simplex (Davies et al. 2022). In this same reference, the case $(a,b)=(4,2)$ is studied. It is a special limit case in which minimizers are not unique (even up to rigid movements): minimizers may be supported on a sphere, on a finite number of points, or may be convex combinations of these possibilities. In the 1D case explicit minimizers are known for some power-law potentials of the form (1.2) (Frank 2022).

Figure 3: Numerical candidate for minimizer of the potential (1.2) with $a=2$, $b=-1$ and $10^{4}$ particles. (We thank David J. Wales for providing this result.) This corresponds to the “strongly repulsive” regime considered in Carrillo et al. 2016. In this regime it is conjectured that a sequence of minimizers $(X_{N})_{N\geq 2}$ must converge to a measure $\nu$ which must have a smooth density with respect to Lebesgue measure. In fact, the most likely candidate measure is given in Carrillo and Huang 2017 and is of the form

$$\nu(x)=C(R^{2}-|x|^{2})_{+}^{1-\frac{b+d}{2}},$$

with $(\cdot)_{+}$ denoting the positive part, with suitable constants $C,R>0$ which ensure in particular that the total mass is $1$. In the case corresponding to this figure we have $d=2$, $b=-1$, so the density of $\nu$ vanishes at the boundary. By contrast, in the case of Figure 2 the exponent $1-(b+d)/2$ becomes negative and the density of $\nu$ should diverge at the boundary.

In addition to this, the question of stability or instability of these minimizers for several dynamical models has been investigated in a few cases; see for example Fellner and Raoul 2010a; Fellner and Raoul 2010b (for 1D stability of sums of Dirac masses for the aggregation equation (1.4)); Balagué et al. 2013b; Kolokolnikov et al. 2011 (stability of spherical shells and rings); Bertozzi et al. 2012 (behavior of patch solutions for $b=2-d$).

We will also mention some results on non-isotropic potentials: see Carrillo et al. 2021 and Carrillo and Shu 2023; Carrillo and Shu 2022.

There are also related minimization problems which would take us too far from the topic; we refer the reader to Blanc and Lewin 2015 and Frank 2023 for an introduction to similar problems and many references on them.

We will illustrate the above properties in 2D with a few numerical results for power-law potentials of the form (1.2). Let us first fix $a=2$ and look at the appearance of numerical minimizers for several values of $b$. As the exponent $b$ increases, the repulsion becomes progressively weaker. It is clearly seen in Figures 4 and 5 that minimizers become more concentrated as $b$ increases, since repulsion becomes locally weaker. Also, they become more concentrated on the boundary of a ball for higher $b$. For $b<0$ we are in the regime studied in Carrillo et al. 2016, so we know minimizers of $E$ should be smooth, but it is not proved that $N$-particle minimizers approach minimizers of $E$. In contrast, for $b>0$ we know from Cañizo and Patacchini 2018 that $N$-particle minimizers must approach minimizers of the continuous energy $E$, but we do not know the specific shape of these minimizers; for example, it is not proved in this range that minimizers must be regular functions on their support, or even $L^{P}$ functions for some $p\geq 1$. Reasonable candidates for minimizers of $E$ are given in Carrillo and Huang 2017, and they are all smooth functions on their support.

Figure 4: Minimizer for the power-law potential (1.2) with $a=2$ and several values of $b$, always with $10^{3}$ particles.

Figure 5: Zoomed-in view of the last 4 minimizers in Figure 4.

For $a=4$ the situation is similar, but now there is a range of exponents with $b\geq 2$ where the potential $V$ is $\mathcal{C}^{2}$. This is the range studied in Davies et al. 2021 and Carrillo et al. 2017. In Figure 6 we observe how the increase of $b$ causes the appearance of empty space in the middle of the particle configuration, and its further increase implies the accumulation of all particles in a ring.

Figure 6: Minimizer for the power-law potential (1.2) with $a=4$ and several values of $b$, always with $10^{3}$ particles.

5 Between crystallization and continuous limits

We will also point out that there is a range of exponents for power-law potentials which has been left almost untouched in the literature. For potentials $V$ which are not locally integrable (that is, with $b\leq-d$) there are ranges of $a$ for which crystallization (in the sense of Q1 to Q3) is not expected to take place, or has not been considered. The main unexplored range is $b\leq-d$, $a>0$, since in the crystallization literature and in problems motivated by statistical mechanics it is always assumed that $V(x)$ decays to $0$ as $|x|\to+\infty$. But there is also a range of exponents $b<a<0$ for which crystallization does not seem to happen, and we are left with some kind of limiting behavior which is not yet understood. We conjecture that in this case there is still a scaling which leads to a well-defined limiting profile along the lines of MO:

MS (macroscopic scaling)

Do there exist $s>d$, a nontrivial measure $\mu$ on $\mathbb{R}^{d}$ and a sequence of minimizers $(X_{N})_{N\geq 2}$ of $E_{N}$ such that (possibly up to a subsequence)

$$\frac{1}{N}\mu_{X_{N}/N^{1/s}}\overset{*}{\rightharpoonup}\mu\qquad\text{as $N\to+\infty$?}$$ (5.1)

Or even stronger: does this take place for all sequences of minimizers, up to rigid movements of these minimizers?

Figure 7: Some minimizers with $b\leq-2$.

There are numerical results in 1D which seem to indicate that this is the case (Carrillo et al. 2014), but this has not been systematically explored, much less the subject of any rigorous proofs. Figure 7 shows 2D configurations of the kind $b<a<0$ and $b<-d$, $a>0$. In some of them, we appreciate certain similarities with crystallization cases (see Figure 1), but we are not aware of any works exploring their behavior. This regime may be related to scaling limits of dynamical problems with non-integrable potentials such as Oelschläger 1990, but the link is far from clear to us.

J. A. Cañizo. Departamento de Matemática Aplicada & IMAG, Universidad de Granada. E-mail address: canizo@ugr.es

A. Ramos-Lora. Departamento de Matemática Aplicada & IMAG, Universidad de Granada. E-mail address: ramoslora@ugr.es