Lattice ground states for Embedded-Atom Models in 2D and 3D

The Embedded-Atom Model (EAM) provides a phenomenological description of atomic arrangements in metallic systems. It consists of a configurational energy depending on atomic positions and featuring the interplay of two-body atomic interactions and nonlocal effects due to the corresponding electronic clouds. The purpose of this paper is to mathematically investigate the minimization of the EAM energy among lattices in two and three dimensions. We present a suite of analytical and numerical results under different reference choices for the underlying interaction potentials. In particular, Gaussian, inverse-power, and Lennard-Jones-type interactions are addressed.


Introduction
Understanding the structure of matter is a central scientific and technological quest, cutting across disciplines and motivating an ever increasing computational effort. First-principles calculations deliver accurate predictions but are often impeded by the inherent quantum complexity, as systems size up [29]. One is hence led to consider a range of approximations. The minimization of empirical atomic pair-potentials represents the simplest of such approximations being able to describe specific properties of large-scaled atomic systems. Still, atomic pair-interactions fall short of describing the basic nature of metallic bonding, which is multi-body by nature, and often deliver inaccurate predictions of metallic systems.
The Embedded-Atom Model (EAM) is a semi-empirical, many-atom potential aiming at describing the atomic structure of metallic systems by including a nonlocal electronic correction. Introduced by Daw and Baskes [16], it has been used to address efficiently different aspects inherent to atomic arrangements including defects, dislocations, fracture, grain boundary structure and energy, surface structure, and epitaxial growth. Proving capable of reproducing experimental observations and being relatively simple to implement, the Embedded-Atom Model is now routinely used in molecular dynamic simulations [18,27]. In particular, it has been applied in a variety of metallic systems [21], including alkali metals Li, Na, K [19,26,37], transition metals Fe, Ni, Cu, Pd, Ag, Pt, Au [13,22,26,27], post-transition metals Al, Pb [13,24,34], the metalloid Si [3], and some of their alloys [13,25].
In the case of a metallic system with a single atomic species, the EAM energy is specified as Here, {x i } indicate atomic positions in R d and the long-range interaction potential φ : R + := (0, ∞) → R + modulates atomic pair-interactions. Atomic positions induce electronic-cloud distributions. The function ρ : R + → R + models the long-range electron-cloud contribution of an atom placed at x j on an atom placed at x i . The sum ρ i describes the cumulative effect on the atom placed at x i of the electronic clouds related to all other atoms. Eventually, the function F : R + → R + describes the energy needed to place (embed) an atom at position x i in the host electron gas created by the other atoms at positions {x j }.
Purely pair-interaction potentials can be re-obtained from the EAM model by choosing F = 0 and have been the subject of intense mathematical research under different choices for φ. The reader is referred to [12] for a survey on the available mathematical results. The setting F = 0 corresponds indeed to the so-called Born-Oppenheimer approximation [29], which is well adapted to the case of very low temperatures and is based on the subsequent solution of the electronic and the atomic problem. As mentioned, this approximation turns out to be not always appropriate for metallic systems at finite temperatures [7,34] and one is asked to tame the quantum nature of the problem. This is however very challenging from the mathematical viewpoint and rigorous optimality results for point configurations in the quantum setting are scarce [10,11]. The EAM model represents hence an intermediate model between zero-temperature phenomenological pair-interaction energies and quantum systems. Electronic effects are still determined by atomic positions, but in a more realistic nonlocal fashion when F is nonlinear, resulting in truly multi-body interaction systems, see [16,20,34] and [18] for a review.
The aim of this paper is to investigate point configurations minimizing the EAM energy. Being interested in periodic arrangements, we restrict our analysis to the class of lattices, namely infinite configurations of the This reduces the optimality problem to finite dimensions, making it analytically and numerically amenable. In particular, the EAM energy-per-atom of the lattice L takes the specific form In the classical pair-interaction case F = 0, the lattice energy E has already received attention and a variety of results are available, see [4,8,9,14,28,31] and the references therein. Such results are of course dependent on the choice of the potential φ. Three reference choices for φ are the Gaussian φ(r) = e −πδr 2 for δ > 0, the inverse-power law φ(r) = r −s for s > d, and the Lennard-Jones-type form φ(r) = ar −α − br −β for d < β < α and a, b > 0. In the Gaussian case, it has been shown by Montgomery [28] that, for all δ > 0, the triangular lattice of unit density is the unique minimizer (up to isometries) of E with F = 0 among unit-density lattices. The same can be checked for the the inverse-power-law case by a Mellin-transform argument. More generally, the minimality of the triangular lattice of unit density is conjectured by Cohn and Kumar in [14,Conjecture 9.4] to hold among all unit-density periodic configurations. This fact is called universal optimality and has been recently proved in dimension 8 and 24 for the lattice E 8 and the Leech lattice Λ 24 , respectively [15]. In the Lennard-Jones case, the minimality in 2d of the triangular lattice at fixed density has been investigated in [4,10], the minimality in 3d of the cubic lattice is proved in [6], and more general properties in arbitrary dimensions have been investigated in [9]. A recap of the main properties of the Lennard-Jones case is presented in Subsection 2.3. These play a relevant role in our analysis.
In this paper, we focus on the general case F = 0, when F is nonlinear. More precisely, we discuss the reference cases of embedding functions F of the form F (r) = r t log(γr) or F (r) = r t for t, γ > 0. The first, logarithmic choice is the classical one chosen to fit with the so-called Universal Binding Curve (see e.g., [30]) and favoring a specific minimizing value r 0 > 0, see [1,13]. The second, power-law form favors on the contrary r 0 = 0 and allows for a particularly effective computational approach. Let us mention that other choices for F could be of interest. In particular, the form F (r) = −c √ r, c > 0, is related to the Finnis-Sinclair model [20] and is discussed in Remark 4.3. Some of our theory holds for general functions F , provided that they are minimized at a sole value r 0 . We call such functions of one-well type.
Our main theoretical results amount at identifying minimizers in the specific reference case of F (r) = r log r and ρ(r) = r −s . More precisely, we find the following: • (Inverse-power law) If φ(r) = r −α , the minimizers of E coincide with those of the Lennard-Jones potential r → r −α − r −s , up to rescaling (Theorem 4.1); • (Lennard-Jones) If φ(r) = ar −α − br −β , under some compatibility assumptions on the parameters, the minimizers of E coincide with those of the Lennard-Jones potential r → r −α −r −s (Theorem 5.2).
Actually, both results hold for more general embedding functions F , see (4.1) and . With this at hand, the problem can be reduced to the pure Lennard-Jones case (i.e., F = 0) which is already well understood. In particular, in the two dimensional case we find that the triangular lattice, up to rescaling and isometries, is the unique minimizer of E in specific parameters regimes. These theoretical findings are illustrated by numerical experiments in two and three dimensions. By alternatively choosing the Gaussian ρ(r) = e −δr 2 , in two dimensions we additionally observe the onset of a phase transition between the triangular and an orthorhombic lattice, as δ decreases. In three dimensions, both in the inverse-powerlaw case ρ(r) = r −s and in the Gaussian case ρ(r) = e −δr 2 , the simple cubic lattice Z 3 is favored against the face-centered and the body-centered cubic lattice for s or δ small, respectively.
In the power-law case F (r) = r t , for ρ of inverse-power-law type and φ of Lennard-Jones type and specific, physically relevant choices of parameters, one can conveniently reduce the complexity of the optimization problem from the analytical standpoint. This reduction allows to explicitly compute the EAM energy for any lattice of unit density, hence allowing to investigate numerically minimality in two and three dimensions. Depending on the parameters, the relative minimality of the triangular, square, and orthorhombic lattices in two dimensions and the simple cubic, body-centered cubic, and face-centered cubic lattices in three dimension is ascertained. This is the plan of the paper: Notation on potentials and energies are introduced in Subsections 2.1 and 2.2. The two subcases F = 0 and φ = 0 are discussed in Subsection 2.3 and in Section 3, respectively. In particular, known results on Lennard-Jones-type interactions are recalled in Subsection 2.3. The inversepower-law case φ(r) = r −α is investigated in Section 4. The Lennard-Jones case φ(r) = ar −α − br −β is addressed theoretically and numerically in Section 5. In particular, Subsection 5.1 contains the classical case F (r) = r log r, and Subsection 5.2 discusses the power-law case F (r) = r t .
is a basis of R d . We write L d (1) ⊂ L d for the set of all lattices with unit density, which corresponds to | det(u 1 , . . . , u d )| = 1.
In dimension two, any lattice L ∈ L 2 (1) can be written as is the so-called (half ) fundamental domain for L 2 (1) (see, e.g., [28,Page 76]). In particular, the square lattice Z 2 and the triangular lattice with unit density, denoted by A 2 ∈ L d (1), are given by the respective choices (x, y) = (0, 1) and (x, y) = 1/2, √ 3/2 , i.e., In dimension three, the fundamental domain of L 3 (1) is much more difficult to describe (see e.g., [35,Section 1.4.3]) and its 5-dimensional nature makes it impossible to plot compared to the 2-dimensional D defined in (2.1). The Face-Centered Cubic (FCC) and Body-Centered Cubic (BCC) lattices with unit density are respectively indicated by D 3 ∈ L 3 (1) and D * 3 ∈ L 3 (1), and are defined as Remark 2.1 (Periodic configurations). All results in this paper are stated in terms of lattices, for the sake of definiteness. Let us however point out that the same statements hold in the more general setting of periodic configurations in dimensions d ∈ {8, 24}, on the basis of the recently proved optimality results from [15]. In dimension d = 2, universal optimality is only known among lattices, see [28]. Still, the validity of the Cohn-Kumar conjecture (see [14,Conjecture 9.4]) would allow us to consider more general periodic configurations as well.

Potentials and energies.
For any dimension d, let S d be the set of all functions f : for some η > 0 as r → ∞. By S + d ⊂ S d we denote the subset of nonnegative functions. We say that a continuous function F : R + → R is a one-well potential if there exists r 0 > 0 such that F is decreasing on (0, r 0 ) and increasing on (r 0 , ∞).
actually corresponds to the Epstein zeta function, which is defined by For any function F : R + → R and for any ρ ∈ S + d , we define the embedding energy E F,ρ : Finally, for any φ ∈ S d , any ρ ∈ S + d , and any F : (2.5) In the following, we investigate E under different choices of the potentials F , ρ, and φ. In some parts, we will require merely abstract conditions on the potentials, such as a monotone decreasing ρ or a one-well potential F . In other parts, we will consider more specific potentials. In particular, we will choose, for γ, δ, t, a, b > 0, s > d, and α > β > d, Note that the choice of s, δ, α, and β implies that φ ∈ S d and ρ ∈ S + d , so that the sums in (2.2) and (2.4) are well defined.
For any L ∈ L d (1), any φ ∈ S d , any ρ ∈ S + d , and any F : R + → R, we define, if they uniquely exist, the following optimal scaling parameters for the energies:

A recap on the Lennard-Jones-type energy.
A classical problem is to study the F = 0 case for a Lennard-Jones-type potential Let us recap some known facts in this case [4,9], which will be used later on. We start by reducing the minimization problem on all lattices to a minimization problem on lattices of unit density only. This is achieved by computing the optimal scaling parameter of the energy λ φ L , see (2.6), for each L ∈ L d (1), which in turn allows to find the minimum of the energy among dilations of L. More precisely, in case (2.7), for all λ > 0 and all lattices L ∈ L d (1), one has where we use (2.3). (This energy was studied first in [4, Section 6.3].) Then, we find the unique minimizer and therefore the energy is given by The latter inequality follows from the fact that α > β. Consequently, for any lattices L, Λ ∈ L d (1), we have that This means that finding the lattice with minimal energy amounts to minimizing the function . This is particularly effective in dimension two where for fixed (α, β) the minimizer can be found numerically by plotting L → min λ E[λL] in the fundamental domain D. Figure 1 shows the case (α, β) = (12, 6), i.e., when φ is the classical Lennard-Jones potential. The global minimum of E φ in L 2 appears to be the triangular lattice For a certain range of parameters (α, β), this observation can be rigorously ascertained. Indeed, for d = 2, it is shown in [4, Theorem 1.2.B.] that the global minimum of E φ is uniquely achieved by a triangular lattice and Γ is the classical Gamma function Γ(r) = ∞ 0 x r−1 e x dx for r > 0. (In the sequel, all statements on uniqueness are intended up to isometries, without further notice.) In fact, under condition (2.10) one has that [4] • A 2 is the unique minimizer in As pointed out in [4,Remark 6.18], it is necessary to choose 2 < β < α < M ≈ 9.2045818 in order to obtain these optimality results by using the method developed there. In particular, this means that the following pairs of integer exponents can be chosen:  ζ L (6) 12 in the fundamental domain D. The triangular lattice A 2 with coordinates (1/2, √ 3/2) appears to be the unique minimizer. Moreover, Z 2 with coordinates (0, 1) appears to be a saddle point.
We now ask ourselves what is the minimal scaling parameter λ and the corresponding lattice L ∈ L d (1) for which E φ [λL] is minimized. Physically, this would correspond to identifying the first minimum of E φ starting from a high-density configuration by progressively decreasing the density. We have the following.
Proposition 2.2 (Smallest volume meeting the global minimum). Let φ be a Lennard-Jones-type potential as in (2.7).
for all L ∈ L d (1). We thus have As we are assuming that where we use that α > β. In view of (2.8), this shows that we have a double equality in (2.12). This implies also equality in (2.11) which is equivalent to e * (L) = e * (L d ). Therefore, it follows that L = L d up to rotation, by uniqueness of the minimizer L d of e * . We refer to Figure 2 for an illustration in the two-dimensional case (α, β) = (12,6). Note that in this case the global minimum is not known. Still, the triangular lattice appears to be the first stable structure reached by increasing the volume (decreasing the density). This is in agreement with Figure 1 and Proposition 2.2. Recall that the triangular lattice also minimizes L → ζ L (β) on L d (1), as required in the statement of Proposition 2.2, see [28].
Notice that in dimension d = 3 there is no rigorous result concerning the minimizer of E φ in L 3 . Only local minimality results for cubic lattices Z 3 , D 3 , D * 3 have been derived in [6]. Numerical investigations suggest that λ φ D3 D 3 is the unique minimizer of E φ in L 3 for any values α > β > d of the exponents, see, e.g., [39, Figure 5], [9, Figures 5 and 6] and [6, Conjecture 1.7]. Therefore, we can conjecture that D 3 is the unique minimizer of L → λ φ L in L 3 (1) by application of Proposition 2.2.

Properties of the embedding energy E F,ρ
In this section we focus on the properties of the embedding energy E F,ρ given in (2.4). Although other choices for the potential F may been considered (see, e.g., [18,20]), we concentrate ourselves on the onewell case (see, e.g., [13] and references therein). In that case, it is clear that the global minimum of E F,ρ in L d can be achieved for any L by simply choosing λ such that E ρ [λL] = r 0 = argmin r>0 F (r). We now ask ourselves what is the minimal scaling parameter λ and the corresponding lattice L ∈ L d (1) for which E F,ρ [λL] achieves min F . In other words, what is the minimizer of L → λ F,ρ L in L d (1) (recall (2.6)). Physically, this would correspond to reach the ground state of the embedding energy min F starting from a high-density configuration by progressively decreasing the density.
Proof. Let r 0 > 0 be the unique minimizer of F , namely, F (r 0 ) = min F . Given any L ∈ L d (1), the fact that ρ ∈ S + d is strictly decreasing implies that λ → E ρ [λL] is strictly decreasing and goes to 0 at infinity and to ∞ at 0. Therefore, there exists a unique λ > 0 such that E ρ [λL] = r 0 . Such λ obviously coincides with λ F,ρ L given in (2.6). This shows the first part of the statement.
We note that Theorem 3.1 can be applied to the choice ρ(r) = r −s , s > d, and the triangular lattice A 2 , the E 8 lattice, or the Leech lattice Λ 24 in dimensions 2, 8, and 24, respectively. In fact, these lattices are the unique minimizer of L → E ρ [λL] for all λ > 0, see [15,28].
Let us mention that, in this setting, asking F to be one-well is not restrictive. In fact, if F is a strictly increasing (resp. decreasing) function, no optimal scaling parameters λ > 0 can be found since, for any L ∈ L d (1), E F,ρ [λL] will be minimized for λ → 0 (resp. λ → ∞).

The EAM energy with inverse-power interaction φ(r) = r −α
In this section, we study the energy E defined in (2.5) when φ is given by the inverse-power interaction φ(r) = r −α . The main result of this section is the following.
Theorem 4.1 (EAM energy for inverse-power interaction). For any α > s > d, let ρ(r) = r −s , let φ(r) = r −α , and let F ∈ C 1 (R + ). We assume that the functions g(r) := r 1−α/s F (r) and h(r) := F (r) − s α rF (r) for r > 0 (4.1) satisfy that g is strictly increasing on I := {F < 0}, that g(I) = (−∞, 0), and that h • g −1 is strictly decreasing on (−∞, 0). (Note that g −1 exists on (−∞, 0) and takes values in R + .) Then, λ E L exists for all L ∈ L d (1) and the following statements are equivalent: In particular, when d = 2 and H(α) < H(s) where H is defined by (2.10), then the unique minimizer of E in L 2 is the triangular lattice The gist of this result is the coincidence of the minimizers of E with those of Eφ forφ(r) = r −α − r −s (up to proper rescaling), under quite general choices of F . This results in a simplification of the minimality problem for E as one reduces to the study of minimality for the Lennard-Jones-type potentialφ, which is already well known, see Subsection 2.3. In particular, in two dimensions and under condition H(α) < H(s), the unique minimizer is a properly rescaled triangular lattice.
Before proving the theorem, let us present some applications to specific choices of F . Remark 4.2 (Application 1 -The classical case F (r) = r log r). We can apply this theorem to F (r) = r t log(γr) for t ∈ (0, α/s) and γ > 0 which is a one-well potential with minimum attained at point r t 0 := 1 γ e −1/t . In particular, the case F (r) = r log r is admissible since s < α. In fact, we have I = (0, r t 0 ) and Since g is strictly increasing on (0, r t 1 ) for r t 1 := 1 γ e 2ts−α t(α−ts) and r t 0 < r t 1 we have that g is strictly increasing on I. Moreover, g(I) = (−∞, 0). On the other hand, h is strictly decreasing on (0, r t 1 ). Therefore, also h • g −1 is strictly decreasing on (−∞, 0). Hence, Theorem 4.1 applies.

Remark 4.3 (Application 2 -Finnis-Sinclair model)
. Theorem 4.1 can also be applied to F (r) = −c √ r for c > 0. This case is known as the long-range Finnis-Sinclair model defined in [34], based on the work of Finnis and Sinclair [20] on the description of cohesion in metals and also used as a model to test the validity of machine-learning algorithms [23]. In this case, we obtain Since s < α < 2α, g is strictly increasing on I = {F < 0} = R + , g(I) = (−∞, 0), and h is strictly decreasing on R + . Therefore, Theorem 4.1 applies.
Remark 4.4 (Application 3 -inverse-power law). Also the inverse-power law F (r) = r −t for t > 0 satisfies the assumption of the theorem. In fact, we have g(r) = −tr −t−α/s and h(r) = 1 + st α r −t .
Proof of Theorem 4.1. In view of (2.3) and (2.5), for any λ > 0 and L ∈ L d (1) we have that The critical points of λ → E[λL] for fixed L are the solutions of This is equivalent to where g is given in (4.1), and e * (L) = ζ L (α) s ζ L (s) α was defined in (2.9). Since g −1 is positive and strictly increasing on (−∞, 0), we have that the unique critical point is given by In view of (4.2), we also have that ∂ λ E[λL] ≥ 0 if and only if g(λ −s ζ L (s)) ≤ − α s e * (L) 1 s , which is equivalent to λ ≥ λ * . In particular, λ → E[λL] is decreasing on (0, λ * ) and increasing on (λ * , ∞). This shows that λ * is a minimizer and thus λ * = λ E L , where λ E L is defined in (2.6).
By using the fact that (4.2) and the identity λ * = λ E L , the minimal energy among dilated copies λL of a given lattice L can be checked to be where h is defined in (4.1). By assumption h • g −1 is strictly decreasing on (−∞, 0). Hence, L d minimizes L → E[λ E L L] in L d (1) (uniquely) if and only if L d minimizes e * (uniquely). This shows the equivalence of the first two items in the statement. The equivalence to the third item has already been addressed in the discussion before (2.9). The two-dimensional case is a simple application of [4, Theorem 1.2.B.] which ensures that A 2 is the unique minimizer of e * in L 2 (1), as it has been already recalled in Subsection 2.3.
To complete the proof, it remains to show the final statement in d dimensions. Assume that L d is the unique minimizer of L → ζ L (s) in L d (1) as well as a minimizer of e * in L d (1). In this case, by using (4.3) and the identity λ * = λ E L , it indeed follows that L d is the unique minimizer of L → λ E L in L d (1), since g −1 is positive and increasing on (−∞, 0).

The EAM energy with Lennard-Jones-type interaction
We now move on to consider the full EAM energy E defined in (2.5) for Lennard-Jones-type potentials φ as in (2.7). We split this section into two parts. At first, we address the classical case F (r) = r log r analytically and numerically. Afterwards, we provide some further numerical studies for the power law case F (r) = r t . 5.1. The classical case F (r) = r log r. We start with two theoretical results and then proceed with several numerical investigations. 5.1.1. Two theoretical results. The following corollary is a straightforward application of Theorem 3.1.

Corollary 5.1 (Existence of parameters for the optimality of
for γ, t > 0, s > 2, α > β > 2, and a, b > 0. Then, given parameters (α, β, γ, s, t) such that H(α) < H(β), where H is defined in (2.10), one can find coefficients a and b such that the unique global minimizer in L 2 of E is the triangular lattice λ A2 A 2 where Moreover, A 2 is the unique minimizer of L → λ E L in L 2 (1).
Proof. We first remark that F and ρ satisfy the assumption of Theorem 3.1. By recalling (2.3), (2.6) and using the fact that argmin F = 1 γ e −1/t , we have is the unique minimizer of the sum of the two energies E F,ρ and E φ . The identity λ F,ρ A2 = λ φ A2 is equivalent to equation For this choice of a and b, we thus get that the unique global minimizer in L 2 of E is the triangular lattice The last statement follows by applying Proposition 2.2 to L d = A 2 .
The drawback of the result is that it is not generic in the sense that it holds only for specific coefficients a and b. We now give a result which holds in any dimension for all coefficients a, b > 0, at the expense of the fact that φ and ρ need to have the same decay O(r −s ). In this regard, the result is in the spirit Theorem 4.1 but under the choice φ(r) = ar −α − br −s .
Theorem 5.2 (EAM energy for Lennard-Jones-type interaction). Let F be as in Theorem 4.1 and additionally suppose that F is convex and in C 2 (R + ). Let for d < s < α and a, b > 0.
Then, λ E L exists for all L ∈ L d (1) and the following statements are equivalent: In particular, when d = 2 and H(α) < H(s) where H is defined by (2.10), then the unique minimizer of E in L 2 is the triangular lattice Proof. In view of (2.3), the energy E can be written as where g and h are defined in (4.1). We first check thatg is strictly increasing onĨ := {F < 0}. Indeed, since F (and henceF ) is convex and α > s, we get that for all r ∈Ĩ. Since by assumption g({F < 0}) = (−∞, 0) andĨ = {F < 0} ⊃ {F < 0}, we find g(Ĩ) = (−∞, 0). Eventually,h •g −1 is strictly decreasing on (−∞, 0), as well. We can hence apply Theorem 4.1 and obtain the assertion. Remark 5.3. As a consequence of Remark 4.2, the previous result can be applied to F (r) = r log r. Already for this F , in the case of a more general Lennard-Jones potential φ(r) = ar −α − br −β , the equation for the critical points of λ → E[λL] for a fixed lattice L is for a = αaζ L (α), b = s 2 ζ L (s), c = sζ L (s)(1 + log ζ L (s)), and d = βbζ L (β). This is generically not solvable in closed form when s = β, and makes the computation of E[λ E L L] more difficult. This is why we choose s = β in the above result.

Numerical investigation in 2d.
We choose s as parameter and fix t = γ = a = b = 1, and α = 12, β = 6, i.e., F (r) = r log r, ρ(r) = r −s , φ(r) = 1 We employ here a gradient descent method, which is rather computationally intensive. Note that a more efficient numerical method will be amenable in Subsection 5.2, as an effect of a different structure of the potentials. Numerically, we observe the following (see Figure 3): • For s > s 1 , s 1 ≈ 5.14, the triangular lattice λ E A2 A 2 is apparently the unique global minimizer of E. • For s < s 1 , the energy does not seem to have a global minimizer.
Furthermore, for s > s 0 , s 0 ≈ 5.09, we have checked (see Figure 4) that whereas the inequality is reversed if s < s 0 .  We now replace ρ by a Gaussian function. Namely, we consider the case In this case, the triangular lattice λ E A2 A 2 still seems to be minimizing E for large δ, see Figure 6. More precisely: • There exists δ 0 ≈ 1.04 such that, for δ > δ 0 , the triangular lattice λ E A2 A 2 is the global minimizer of E in L 2 .
• For δ < δ 0 , the global minimizer of E seems to move (continuously) in D increasingly following the y-axis as δ decreases to 0. For instance, -If δ = 1, then the minimizer is (0, y 1 ) where y 1 ≈ 1.014.
We have numerically studied the following function for s > 3, see Figure 7. We have found that there exist s 0 < s 1 < s 2 where s 0 ≈ 5.4985, s 1 ≈ 5.576, and s 2 ≈ 5.584 such that It is remarkable that for small values of s the simple cubic lattice Z 3 has lower energy with respect to the usually energetically favored D 3 and D * 3 . In the following, we will call θ L (δ) the lattice theta function with parameter δ > 0. Note however that under this name one usually refers to such sum including the term for p = 0 and with weight e −δπ|p| 2 .
In Figure 8 we plot the functions δ → min λ>0 E[λL] for L ∈ {D 3 , D * 3 , Z 3 }. We numerically observe that there exist 0 < δ 1 < δ 2 < δ 3 , where δ 1 ≈ 1.13, δ 2 ≈ 1.21, and δ 3 ≈ 1.223 such that It is indeed important that the EAM energy favors D 3 or D * 3 for some specific choice of parameters. In fact, FCC and BCC lattices are commonly emerging in metals. It is also remarkable that the simple cubic lattice Z 3 (up to rescaling) is favored with respect to D 3 or D * 3 for some other choice of parameters. In [6], we were able to identify a range of densities such that cubic lattices are locally optimal at fixed density, but it is the first time -according to our knowledge -that such phenomenon is observed at the level of the global minimizer.

5.2.
The power-law case F (r) = r t . In this subsection, we study the case where F (r) = r t , t > 0. Although F is not a one-well potential, this case turns out to be mathematically interesting. Indeed, we are able to present a special case where we can explicitly compute min λ E[λL] for any L ∈ L d (1). As we have seen above, this dimension reduction is extremely helpful when one looks for the ground state of E in L d , especially for d = 2, since we can plot L → min λ E[λL] in the fundamental domain D.

5.2.1.
A special power-law case. Let us now assume that for t > 0, s > d, α > β > d, and a, b > 0. Therefore, by (2.3) we have, for any λ > 0 and any L ∈ L d (1), that . For a fixed lattice L, the critical points of λ → E[λL] are the solutions of the following equation Solving this equation in closed form is impracticable out of a discrete set of parameter values. Correspondingly, comparing energy values is even more complicated than in the pure Lennard-Jones-type case, which is already challenging when treated in whole generality.
Having pointed out this difficulty, we now focus on some additional specifications of the parameters, allowing to proceed further with the analysis. We have the following.

5.2.2.
Numerical investigations of the special power-law case in 2d and 3d. We let t ∈ (0, 9/d) vary and fix a = b = 1, α = 12, β = 6, s = 9/t, so that F (r) = r t , ρ(r) = r −9/t , φ(r) = 1 r 12 − 1 r 6 . (5.6) Note that (5.4) holds under these assumptions. In two dimensions, by testing as t ∈ (0, 4.5) increases, we observe numerically the following: • If t ∈ (0, t 1 ), t 1 ≈ 1.605, then A 2 minimizes e * (see Figures 9 and 10); • If t ∈ (t 1 , t 2 ), where t 2 ≈ 1.633, then Z 2 is a local minimizer of e * but there seems to be no minimizer for e * (see Figure 11); • if t ∈ (t 2 , 4.5), there seems to be no minimizer for e * , and Z 2 is a saddle point (see Figure 12).  Similarly to the discussion of Subsection 5.1, for some choice of parameters, a square lattice seems to be locally minimizing the EAM energy, at least within the range of our numerical testing. In [5], we have identified a range of densities for which a square lattice is optimal at fixed density. This seems however to be the first occurrence of such minimality among all possible lattices, without a density constraint. Indeed, when minimizing among all lattices, the square lattice Z 2 usually happens to be a saddle point, see, e.g., Figure 1 for the Lennard-Jones case.  For t = 1.632 (left) the square lattice is a local minimizer of e * whereas A 2 is a local maximizer. For t = 2 (right) it seems that e * does not have any local minimizer and A 2 stays a local maximizer. In both cases, there is no global minimum.
When t → 0, since s = 9/t → ∞ and r t → 1 for fixed r > 0, it is expected that the global minimizer of E in L 3 converges to the one of E φ , which in turn is expected to be a FCC lattice. This is supported by our numerics for t < t 1 . Figure 13. Special power-law case (5.6) in three dimensions. Plot of t → e * (L) for L = D 3 (red), L = D * 3 (blue) and L = Z 3 (black) for t ∈ (0, 3). The graph on the right is a close-up of the two transitions at t 1 and t 2 .