# Optimality of the triangular lattice for Lennard-Jones type lattice energies: a computer-assisted method

Laurent Bétermin

<figure><img src="image_1.png" /></figure>

To cite this version:

Laurent Bétermin. Optimality of the triangular lattice for Lennard-Jones type lattice energies: a computer-assisted method. Journal of Physics A: Mathematical and Theoretical, 2023, 56 (14), pp.145204. 10.1088/1751-8121/acc21d hal-04454224

HAL Id: hal-04454224

https://hal.science/hal-04454224v1

Submitted on 13 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimality of the triangular lattice for Lennard-Jones type lattice energies: a computer-assisted method

Laurent Bétermin

Institut Camille Jordan, Université Claude Bernard Lyon 1
43 Bd du 11 novembre 1918, 69622 Villeurbanne Cedex, France
betermin@math.univ-lyon1.fr. ORCID id: 0000-0003-4070-3344

February 9, 2023

## Abstract

It is well-known that any Lennard-Jones type potential energy must have a periodic ground state given by a triangular lattice in dimension 2. In this paper, we describe a computer-assisted method that rigorously shows such global minimality result among 2-dimensional lattices once the exponents of the potential have been fixed. The method is applied to the widely used classical (12, 6) Lennard-Jones potential, which is the main result of this work. Furthermore, a new bound on the inverse density (i.e. the co-volume) for which the triangular lattice is minimal is derived, improving those found in [L. Bétermin and P. Zhang, Commun. Contemp. Math., 17 (2015), 1450049] and [L. Bétermin, SIAM J. Math. Anal., 48 (2016), 3236-3269]. The same results are also shown to hold for other exponents as additional examples and a new conjecture implying the global optimality of a triangular lattice for any parameters is stated.

**AMS Classification:** Primary 74G65, Secondary 82B20, 52C15.

**Keywords:** Lennard-Jones potential, Crystallization, Lattices, Triangular lattice, Epstein zeta function, Ground state.

## 1 Introduction and main results

The question of existence and uniqueness of periodic ground states for interacting particles systems has been addressed in many contexts (see [13]). The crucial question “Why are solids crystalline?” asked for instance by Radin [37] in the context of low-temperature matter is still under investigation and only few rigorous answers have emerged, leaving the “Crystallization Conjecture” largely open. Even when particles are interacting through a pairwise potential, finding theoretically or numerically the optimal configuration of points minimizing the potential energy of such system is extremely challenging due to the huge number of parameters and the presence of lots of critical points.

We usually expect for a pairwise interaction potential between radially symmetric neutral particles at very low temperature to be repulsive at small distances (due to Pauli Principle) and attractive at large distances (due to Van der Waals forces). A prototypical example of such interaction is the classical Lennard-Jones potential

$$f(r) = \frac{a}{r^{12}} - \frac{b}{r^6}, \quad \text{where} \quad (a, b) \in (0, \infty)^2, \qquad (1.1)$$

popularized by Lennard-Jones in [28] to study the thermodynamical properties of rare gases (e.g. liquid argon), but originally proposed by Mie in [33] in the general form $f(r) = ar^{-\alpha} - br^{-\beta}$ for $\alpha > \beta$ (we will however call such potential “Lennard-Jones type potentials”). The repulsion is therefore given by the
inverse power law $ar^{-12}$ and the behavior at infinity in $-br^{-6}$ mimics the Van der Waals attraction. This potential is widely used in Molecular Dynamics (see e.g. [24, 29]). In dimension 2, particles interacting through such potential are expected to crystallize, i.e. to find their ground state, on a perfect triangular lattice (see Figure 1 and (1.4)) as their number goes to infinity (see e.g. [12] and [13, Fig. 1]). In dimension 3, the same is expected to happen on a Hexagonal Close Packing structure (see e.g. [1, 4]). The rare attempts to show the two-dimensional case have lead to optimality results for the triangular lattice in the case of hard-sphere potentials [26, 32] and approximations of them [25, 36, 42]. Notice also that crystallization results for long-range pairwise potentials have only been proven in [42] (on a triangular lattice) and recently in [8] (on a square lattice) where the potentials are strongly repulsive at the origin and have very narrow wells (for the triangular lattice, in order to catch only the nearest-neighbors) or sufficiently large well one (for the square lattice, in order to also catch the next-nearest neighbors) before converging rapidly to zero.

<figure><img src="image_3.png" /><figcaption>Figure 1: Patch of the triangular lattice of unit density $A_2 \subset \mathbb{R}^2$ defined by (1.4) where all primitive triangles are equilateral with side length $\sqrt{2/\sqrt{3}}$.</figcaption></figure>

When crystallization is assumed, finding the most energetically favorable periodic lattice structure (see Definition 2.1) is an interesting and much simpler problem due to the small number of parameters in action. However, the apparent simplicity of the problem hides a lot of technical issues, even in dimension 2. Indeed, the usual way to minimize energies of type

$$E_f[L] := \sum'_{p \in L} f(|p|), \quad L \in \mathcal{L}_2 := \{Zu \oplus Zv \subset \mathbb{R}^2 : (u, v) \text{ is a basis of } \mathbb{R}^2\}, \quad (1.2)$$

where $|\cdot|$ is the Euclidean norm on $\mathbb{R}^2$ and $\sum'$ denotes the sum without the $p=0$ term, among lattices $L \in \mathcal{L}_2$ is to parametrize the lattices with 3 real numbers and to optimize this 3-variable function by computing its derivatives. When $f(r) = r^{-s}$ is an inverse power-law, $E_f$ is called the Epstein zeta function defined, for $s > 2$ and $L \in \mathcal{L}_2$, by

$$\zeta_L(s) := \sum'_{p \in L} \frac{1}{|p|^s}, \quad (1.3)$$

and has been studied by Rankin [38], Ennola [22], Cassels [14] and Diananda [21] (see also [27] for a review of these results). They have shown, for all fixed $s > 2$, the minimality of the triangular lattice $A_2 \subset \mathbb{R}^2$ defined by

$$A_2 := \sqrt{\frac{2}{\sqrt{3}}} \left[ \mathbb{Z}(1, 0) \oplus \mathbb{Z} \left( \frac{1}{2}, \frac{\sqrt{3}}{2} \right) \right] \quad (1.4)$$

among two-dimensional lattices of fixed unit density. This result have been generalized to all $f(r) = F(r^2)$, where $F$ is a completely monotone function (i.e. the Laplace transform of a nonnegative Radon
measure), by Montgomery in [34] where he showed the optimality of the triangular lattice for any Gaussian function $f(r) = e^{-\alpha r^2}$, $\alpha > 0$ (this is also called “universal optimality among lattices”). Concerning the Lennard-Jones type potentials $f(r) = ar^{-\alpha} - br^{-\beta}$, several results have been derived in [2, 3, 5, 9, 11], especially the optimality of the triangular lattice among lattices with fixed high density as well as the minimality of a triangular lattice for $E_f$, with small exponents, among all possible two-dimensional lattices. We have also shown in [9] that the type (or “shape”) of the global minimizer of $E_f$ is independent of $(a, b) \in (0, \infty)^2$ (see also Theorem 2.6). Unfortunately, the method of [2, 11] based on Montgomery’s result and a Riemann’s splitting of $E_f$ was not adapted to show that a triangular lattice is the minimizer of $E_f$ among all lattices when $f$ is the classical Lennard-Jones potential (1.1), but only for certain pairs of small exponents where $\alpha \le 8$ (see [2, Thm. 1.2.B.2 and Rmk. 6.18] for more details).

The present work aims to correct this issue by deriving a rigorous computer-assisted method designed to show that, given the exponents $\alpha > \beta > 2$, the minimizer of $E_f$ is a triangular lattice. This general optimality result has been conjectured in [2, 3, 9, 11] and numerically checked several times (see also [43]). Even though we have checked our proof for many different pairs of exponents (see Section 3.3 and in particular Equation (3.7)), we state our main theorem for the classical Lennard-Jones potential (1.1) since this is the most popular and interesting one among them in mathematical physics and molecular dynamics.

**Theorem 1.1 (Optimality of a triangular lattice for the classical Lennard-Jones energy).**

For any $(a, b) \in (0, \infty)^2$, the triangular lattice $\sqrt{V_{A_2}}A_2$ of co-volume $V_{A_2} = \left(\frac{2a\zeta_{A_2}(12)}{b\zeta_{A_2}(6)}\right)^{\frac{1}{3}}$ is the unique minimizer in $\mathcal{L}_2$ defined by (1.2), up to rotation, of

$$L \mapsto E_f[L] := \sum'_{p \in L} \left[ \frac{a}{|p|^{12}} - \frac{b}{|p|^6} \right] = a\zeta_L(12) - b\zeta_L(6), \quad (1.5)$$

and where the Epstein zeta function $\zeta_L(s)$ is defined by (1.3).

The proof’s strategy is precisely explained in Section 3.1 but we briefly sketch it here. This is a computer-assisted proof since it needs a computer to check a finite number of values in order to conclude. The method is also designed to work for a general pair of exponents $(\alpha, \beta)$ replacing the parameters (12, 6) in Theorem 1.1, i.e. when $f$ is a Lennard-Jones type potential $f(r) = ar^{-\alpha} - br^{-\beta}$. First, our minimization problem among all lattices is transformed into the new problem $Q_{\alpha,\beta}(L) > \alpha/\beta$, where, for all $L \neq A_2$,

$$Q_{\alpha,\beta}(L) := \frac{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)},$$

among unit density lattices. We therefore restrict our problem to a compact set $\overline{K}_{\alpha,\beta}$ (see (3.4) and Figure 2) of unit density lattices, since the above inequality appears to be true outside of it. This is hence possible to design a square grid of points $G_{\alpha,\beta}^{\delta} \subset \overline{K}_{\alpha,\beta}$ (see Figure 2) with appropriate small size $\delta$ such that the minimum value of $Q_{\alpha,\beta}$ in $G_{\alpha,\beta}^{\delta}$ gives a sufficiently precise value of the minimum of $Q_{\alpha,\beta}$ in $\overline{K}_{\alpha,\beta}$. Finally, if $Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta}$ among these values with the appropriate degree of accuracy, then Theorem 1.1 is proved for the pair of exponents $(\alpha, \beta)$. This is in particular true for $(\alpha, \beta) = (12, 6)$ as precisely stated in Theorem 1.1 (see Section 3.2) and for other exponents as shown in Section 3.3.

This type of computer-assisted proof applied to lattice energy minimization problems has been also used by Sarnak and Strombergsson in [39] for showing that the FCC lattice is the minimizer of the three-dimensional height among flat tori (i.e. basically the derivative of the Epstein zeta function at $s = 0$). However, our problem concerning $Q_{\alpha,\beta}$ is really less complicated than the latter since we are in dimension 2 (i.e. the computational time is rather short) and we do not need extremely precise estimates to conclude (the inequality $Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta}$ is actually far from being sharp, see Sections 3.2
and 3.3). For these reasons, we do not provide any code in this paper since we only need to compute a good approximation of Epstein zeta functions (i.e. with enough terms) as well as enough values on the grid $G_{\alpha, \beta}^{\delta}$ for showing our result, which are rather simple tasks. Also, our method can be applied to any chosen parameters $(\alpha, \beta)$ (see Section 3.3) as it is expected that the minimality of the triangular lattice holds for all exponents (see e.g. [9, 43]). However, our proof is specific to Lennard-Jones type potentials since it strongly relies on homogeneity of the Epstein zeta function (see Lemma 2.4) as well as on the fact that $f$ is a one-well potential. Therefore, we do not expect our method to be adaptable to other types of potentials. Furthermore, since no optimality result is available for the three-dimensional Epstein zeta function (for which the Face-Centred-Cubic lattice is expected to be a minimizer when $s > 3$) and since the Hexagonal-Close-Packing structure is not a lattice in the sense of our Definition 2.1, our method is not directly applicable in dimension $d = 3$ (see also Remark 3.5).

It is interesting to notice that Theorem 1.1 – and its analogue for any $(\alpha, \beta)$ where the computer-assisted proof is checked – directly implies optimality results of a triangular lattice (with the corresponding optimal scaling) for Embedded-Atom Models of the form

$$ \mathcal{E}(L) := F(E_{\rho}[L]) + E_{f}[L], \quad L \in \mathcal{L}_2, \quad E_{f} \text{ and } E_{\rho} \text{ defined by (1.2)}, $$

and where, for all $r > 0$, $\rho(r) = r^{-\beta}$, $f(r) \in \{r^{-\alpha}, ar^{-\alpha} - br^{-\beta}\}$ and $F(r) \in \{r^t, -\log r, -c\sqrt{r}\}$, with $(c, t) \in (0, \infty)^2$ and $\alpha > \beta > 2$, as recently shown in [6, Thm. 4.1 and 5.2]. The function $\rho$ (resp. $f$) corresponds to the electrons (resp. nuclei) energy contribution of atoms located at lattice sites. This type of semiempirical model, initially introduced by Daw and Baskes [19] and which can be viewed as an intermediate model between zero-temperature phenomenological pair-interaction energies and quantum systems, is routinely used in Molecular Simulation for Material Sciences (see e.g. [20, 31]).

Additionally, another consequence of our proof is the following new bound for the co-volume $V$ (i.e. the inverse density) for which $\sqrt{V}A_2$ is the unique minimizer of $E_f$ in $\mathcal{L}_2(V)$, i.e. among lattices with fixed co-volume $V$ (see Definition 2.1 for more details).

**Proposition 1.2 (Optimality at fixed high density of the triangular lattice).** For any $(a, b) \in (0, \infty)^2$, if $0 < V \le (\frac{2a}{b})^{\frac{1}{3}}$ then the triangular lattice $\sqrt{V}A_2$ is the unique minimizer of $E_f$ defined by (1.5) in $\mathcal{L}_2(V)$ up to rotation.

This result improves the previous bound we found in [2, 11], which was $(\frac{\pi^3 a}{60b})^{\frac{1}{3}}$. Furthermore, our proof is actually designed to show that the general bound for fixed $(\alpha, \beta)$ and any $(a, b)$ is actually $(\frac{a\alpha}{b\beta})^{\frac{2}{\alpha-\beta}}$ (see Section 3.1 and Section 3.3 for applications to other exponents), which again improves the bound of [2, Prop. 6.11] obtained by using a Riemann’s splitting argument. Notice that this bound is still not the sharpest one since we expect the triangular lattice to be minimal in $\mathcal{L}_2(V)$ for slightly larger values of $V$ (see [3, 43]). It is also important to notice that Proposition 1.2 actually implies Theorem 1.1. However, we have chosen to write the results in this order for highlighting our main global optimality result.

We would like to add to this introduction the statement of an important conjecture related to our problem, which implies Theorem 1.1 and Proposition 1.2 for any pair of exponents. Indeed, the present work shows a method to prove the optimality of the triangular lattice for any Lennard-Jones type potentials, but only once the exponents $(\alpha, \beta)$ are fixed. It would be obviously very interesting to find a more general proof since it is expected that these results hold for all exponents (see [2, 3, 9]). Theorem 1.1 and Proposition 1.2 would actually hold for any parameters $2 < \beta < \alpha$ if we were able to show the following conjecture.

**Conjecture 1.3.** For all $s > 2$, $A_2$ is the unique minimizer, among the two-dimensional lattices of unit density, of

$$ \mathcal{F}_s(L) := - \sum'_{p \in L} \frac{s \log |p| + 1}{|p|^s}. $$
Indeed, this conjecture is equivalent to the fact that $Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta}$ for all unit density lattice L and all $2 < \beta < \alpha$ which implies Proposition 1.2 and then Theorem 1.1 for all the exponents (see Section 3.3 for details). We have actually already shown this result in [2, Thm. 1.2.B.2] for $2 < s < 2(\psi^{-1}(\log \pi) - 1) \approx 5.25$ (see also [2, Rmk. 6.18] for details) but the conjecture stays totally open for larger s. Furthermore, since $r \mapsto -(s \log r + 1)r^{-s}$ is a one-well potential with minimum at $r = 1$, Conjecture 1.3 stays however difficult to solve and is therefore an interesting challenging problem.

Finally, it has to be noticed that our study's setting does not cover other important physical cases such as periodic topological defect pattern in the background of the lattice structures or various symmetry-breaking amorphous lattices. For example, dislocations may be periodically distributed over the triangular lattice. Considering the long tail of the Lennard-Jones potential, the above-mentioned two kinds of cases could not be readily excluded as the energy minima. The reader can for example refer to [35] and references therein.

**Plan of the paper.** We first recall in Section 2 the notions of lattices and energies (see Section 2.1), as well previous results on Epstein zeta function and Lennard-Jones type energies (see Section 2.2). Section 3 is devoted to the proof of Theorem 1.1 and Proposition 1.2 where the general method is explained in Section 3.1, the application to the parameters $(\alpha, \beta) = (12, 6)$ is shown in Section 3.2 and to others exponents in Section 3.3. Finally, a central technical auxiliary lemma is proved in Section 4.

# 2 Definitions and previous results on lattice energies

## 2.1 Lattices and Energies

We first define the set of periodic discrete configurations of points we are interested in, as well as their associated quadratic forms.

**Definition 2.1 (Lattices and associated quadratic form).** We call $\mathcal{L}_2$ the set of two-dimensional lattices of the form $L = \mathbb{Z}u \oplus \mathbb{Z}v$ where $\{u, v\} \subset \mathbb{R}^2$ is a basis of $\mathbb{R}^2$. Furthermore, we write $\mathcal{L}_2(V) \subset \mathcal{L}_2$ the set of lattices $L = \mathbb{Z}u \oplus \mathbb{Z}v$ with co-volume $|\det(u, v)| = V$.

Moreover, to any lattice $L = \mathbb{Z}u \oplus \mathbb{Z}v \in \mathcal{L}_2$ is associated a quadratic form $Q_L$ defined by

$$Q_L(m, n) := |mu + nv|^2, \quad \forall (m, n) \in \mathbb{Z}^2.$$

Furthermore, the theory of quadratic forms' reduction (see e.g. [41] or [34]) allows to restrict our set of lattices to a fundamental domain where each of them appears only once and can be parametrized by a point $(x, y) \in \mathbb{R}^2$.

**Proposition 2.1 (Parametrization and fundamental domain).** Any lattice $L \in \mathcal{L}_2(1)$ can be uniquely parametrized by $(x, y) \in \mathcal{D} \subset \mathbb{R}^2$ such that

$$\mathcal{D} := \{(x, y) \in \mathbb{R}^2 : y > 0, x \in [0, 1/2], x^2 + y^2 \ge 1\}. \quad (2.1)$$

The set $\mathcal{D}$, depicted in Figure 2, is called the half-fundamental domain, the parametrization can be written

$$L = \mathbb{Z}u \oplus \mathbb{Z}v, \quad \text{where} \quad u = \left(\frac{1}{\sqrt{y}}, 0\right) \quad \text{and} \quad v = \left(\frac{x}{\sqrt{y}}, \sqrt{y}\right),$$

and the quadratic form of such L is therefore

$$Q_L(m, n) = \frac{1}{y}(m + xn)^2 + yn^2, \quad \forall (m, n) \in \mathbb{Z}^2.$$
**Examples 2.2** (The triangular and square lattices). The square lattice $\mathbb{Z}^2$ and the triangular lattice $A_2$, respectively defined by

$$ \mathbb{Z}^2 = \mathbb{Z}(1,0) \oplus \mathbb{Z}(0,1) \quad \text{and} \quad A_2 := \sqrt{\frac{2}{\sqrt{3}}} \left[ \mathbb{Z}(1,0) \oplus \mathbb{Z} \left( \frac{1}{2}, \frac{\sqrt{3}}{2} \right) \right], $$

belong to $\mathcal{L}_2(1)$ and are respectively parametrized in $\mathcal{D}$ by $(0,1)$ and $(\frac{1}{2}, \frac{\sqrt{3}}{2})$.

We now define our main lattice energies, namely the Epstein zeta function and the Lennard-Jones type energy which is actually a linear combination of the first one. The Epstein zeta function was originally introduced by Epstein [23] and appears in many contexts (see e.g. [7, 27, 40]).

**Definition 2.2** (Epstein zeta function and Lennard-Jones type energy). For any $s > 2$, the Epstein zeta function $\zeta_L(s)$ of $L \in \mathcal{L}_2$ and its $N^{th}$ partial sum $\zeta_L^N(s)$, for some $N \in \mathbb{N}$, are defined by

$$ \zeta_L(s) := \sum'_{p \in L} \frac{1}{|p|^s} = \sum'_{(m,n) \in \mathbb{Z}^2} \frac{1}{Q_L(m,n)^{\frac{s}{2}}} \quad \text{and} \quad \zeta_L^N(s) := \sum'_{\substack{|m| \le N, \\ (m,n) \in \mathbb{Z}^2}} \frac{1}{Q_L(m,n)^{\frac{s}{2}}}, \quad (2.2) $$

where $\sum'$ is the sum where the origin is omitted.

Furthermore, for any $\alpha > \beta > 2$ and $(a,b) \in (0, \infty)^2$, the Lennard-Jones type potential is defined by

$$ f(r) := \frac{a}{r^{\alpha}} - \frac{b}{r^{\beta}}, $$

and the corresponding Lennard-Jones type energy of $L \in \mathcal{L}_2$ is given by

$$ E_f[L] := \sum'_{p \in L} f(|p|) = a\zeta_L(\alpha) - b\zeta_L(\beta) = \sum'_{p \in L} \left[ \frac{a}{|p|^\alpha} - \frac{b}{|p|^\beta} \right] = \sum'_{(m,n) \in \mathbb{Z}^2} \left[ \frac{a}{Q_L(m,n)^{\frac{\alpha}{2}}} - \frac{b}{Q_L(m,n)^{\frac{\beta}{2}}} \right]. $$

**Remark 2.3.** In the proof of Theorem 1.1 and Proposition 1.2, we will use a computer to check the values of a quotient $Q_{\alpha,\beta}(L)$ among unit density lattices $L$. For that, we will obviously need to work with $\zeta_L^N(s)$ defined by (2.2) for which a good approximation is known (see e.g. [18]). Since our goal is to work with exponents that are rather large (typically $(\alpha, \beta) = (12, 6)$), only few terms (i.e. small $N$, which actually depends on $y$ in the parametrization in $\mathcal{D}$) are sufficient to accurately approximate $\zeta_L(s)$ by $\zeta_L^N(s)$.

## 2.2 Results on Epstein zeta functions and Lennard-Jones type energies

In this section, we summarize important results concerning the Lennard-Jones type energies that we will need for our proof. The reader can refer to [2, 3, 9, 43] for more details, proofs and results, including the optimality of the triangular lattice for $E_f$ for small exponents shown in [2, Theorem 1.2.B.2]. We first start by a straightforward scaling formula.

**Lemma 2.4** (Scaling formulas). For any $L \in \mathcal{L}_2(1)$ and any $V > 0$, $\sqrt{V}L \in \mathcal{L}_2(V)$ and, furthermore, for any $s > 2$ and any $\alpha > \beta > 2$, we have

$$ \zeta_{\sqrt{V}L}(s) = V^{-\frac{s}{2}}\zeta_L(s), \quad \text{and} \quad E_f[\sqrt{V}L] = aV^{-\frac{\alpha}{2}}\zeta_L(\alpha) - bV^{-\frac{\beta}{2}}\zeta_L(\beta). \quad (2.3) $$

We also recall the famous optimality result found by Rankin-Ennola-Cassels-Diananda in [14, 21, 22, 38] and generalized by Montgomery in [34].

**Theorem 2.5** (Optimality of the triangular lattice for the Epstein zeta function). For any $s > 2$, $A_2$ is the unique minimizer of $L \mapsto \zeta_L(s)$ in $\mathcal{L}_2(1)$, up to rotation.
Finally, we state the variational properties of $E_f$ regarding the minimizer among the dilated versions of a lattice $L$ as well as an upper bound for the co-volume of a global minimizer of $E_f$ (see [2, 9]).

**Theorem 2.6** (Minimality properties of the Lennard-Jones type energy). Let $\alpha > \beta > 2$. For any $L \in \mathcal{L}_2(1)$, the unique minimizer of $V \mapsto E_f[\sqrt{V}L]$ is

$$V_L := \left( \frac{a\alpha\zeta_L(\alpha)}{b\beta\zeta_L(\beta)} \right)^{\frac{2}{\alpha-\beta}}, \qquad (2.4)$$

and the minimal energies among the dilated versions of $L$ is therefore

$$\min_{V>0} E_f[\sqrt{V}L] = E_f[\sqrt{V_L}L] = \frac{b^{\frac{\alpha}{\alpha-\beta}} \zeta_L(\beta)^{\frac{\alpha}{\alpha-\beta}}}{a^{\frac{\beta}{\alpha-\beta}} \zeta_L(\alpha)^{\frac{\beta}{\alpha-\beta}}} \left(\frac{\beta}{\alpha}\right)^{\frac{\beta}{\alpha-\beta}} \left(\frac{\beta}{\alpha} - 1\right) < 0.$$

In particular, the minimizer of $L \mapsto E_f[\sqrt{V_L}L]$ in $\mathcal{L}_2(1)$ only depends on $(\alpha, \beta)$.

Furthermore, if $L_0 \in \mathcal{L}_2(V_{L_0})$ is a global minimizer of $E_f$ in $\mathcal{L}_2$, then its volume $V_{L_0}$ is bounded as follows for any set of parameters $(\alpha, \beta, a, b)$:

$$0 < V_{L_0} \le \left( \frac{a\alpha}{b\beta} \right)^{\frac{2}{\alpha-\beta}}. \qquad (2.5)$$

**Remark 2.7.** In the above theorem, notice that the minimizer of $L \mapsto E_f[\sqrt{V_L}L]$ (i.e. the “shape” of the minimizer of $E_f$, see [9, Sec. 1.1]) does not depend on $(a, b) \in (0, +\infty)$ but the minimal value of the energy $E_f$ does indeed depend on these parameters.

# 3 Proof of our results and generalization

## 3.1 Description of our general computer-assisted method

We aim to present a proof of Theorem 1.1 and Proposition 1.2 as general as possible in order to be able to use our computer-assisted method for any pair of exponents $(\alpha, \beta)$. By Theorem 2.6, we know that the “shape” of the global minimizer of $E_f$, i.e. the minimizer of $L \mapsto E_f[\sqrt{V_L}L]$ in $\mathcal{L}_2(1)$ where $V_L$ is given by (2.4), is independent of $(a, b)$. Therefore, we choose $a = 1$ and $b = \frac{\alpha}{\beta}$ in such a way that the Lennard-Jones potential

$$f(r) = \frac{1}{r^{\alpha}} - \frac{\alpha}{\beta} \frac{1}{r^{\beta}}$$

reaches its minimum for $r = 1$. Furthermore, it also follows that any global minimizer of $E_f$ has its volume smaller than 1 according to (2.5) since $\frac{a\alpha}{b\beta} = 1$. Our goal is therefore to show that the triangular lattice $\sqrt{V}A_2$ is the unique minimizer of $E_f$ in $\mathcal{L}_2(V)$ for all $V \le 1$ and Theorem 1.1 will be proved. This is actually equivalent, by a simple scaling argument, with the following statement generalizing Proposition 1.2:

* For any $(\alpha, \beta, a, b) \in (0, \infty)^4$, if $V \le \left(\frac{a\alpha}{b\beta}\right)^{\frac{2}{\alpha-\beta}}$ then the triangular lattice $\sqrt{V}A_2$ is the unique minimizer of $E_f$ in $\mathcal{L}_2(V)$, up to rotation.

This automatically holds if the proof we are presenting in this section is applicable. In particular, we will show this result for $(\alpha, \beta) = (12, 6)$ in Section 3.2 and for other exponents in Section 3.3. To summarize, Theorem 1.1 is actually a corollary of Proposition 1.2 as a simple consequence of Theorem 2.6. However, we have chosen to present Theorem 1.1 as the main result of this paper.

We remark that, using the homogeneity formulas stated in (2.3), for $V > 0$, we have

$$E_f[\sqrt{V}L] - E_f[\sqrt{V}A_2] \ge 0, \forall L \in \mathcal{L}_2(1)$$
$$ \iff V^{-\frac{\alpha}{2}} (\zeta_L(\alpha) - \zeta_{A_2}(\alpha)) \geq \frac{\alpha}{\beta} V^{-\frac{\beta}{2}} (\zeta_L(\beta) - \zeta_{A_2}(\beta)), \forall L \in \mathcal{L}_2(1) $$

$$ \iff V \leq V_{\alpha, \beta} := \left(\frac{\beta}{\alpha}\right)^{\frac{2}{\alpha-\beta}} \inf_{\substack{L \in \mathcal{L}_2(1) \\ L \neq A_2}} \left(\frac{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)}\right)^{\frac{2}{\alpha-\beta}}, $$

where the inequality in the second equivalence follows from Theorem 2.5, i.e. the optimality of $A_2$ in $\mathcal{L}_2(1)$ for $L \mapsto \zeta_L(s)$ for all $s > 2$. If we show that $V_{\alpha, \beta} \ge 1$, then the theorem is proved. We therefore call, for any $L \in \mathcal{L}_2(1) \setminus \{A_2\}$,

$$ Q_{\alpha, \beta}(L) := \frac{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)}, \qquad (3.1) $$

and $V_{\alpha, \beta} \ge 1$ is equivalent with

$$ Q_{\alpha, \beta}(L) > \frac{\alpha}{\beta}, \quad \forall L \in \mathcal{L}_2(1) \setminus \{A_2\}. $$

More precisely, using the parametrization in the half-fundamental domain (see Proposition 2.1), we naturally restrict our study to $L$ parametrized by $(x, y) \in \mathcal{D}$ defined by (2.1). We therefore write $Q_{\alpha, \beta}(x, y)$ instead of $Q_{\alpha, \beta}(L)$ when $L$ has such parametrization and we aim to show, in three steps presented below, that

$$ Q_{\alpha, \beta}(x, y) > \frac{\alpha}{\beta}, \quad \forall (x, y) \in \mathcal{D} \setminus \left\{ \left( \frac{1}{2}, \frac{\sqrt{3}}{2} \right) \right\}. \qquad (3.2) $$

**Step 1. Smoothness at $(x, y) = (1/2, \sqrt{3}/2)$.** We first show that $Q_{\alpha, \beta} \in C^1(\mathcal{D})$. The smoothness on $\mathcal{D}\{(1/2, \sqrt{3}/2)\}$ is clear since $L \mapsto \zeta_L(s)$ is a smooth function on $\mathcal{D}$ when $s > 2$. The continuity at $(1/2, \sqrt{3}/2)$ follows from the fact that $Q_{\alpha, \beta}(x, y)$ converges to a positive constant as $(x, y) \to (1/2, \sqrt{3}/2)$ when $(x, y) \in \mathcal{D}\{(1/2, \sqrt{3}/2)\}$. Indeed, it has been proven by Coulangeon and Schürmann in [16] that $A_2$ is a strict local minimum in $\mathcal{L}_2(1)$ of the Epstein zeta function $L \mapsto \zeta_L(s)$ for all $s > 2$, which means by L'Hospital's Rule that

$$ \lim_{\substack{(x,y) \to (1/2, \sqrt{3}/2) \\ (x,y) \in \mathcal{D}\{(1/2, \sqrt{3}/2)\}}} Q_{\alpha, \beta}(x,y) = \lim_{\substack{L \to A_2 \\ L \in \mathcal{L}_2(1) \setminus \{A_2\}}} \frac{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)} > 0, $$

since both numerator and denominator are positive. Furthermore, for any $z \in \{x, y\}$ and any $L$ represented by $(x, y) \in \mathcal{D}$, we have that

$$ \partial_z Q_{\alpha, \beta}(x, y) = \left( \frac{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)} \right) \left( \frac{\partial_z \zeta_L(\alpha)}{\zeta_L(\alpha) - \zeta_{A_2}(\alpha)} - \frac{\partial_z \zeta_L(\beta)}{\zeta_L(\beta) - \zeta_{A_2}(\beta)} \right) \qquad (3.3) $$

As seen above, the first factor of this expression is going to a positive constant. By L'Hospital's Rule again, the second factor goes to the difference of quotients of the z-derivatives of the quadratic form associate to the 2nd order derivative of $L \mapsto \zeta_L(s)$ at $L = A_2$, $s \in \{\alpha, \beta\}$, divided by the same quadratic form at the same point. These quotients are both finite by smoothness of $L \mapsto \zeta_L(s)$ and by the strict minimality of $A_2$ for these Epstein zeta functions (see [17] and [9, Eq. (6.6)] for an analogous computation involving Epstein zeta functions). We therefore obtain that

$$ \forall z \in \{x, y\}, \quad \lim_{\substack{(x,y) \to (1/2, \sqrt{3}/2) \\ (x,y) \in \mathcal{D}\{(1/2, \sqrt{3}/2)\}}} \partial_z Q_{\alpha, \beta}(x, y) < \infty. $$

We therefore obtain that $\nabla_{(x,y)} Q_{\alpha, \beta}(x, y)$ is continuous on $\mathcal{D}$, which means that $Q_{\alpha, \beta} \in C^1(\mathcal{D})$.
**Step 2. Restriction to a compact set.** Since $\mathcal{D}$ is an infinite set (in the $y$-direction), we cannot check (3.2) numerically. However, since, for any $(x, y) \in \mathcal{D}$,

$$Q_{\alpha,\beta}(x, y) = \frac{y^{\frac{\alpha}{2}} \sum'_{m,n} \frac{1}{((m+nx)^2 + y^2n^2)^{\frac{\alpha}{2}}} - \zeta_{A_2}(\alpha)}{y^{\frac{\beta}{2}} \sum'_{m,n} \frac{1}{((m+nx)^2 + y^2n^2)^{\frac{\beta}{2}}} - \zeta_{A_2}(\beta)},$$

it is straightforward that $\lim_{y \to +\infty} Q_{\alpha,\beta}(x, y) = +\infty$, which means that (3.2) is satisfied at least for $y$ large enough and any $x \in [0, 1/2]$. More precisely, we have the explicit value $y_{\alpha,\beta}$ for such $y$ given by the following lemma whose proof is postponed to Section 4 to increase the readability of this section.

**Lemma 3.1.** For any $\alpha > \beta > 2$, we have that, for all $(x, y) \in \mathcal{D}$,

$$y > y_{\alpha,\beta} \Rightarrow Q_{\alpha,\beta}(x, y) > \frac{\alpha}{\beta},$$

where

$$y_{\alpha,\beta} := \begin{cases} \frac{9}{2^{1+\frac{2}{\beta}} \zeta(2\beta)^{\frac{2}{\beta}}} \left( \zeta_{\mathbb{Z}^2}(\beta) + \sqrt{\zeta_{\mathbb{Z}^2}(\beta)^2 - \frac{2}{3^\beta} \zeta(2\beta) (2\zeta_{A_2}(\beta) - \zeta_{A_2}(2\beta))} \right)^{\frac{2}{\beta}} & \text{if } \alpha = 2\beta, \\ \frac{3^{\frac{\alpha}{\alpha-\beta}}}{2^{1+\frac{2}{\alpha-\beta}}} \left( \frac{\alpha \zeta_{\mathbb{Z}^2}(\beta)}{\beta \zeta(\alpha)} \right)^{\frac{2}{\alpha-\beta}} & \text{otherwise,} \end{cases}$$

and where $\zeta(s) := \sum_{m \in \mathbb{N}} \frac{1}{m^s}$ is the Riemann zeta function.

**Remark 3.2** (Difference between the two values of $y_{\alpha,\beta}$). The second expression giving $y_{\alpha,\beta}$ still holds for $\alpha = 2\beta$, but is larger. Actually, for $(\alpha, \beta) = (12, 6)$, the two above expressions are quite similar (the difference is of order $10^{-3}$), but it has to be noticed that the difference increases with $\beta$. That is why we have derived a slightly better formula for $\alpha = 2\beta$, more precisely for saving some computational time in Step 3.

We therefore write

$$\mathcal{K}_{\alpha,\beta} := \{(x, y) \in \mathbb{R}^2 : x \in [0, 1/2], y \in [\sqrt{3}/2, y_{\alpha,\beta}], x^2 + y^2 \ge 1\},$$

in such a way that (3.2) holds in $\mathcal{D} \setminus \mathcal{K}_{\alpha,\beta}$.

**Step 3. Numerical checking of (3.2) in $\mathcal{K}_{\alpha,\beta}$.** It remains to show that (3.2) holds in $\mathcal{K}_{\alpha,\beta}$. For that, since we are only summing a finite number of terms $\zeta_L^N(s)$ (see (2.2)) of the Epstein zeta functions where $s \in \{\alpha, \beta\}$, we need to know how to control the value of the rest. The quantities summed in the Epstein zeta functions we are considering are converging fast to zero and the error terms $\zeta_L(s) - \zeta_L^N(s)$ have been studied for instance in [18] where an algorithm is written to compute $\zeta_L(s)$ with any degree of accuracy. It is also clear that the larger $y$ is, the more terms we need (i.e. more $N$ should be large) since there are more points at small distance from each others in $L$. The same can be done for the Riemann zeta function $\zeta(\alpha)$, for which it is well-known that the error term is bounded by $\frac{\alpha}{\alpha-1}N^{1-\alpha}$ (see e.g. [30, p. 3]). Since $y_{\alpha,\beta}$ is written in terms of Epstein and Riemann zeta functions, it is therefore easy to compute a new approximating value $\overline{y_{\alpha,\beta}}$ such that, for chosen $k \in \mathbb{N}$,

$$y_{\alpha,\beta} \le \overline{y_{\alpha,\beta}} \quad \text{and} \quad 10^k \overline{y_{\alpha,\beta}} \in \mathbb{N}.$$
This upper bound $\overline{y_{\alpha,\beta}}$, and more precisely the integer $k$, can be changed according to the accuracy we use in the sequel. We now write

$$ \overline{\mathcal{K}}_{\alpha,\beta} = \{(x, y) \in \mathbb{R}^2 : x \in [0, 1/2], y \in [\sqrt{3}/2, \overline{y_{\alpha,\beta}}], x^2 + y^2 \ge 1\} \quad (3.4) $$

and we want to show that (3.2) holds in $\overline{\mathcal{K}}_{\alpha,\beta} \supset \mathcal{K}_{\alpha,\beta}$.

Our method relies on the fact that one can only compute a finite number of values $\{Q_{\alpha,\beta}(x_i, y_j)\}_{i,j}$ where $(x_i, y_j) \in \overline{\mathcal{K}}_{\alpha,\beta}$ in order to approximate $\min_{(x,y) \in \overline{\mathcal{K}}_{\alpha,\beta}} Q_{\alpha,\beta}(x, y)$ with a sufficient degree of accuracy. Therefore, the fact that $\min_{i,j} Q_{\alpha,\beta}(x_i, y_j) > \frac{\alpha}{\beta}$ will imply (3.2). We actually create a grid of $G_{\alpha,\beta}^{\delta}$ of $I \times J$ points $(x_i, y_j) \in \overline{\mathcal{K}}_{\alpha,\beta}$, as illustrated in Figure 2, such that:

$$ \bullet \ x_1 = 0 < x_2 < \dots < x_I = 1/2 \text{ and } y_1 \approx \frac{\sqrt{3}}{2} < y_2 < \dots < y_J = \overline{y_{\alpha,\beta}}, $$

$$ \bullet \ x_{i+1} - x_i = y_{j+1} - y_j = \delta = \frac{1}{2I} = \frac{\overline{y_{\alpha,\beta}} - y_1}{J} \text{ for all } i \in \{1, \dots, I-1\}, \text{ for all } j \in \{1, \dots, J-1\}, $$

and where $y_1$ is an approximation of $\frac{\sqrt{3}}{2}$ by above with the accuracy given by $\delta$, i.e. $|y_1 - \sqrt{3}/2| < \delta$. We also choose $N$ such that $\zeta_L^N(s)$ is well approximated where $L$ is parametrized by any $(x, \overline{y_{\alpha,\beta}})$, $x \in [0, 1/2]$, in such a way that any other lattice parametrized by $(x, y) \in \overline{\mathcal{K}}_{\alpha,\beta}$ is also well approximated.

<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of the fundamental domain $\mathcal{D}$ given by (2.1), the compact set $\overline{\mathcal{K}}_{\alpha,\beta} \subset \mathcal{D}$ (in grey), the threshold value $\overline{y_{\alpha,\beta}}$ approximating by above $y_{\alpha,\beta}$ given by Lemma 3.1, the point $(\frac{1}{2}, \frac{\sqrt{3}}{2})$ representing the triangular lattice $\mathbb{A}_2$ as well as a possible square grid $G_{\alpha,\beta}^{\delta} = \{(x_i, y_j)\}_{i,j}$ with increment $\delta$ where we compute the values of $Q_{\alpha,\beta}$.</figcaption></figure>

Since $Q_{\alpha,\beta} \in C^1(\overline{\mathcal{K}}_{\alpha,\beta})$, it follows that $Q_{\alpha,\beta}$ is a Lipschitz function in $\overline{\mathcal{K}}_{\alpha,\beta}$, i.e. there exists a positive constant $M_{\alpha,\beta}$ such that, for all $(x, y)$, $(\bar{x}, \bar{y}) \in \overline{\mathcal{K}}_{\alpha,\beta}$,

$$ |Q_{\alpha,\beta}(x, y) - Q_{\alpha,\beta}(\bar{x}, \bar{y})| \le M_{\alpha,\beta}|(x, y) - (\bar{x}, \bar{y})|, \quad (3.5) $$

which implies in particular that, for all $(x, y) \in \overline{\mathcal{K}}_{\alpha,\beta}$ and all $(x_i, y_j) \in G_{\alpha,\beta}^{\delta}$,

$$ Q_{\alpha,\beta}(x, y) \ge Q_{\alpha,\beta}(x_i, y_j) - M_{\alpha,\beta}|(x, y) - (x_i, y_j)|. $$
This means that the latter inequality holds for $(x, y)$ in any square centred at $(x_i, y_j)$ with sidelength $\delta$ parallel to the axis. In that case, we have $|(x, y) - (x_i, y_j)| < \frac{\sqrt{2}}{2}\delta$. We now choose $\delta$ such that $M_{\alpha,\beta}\delta\frac{\sqrt{2}}{2}$ is small enough in such a way that we can numerically check that

$$ \min_{(x_i, y_j) \in G_{\alpha, \beta}^{\delta}} Q_{\alpha, \beta}(x_i, y_j) - M_{\alpha, \beta}\delta \frac{\sqrt{2}}{2} > \frac{\alpha}{\beta}. $$

Therefore, one can be certain that (3.2) holds in $\overline{\mathcal{K}}_{\alpha,\beta}$ and the proof is completed.

*Computation of $M_{\alpha,\beta}$.* We do not give here any explicit formula for such $M_{\alpha,\beta}$ which covers all the possible case, but the following computation yields to a simple estimate given below in (3.6). We begins by writing, for $s \in \{\alpha, \beta\}$,

$$ \Delta_s(x, y) := \zeta_L(s) - \zeta_{A_2}(s). $$

Therefore $Q_{\alpha,\beta}(x, y) = \frac{\Delta_{\alpha}(x,y)}{\Delta_{\beta}(x,y)}$ and we have $\forall (x, y) \in \overline{\mathcal{K}}_{\alpha,\beta}$, using (3.3),

$$ \| \nabla_{(x,y)} Q_{\alpha,\beta} \| \le Q_{\alpha,\beta}(x,y) \left\| \frac{\nabla_L \zeta_L(\alpha)}{\Delta_{\alpha}(x,y)} \right\| + Q_{\alpha,\beta}(x,y) \left\| \frac{\nabla_L \zeta_L(\beta)}{\Delta_{\beta}(x,y)} \right\|. $$

Moreover, we recall that $\Delta_s(x, y) \ge 0$ vanishes only for $(x, y) = (1/2, \sqrt{3}/2)$ (by uniqueness of the triangular minimizer) and $\nabla_L \zeta_L(s)$ vanishes only for $(x, y) \in \{(0, 1), (1/2, \sqrt{3}/2)\}$ (the only critical points of the Epstein zeta function are the square and the triangular lattices).

We can now bound above $Q_{\alpha,\beta}(x, y)$ and $\left\| \frac{\nabla_L \zeta_L(s)}{\Delta_s(x,y)} \right\|$ for $(x, y) \in \overline{\mathcal{K}}_{\alpha,\beta}$ in the same way by splitting $\overline{\mathcal{K}}_{\alpha,\beta}$ into two parts:

* in a ball $B_{\varepsilon}(A_2)$ (for the Euclidean norm) centred in $A_2$ and with radius $\varepsilon$ (that we will fix afterwards). A Taylor expansion gives us here a simple upper bound for $\left\| \frac{\nabla_L \zeta_L(s)}{\Delta_s(x,y)} \right\|$ in terms of the third (for the numerator) and the second derivative (for the denominator) of $L \mapsto \zeta_L(s)$ that can be easily bounded in the ball by iteration of the formula given below for $\|\nabla_L \zeta_L(s)\|$. The same Taylor expansion approach is used to bound above $Q_{\alpha,\beta}$.

* outside $B_{\varepsilon}(A_2)$. We roughly estimate the quotients by using the fact that, for $s > 2$, we have, using the fact that $x \in [0, 1/2]$ and $y \le \overline{y_{\alpha,\beta}}$,

$$ \begin{aligned} \|\nabla_L \zeta_L(s)\|^2 &= (\partial_x \zeta_L(s))^2 + (\partial_y \zeta_L(s))^2 \\ &= \left( \partial_x \left\{ y^{\frac{s}{2}} \sum_{m,n} \frac{1}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}}} \right\} \right)^2 + \left( \partial_y \left\{ y^{\frac{s}{2}} \sum_{m,n} \frac{1}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}+1}} \right\} \right)^2 \\ &= \left( -\frac{s}{2} y^{\frac{s}{2}} \sum_{m,n} \frac{2n(m+nx)}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}+1}} \right)^2 + \left( \frac{s}{2} y^{\frac{s}{2}-1} \sum_{m,n} \frac{(m+nx)^2 - y^2n^2}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}+1}} \right)^2 \\ &= \frac{s^2}{4} y^{s-2} \left\{ y^2 \left( \sum_{m,n} \frac{2n(m+nx)}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}+1}} \right)^2 + \left( \sum_{m,n} \frac{(m+nx)^2 - y^2n^2}{((m+nx)^2 + y^2n^2)^{\frac{s}{2}+1}} \right)^2 \right\} \\ &\le 2^s s^2 \overline{y_{\alpha,\beta}}^{s-2} \left\{ \overline{y_{\alpha,\beta}}^2 \left( \sum_{m,n} \frac{2|n|(|m| + 0.5|n|)}{(m^2 + n^2)^{\frac{s}{2}+1}} \right)^2 + \left( \sum_{m,n} \frac{m^2 + 0.25n^2 + |m||n|}{(m^2 + n^2)^{\frac{s}{2}+1}} \right)^2 \right\} \\ &=: S_{\alpha,\beta}(s), \end{aligned} $$
and using a simple lower/upper bounds for $\Delta_s(x, y)$ – depending on $\varepsilon$ for the lower bound. In particular, we easily have, by growth of $L \mapsto \zeta_L(s)$ in the $y$-direction and its degrowth in the $x$-direction (see [34]), for any $(x, y) \in \overline{\mathcal{K}}_{\alpha, \beta}$,

$$ \Delta_{\alpha}(x, y) = \zeta_{L}(\alpha) - \zeta_{A_{2}}(\alpha) \leq \zeta_{\overline{L}}(\alpha), $$

where $\overline{L}$ is parametrized by the point $(0, \overline{y_{\alpha, \beta}}) \in \mathcal{D}$.

Finally, choosing $\varepsilon$ such that $\Delta_s(x, y) > 1$ – which is straightforward from the variation of $L \mapsto \zeta_L(s)$ (decreasing in the $x$-direction and increasing in the $y$-direction as shown in [34]) – allows, using estimates on the second derivative of $L \mapsto \zeta_L(s)$, to find a constant $M_{\alpha, \beta}$ that we approximate in order to make it numerically tractable. We observe that this constant is given by the maximum between the values we found in the neighborhood of $A_2$ and basically the values in the upper boundary of $\overline{\mathcal{K}}_{\alpha, \beta}$. For the examples we cover in this article, we found that

$$ M_{\alpha, \beta} = \zeta_{\overline{L}}(\alpha) (S_{\alpha, \beta}(\alpha) + S_{\alpha, \beta}(\beta)). \qquad (3.6) $$

**Remark 3.3 (List of the needed approximations - A summary).** To summarize, we actually need to estimate the following quantities in order to complete our Step 3:

* the Epstein zeta functions $\zeta_L(\alpha)$, $\zeta_L(\beta)$, that are approximated by $\zeta_L^N(\alpha)$, $\zeta_L^N(\alpha)$ given in (2.2), and for which a bound of the error terms is known (see e.g. [18]). Notice that, since we are considering only sufficiently large exponents that are not treated in [2] (for which we already know that Theorem 1.1 holds), all our sums are converging very fast to zero.

* the value of $\overline{y_{\alpha, \beta}}$ based on the real $y_{\alpha, \beta}$ given by Lemma 3.1 which is written in terms of

- the Epstein zeta function (see previous point),

- the Riemann zeta function $\zeta(\alpha)$ for which the error term is bounded by $\frac{\alpha}{\alpha-1}N^{1-\alpha}$ is well-known (see e.g. [30, p. 3]).

* the upper bound $M_{\alpha, \beta}$ of $|\nabla_{(x,y)} Q_{\alpha, \beta}(x, y)|$ in the compact set $\overline{\mathcal{K}}_{\alpha, \beta}$ – for which we need a very rough estimate as given by (3.6) – that depends on

- the values of $S_{\alpha, \beta}(s)$ for $s \in \{\alpha, \beta\}$ that we estimate by its first terms when $|m| \le N$ and $|n| \le N$, as we do for $\zeta_L(\alpha)$ by $\zeta_L^N(\alpha)$.

- the value of $\overline{y_{\alpha, \beta}}$ (see previous point).

- in the following computations, we will round the approximation to the upper integer part. It does not really matter since we therefore choose $\delta$ with respect to $M_{\alpha, \beta}$ and the computation time stays reasonable.

Therefore, the fact that all the expressions we are summing or integrating are going to zero very fast allows us to choose some reasonable numbers $N$ (in the sense of “not too large”) and $\delta$ (in the sense of “not too small”) for computing our quantities and be sure that the inequality $Q_{\alpha, \beta}(x, y) > \frac{\alpha}{\beta}$ holds in $\overline{\mathcal{K}}_{\alpha, \beta}$ whereas it is established in Lemma 3.1 that the inequality holds in $\mathcal{D} \setminus \overline{\mathcal{K}}_{\alpha, \beta}$.

**Remark 3.4 (Possible improvements of our method).** We obviously do not claim that our estimates are the best one. The bound $y_{\alpha, \beta}$ given in Lemma 3.1 might be improved and a precise value of $M_{\alpha, \beta}$ could be computed by using more optimal bounds on $Q_L(m, n)$. However, in dimension $d = 2$, we can easily compute a large number of values such that the total computational time is less than an hour. This is already quite long, but rather far from the 19 hours needed in Sarnak-Strombergsson’s method [39]. Of course, a three-dimensional adaption (if any) of such method should lead to better estimates in order to considerably reduce this computational time.
**Remark 3.5 (Non-adaptability to other potentials and in other dimensions).** It has to be noticed that our method strongly relies on the homogeneity of the Epstein zeta function (Lemma 2.4) as well as the fact that $f$ is a one-well potential which implies the bound (2.5). It is therefore unclear whether our method can be adapted to other interactions such as Morse type potentials. Furthermore, the optimality of the triangular lattice for the Epstein zeta function given in Theorem 2.5 is also a key point in our method. Such optimality result is only available in dimensions $d \in \{2, 8, 24\}$ (see [15]) but nothing has been shown in dimension $d = 3$ so far (see e.g. [4, 39] or the numerical study in [10]). This makes our method impossible to apply in dimension $d = 3$ yet, whereas computational time is generally extremely long for such computer-assisted method in dimensions 8 and 24.

## 3.2 Application to the (12, 6) Lennard-Jones potential

We now check these three steps for $(\alpha, \beta) = (12, 6)$ in order to show Proposition 1.2 which implies Theorem 1.1. We use the software Scilab to perform our computations. Choosing $N = 40$, we therefore find

$$y_{12,6} \le \overline{y_{12,6}} = 7.52,$$

and $\zeta_L(6)$, $\zeta_L(12)$ are accurate with an error of order $10^{-4}$ for $L$ parametrized by $(x, \overline{y_{12,6}})$. Using (3.6), we compute our Lipschitz constant (approximated by its upper integer part)

$$M_{12,6} = 181,$$

and it is actually enough to compute our values with an increment of $\delta = 10^{-2}$ (in that case we have $M_{12,6}\delta\sqrt{2}/2 \approx 1.28$). Therefore, we split $[0, 1/2]$ into $I = 50$ equidistant values and $[0.87, 7.52]$, where $y_1 = 0.87 \approx \frac{\sqrt{3}}{2}$ (with an accuracy of $\delta$), into $J = 666$ equidistant values (both with increment $\delta$). We numerically checked that

$$\min_{(x_i, y_j) \in G_{12,6}^{\delta}} Q_{12,6}(x_i, y_j) - M_{12,6}\delta \frac{\sqrt{2}}{2} > \frac{12}{6} = 2$$

which completes the proof of Proposition 1.2 and Theorem 1.1.

## 3.3 Application to other exponents and justification of Conjecture 1.3

The global optimality of a triangular lattice for $E_f$ has been conjectured to hold for any pair of exponents $(\alpha, \beta)$ (see [2, 3, 9, 43]) and we have used the above strategy for a lot of them, which supports this conjecture. In particular, we have checked that our proof works for

$$(\alpha, \beta) \in \{(k, 6) : k \in \{14, 16, 18, 20, 22, 24\}\}, \quad (3.7)$$

i.e. when we keep the Van der Waals attraction $-br^{-6}$ whereas the repulsion parameters $\alpha$ varies. The range of exponents $\alpha$ is actually inspired by Kaplan's remark in [29, p. 184] about applicability of Lennard-Jones type potentials in physics and the fact that we kept the same value for $\beta$ is therefore only motivated by Physical arguments. The reader can also easily check the applicability of our method for her/his favorite pair of exponents $(\alpha, \beta)$.

In order to complete Step 3 for the exponents given in (3.7), as in the (12, 6) classical case we have chosen to keep $N = 40$, $\delta = 10^{-2}$, $I = 50$, $J = 666$ since $\alpha \ge 6$. All our numerical findings are summarized in Table 1. Thus, for these exponents $(\alpha, \beta)$ and all $(a, b) \in (0, \infty)^2$, the triangular lattice $\sqrt{V_{A_2}}A_2$ where

$$V_{A_2} = \left( \frac{a\alpha\zeta_{A_2}(\alpha)}{b\beta\zeta_{A_2}(\beta)} \right)^{\frac{2}{\alpha - \beta}}$$
<table><caption>Table 1: For the exponents $(\alpha, \beta)$ given in (3.7): values of $\overline{y_{\alpha,\beta}}$ and $M_{\alpha,\beta}$ (for $N = 40$ in the approximations).</caption><thead><tr><th>α</th><th>β</th><th>y&#773;<sub>α,β</sub></th><th>M<sub>α,β</sub></th></tr></thead><tbody><tr><td>14</td><td>6</td><td>5.23</td><td>95</td></tr><tr><td>16</td><td>6</td><td>4.18</td><td>72</td></tr><tr><td>18</td><td>6</td><td>3.60</td><td>50</td></tr><tr><td>20</td><td>6</td><td>3.22</td><td>39</td></tr><tr><td>22</td><td>6</td><td>2.97</td><td>34</td></tr><tr><td>24</td><td>6</td><td>2.77</td><td>33</td></tr></tbody></table>

is the unique minimizer of $E_f$ in $\mathcal{L}_2$, up to rotation, i.e. the analogue of Theorem 1.1 is shown for these parameters. Furthermore, for the same exponents and for all $(a, b) \in (0, \infty)^2$, $\sqrt{V}A_2$ is the unique minimizer of $E_f$ in $\mathcal{L}_2(V)$, up to rotation, if

$$0 < V < \left( \frac{a\alpha}{b\beta} \right)^{\frac{2}{\alpha-\beta}},$$

which shows Proposition 1.2 for the exponents given in (3.7).

Let us now explain why Conjecture 1.3 should hold. Indeed, showing that $Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta}$ for all $2 < \beta < \alpha$ and all $L \in \mathcal{L}_2(1)$ is equivalent with showing that

$$s \mapsto \frac{\zeta_L(s) - \zeta_{A_2}(s)}{s}$$

is strictly increasing on $(2, +\infty)$ for all $L \in \mathcal{L}_2(1)$, i.e.

$$\partial_s \left( \frac{\zeta_L(s) - \zeta_{A_2}(s)}{s} \right) > 0, \quad \forall L \in \mathcal{L}_2(1), \quad \forall s > 2.$$

Computing the above derivative, we find

$$\begin{gathered} \partial_s \left( \frac{\zeta_L(s) - \zeta_{A_2}(s)}{s} \right) = \frac{s(\partial_s \zeta_L(s) - \partial_s \zeta_{A_2}(s)) - \zeta_L(s) - \zeta_{A_2}(s)}{s^2} > 0 \quad \forall L \in \mathcal{L}_2(1), \quad \forall s > 2 \\ \iff s \partial_s \zeta_L(s) - \zeta_L(s) > s \partial_s \zeta_{A_2}(s) - \zeta_{A_2}(s) \quad \forall L \in \mathcal{L}_2(1), \quad \forall s > 2. \end{gathered}$$

Since

$$s \partial_s \zeta_L(s) - \zeta_L(s) = s \sum'_{p \in L} \partial_s \left( \frac{1}{|p|^s} \right) - \sum'_{p \in L} \frac{1}{|p|^s} = - \sum'_{p \in L} \frac{s \log |p| + 1}{|p|^s},$$

Conjecture 1.3 implies that $Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta}$ for all $2 < \beta < \alpha$ and all $L \in \mathcal{L}_2(1)$ and therefore implies the optimality of the triangular lattice $\sqrt{V}A_2$ for $E_f$ for all $2 < \beta < \alpha$.

# 4 Appendix: Proof of Lemma 3.1

Let us show how we obtain our threshold value $y_{\alpha,\beta}$.

*Proof of Lemma 3.1.* We start the proof by showing that, for all $(m, n) \in \mathbb{Z}$ and all $(x, y) \in \mathcal{D}$,

$$\frac{m^2 + n^2}{2} \le (m + xn)^2 + n^2 y^2 \le \frac{3}{2}m^2 + \left(\frac{3}{4} + y^2\right) n^2,$$
where we have combined Young's inequality with the fact that $x^2 + y^2 \ge 1$ and $x \in [0, 1/2]$. It follows that, for all $L \in \mathcal{D}$,

$$Q_{\alpha,\beta}(L) = \frac{y^{\frac{\alpha}{2}} \sum_{m,n}^{\prime} \frac{1}{((m+nx)^2 + y^2 n^2)^{\frac{\alpha}{2}}} - \zeta_{A_2}(\alpha)}{y^{\frac{\beta}{2}} \sum_{m,n}^{\prime} \frac{1}{((m+nx)^2 + y^2 n^2)^{\frac{\beta}{2}}} - \zeta_{A_2}(\beta)} \ge \frac{y^{\frac{\alpha}{2}} \sum_{m,n}^{\prime} \frac{1}{(\frac{3}{2}m^2 + (\frac{3}{4} + y^2)n^2)^{\frac{\alpha}{2}}} - \zeta_{A_2}(\alpha)}{y^{\frac{\beta}{2}} \sum_{m,n}^{\prime} \frac{2^{\frac{\beta}{2}}}{(m^2 + n^2)^{\frac{\beta}{2}}} - \zeta_{A_2}(\beta)} \\ \ge \frac{\frac{2^{1+\frac{\alpha}{2}} y^{\frac{\alpha}{2}}}{3^{\frac{\alpha}{2}}} \zeta(\alpha) - \zeta_{A_2}(\alpha)}{2^{\frac{\beta}{2}} y^{\frac{\beta}{2}} \zeta_{\mathbb{Z}^2}(\beta) - \zeta_{A_2}(\beta)},$$

where $2^{\frac{\beta}{2}}y^{\frac{\beta}{2}}\zeta_{\mathbb{Z}^2}(\beta) - \zeta_{A_2}(\beta) > 0$ is ensured by the fact that $\zeta_{\mathbb{Z}^2}(\beta) > \zeta_{A_2}(\beta)$ for all $\beta > 2$ by Theorem 2.5, and $(2y)^{\frac{\beta}{2}} \ge (\sqrt{3})^{\frac{\beta}{2}} > 1$ since $y \ge \sqrt{3}/2$. We therefore have

$$\frac{2^{1+\frac{\alpha}{2}} y^{\frac{\alpha}{2}}}{3^{\frac{\alpha}{2}}} \zeta(\alpha) - \zeta_{A_2}(\alpha) > \frac{\alpha}{\beta} \iff \eta_1 y^{\frac{\alpha}{2}} - \eta_2 y^{\frac{\beta}{2}} + \eta_3 > 0,$$

where $\{\eta_1, \eta_2, \eta_3\}$ are defined by

$$\eta_1 := \frac{2^{1+\frac{\alpha}{2}}}{3^{\frac{\alpha}{2}}} \zeta(\alpha) > 0, \quad \eta_2 := 2^{\frac{\beta}{2}} \frac{\alpha}{\beta} \zeta_{\mathbb{Z}^2}(\beta) > 0 \quad \text{and} \quad \eta_3 := \frac{\alpha}{\beta} \zeta_{A_2}(\beta) - \zeta_{A_2}(\alpha).$$

Let us now show that $\alpha_3 > 0$ by showing that the function $f$ defined by

$$f(s) := \frac{\zeta_{A_2}(s)}{s}$$

is decreasing on $(2, \infty)$. This will indeed imply that $\eta_3 = \alpha(f(\beta) - f(\alpha)) > 0$ for all $\alpha > \beta > 2$. Differentiating $f$ gives

$$f'(s) = -\frac{\zeta_{A_2}(s)}{s^2} - \frac{1}{s} \sum_{p \in A_2} \frac{\log |p|}{|p|^s}$$

and noticing that $\log|p| \ge \log\sqrt{2/\sqrt{3}} > 0$ for all $p \in A_2 \setminus \{0\}$ implies that $f'(s) < 0$ for all $s > 2$ and $f$ is therefore decreasing. In the general case, we remark that

$$\eta_1 y^{\frac{\alpha}{2}} - \eta_2 y^{\frac{\beta}{2}} + \eta_3 > \eta_1 y^{\frac{\alpha}{2}} - \eta_2 y^{\frac{\beta}{2}}$$

and $\eta_1 y^{\frac{\alpha}{2}} - \eta_2 y^{\frac{\beta}{2}} > 0$ if and only if $y > y_{\alpha,\beta} := (\frac{\eta_2}{\eta_1})^{\frac{2}{\alpha-\beta}}$. Therefore, the case $\alpha \ne 2\beta$ is proved. Furthermore, if $\alpha = 2\beta$, we can easily solve the inequality

$$\eta_1 y^{\frac{\alpha}{2}} - \eta_2 y^{\frac{\beta}{2}} + \eta_3 = \eta_1 X^2 - \eta_2 X + \eta_3 > 0, \quad \text{where} \quad X = y^{\frac{\beta}{2}}.$$

The discriminant of the polynomial $P(X) = \eta_1 X^2 - \eta_2 X + \eta_3$ is $\Delta = \eta_2^2 - 4\eta_1\eta_3 > 0$ and therefore

$$X > \frac{\eta_2 + \sqrt{\eta_2^2 - 4\eta_1\eta_3}}{2\eta_1} \Rightarrow P(X) > 0,$$

which means that

$$y > y_{\alpha,\beta} := \left( \frac{\eta_2 + \sqrt{\eta_2^2 - 4\eta_1\eta_3}}{2\eta_1} \right)^{\frac{2}{\beta}} \Rightarrow Q_{\alpha,\beta}(L) > \frac{\alpha}{\beta},$$

and proof is completed.

□

**Acknowledgement:** I would like to thank the Austrian Science Fund (FWF) and the German Research Foundation (DFG) for their financial support through the joint project FR 4083/3-1/I 4354 during my stay in Vienna. I am also grateful to the anonymous referees for their helpful suggestions.
References

[1] M. Zschornak ans T. Leisegang, F. Meutzner, H. Stöcker, T. Lemser, T. Tauscher, C. Funke, C. Cherkouk, and D. C. Meyer. Harmonic Principles of Elemental Crystals — From Atomic Interaction to Fundamental Symmetry. *Symmetry*, 10(6):228, 2018.

[2] L. Bétermin. Two-dimensional Theta Functions and Crystallization among Bravais Lattices. *SIAM J. Math. Anal.*, 48(5):3236–3269, 2016.

[3] L. Bétermin. Local variational study of 2d lattice energies and application to Lennard-Jones type interactions. *Nonlinearity*, 31(9):3973–4005, 2018.

[4] L. Bétermin. Local optimality of cubic lattices for interaction energies. *Anal. Math. Phys.*, 9(1):403–426, 2019.

[5] L. Bétermin. On energy ground states among crystal lattice structures with prescribed bonds. *J. Phys. A*, 54(24):245202, 2021.

[6] L. Bétermin, M. Friedrich, and U. Stefanelli. Lattice ground states for Embedded-Atom Models in 2D and 3D. *Lett. Math. Phys.*, 111(107), 2021.

[7] L. Bétermin and H. Knüpfer. On Born's conjecture about optimal distribution of charges for an infinite ionic crystal. *J. Nonlinear Sci.*, 28(5):1629–1656, 2018.

[8] L. Bétermin, L. De Luca, and M. Petrache. Crystallization to the square lattice for a two-body potential. *Arch. Ration. Mech. Anal.*, 240:987–1053, 2021.

[9] L. Bétermin and M. Petrache. Optimal and non-optimal lattices for non-completely monotone interaction potentials. *Anal. Math. Phys.*, 9(4):2033–2073, 2019.

[10] L. Bétermin, L. Šamaj, and I. Travènec. Three-dimensional lattice ground states for Riesz and Lennard-Jones type energies. *Studies in Applied Mathematics*, 2022. online first. DOI:10.1111/sapm.12533.

[11] L. Bétermin and P. Zhang. Minimization of energy per particle among Bravais lattices in $R^2$: Lennard-Jones and Thomas-Fermi cases. *Commun. Contemp. Math.*, 17(6):1450049, 2015.

[12] X. Blanc, C. Le Bris, and B. H. Yedder. A Numerical Investigation of the 2-Dimensional Crystal Problem. Preprint du laboratoire J.-L. Lions, Université de Paris 6, 2003.

[13] X. Blanc and M. Lewin. The Crystallization Conjecture: A Review. *EMS Surv. in Math. Sci.*, 2:255–306, 2015.

[14] J.W.S. Cassels. On a Problem of Rankin about the Epstein Zeta-Function. *Proceedings of the Glasgow Mathematical Association*, 4:73–80, 7 1959.

[15] H. Cohn, A. Kumar, S. D. Miller, D. Radchenko, and M. Viazovska. Universal optimality of the $E_8$ and Leech lattices and interpolation formulas. *to appear in Annals of Math.*, 2019. arXiv:1902.05438.

[16] R. Coul angeon and A. Schürmann. Energy Minimization, Periodic Sets and Spherical Designs. *Int. Math. Res. Not. IMRN*, pages 829–848, 2012.

[17] R. Coul angeon and A. Schürmann. Local energy optimality of periodic sets. *Communications in Number Theory and Physics*, 15(3):489–516, 2021.

[18] R. E. Crandall. Fast evaluation of the Epstein zeta functions. http://perfsci.com/free/techpapers/index.html, 1998.

[19] M. S. Daw and M. I. Baskes. Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals. *Phys. Rev. Lett.*, 50(17):1285–1288, 1983.

[20] M. S. Daw, S. M. Foiles, and M. I. Baskes. The embedded-atom method: a review of theory and applications. *Materials Science Reports*, 9(7-8):251–310, 1993.

[21] P. H. Diananda. Notes on Two Lemmas concerning the Epstein Zeta-Function. *Proceedings of the Glasgow Mathematical Association*, 6:202–204, 7 1964.

[22] V. Ennola. A Lemma about the Epstein Zeta-Function. *Proceedings of The Glasgow Mathematical Association*, 6:198–201, 1964.

[23] P. Epstein. Zur Theorie allgemeiner Zetafunctionen. *Mathematische Annalen*, 56(4):615–644, 1903.

[24] D. Frenkel and B. Smit. *Understanding Molecular Simulation: from Algorithms to Applications*. Academic Press, second edition edition, 2002.

[25] G. Friesecke, B. Schmidt, and Y. A. Yeung. Minimizing atomic configurations of short range pair potentials in two dimensions: crystallization in the Wulff shape. *Calc. Var. Partial Differential Equations*, 44(1-2):81–100, 2012.

[26] R. C. Heitmann and C. Radin. The Ground State for Sticky Disks. *J. Stat. Phys.*, 22:281–287, 1980.
[27] A. Henn. The Hexagonal Lattice and the Epstein Zeta Function. Dynamical Systems, Number Theory and Applications A Festschrift in Honor of Armin Leutbecher's 80th Birthday (Chapter 7), 2016.

[28] J. E. Jones. On the determination of molecular fields II. From the equation of state of a gas. Proc. R. Soc. London, Ser. A, 106:463, 1924.

[29] I. G. Kaplan. *Intermolecular Interactions: Physical Picture, Computational Methods, Model Potentials*. John Wiley and Sons Ltd, 2006.

[30] A. A. Karatsuba and S. M. Voronin. *The Riemann Zeta-Function*. De Gruyter, 1992.

[31] R. LeSar. *Introduction to Computational Materials Science. Fundamentals to Applications*. Cambridge University Press, 2013.

[32] L. De Luca and G. Friesecke. Crystallization in two dimensions and a discrete Gauss-Bonnet theorem. *J. Nonlinear Sci.*, 28(1):69–90, 2018.

[33] G. Mie. Zur kinetischen Theorie der einatomigen Körper. *Annalen der Physik*, 316(8):657–697, 1903.

[34] H. L. Montgomery. Minimal Theta Functions. *Glasg. Math. J.*, 30(1):75–85, 1988.

[35] Lívia B. Pártay, Christoph Ortner, Albert P. Bartók, Chris J. Pickard, and Gábor Csányi. Polytypism in the ground state structure of the Lennard-Jonesium. *Physical Chemistry Chemical Physics*, 19:19369–19376, 2017.

[36] C. Radin. The Ground State for Soft Disks. *J. Stat. Phys.*, 26(2):365–373, 1981.

[37] C. Radin. Low temperature and the origin of crystalline symmetry. *International Journal of Modern Physics B*, 1(5 and 6):1157–1191, 1987.

[38] R. A. Rankin. A Minimum Problem for the Epstein Zeta-Function. *Proceedings of The Glasgow Mathematical Association*, 1:149–158, 1953.

[39] P. Sarnak and A. Strömbergsson. Minima of Epstein's Zeta Function and Heights of Flat Tori. *Invent. Math.*, 165:115–151, 2006.

[40] S. L. Sobolev and V. Vaskevitch. *The Theory of Cubature Formulas*. Springer Netherlands, 1997.

[41] A. Terras. *Harmonic Analysis on Symmetric Spaces and Applications II*. Springer New York, 1988.

[42] F. Theil. A Proof of Crystallization in Two Dimensions. *Comm. Math. Phys.*, 262(1):209–236, 2006.

[43] I. Travěnec and L. Šamaj. Two-dimensional Wigner crystals of classical Lennard-Jones particles. *J. Phys. A: Math. Theor.*, 52(20):205002, 2019.