**1. Introduction**

Consider a system of *N* molecules, modeled as identical spherical particles, enclosed in a bounded region B of R3. At any given instant (or in an equilibrium configuration), the total potential energy of the molecular array is given by:

$$E = \sum\_{i$$

where *φ* is the intermolecular potential energy with · the standard Euclidean or two-vector norm, and**r**1, ... ,**r***<sup>N</sup>* ∈ R<sup>3</sup> are the positions of the particles. More general energy potentials have been considered of which (1) is a special case (cf. [1,2]), or those based on the eigenvalues of adjacency matrices like in [3]. The problem of minimizing (1) subject to certain type of global or local conditions have been studied extensively (see e.g., [4–6] and the references there in). In these models either the array is infinite, with some local repeating structure, or finite but with *N* → ∞. In this paper, none of these conditions are required but we expect that our results can be extrapolated to such more general scenarios. Also, we do not commit to any particular intermolecular potential *φ* (but give examples for instance for Lennard–Jones type potentials) so that our results are applicable to any such smooth potential.

The particular problems that we consider in this paper are those of characterizing the minimum energy configurations of (1) in the case of three particles ( *N* = 3) arranged in a triangle and that of four particles (*N* = 4) in a tetrahedral array. The minimization problem is subject to the constraint of fixed area for the triangular array and of fixed volume in the tetrahedral case. We are particularly interested on the dependence of the minimizing states on the parameter of area or volume in the constraint. Both problems have the particularity that they can be formulated in terms of the intermolecular distances only, that is, without specifying the coordinates corresponding to the positions of the particles, thus substantially reducing the number of unknowns in each problem.

The motivation to study these problems comes from the following phenomena observed both in laboratory experiments and molecular dynamics simulations (see e.g., [7,8]). As the density of a fluid is progressively lowered (keeping the temperature constant), there is a certain "critical" density such that if the density of the fluid is lower than this critical value, then bubbles or regions with very low density appear within the fluid. This phenomenon is usually called "cavitation" and it has been extensively studied as well in solids. (See for instance [9,10] for discussions and further references.) When using discrete models of materials like (1), distinguishing between regions of low vs. high density, or whether bubbles or holes have form within the array, is not obvious since one is dealing essentially with a set of points. (See for instance [11,12] for the use of Voronoi polyhedra to study such arrays.) Thus, to study this phenomenon within this discrete model, one is naturally led to study or characterize the stability of homogeneous energy minimizing configurations of such an array, as the density of the array changes. The problems considered in this paper are the simplest problems within such a model. Our main contribution is on the application of global bifurcation theory (as opposed to just local) to study the set of equilibrium configurations for (1) under the stated constraints. In particular, to give specific conditions in terms of the intermolecular potential *φ* for the existence of nontrivial states.

In Section 3 we consider the problem of three particles. By Heron's formula for the area of a triangle, any three numbers (representing the intermolecular distances) that yield a positive value for the area formula, represent a triangle. In this case, we show that the functional (1) subject to the constraint of fixed area *A*, has for any value of *A*, a critical point representing an equilateral triangle. Moreover, in Theorem 2 we give a necessary and sufficient condition (cf. (22)), in terms of the intermolecular potential, for this equilibrium point to be a (local) minimizer of the energy functional. This condition leads to a set of values A for the area parameter *A* for which the equilateral triangle is a stable configuration. We give examples of how this set looks for various intermolecular potentials including the classical Lennard–Jones [13] and Buckingham [14] potentials, and those that model hard and soft springs including the usual Hooke's law.

Next in Section 3.2 we turn to the question of whether there exists other (not equilateral) equilibrium configurations for those values of the area parameter *A* for which the equilateral triangle becomes unstable, that is, when it ceases to be a local minimizer. To answer this question, we treat the system of equations characterizing the equilibrium points (cf. (13)) as a bifurcation problem with the parameter *A* as a bifurcation parameter, and the set of equilateral equilibrium configurations as the trivial solution branch. We find that the necessary condition for bifurcation from the trivial branch for this system occurs exactly at the boundary points of the set A given by the stability condition (22). To check the sufficiency condition for bifurcation, one must consider the linearization of the system (13) about the trivial branch at a boundary point *A*0 of A. However, since the kernel of this linearization is two-dimensional we cannot immediately apply the usual or standard results from bifurcation theory (cf. [15–18]). Because of the symmetries present in this problem (cf (32)), we can apply bifurcation equivariant theory (cf. [16,19]) to construct a suitable reduced problem corresponding to isosceles triangular equilibrium configurations. The linearization of the reduced problem at the point where *A* = *A*0 has now a one-dimensional kernel and provided that a certain transversality condition is satisfied (cf. (34)), we can show that there are three branches corresponding to isosceles triangles bifurcating from the trivial branch at the point where *A* = *A*0. Since the stability of these bifurcating branches can only be determined numerically (because one must linearize about an unknown solution), in Section 5 we construct numerically the bifurcation diagrams, with their respective stability patterns, for instances of the Lennard–Jones and Buckingham potentials. These examples show that the primary bifurcations off the trivial branch are of trans-critical type, and that at least for the Lennard–Jones potential, there are secondary bifurcations corresponding to stable scalene triangles. Moreover, there are intervals of values of the area parameter, for which there exists multiple stable states of the system for each value of *A* in such an interval.

In Section 4 we consider an array of four molecules in a tetrahedron. The general treatment in this case is similar to the three particle case but with two main differences. First the characterization of when six numbers (representing the lengths of the sides of the tetrahedron) determine a tetrahedron, is given in terms of the Cayley-Menger determinant and the triangle inequalities of one of its faces (cf. [20]). The next complication arises from the fact that the tetrahedron has 24 symmetries as compared to only six for the triangle! To deal with this many possibilities, once again we make use of the basic techniques of equivariant theory to ge<sup>t</sup> suitable reduced problems to work with. As the Cayley–Menger determinant is proportional to the volume of the corresponding tetrahedron, the volume constraint in our problem is basically that of setting this determinant to a given value *V* for the volume. In this case we show in Section 4.1 that the functional (1) subject to the constraint of fixed volume *V*, has for any value of *V* a critical point representing a regular or equilateral tetrahedron. Moreover, in Theorem 4 we give necessary and sufficient conditions (cf. (43)), in terms of the intermolecular potential, for this equilibrium point to be a (local) minimizer of the energy functional. As for the triangular case, these conditions determine a set of values V for the volume parameter *V* for which the equilateral tetrahedron is a stable configuration.

In Section 4.2 we consider the question of the existence of non-equilateral equilibrium configurations. The equilibrium configurations in this case are given as solutions to a nonlinear system of seven equations in eight unknowns (cf. (40)). We treat this system as a bifurcation problem with the parameter *V* as a bifurcation parameter, and the set of equilateral tetrahedrons as the trivial solution branch. The necessary condition for bifurcation from the trivial branch for this system occurs exactly at the boundary points of the set V given by the stability conditions (43). At a boundary point *V*0 of V there are two possibilities: the kernel of the linearization has dimension two or three. Using some of the machinery of equivariance theory as in [16], we can construct suitable reduced problems in each of these two cases, which enables us to establish the existence of non-equilateral equilibrium configurations and to ge<sup>t</sup> a full description of the symmetries of the bifurcating branches (cf. Theorems 5–7). As in the triangular case, the stability of this bifurcating branches can only be established numerically because one must linearize about the unknown bifurcating branch.

**Notation:** We let R*n* denote the *n* dimensional space of column vectors with elements denoted by **x**, **y**,... The inner product of **<sup>x</sup>**,**y** ∈ R*n* is denoted either by **x**,**<sup>y</sup>** or **x***t* **y**, where the superscript "*t*" denotes transpose. We denote the set of *n* × *m* matrices by R*<sup>n</sup>*×*m*. For *L* ∈ R*<sup>n</sup>*×*m*, we let ker(*L*) = {**x** ∈ R*n* : *L***x** = **0**} and Range(*L*) = {*L***x** : **x** ∈ <sup>R</sup>*<sup>n</sup>*}. For a function **F** : R*n* → R*<sup>m</sup>*, we denote its Fréchet derivative by D **F** which is given by the *m* × *n* matrix of partial derivatives of the components of **F**. If the variables in **F** are given by (**x**,**y**), then D**x F** denotes the derivative of **F** with respect to the vector of variables **x**, i.e., the matrix of the partial derivatives of the components of **F** with respect to the variables corresponding to**x**.

#### **2. Equivariant Bifurcation from a Simple Eigenvalue**

In this section, we provide an overview on some of the basic results on bifurcation theory from a simple eigenvalue for mappings between finite dimensional spaces, where the maps possess certain symmetries. The literature on this subject is extensive but we refer to [6,15,17] for details on the material presented in this section and further developments like for instance, the infinite dimensional case.

Let **F** : U × R → R*n* where U is an open subset of R*<sup>n</sup>*, be a *C*<sup>2</sup> function, and consider the problem of characterizing the solution set of:

$$\vec{\mathbf{F}}(\vec{\mathbf{x}},A) = \vec{\mathbf{0}}, \quad (\vec{\mathbf{x}},A) \in \mathcal{U} \times \mathbb{R}. \tag{2}$$

We assume that there exists a (known) smooth function**g**(·) such that:

$$\vec{\mathbf{F}}(\vec{\mathbf{g}}(A), A) = \vec{\mathbf{0}}, \quad \forall \ A.$$

The set T = {(**g**(*A*), *A*) : *A* ∈ R} is called the *trivial branch* of solutions of (2). We say that (**x**0, *<sup>A</sup>*0) ∈ T is a *bifurcation point* off the trivial branch T , if every neighborhood of (**x**0, *<sup>A</sup>*0) contains solutions of (2) not in T . If we let

$$L(A) = \mathbf{D}\_{\vec{\mathbf{x}}} \mathsf{F}(\vec{\mathbf{g}}(A), A)\_{\prime}$$

then by the Implicit Function Theorem, a necessary for (**x**0, *<sup>A</sup>*0) to be a bifurcation point is that *<sup>L</sup>*(*<sup>A</sup>*0) must be singular, a condition well known to be not sufficient.

In many applications of bifurcation theory and for the problems considered in this paper, the mapping **F** possesses symmetries due to the geometry of the underlying physical problem. The use of these symmetries in the analysis is useful for example to deal with problems in which dim ker(*L*(*<sup>A</sup>*0)) > 1. Thus, we assume that for a proper subgroup G of R*<sup>n</sup>*×*n*, characterizing the symmetries in the problem, the mapping **F** satisfies:

$$\vec{\mathbf{F}}(P\vec{\mathbf{x}},A) = P\vec{\mathbf{F}}(\vec{\mathbf{x}},A), \quad \forall P \in \mathcal{G}.\tag{3}$$

Let**v** ∈ ker(*L*(*<sup>A</sup>*0)) and define the *isotropy subgroup* of G at**v** by

$$\mathcal{H} = \{ P \in \mathcal{G} \; : \; P\vec{\mathbf{v}} = \vec{\mathbf{v}} \} \; , \tag{4}$$

and the H–*fixed point set* by

$$\mathbb{R}\_{\mathcal{H}}^{\boldsymbol{\mu}} = \left\{ \vec{\mathbf{x}} \in \mathbb{R}^{\boldsymbol{n}} \; : \; P\vec{\mathbf{x}} = \vec{\mathbf{x}} \; \forall P \in \mathcal{H} \right\}. \tag{5}$$

Clearly**v** ∈ <sup>R</sup>*<sup>n</sup>*H.

Let PH : R*n* → R*n* be a linear map that projects onto <sup>R</sup>*<sup>n</sup>*H that is Range(<sup>P</sup>H) = <sup>R</sup>*<sup>n</sup>*H and <sup>P</sup>H(R*<sup>n</sup>*H) = <sup>R</sup>*<sup>n</sup>*H. With UH = <sup>P</sup>H(U) = U ∩ <sup>R</sup>*<sup>n</sup>*H, we define**F**H : UH × R → <sup>R</sup>*<sup>n</sup>*H by:

$$\P\_{\mathcal{H}}(\vec{\mathbf{u}}, A) = \mathbb{P}\_{\mathcal{H}}\vec{\mathbf{F}}(\vec{\mathbf{u}}, A), \quad (\vec{\mathbf{u}}, A) \in \mathcal{U}\_{\mathcal{H}} \times \mathbb{R}. \tag{6}$$

An easy calculation now gives that

$$D\_{\vec{\mathbf{u}}}\mathbb{F}\_{\mathcal{H}}(\vec{\mathbf{u}},A) = \mathbb{P}\_{\mathcal{H}}\mathbb{F}\_{\vec{\mathbf{x}}}(\vec{\mathbf{u}},A)\mathbb{P}\_{\mathcal{H}}.$$

We assume that**g**(*A*) ∈ <sup>R</sup>*<sup>n</sup>*H for all *A*, so that

$$\vec{\mathbf{F}}\_{\mathcal{H}}(\vec{\mathbf{g}}(A)\_{\prime}A) = \vec{\mathbf{0}}\_{\prime} \quad \forall \ A.$$

It follows now that *<sup>L</sup>*H(*A*) : <sup>R</sup>*<sup>n</sup>*H → <sup>R</sup>*<sup>n</sup>*H is given by:

$$L\_{\mathcal{H}}(A) = \mathrm{D}\_{\overline{\mathfrak{u}}}\mathbb{H}\_{\mathcal{H}}(\overline{\mathfrak{g}}(A), A) = \mathbb{P}\_{\mathcal{H}}L(A)\mathbb{P}\_{\mathcal{H}}.$$

Clearly**v** ∈ ker(*<sup>L</sup>*H(*<sup>A</sup>*0)). The H-*reduced problem* is now given by:

$$
\vec{\mathbb{F}}\_{\mathcal{H}}(\vec{\mathbf{u}}, A) = \vec{\mathbf{0}}, \quad (\vec{\mathbf{u}}, A) \in \mathcal{U}\_{\mathcal{H}} \times \mathbb{R}. \tag{7}
$$

An important property relating (2) and (7) is that (**x**, *A*) ∈ UH × R is a solution of (2) if and only if (**x**, *A*) is a solution of (7). The following result provides the required sufficient conditions for (**x**0, *<sup>A</sup>*0) to be a bifurcation point of the H-reduced problem.

**Theorem 1** (Equivariant Bifurcation Theorem [19])**.** *Assume that for A* = *A*0 *there exists***v** ∈ ker(*L*(*<sup>A</sup>*0)) *that defines a proper isotropy subgroup* H *such that:*

$$\ker(L\_{\mathcal{H}}(A\_0)) = \text{span}\left\{\vec{\mathbf{v}}\right\}, \quad L\_{\mathcal{H}}'(A\_0)\vec{\mathbf{v}} \notin \text{Range}(L\_{\mathcal{H}}(A\_0)).$$

*Then there exists a branch* CH *of nontrivial solutions of* **<sup>F</sup>**H(**u**, *A*) =**0** *bifurcating from the trivial branch* T *at the point where A* = *A*0 *and such that either:*


The proof of this theorem is basically an application of a result from Krasnoselski [18] that uses the homotopy invariance of the topological degree. The three alternatives in the statement of the theorem are usually referred to as the Crandall and Rabinowitz alternatives. The local version of this result (cf. [6]) that is, without the Crandall and Rabinowitz alternatives, can be obtained via the Lyapunov–Schmidt reduction method. A useful consequence of this reduction is an approximate formula for the bifurcating branch in a neighborhood of the bifurcation point. Let ker(*<sup>L</sup>*H(*<sup>A</sup>*0)*<sup>t</sup>*) = span {**v**∗}, where **v**∗,**<sup>v</sup>** = 1, so that Range(*<sup>L</sup>*H(*<sup>A</sup>*0)) = **y** ∈ <sup>R</sup>*<sup>n</sup>*H : **v**∗,**<sup>y</sup>** = <sup>0</sup>. Now if we define

$$A^0 = \langle \vec{\mathbf{v}}^\*, (\mathcal{D}\_{\vec{\mathbf{u}}\vec{\mathbf{u}}} \mathcal{F}\_{\mathcal{H}}^0 \vec{\mathbf{v}}) \vec{\mathbf{v}} \rangle, \quad B^0 = \langle \vec{\mathbf{v}}^\*, L\_{\mathcal{H}}'(A\_0) \vec{\mathbf{v}} \rangle. \tag{8}$$

(here the zero superscripts mean evaluated at (**x**0, *<sup>A</sup>*0)), then the bifurcating branch have the following asymptotic expansion (cf. [16]):

$$\left( (\vec{\mathbf{x}}, A) = \left( \vec{\mathbf{g}} (A\_0 + \varepsilon) + \varepsilon m \vec{\mathbf{v}} + O(\varepsilon^2), A\_0 + \varepsilon \right), \tag{9}$$

where

$$m = -\frac{2B^0}{A^0}, \quad A^0 \neq 0. \tag{10}$$

#### **3. The Three Particle Case**

In this section, we consider the case in which the molecular array consists of three particles. The intermolecular energy is given by a smooth function *φ* : (0, ∞) → R called the *potential*. If (*a*, *b*, *c*) are the distances between the particles in the array, the total energy of the system is given by:

$$E(a,b,c) = \phi(a) + \phi(b) + \phi(c), \quad a, b, c \gg 0. \tag{11}$$

Also, the square of the area of the triangular array is given by Heron's formula:

$$s(a,b,c) \equiv s(s-a)(s-b)(s-c),$$

where *s* = (*a* + *b* + *c*)/2.

> For any given number *A* > 0, we consider the following constrained minimization problem:

$$\begin{cases} \min\_{a,b,c>0} E(a,b,c) \\ \text{subject to } \gcd(a,b,c) = A^2. \end{cases} \tag{12}$$

Thus, we are looking for minimizers of the energy functional *E* subject to the constraint that the area of the array is *A*. The first-order necessary conditions for a solution of this problem are given by:

$$\begin{cases} \begin{array}{rcl} \mathcal{g}(a,b,c) - A^2 &=& 0, \\ \nabla E(a,b,c) + \lambda \nabla \mathcal{g}(a,b,c) &=& \vec{\mathbf{0}}, \end{array} \end{cases} \tag{13}$$

where *λ* ∈ R is the Lagrange multiplier corresponding to the restriction *g*(*<sup>a</sup>*, *b*, *c*) = *A*2. For any given value of *A* > 0, this is a nonlinear system of equations for the unknowns (*<sup>λ</sup>*, *a*, *b*, *c*) in terms of *A*. In general this system can have multiple solutions depending on the characteristics of the potential *φ* and the value of *A*.

#### *3.1. Existence and Stability of Trivial States*

An easy calculation shows that:

$$
\nabla E = \begin{bmatrix}
\phi'(a) \\
\phi'(b) \\
\phi'(c)
\end{bmatrix}, \quad \nabla g = \frac{1}{4} \begin{bmatrix}
a(b^2 + c^2 - a^2) \\
b(a^2 + c^2 - b^2) \\
c(a^2 + b^2 - c^2)
\end{bmatrix}. \tag{14}
$$

Thus, the system (13) is equivalent to:

$$\begin{cases} \frac{1}{8}(a^2b^2 + a^2c^2 + b^2c^2) - \frac{1}{16}(a^4 + b^4 + c^4) - A^2 &=& 0, \\ \phi'(a) + \frac{\lambda}{4}a(b^2 + c^2 - a^2) &=& 0, \\ \phi'(b) + \frac{\lambda}{4}b(a^2 + c^2 - b^2) &=& 0, \\ \phi'(c) + \frac{\lambda}{4}c(a^2 + b^2 - c^2) &=& 0. \end{cases} \tag{15}$$

This system always has a solution with *a* = *b* = *c*. In fact, upon setting *a* = *b* = *c* in (15), this system reduces to:

$$\frac{3}{16}a^4 = A^2, \quad \phi'(a) = -\frac{\lambda a^3}{4}.\tag{16}$$

Thus, we have the following result:

**Lemma 1.** *For any value of A* > 0*, the system* (15) *has a solution of the form* (*<sup>λ</sup>A*, *aA*, *aA*, *aA*, *A*) *where:*

$$a\_A = \frac{2\sqrt{A}}{\sqrt[4]{3}}, \quad \lambda\_A = -\frac{4\phi'(a\_A)}{a\_A^3}.\tag{17}$$

We now characterize for which values of *A*, the solution provided in Lemma 16 is actually a minimizer, i.e., a solution of (12). To do this we need to examine the matrix 6∇2*E* + *<sup>λ</sup>*∇<sup>2</sup>*g*7 (**v***A*) where **v***A* = (*<sup>λ</sup>A*, *aA*, *aA*, *aA*), over the subspace:

$$\mathcal{M} = \left\{ (\mathbf{x}, \mathbf{y}, \mathbf{z}) \, : \, \vec{\nabla} \mathcal{g}(\vec{\mathbf{v}}\_A) \cdot (\mathbf{x}, \mathbf{y}, \mathbf{z}) = 0 \right\}. \tag{18}$$

A straightforward calculation gives that:

$$
\nabla^2 E \quad = \begin{bmatrix}
\phi^{\prime\prime}(a) & 0 & 0 \\
0 & \phi^{\prime\prime}(b) & 0 \\
0 & 0 & \phi^{\prime\prime}(c)
\end{bmatrix},
\tag{19}
$$

$$\nabla^2 \mathbf{g} \quad = \frac{1}{4} \begin{bmatrix} b^2 + c^2 - 3a^2 & 2ab & 2ac \\ 2ab & a^2 + c^2 - 3b^2 & 2bc \\ 2ac & 2bc & a^2 + b^2 - 3c^2 \end{bmatrix} . \tag{20}$$

*Symmetry* **2019**, *11*, 158

> It follows now that

$$
\begin{bmatrix}
\nabla^2 E + \lambda \nabla^2 \mathcal{g}
\end{bmatrix}
\begin{pmatrix}
\vec{\mathbf{v}}\_A \\
\vec{\mathbf{v}}\_A
\end{pmatrix} = 
\begin{bmatrix}
\phi''(a\_A) - \frac{\lambda\_A a\_A^2}{4} & \frac{\lambda\_A a\_A^2}{2} & \frac{\lambda\_A a\_A^2}{2} \\
\frac{\lambda\_A a\_A^2}{2} & \phi''(a\_A) - \frac{\lambda\_A a\_A^2}{4} & \frac{\lambda\_A a\_A^2}{2} \\
\frac{\lambda\_A a\_A^2}{2} & \frac{\lambda\_A a\_A^2}{2} & \phi''(a\_A) - \frac{\lambda\_A a\_A^2}{4}
\end{bmatrix},
\tag{21}
$$

and that

$$\mathcal{M} = \{ (\mathbf{x}, y, z) \; : \; \mathbf{x} + y + z = 0 \}$$

Using (17) one can show now that the matrix (21) is positive definite over M if and only if:

$$
\phi''(a\_A) + \frac{3}{a\_A} \phi'(a\_A) > 0. \tag{22}
$$

.

Thus, we have the following result:

**Theorem 2.** *Let φ* : (0, ∞) → R *be twice continuously differentiable function. Then the uniform array* (*a*, *b*, *<sup>c</sup>*)=(*aA*, *aA*, *aA*) *in Lemma 1 is a relative minimizer for the problem* (12) *for those values of A for which* (22) *holds.*

**Example 1.** *Consider the case of a potential that has the following form:*

$$\phi(r) = \frac{c\_1}{r^{\delta\_1}} - \frac{c\_2}{r^{\delta\_2}},\tag{23}$$

*where c*1, *c*2 *are positive constants and δ*1 > *δ*2 > 2*. (These constants determine the physical properties of the particle or molecule in question. The classical Lennard–Jones [13] potential is obtained upon setting δ*1 = 12 *and δ*2 = 6*.) For this function:*

$$\phi'(r) = -\frac{c\_1 \delta\_1}{r^{\delta\_1 + 1}} + \frac{c\_2 \delta\_2}{r^{\delta\_2 + 1}}, \quad \phi''(r) = \frac{c\_1 \delta\_1 (\delta\_1 + 1)}{r^{\delta\_1 + 2}} - \frac{c\_2 \delta\_2 (\delta\_2 + 1)}{r^{\delta\_2 + 2}},$$

*so that:*

$$
\phi''(r) + \frac{3}{r}\phi'(r) = \frac{c\_1\delta\_1(\delta\_1 - 2)}{r^{\delta\_1 + 2}} - \frac{c\_2\delta\_2(\delta\_2 - 2)}{r^{\delta\_2 + 2}}.
$$

*Since aA is directly proportional to* √*A (see* (17)*), we have that for* (23)*, the stability condition* (22) *holds if and only if A* < *A*0*, where A*0 *is determined from the condition:*

$$\frac{c\_1 \delta\_1(\delta\_1 - 2)}{a\_A^{\delta\_1 + 2}} - \frac{c\_2 \delta\_2(\delta\_2 - 2)}{a\_A^{\delta\_2 + 2}} = 0\_\prime$$

*from which it follows that:*

$$A\_0 = \frac{\sqrt{3}}{4} \left[ \frac{c\_1 \delta\_1 (\delta\_1 - 2)}{c\_2 \delta\_2 (\delta\_2 - 2)} \right]^{\frac{2}{\delta\_1 - \delta\_2}}.\tag{24}$$

*Thus,* (17) *is a (local) solution of* (12) *if and only if A* < *A*0*. We will show that for A* > *A*0 *there exist solutions that break the symmetry a* = *b* = *c.*

**Example 2.** *A Buckingham potential has the form ([14,21]):*

$$\Phi(r) = a \mathbf{e}^{-\beta r} - \frac{\gamma}{r^{\eta}},\tag{25}$$

*where α*, *β*, *γ*, *η are positive constants. Thus*

$$\phi'(r) = -a\beta \mathbf{e}^{-\beta r} + \frac{\gamma \eta}{r^{\eta + 1}}, \quad \phi''(r) = a\beta^2 \mathbf{e}^{-\beta r} - \frac{\gamma \eta (\eta + 1)}{r^{\eta + 2}}, \dots$$

*from which it follows that:*

$$
\phi''(r) + \frac{3}{r} \phi'(r) = \alpha \beta \left[\beta - \frac{3}{r}\right] \mathbf{e}^{-\beta r} - \frac{\gamma \eta \left(\eta - 2\right)}{r^{\eta + 2}}.
$$

*After rearrangement, the stability condition* (22) *is equivalent to:*

$$F(a\_A) > G(a\_A),\tag{26}$$

*where*

$$F(r) = \alpha \beta (\beta r - 3) e^{-\beta r}, \quad G(r) = \frac{\gamma \eta (\eta - 2)}{r^{\eta + 1}}.$$

*These functions generically look as in Figure 1 where for G we assumed that η* > 2*. Now clearly <sup>F</sup>*(*r*) < *<sup>G</sup>*(*r*) *for r sufficiently large. Thus generically we expect the set of values of A for which* (26) *is satisfied to be of the form* (*<sup>A</sup>*0, *<sup>A</sup>*1)*. Since F has a maximum at rm* = 4*β , a sufficient condition for this is that <sup>F</sup>*(*rm*) > *<sup>G</sup>*(*rm*)*, or after rearrangement that the coefficients and exponents in* (25) *satisfy:*

$$
\alpha \beta \mathbf{e}^{-4} > \gamma \eta (\eta - 2) \left(\frac{\beta}{4}\right)^{\eta + 1} . \tag{27}
$$

*To check this condition against the results in [21], we let D*, *R*, *ξ* > 0 *with ξ* > *η, and define*

$$
\alpha = \frac{D\eta \mathbf{e}^{\overline{\xi}}}{\overline{\xi} - \eta'} \quad \beta = \frac{\overline{\xi}}{R}, \quad \gamma = \frac{D\xi R^{\eta}}{\overline{\xi} - \eta}. \tag{28}
$$

*It follows that* (25) *is now given in terms of D*, *R*, *ξ by:*

$$\phi(r) = D \left[ \frac{\eta}{\xi - \eta} \mathbf{e}^{\xi \left(1 - \frac{r}{R}\right)} - \frac{\xi}{\xi - \eta} \left(\frac{R}{r}\right)^{\eta} \right].$$

*It is easy to check now that provided ξ* > *η* + 1*, then φ has negative minimum value at r* = *R. The results in ([21], Table 1, Page 202) show that the best fit of a normalized Buckingham potential to a normalized Lennard–Jones (12-6) potential (δ*1 = 12 *and δ*2 = 6 *in* (23)*) is achieved for ξ* = 14.3863 *and η* = 5.6518*. For these values one can check that* (28) *satisfy the inequality* (27) *independent of the values of D and R.*

**Example 3.** *Consider a potential of the form*

$$
\phi(r) = \frac{1}{2}kr^2 + \frac{1}{4}\,\beta r^4,\tag{29}
$$

*with k* > 0 *and β* ∈ R*. This potential corresponds to a Hook-type spring when β* = 0*, a hard spring if β* > 0*, and a soft spring if β* < 0*. (More general versions of* (29) *have been used in the study of the control of multi-agent systems, e.g., [22,23].) For this potential*

$$
\phi''(r) + \frac{3}{r}\phi'(r) = 4k + 6\beta r^2.
$$

*Please note that the stability condition* (22) *holds when β* ≥ 0 *independent of the value of A! That is, the symmetric state* (17) *is a minimizer for all values of A. In the case β* = 0 *is easy to show that this state is a global minimizer.*

*On the other hand, if β* < 0*, the stability condition holds if and only if A* < *A*0 *where*

*A*0= − *k* 2√3*β*.

**Figure 1.** Generic graphs of the functions *F* and *G* appearing in the stability condition (26) for a Buckingham potential.

#### *3.2. Existence and Stability of Nontrivial Solutions*

We say that solutions of (13) are *trivial* if *a* = *b* = *c* and call the set

$$\mathcal{T} = \left\{ (\lambda\_A, a\_A, a\_A, a\_A, A) \; : \; \lambda\_A, a\_A \text{ given by (17), } A > 0 \right\}, \tag{30}$$

the *trivial branch* parametrized by *A*. In this section, we show that there exist nontrivial solutions of (13) that bifurcate from the trivial branch.

If we let**x** = (*<sup>λ</sup>*, *a*, *b*, *c*) and **G** : R × (0, ∞)<sup>4</sup> → R<sup>4</sup> be the left-hand side of (15), then this system is equivalent to **G** (**x**, *A*) <sup>=</sup>**0**. An easy calculation gives that

$$D\_{\vec{\mathbf{x}}}\vec{\mathbf{G}}(\vec{\mathbf{x}},A) = \begin{bmatrix} 0 & (\vec{\nabla}\mathcal{g})^t \\ \vec{\nabla}\mathcal{g} & \nabla^2E + \lambda\nabla^2\mathcal{g} \end{bmatrix}.\tag{31}$$

If we evaluate now at the trivial branch (17), we ge<sup>t</sup> that

$$\mathbf{D}\_{\overrightarrow{\mathbf{x}}} \overrightarrow{\mathbf{G}}(\overrightarrow{\mathbf{v}}\_{A'} A) = \begin{bmatrix} 0 & \gamma & \gamma & \gamma \\ \gamma & \alpha & \beta & \beta \\ \gamma & \beta & \alpha & \beta \\ \gamma & \beta & \beta & \alpha \end{bmatrix} \prime$$

where**v***A* = (*<sup>λ</sup>A*, *aA*, *aA*, *aA*) and

$$\kappa = \phi^{\prime\prime}(a\_A) - \frac{\lambda\_A a\_A^2}{4}, \quad \beta = \frac{\lambda\_A a\_A^2}{2}, \quad \gamma = \frac{a\_A^3}{4}.$$

The eigenvalues of this matrix are *α* − *β*, which is a double eigenvalue with geometric multiplicity two, and the simple eigenvalues

$$\frac{1}{2}\left[\alpha + 2\beta \pm \sqrt{(\alpha + 2\beta)^2 + 12\gamma^2}\right] \dots$$

which are always nonzero. A pair of linearly independent eigenvectors corresponding to *α* − *β* is (0, −1, 1, <sup>0</sup>)*<sup>t</sup>*,(0, −1, 0, <sup>1</sup>)*<sup>t</sup>*. Since

$$
\alpha - \beta = \phi^{\prime\prime}(a\_A) - \frac{3}{4} \lambda\_A a\_A^2 = \phi^{\prime\prime}(a\_A) + \frac{3}{a\_A} \phi^{\prime}(a\_A) \lambda\_A
$$

the double eigenvalue *α* − *β* becomes zero exactly at the value *A*0 where the stability condition (22) fails by becoming zero. Thus according to standard theory of bifurcation theory, we can have either none, two or four branches of solutions of (15) bifurcating at the point where *A* = *A*0. We now show that there are exactly four branches bifurcating from such a point: the trivial branch and three branches corresponding to isosceles triangles.

To avoid the complications of dealing with the two-dimensional kernel of (31) when evaluated at the trivial branch at *A* = *A*0, we make use of the symmetries possessed by the mapping **G** . In particular, if we denote by G the subgroup of R4×<sup>4</sup> of permutations that permute just the *a*, *b*, *c* components of any **x** = (*<sup>λ</sup>*, *a*, *b*, *c*) ∈ R4, then

$$
\vec{\mathbf{G}}(P\vec{\mathbf{x}},A) = P\vec{\mathbf{G}}(\vec{\mathbf{x}},A), \quad \forall P \in \mathcal{G}.\tag{32}
$$

Please note that every permutation in G changes the eigenvectors (0, −1, 1, <sup>0</sup>)*<sup>t</sup>*,(0, −1, 0, 1)*<sup>t</sup>* of *α* − *β*. However, the eigenvector

$$\vec{\mathbf{v}} \equiv (0, -2, 1, 1)^t = (0, -1, 1, 0)^t + (0, -1, 0, 1)^t,$$

is unchanged by the proper subgroup of permutations H of G that permutes just the *b*, *c* components of any **x** = (*<sup>λ</sup>*, *a*, *b*, *c*) ∈ R4. Thus, H is the isotropy subgroup of G at**v**. The H fixed point set is given by:

$$\mathbb{R}^4\_{\mathcal{H}} = \left\{ (\lambda, a, b, b)^t \,:\, \lambda, a, b \in \mathbb{R} \right\}.$$

The projection PH : R<sup>4</sup> → <sup>R</sup>4H has matrix representation:

$$\mathbb{P}\_{\mathcal{H}} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} \end{bmatrix} / $$

and the H reduced problem is now:

$$
\check{\mathbf{G}}\_{\mathcal{H}}(\vec{\mathbf{u}}, A) \equiv \mathbb{P}\_{\mathcal{H}} \check{\mathbf{G}}(\vec{\mathbf{u}}, A) = \vec{\mathbf{0}}, \quad (\vec{\mathbf{u}}, A) \in \mathbb{R}\_{\mathcal{H}}^{4} \times (0, \infty).
$$

Since T ⊂ <sup>R</sup>4H × (0, <sup>∞</sup>), it follows that T is a branch of solutions for the H reduced problem. Also, since

$$\mathrm{D}\_{\vec{\mathbf{u}}} \vec{\mathbf{G}}\_{\mathcal{H}}(\vec{\mathbf{u}}, A) = \mathbb{P}\_{\mathcal{H}} \mathrm{D}\_{\vec{\mathbf{x}}} \vec{\mathbf{G}}(\vec{\mathbf{u}}, A) \mathbb{P}\_{\mathcal{H}}.$$

we have that *<sup>L</sup>*H(*A*) : <sup>R</sup>4H → <sup>R</sup>4H is given by:

$$L\_{\mathcal{H}}(A) = \mathrm{D}\_{\vec{\mathbf{u}}} \vec{\mathbf{G}}\_{\mathcal{H}}(\vec{\mathbf{v}}\_{A\prime} A) = \mathbb{P}\_{\mathcal{H}} \mathrm{D}\_{\vec{\mathbf{x}}} \vec{\mathbf{G}}(\vec{\mathbf{v}}\_{A\prime} A) \mathbb{P}\_{\mathcal{H}}.$$

Let

$$
\mu(A) = a - \beta = \phi^{\prime\prime}(a\_A) + \frac{3}{a\_A} \phi^{\prime}(a\_A). \tag{33}
$$

We now establish a result for the existence of bifurcating branches for the reduced problem.

**Theorem 3.** *Let μ*(*<sup>A</sup>*0) = 0 *and assume that*

$$\frac{\mathrm{d}\mu}{\mathrm{d}A}(A\_0) \neq 0.\tag{34}$$

*Consider the system* (13) *and its trivial branch of solutions* (30)*. Then from the point* (*<sup>λ</sup>*0, *a*0, *a*0, *<sup>a</sup>*0) ∈ T *bifurcate three branches of nontrivial solutions of* (13) *each corresponding to isosceles triangles.*

**Proof.** With the definitions and notation as above, a lengthy but otherwise elementary calculation shows that for any *A* > 0, *μ*(*A*) is a simple eigenvalue of *<sup>L</sup>*H(*A*) restricted to <sup>R</sup>4H with corresponding eigenvector**v** = (0, −2, 1, <sup>1</sup>)*<sup>t</sup>*. Thus

$$L\_{\mathcal{H}}(A)\vec{\mathbf{v}} = \mu(A)\vec{\mathbf{v}}, \quad \forall A > 0.$$

In particular ker(*<sup>L</sup>*H(*<sup>A</sup>*0)) = span{**v**}. If we differentiate with respect to *A* in the equation above and set *A* = *A*0, we ge<sup>t</sup> that

$$L'\_{\mathcal{H}}(A\_0)\vec{\mathbf{v}} = \mu'(A\_0)\vec{\mathbf{v}}.$$

Since *<sup>L</sup>*H(*<sup>A</sup>*0) is symmetric, we have that Range(*<sup>L</sup>*H(*<sup>A</sup>*0)) = **y** ∈ <sup>R</sup>*<sup>n</sup>*H : **v**,**<sup>y</sup>** = <sup>0</sup>. Thus, the hypotheses in Theorem 1 are satisfied if and only if *μ* (*<sup>A</sup>*0) = 0. Thus, we ge<sup>t</sup> a branch of solutions of the reduced problem, equivalently (15), bifurcating from the trivial branch at the point where *A* = *A*0. Since this branch belongs to <sup>R</sup>4H × R, we can use (32) to ge<sup>t</sup> that there exist two additional branches of solutions, one belonging to (*<sup>λ</sup>*, *b*, *a*, *b*, *A*)*<sup>t</sup>* : *λ*, *a*, *b*, *A* ∈ R and the other in (*<sup>λ</sup>*, *b*, *b*, *a*, *A*)*<sup>t</sup>* : *λ*, *a*, *b*, *A* ∈ <sup>R</sup>.

#### **4. Four Particles in a Tetrahedron**

We now consider the case of four particles arranged in a tetrahedron *T*. Let *a*, *b*, *c*, *A*, *B*, *C* be the distances between the particles where *a*, *b*, *c* denote the lengths of the edges joining a vertex of *T*, *A* the length of the edge opposite to *a*, *B* the length of the edge opposite to *b*, and *C* the length of the edge opposite to *c*. The six-tuple**a** = (*a*, *b*, *c*, *A*, *B*, *C*)*<sup>t</sup>* generates a tetrahedron ([20]) if and only if

$$\log(\vec{\mathbf{a}}) > 0, \quad A < B + \mathbb{C}, \quad B < A + \mathbb{C}, \quad \mathbb{C} < A + B,\tag{35}$$

where *g*(**a**) is given by the Cayley–Menger determinant:

$$g(\vec{\mathbf{a}}) = \begin{vmatrix} 0 & a^2 & b^2 & c^2 & 1 \\ a^2 & 0 & C^2 & B^2 & 1 \\ b^2 & C^2 & 0 & A^2 & 1 \\ c^2 & B^2 & A^2 & 0 & 1 \\ 1 & 1 & 1 & 1 & 0 \end{vmatrix}.\tag{36}$$

If we let R+ denote the set of positive real numbers, then we define

$$\mathcal{S} = \left\{ \vec{\mathbf{a}} = (a, b, c, A, B, \mathbf{C})^t \in \mathbb{R}\_+^6 \, : \, \text{(35) holds} \right\} \dots$$

Please note that S is open in R6. Moreover, any *regular* tetrahedron in which *a* = *b* = *c* = *A* = *B* = *C* > 0, is contained in S.


Since there are six permutations of the type *R* and four of the type *Q*, we ge<sup>t</sup> that there are 24 permutations of the form *P* = *RQ*. These 24 permutations form a subgroup R of the group of permutations of six letters. Also, it is easy to show that

$$\mathfrak{g}(P\vec{\mathfrak{a}}) = \mathfrak{g}(\vec{\mathfrak{a}}), \quad \forall P \in \mathcal{R}. \tag{37}$$

As the Cayley–Menger determinant is directly proportional to the square of the volume of the tetrahedron (cf. (39)), this identity simply states that the volume of the tetrahedron remains the same after rotations of the base and independent of which face we use as the base.

The total energy of the system of four particles is given now by:

$$E(\vec{\mathbf{a}}) = \phi(a) + \phi(b) + \phi(c) + \phi(A) + \phi(B) + \phi(C),\tag{38}$$

where the intermolecular potential *φ* is as before. For any *V* > 0 we consider the constrained minimization problem:

$$\begin{cases} \min\_{\vec{S}} E(\vec{\mathbf{a}})\\ \text{subject to } g(\vec{\mathbf{a}}) = 288V^2. \end{cases} \tag{39}$$

The constraint here specifies that the tetrahedron determined by**a** has volume *V* (cf. [20]). The first-order necessary conditions for a solution of this problem are given by (Since the inequality constraints in the definition of the set S are strict (non-active), the multipliers corresponding to these constraints are zero.): confirm. 

$$\begin{cases} \begin{array}{rcl} \mathcal{G}(\vec{\bf a}) - 288V^2 &=& 0, \\ \nabla E(\vec{\bf a}) + \lambda \, \nabla \mathcal{g}(\vec{\bf a}) &=& \vec{\bf 0}, \end{array} \end{cases} \tag{40}$$

which is now a nonlinear system for the seven unknowns (*<sup>λ</sup>*,**a**) in terms of the parameter *V*.

#### *4.1. Existence and Stability of Trivial States*

Expanding the determinant in (36) and computing its partial derivatives, we ge<sup>t</sup> that

$$\begin{array}{rcl}\nabla\_{\mathcal{S}}(\vec{a}) &=& 4\left[a(A^2(b^2+c^2+B^2+\mathcal{C}^2-2a^2-A^2)+(b^2-c^2)(B^2-\mathcal{C}^2)),b\right] \\ & b(B^2(a^2+c^2+A^2+\mathcal{C}^2-2b^2-B^2)+(a^2-c^2)(A^2-\mathcal{C}^2)), \\ & c(\mathcal{C}^2(a^2+b^2+A^2+B^2-2c^2-\mathcal{C}^2)+(a^2-b^2)(A^2-B^2)), \\ & A(a^2(b^2+c^2+B^2+\mathcal{C}^2-2A^2-a^2)-(b^2-\mathcal{C}^2)(c^2-B^2)), \\ & B(b^2(a^2+c^2+A^2+\mathcal{C}^2-2B^2-b^2)-(a^2-\mathcal{C}^2)(c^2-A^2)), \\ & C(c^2(a^2+b^2+A^2+B^2-2\mathcal{C}^2-c^2)-(a^2-B^2)(b^2-A^2))). \end{array}$$

Since ∇ *<sup>E</sup>*(**a**)=(*φ* (*a*), *φ* (*b*), *φ* (*c*), *φ* (*A*), *φ* (*B*), *φ* (*C*))*<sup>t</sup>*, the system (40) when evaluated at the regular tetrahedron**a** = (*a*, *a*, *a*, *a*, *a*, *<sup>a</sup>*), reduces to

$$4a^6 = 288V^2, \quad \phi'(a) + 4\lambda a^5 = 0.$$

Thus, we have the following result. **Lemma 2.** *For any V* > 0 *the system* (40) *has the solution* (*<sup>λ</sup>V*,**a***V*, *V*) *where***a***V* = (*aV*, *aV*, *aV*, *aV*, *aV*, *aV*) *and*

$$a\_V^3 = 6\sqrt{2}\,\mathrm{V}, \quad \lambda\_V = -\frac{\phi'(a\_V)}{4a\_V^5}.\tag{41}$$

We now examine the stability of the trivial state (41). A lengthy but otherwise elementary calculation shows that

$$H\_V \equiv \left[\nabla^2 E + \lambda\_V \nabla^2 \mathfrak{g}\right](\tilde{\mathfrak{a}}\_V) = \begin{bmatrix} a & \beta & \beta & 0 & \beta & \beta \\ \beta & a & \beta & \beta & 0 & \beta \\ \beta & \beta & a & \beta & \beta & 0 \\ 0 & \beta & \beta & a & \beta & \beta \\ \beta & 0 & \beta & \beta & a & \beta \\ \beta & \beta & 0 & \beta & \beta & a \end{bmatrix} \tag{42}$$

where

$$
\alpha = \phi''(a\_V) + \frac{3}{a\_V} \phi'(a\_V), \quad \beta = -\frac{2}{a\_V} \phi'(a\_V).
$$

Since ∇ *g*(**a***V*) = <sup>4</sup>*a*5*V*(1, 1, 1, 1, 1, <sup>1</sup>)*<sup>t</sup>*, we need examine the structure of *HV* on the subspace of R<sup>6</sup> given by

$$\mathcal{M} = \left\{ \vec{\mathbf{y}} \in \mathbb{R}^6 \, : \, \mathcal{Y}\_1 + \mathcal{Y}\_2 + \mathcal{Y}\_3 + \mathcal{Y}\_4 + \mathcal{Y}\_5 + \mathcal{Y}\_6 = 0 \right\}.$$

We have now the following result:

**Theorem 4.** *Let φ* : (0, ∞) → R *be twice continuously differentiable function. Then the matrix* (42) *is positive definite over* M *if and only if α* > 0 *and α* − 2*β* > 0*, which in turn are equivalent to*

$$
\phi''(a\_V) + \frac{3}{a\_V} \phi'(a\_V) > 0, \quad \phi''(a\_V) + \frac{7}{a\_V} \phi'(a\_V) > 0. \tag{43}
$$

.

*Thus, the regular tetrahedron***a***V is a relative minimizer for the problem* (39) *for those values of V for which conditions* (43) *hold.*

**Proof.** It is easy to check that M = range(*M*) where

$$M = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ -1 & -1 & -1 & -1 & -1 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}.$$

The matrix corresponding to the quadratic form of *HV* restricted to M is now given by *UV* = *<sup>M</sup>tHVM*. The eigenvalues of *UV* are:

$$\text{a (double)}, \quad a - 2\beta, \quad \frac{1}{2} \left[ 7a - 6\beta \pm \sqrt{16a^2 + 9(a - 2\beta)^2} \right],$$

Since the product of the last two of these eigenvalues is <sup>6</sup>*α*(*α* − <sup>2</sup>*β*), and 7*α* − 6*β* = <sup>3</sup>(*α* − 2*β*) + 4*α*, we can conclude now that they are all positive if and only if *α* > 0 and *α* − 2*β* > 0. Thus, *HV* restricted to M is positive definite provided these two conditions hold, which in turn implies that**a***V* is a relative minimizer for problem (39). That *α* > 0 and *α* − 2*β* > 0 are equivalent to (43) follows from the definitions of *α* and *β*.

**Example 4.** *For the Lennard–Jones potential* (23) *we have that:*

$$\begin{array}{rcl} \phi''(r) + \frac{3}{r}\phi'(r) &=& \frac{c\_1\delta\_1(\delta\_1 - 2)}{r^{\delta\_1 + 2}} - \frac{c\_2\delta\_2(\delta\_2 - 2)}{r^{\delta\_2 + 2}},\\ \phi''(r) + \frac{7}{r}\phi'(r) &=& \frac{c\_1\delta\_1(\delta\_1 - 6)}{r^{\delta\_1 + 2}} - \frac{c\_2\delta\_2(\delta\_2 - 6)}{r^{\delta\_2 + 2}}.\end{array}$$

*For simplicity, we assume δ*1 > 6*. We now have two cases:*

*1. Assume that δ*2 ∈ (2, 6]*. Then the second condition in* (43) *is automatically satisfied and the first condition holds if and only if V* < *V*0*, where V*0 *is determined from the condition (cf.* (41)*):*

$$\frac{c\_1 \delta\_1(\delta\_1 - 2)}{r\_0^{\delta\_1 + 2}} - \frac{c\_2 \delta\_2(\delta\_2 - 2)}{r\_0^{\delta\_2 + 2}} = 0, \quad r\_0 = a\_{V\_0}.$$

*from which it follows that:*

$$V\_0 = \frac{\sqrt{2}}{12} \left[ \frac{c\_1 \delta\_1 (\delta\_1 - 2)}{c\_2 \delta\_2 (\delta\_2 - 2)} \right]^{\frac{3}{\delta\_1 - \delta\_2}}.$$

.

.

*Thus, in this case the regular tetrahedron***a***V is a (local) solution of* (39) *if and only if V* < *V*0*.*

*2. If δ*2 > 6*, then the second condition in* (43) *holds if and only if V* < *V*1*, where V*1 *is determined from the condition:*

$$\frac{c\_1 \delta\_1(\delta\_1 - \theta)}{r\_1^{\delta\_1 + 2}} - \frac{c\_2 \delta\_2(\delta\_2 - \theta)}{r\_1^{\delta\_2 + 2}} = 0, \quad r\_1 = a\_{V\_1 \prime}$$

*from which it follows that:*

$$V\_1 = \frac{\sqrt{2}}{12} \left[ \frac{c\_1 \delta\_1 (\delta\_1 - 6)}{c\_2 \delta\_2 (\delta\_2 - 6)} \right]^{\frac{3}{\delta\_1 - \delta\_2}}$$

*Since δ*1 > *δ*2 > 6*, it follows that V*1 < *V*0*. Thus, in this case the regular tetrahedron* **a***V is a (local) solution of* (39) *if and only if V* < *V*1*.*

**Example 5.** *For the Buckingham* (25)*, we have that*

$$\begin{aligned} \boldsymbol{\phi}^{\prime\prime}(\boldsymbol{r}) + \frac{3}{r} \boldsymbol{\phi}^{\prime}(\boldsymbol{r}) &= \boldsymbol{\alpha} \boldsymbol{\beta} \left[ \boldsymbol{\beta} - \frac{3}{r} \right] \mathbf{e}^{-\beta \boldsymbol{r}} - \frac{\gamma \eta (\eta - 2)}{r^{\eta + 2}} \boldsymbol{\kappa} \\ \boldsymbol{\phi}^{\prime\prime}(\boldsymbol{r}) + \frac{7}{r} \boldsymbol{\phi}^{\prime}(\boldsymbol{r}) &= \boldsymbol{\alpha} \boldsymbol{\beta} \left[ \boldsymbol{\beta} - \frac{7}{r} \right] \mathbf{e}^{-\beta \boldsymbol{r}} - \frac{\gamma \eta (\eta - 6)}{r^{\eta + 2}} \boldsymbol{\kappa} \end{aligned}$$

*Please note that the first condition in* (43) *holds for an interval* (*<sup>V</sup>*0, *<sup>V</sup>*1) *of volume values under the conditions* (27) *in Example 2. The analysis now becomes rather complicated and we just describe it qualitatively. If η* > 6*, then the second condition in* (43) *would hold as well for values of V in an interval of the form* (*<sup>V</sup>*2, *<sup>V</sup>*3) *provided some condition similar to* (27) *holds. Depending as to whether or not the intersection* (*<sup>V</sup>*0, *<sup>V</sup>*1) ∩ (*<sup>V</sup>*2, *<sup>V</sup>*3) *is non-empty, we might get stable regular tetrahedrons. On the other hand, if η* ∈ (2, 6]*, then the second condition in* (43) *would hold as well for values of V in an interval of the form* (*<sup>V</sup>*4, ∞) *and again the existence of trivial states will depend on whether the corresponding intersection is non-empty.*

**Example 6.** *For the potential* (29)

$$
\phi''(r) + \frac{3}{r}\phi'(r) = 4k + 6\beta r^2, \quad \phi''(r) + \frac{7}{r}\phi'(r) = 8k + 10\beta r^2.
$$

*Please note that the stability conditions* (43) *holds when β* ≥ 0 *independent of the value of V! That is, the regular tetrahedron***a***V is a relative minimizer for the problem* (39) *for all values of V. In the case β* = 0*, since the functional* (38) *is convex, this state is a global minimizer.*

*On the other hand, if β* < 0*, the first condition in* (43) *holds if V* < *V*1 *and the second condition if V* < *V*2 *where*

$$V\_1^2 = -\frac{k^3}{243\,\beta^3}, \quad V\_2^2 = -\frac{8k^3}{1125\,\beta^3}.$$

*Since V*1 < *V*2 *we get that conditions* (43) *hold both for V* < *V*1*. For V* > *V*1 *either one or both conditions fail.*

#### *4.2. Existence and Stability of Nontrivial States*

Let **G** : R ×S× (0, ∞) → R<sup>7</sup> be given by the left-hand side of (40):

$$
\vec{\mathbf{G}}(\vec{\mathbf{x}}, V) = \begin{bmatrix}
\mathcal{g}(\vec{\mathbf{a}}) - 288V^2 \\
\nabla E(\vec{\mathbf{a}}) + \lambda \vec{\nabla} \mathcal{g}(\vec{\mathbf{a}})
\end{bmatrix} / 2
$$

where**x** = (*<sup>λ</sup>*,**a**). We have now that

$$D\_{\vec{\mathbf{x}}}\vec{\mathbf{G}}(\vec{\mathbf{x}},V) = \begin{bmatrix} 0 & (\vec{\nabla}\mathcal{g}(\vec{\mathbf{a}}))^t \\ \vec{\nabla}\mathcal{g}(\vec{\mathbf{a}}) & \nabla^2 E(\vec{\mathbf{a}}) + \lambda \nabla^2 \mathcal{g}(\vec{\mathbf{a}}) \end{bmatrix}. \tag{44}$$

It follows from (37) that

$$\begin{array}{rcl} \vec{\nabla} \mathcal{g} (P \vec{\mathbf{a}}) & = & P \vec{\nabla} \mathcal{g} (\vec{\mathbf{a}}), \end{array} \tag{45a}$$

$$\nabla^2 \mathcal{g}(P\vec{\mathfrak{a}}) \quad = \quad P\nabla^2 \mathcal{g}(\vec{\mathfrak{a}})P^t, \quad P \in \mathcal{R}\_\prime \tag{45b}$$

with similar relations for the total energy *E*. It follows now from (45a) that

$$\mathbf{G}(Q\vec{\mathbf{x}}, V) = Q\mathbf{G}(\vec{\mathbf{x}}, V), \quad Q \in \mathcal{G}, \tag{46}$$

.

.

where

$$\mathcal{G} = \left\{ Q = \left[ \begin{array}{cc} 1 & \vec{\mathbf{0}}^t \\ \mathbf{0} & P \end{array} \right] \, : \, P \in \mathcal{R} \right\}.$$

Thus, the system (40) remains the same, up to reordering of the equations, when**x** = (*<sup>λ</sup>*,**a**) is replaced by *Q***x**.

We now begin the analysis of the existence of solutions of the system (40) bifurcating from the trivial branch:

$$\mathcal{T} = \{ (\lambda\_{V'} \vec{\mathbf{a}}\_{V'} V) \; : \; \lambda\_{V'} \ \vec{\mathbf{a}}\_{V} \text{ given by (41), } V > 0 \}$$

If we evaluate (44) at the trivial state (*<sup>λ</sup>V*,**a***V*, *<sup>V</sup>*), then this matrix reduces to:

$$\mathbf{D}\_{\vec{\mathbf{x}}} \vec{\mathbf{G}}(\lambda\_{V}, \vec{\mathbf{a}}\_{V}, V) = \begin{bmatrix} 0 & (\vec{\nabla} \mathcal{g}(\vec{\mathbf{a}}\_{V}))^{t} \\ \vec{\nabla} \mathcal{g}(\vec{\mathbf{a}}\_{V}) & H\_{V} \end{bmatrix}^{t} \tag{47}$$

where ∇ *g*(**a***V*) = <sup>4</sup>*a*5*V*(1, 1, 1, 1, 1, 1)*<sup>t</sup>* and *HV* is given by (42). The matrix (47) has two eigenvalues which are nonzero for every value of *V* > 0, with the remaining eigenvalues given by:

1. *μ*1(*V*) = *φ* (*aV*) + 3*aV φ* (*aV*) with algebraic and geometric multiplicity three, and corresponding eigenvectors:

$$(0, -1, 0, 0, 1, 0, 0)^t, \quad (0, 0, -1, 0, 0, 1, 0)^t, \quad (0, 0, 0, -1, 0, 0, 1)^t. \tag{48}$$

2. *μ*2(*V*) = *φ* (*aV*) + 7*aV φ* (*aV*) with algebraic and geometric multiplicity two, and corresponding eigenvectors:

$$(0, -1, 1, 0, -1, 1, 0)^t, \quad (0, -1, 0, 1, -1, 0, 1)^t. \tag{49}$$

**Remark 1.** *Please note that the expressions for these eigenvalues are the ones that appear in Theorem 4 characterizing the stability of the trivial state* (*<sup>λ</sup>V*,**a***V*, *<sup>V</sup>*)*. Thus, the trivial state can change stability exactly when one of these two eigenvalues becomes zero.*

To deal with these kernels with dimension greater than one, we proceed as in the previous section by considering a suitable reduced problem in each case. These reductions are determined by the symmetries present in this problem which are embodied in (46).

4.2.1. The Eigenvalue *μ*1(*V*)

Let us take the eigenvector**g** = (0, 0, 0, −1, 0, 0, 1)*<sup>t</sup>* of the eigenvalue *μ*1(*V*) above. (The analysis for the other two eigenvectors is similar.) By inspection it is easy to ge<sup>t</sup> that the isotropy subgroup H of G at**g** is given by:

$$
\begin{array}{rcl}
\mathcal{H} &=& \left\{ \left( \begin{array}{cccc} \lambda & a & b & c & A & B & \mathbb{C} \\ \lambda & a & b & c & A & B & \mathbb{C} \end{array} \right) , \left( \begin{array}{cccc} \lambda & a & b & c & A & B & \mathbb{C} \\ \lambda & b & a & c & B & A & \mathbb{C} \end{array} \right) , \left( \begin{array}{cccc} \lambda & a & b & c & A & B & \mathbb{C} \\ \lambda & b & a & c & B & A & \mathbb{C} \end{array} \right) , \left( \begin{array}{cccc} \lambda & a & b & c & A & B & \mathbb{C} \\ \lambda & a & b & c & A & B & \mathbb{C} \\ \lambda & a & b & c & a & b & \mathbb{C} \end{array} \right) \right\}
\end{array}
$$

The H-*fixed point set* is now given by:

$$\mathbb{R}\_{\mathcal{H}}^{\mathcal{T}} = \left\{ (\lambda, a, a, c, a, a, \mathcal{C})^{t} \, : \, \lambda, a, c, \mathcal{C} \in \mathbb{R} \right\}. \tag{50}$$

.

The projection PH : R<sup>7</sup> → <sup>R</sup>7H has matrix representation:

$$\mathbb{P}\_{\mathcal{H}} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \sim$$

and the H reduced problem is now:

$$
\vec{\mathbf{G}}\_{\mathcal{H}}(\vec{\mathbf{u}}, V) \equiv \mathbb{P}\_{\mathcal{H}} \vec{\mathbf{G}}(\vec{\mathbf{u}}, V) = \vec{\mathbf{0}}, \quad (\vec{\mathbf{u}}, V) \in \mathbb{R}\_{\mathcal{H}}^{\mathcal{T}} \times (0, \infty).
$$

Since T ⊂ <sup>R</sup>7H × (0, <sup>∞</sup>), it follows that T is a branch of solutions for the H reduced problem. Also, since

$$\mathrm{D}\_{\vec{\mathbf{u}}} \vec{\mathbf{G}}\_{\mathcal{H}}(\vec{\mathbf{u}}, V) = \mathbb{P}\_{\mathcal{H}} \mathrm{D}\_{\vec{\mathbf{x}}} \vec{\mathbf{G}}(\vec{\mathbf{u}}, V) \mathbb{P}\_{\mathcal{H}}.$$

we have that *<sup>L</sup>*H(*V*) : <sup>R</sup>7H → <sup>R</sup>7H is given by:

$$L\_{\mathcal{H}}(V) = \mathbf{D}\_{\overline{\mathbf{u}}} \vec{\mathbf{G}}\_{\mathcal{H}}(\lambda\_V, \vec{\mathbf{a}}\_V, V) = \mathbb{P}\_{\mathcal{H}} \mathbf{D}\_{\overline{\mathbf{x}}} \vec{\mathbf{G}}(\lambda\_{V'}, \vec{\mathbf{a}}\_{V'}, V) \mathbb{P}\_{\mathcal{H}}.$$

It easy to check now that *μ*1(*V*) is a simple eigenvalue of *<sup>L</sup>*H(*V*) restricted to <sup>R</sup>7H with corresponding eigenvector**g** that is

$$L\_{\mathcal{H}}(V)\vec{\mathbf{g}} = \mu\_1(V)\vec{\mathbf{g}}\_{\mathcal{V}} \quad \forall \; V > 0.$$

We now have the result for the existence of bifurcating branches for the reduced problem. We omit the proof as it is similar to that of Theorem 3.

**Theorem 5.** *Let μ*1(*<sup>V</sup>*1) = 0 *and assume that* d*μ*1 d*V* (*<sup>V</sup>*1) = 0*. Then the system* (40) *has a branch of nontrivial solutions in* <sup>R</sup>7H × (0, ∞) *bifurcating from the trivial branch* T *at the point where V* = *V*1*, where* <sup>R</sup>7H *is given by* (50)*.*

**Remark 2.** *It follows from* (46) *that there are two additional branches of nontrivial solutions of the system* (40) *of the forms:*

$$\begin{aligned} \left\{ (\lambda, a, c, a, a, c, \mathsf{C}, a) \; : \; \lambda \in \mathbb{R}, \; a, c, \mathsf{C} > 0 \right\}, \\ \left\{ (\lambda, c, a, a, a, \mathsf{C}, a, a) \; : \; \lambda \in \mathbb{R}, \; a, c, \mathsf{C} > 0 \right\}. \end{aligned}$$

We now consider the case of the eigenvector **g** = (0, −1, −1, −1, 1, 1, 1)*<sup>T</sup>* of *μ*1(*V*). This eigenvector is obtained by adding the three eigenvectors in (48). By inspection, the isotropy subgroup H of G at **g** is given by those permutations in G that permute the symbols (*a*, *b*, *c*) and (*<sup>A</sup>*, *B*, *C*) in (*<sup>λ</sup>*, *a*, *b*, *c*, *A*, *B*, *C*) with the same permutation. (Thus, H has six elements.) The H–fixed point set is now given by:

$$\mathbb{R}\_{\mathcal{H}}^{\mathcal{T}} = \left\{ (\lambda, a, a, a, A, A, A)^t \, : \, \lambda, a, A \in \mathbb{R} \right\}. \tag{51}$$

The projection PH : R<sup>7</sup> → <sup>R</sup>7H has matrix representation:


It follows now that for the H reduced problem, *μ*1(*V*) is a simple eigenvalue with corresponding eigenvector**g**. The proof of the following result is as that of Theorem 3.

**Theorem 6.** *Let μ*1(*<sup>V</sup>*1) = 0 *and assume that* d*μ*1 d*V* (*<sup>V</sup>*1) = 0*. Then the system* (40) *has a branch of nontrivial solutions in* <sup>R</sup>7H × (0, ∞) *bifurcating from the trivial branch* T *at the point where V* = *V*1*, where* <sup>R</sup>7H *is given by* (51)*.*

**Remark 3.** *By applying all the transformations in* G*, it follows from* (46) *that there are three additional branches of solutions of the system* (40) *of the forms:*

$$\begin{aligned} & \left\{ (\lambda, A, A, a, a, a, A) \;:\; \lambda \in \mathbb{R}, \; a, A > 0 \right\}, \\ & \left\{ (\lambda, a, A, A, A, a, a) \;:\; \lambda \in \mathbb{R}, \; a, A > 0 \right\}, \\ & \left\{ (\lambda, A, a, A, a, A, a) \;:\; \lambda \in \mathbb{R}, \; a, A > 0 \right\}. \end{aligned}$$

*Thus, combining both theorems, we get that there are seven branches of nontrivial solutions bifurcating from the trivial branch* {(*<sup>λ</sup>V*,**a***V*, *V*) : *V* > 0} *at the value of V* = *V*1 *where μ*1(*<sup>V</sup>*1) = 0 *and <sup>μ</sup>* 1(*<sup>V</sup>*1) = 0*.*

4.2.2. The Eigenvalue *μ*2(*V*)

We now consider the case of the eigenvector**g** = (0, −2, 1, 1, −2, 1, 1)*<sup>T</sup>* of *μ*2(*V*). This eigenvector is obtained by adding the eigenvectors in (49). By inspection, the isotropy subgroup H of G at **g** is given by

$$
\begin{array}{ccl}
\mathcal{H} &=& \left\{ \left( \begin{array}{cccccccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & a & b & c & A & B & \mathbb{C} \\
\end{array} \right) ' \left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & a & c & b & A & \mathbb{C} \\
\end{array} \right) ' \\
\left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & A & \mathbb{C} & b & a & c & B \\
\end{array} \right) ' \left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & A & b & \mathbb{C} & a & B & \mathbb{C} \\
\end{array} \right) ' \\
\left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & A & B & c & A & B & \mathbb{C} \\
\end{array} \right) ' \left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & A & c & B & a & \mathbb{C} \\
\end{array} \right) ' \left( \begin{array}{cccc}
\lambda & a & b & c & A & B & \mathbb{C} \\
\lambda & A & c & B & a & \mathbb{C} \\
\end{array} \right) ' \right)
$$

with H–fixed point set given by:

$$\mathbb{R}\_{\mathcal{H}}^{\mathbb{Z}} = \left\{ (\lambda, a, b, b, a, b, b) \; : \; \lambda, a, b \in \mathbb{R} \right\}.\tag{52}$$

,

The projection PH : R<sup>7</sup> → <sup>R</sup>7H has matrix representation:

$$\mathbb{P}\_{\mathcal{H}} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} \\ 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} \\ 0 & \frac{1}{2} & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} \\ 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} \end{bmatrix}.$$

It follows now that for the H-reduced problem, *μ*2(*V*) is a simple eigenvalue with corresponding eigenvector**g**. The proof of the following result is as that of Theorem 3.

**Theorem 7.** *Let μ*2(*<sup>V</sup>*2) = 0 *and assume that* d*μ*2 d*V* (*<sup>V</sup>*2) = 0*. Then the system* (40) *has a branch of nontrivial solutions in* <sup>R</sup>7H × (0, ∞) *bifurcating from the trivial branch* T *at the point where V* = *V*2*, where* <sup>R</sup>7H *is given by* (52)*.*

**Remark 4.** *By applying all the transformations in* G*, it follows from* (46) *that there are two additional branches of solutions of the system* (40) *of the forms:*

$$\begin{aligned} \{ (\lambda, b, a, b, b, a, b) \;:\; \lambda \in \mathbb{R}, \; a, b > 0 \} \;. \\ \{ (\lambda, b, b, a, b, b, b, a) \;:\; \lambda \in \mathbb{R}, \; a, b > 0 \} \;. \end{aligned}$$
