**5. Numerical Examples**

In this section, we present some numerical examples illustrating the results of the previous sections. For simplicity we limit ourselves to the three particle problem. The examples show that the structure of the bifurcation diagrams is quite rich and complex. To construct the pictures in this section, we use the results of Theorem 3, in particular the symmetries given by (32), together with various numerical techniques to ge<sup>t</sup> full or detailed descriptions of the corresponding bifurcation diagram.

To compute approximations of the bifurcating solutions predicted by Theorem 3, one employs a predictor-corrector continuation method (cf. [24,25]). The bifurcation points off the trivial branch can be determined, by Theorem 3, from the solutions of the equation *μ*(*A*) = 0 (cf. (17), (33)). Secondary bifurcation points off nontrivial branches can be detected by monitoring the sign of a certain determinant. Once a sign change in this determinant is detected, the bifurcation point can be computed by a bisection or secant type iteration. After detection and computation of a bifurcation point, then one can use formulas (8)–(10) to ge<sup>t</sup> an approximate point on the solution curve from which the continuation of the bifurcating branch can proceed.

Any trivial or nontrivial computed solution (**x**∗, *A*∗) will be called *stable*, if the matrix (∇2*E* + *<sup>λ</sup>*∇<sup>2</sup>*g*)(**x**∗, *A*∗) (cf. (31)) is positive definite when restricted to the tangent space at**x**<sup>∗</sup> of the constraint of fixed area. Otherwise the point (**x**∗, *A*∗) will be called *unstable*. We recall that the tangent space at**x**<sup>∗</sup> of the constraint of fixed area is given by

$$\mathcal{M} = \left\{ (\mathfrak{x}, \mathfrak{y}, \mathfrak{z}) \,:\, \nabla \mathfrak{g}(\vec{\mathfrak{x}}^\*) \cdot (\mathfrak{x}, \mathfrak{y}, \mathfrak{z}) = 0 \right\} \dots$$

Please note that since this space depends on the point**x**∗, then except for the trivial branch where the solution is known explicitly, the stability of a solution can only be determined numerically.

In our first example we consider the Lennard–Jones potential (23) with *c*1 = 1, *c*2 = 2, *δ*1 = 12, and *δ*2 = 6. (We obtained similar results for other values of *c*1, *c*2 like those for argon in which *<sup>c</sup>*1/*<sup>c</sup>*2 = 3.4<sup>6</sup> Å6.) From equation (24) we ge<sup>t</sup> that the bifurcation point off the trivial branch is given approximately by *A*0 = 0.5877. For the case of Theorem 3 in which *a* = *b*, we show in Figure 2 a close-up of the bifurcating branch near this bifurcation point, for the projection onto the *A*–*<sup>a</sup>* plane. In this figure and the others, the color red indicates unstable solutions while the stable ones are marked in green. Please note that the bifurcation is of the trans-critical type. It is interesting to note that for an interval of values of the parameter *A* to the left of *A*0 in the figure (approximately (0.5855, 0.5877)), there are multiple states (trivial and nontrivial) which are stable, the trivial one with an energy less than the nontrivial state in this case. In Figure 3 we look at the same branches of solutions, again the projection onto the *A*–*<sup>a</sup>* plane, but for a larger interval of values of *A*. We now discover that there are two secondary bifurcation points (In Figure 3 there are bifurcations only corresponding to the values of *A*0 = 0.5877, *A*1 = 0.6251 and *A*2 = 0.6670. The apparent crossing of a branch of scalene triangles and the trivial branch is just an artifact of the projection onto the *A*–*<sup>a</sup>* plane.) at approximately *A*1 = 0.6251 and *A*2 = 0.6670, and once again we have multiple stable states (with different symmetries) existing for an interval of values of the parameter *A*. The branches of solutions bifurcating at these values of *A* correspond to stable scalene triangles. Once a branch of solutions is computed, we can use the symmetries (32) to generate other branches of solutions. In Figure 4 we show all the solutions obtained via this process, projected to the *abc* space (no *A* dependence). Figure 5 show the same set of solution but with the branch or axis of trivial solutions coming out of the page. The figures clearly show the variety of solutions (stable and unstable) for the problem (12) as well as the rather complexity of the corresponding solution set.

**Figure 2.** Bifurcation diagram for the *a* component vs. *A* for the system (13) in the case *a* = *b* and a Lennard–Jones potential. The points in green represent local minima of (12) while those in red are either maxima or none.

**Figure 3.** Bifurcation diagram for the system (13) in the case of a Lennard–Jones potential for a larger interval of values of *A*. There are secondary bifurcations into stable scalene triangles.

**Figure 4.** Solution set for the system (13) in the case of a Lennard–Jones potential without the *A* dependence.

For our next numerical example, we consider the Buckingham potential (25) with parameter values *α* = *β* = *γ* = 1 and *η* = 4, which satisfy (27). In this case, we have two bifurcation points off the trivial branch (which correspond to solutions of *μ*(*A*) = 0) at approximately *A*0 = 5.3154 and *A*1 = 74.2253. The trivial branch is stable for *A* ∈ (*<sup>A</sup>*0, *<sup>A</sup>*1) and unstable otherwise. Both bifurcations are into isosceles triangles, and both are of trans-critical type but with different stability patterns. In Figure 6 we show the solution set for the case *a* = *b*. The plot shows the dependence of the *a* and *c* components on the area parameter *A*. In Figure 7 we show the projection of this set onto the *c* vs. *A* plane where one can appreciate somewhat better the stability patterns at each bifurcation point, and that there is a turning point for *A* ≈ 46 on the branch bifurcating from *A*0. Please note that once again we have multiple stable states existing for an interval of values of the parameter *A*. No secondary bifurcations were detected in this case.

**Figure 5.** Solution set for the system (13) in the case of a Lennard–Jones potential without the *A* dependence with the branch of trivial solutions coming out of the page.

**Figure 6.** Dependence of the *a* and *c* components on the parameter *A* for the system (13) in the case *a* = *b* and for a Buckingham potential.

**Figure 7.** Projection onto the *c* vs. *A* plane of the set in Figure 6.

## **6. Final Comments**

The variety or type of solutions obtained from Theorems 3, 5–7, could be predicted generically from an analysis of the symmetries present in our problem as given by (32) and (46). However, such an analysis does not guaranty the existence of solutions with the predicted symmetries. It is the application of the Equivariant Bifurcation Theorem 1 that actually yields the result that such solutions exist. The generic analysis however is a preliminary step in identifying the spaces in which Theorem 1 can be applied. We should also point out that the results on the bifurcating branches in Theorems 3, 5–7 are global in the sense that the so-called Crandall and Rabinowitz alternatives in Theorem 1 hold. That is, any bifurcating branch is either unbounded, or it intersects the boundary of the domain of definition of the operator in the equilibrium conditions, or it intersects the trivial branch at another eigenvalue.

The results of this paper might be useful in the analysis of the more general and complex problem of arrays with many particles. As the total area or volume of such an array is increased, thus reducing its density, one might expect that locally situations similar to the ones discussed in this paper might be taking place in different parts of the array. It is interesting to note here that the existence of multiple stable configurations detected in the numerical examples of Section 5, opens the possibility for the existence of multiple equilibrium (local) states in a large molecular array, reminiscent of the bubble formation phenomena mentioned in the introduction. Thus, as further analysis either via molecular dynamics simulations or theoretically would be the questions as to whether the local results in this paper are related or can predict the initiation of cavitation bubbles of different sizes in actual fluid or gases.

**Author Contributions:** Conceptualization, P.N.M., M.L.S.; methodology, P.N.M., M.L.S.; validation, M.L.S.; simulation, M.L.S.; formal analysis and investigation, P.N.M., M.L.S.; writing original draft preparation, writing and editing, P.N.M.; funding acquisition, P.N.M., M.L.S.

**Funding:** This research was funded by the National Science Foundation of USA under the Partnership for Research and Education in Materials Program of the University of Puerto Rico at Humacao (Grant No. NSF–DMR–1523463).

**Acknowledgments:** The author would like to thank the academic editor and the three anonymous referees for their constructive comments that have improved the presentation of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
