Theory of dynamical stability for two- and three-dimensional Lennard-Jones crystals

Lattice vibration is of importance to understand the dynamical properties of solids. The crystal lattice is dynamically stable if the vibrational frequency is real over the entire Brillouin zone, while if a specific vibrational mode has an imaginary frequency, the crystal lattice is unstable against such a mode grimvall . Due to the development of first-principles approach such as density-functional theory (DFT) KS , the dynamical stability has been studied in a variety of materials kawazoe ; dfpt ; togo . For simple metals, the face-centered cubic (FCC) metals in the hexagonal close packed (HCP) structure and the HCP metals in the FCC structure are dynamically stable, i.e., metastable state. However, the FCC metals in the body-centered cubic (BCC) structure and vice versa are, in general, dynamically unstable grimvall . Also metastable phases have attracted significant interest in recent years due to their unconventional properties huang ; schonecker .

Recently, one of the authors has studied the dynamical stability of two-dimensional (2D) layers of simple metals using the first-principles approach ono2020_1 ; ono2020_2 . We have considered three-types of 2D lattice structure: The hexagonal (HX), the buckled honeycomb (bHC), and the buckled square (bSQ) structures, as shown in Fig. 1 ono2020_2 . It has been proposed that (i) If the HX structure is dynamically stable, the FCC and/or HCP structures are also dynamically stable; (ii) If the bHC structure is dynamically stable only, the BCC structure is dynamically stable; (iii) If the bSQ structure is dynamically stable only, the HCP structure is dynamically stable; and (iv) If the bHC and bSQ structures are dynamically stable, the HCP or FCC structures are stable depending on the group in the periodic table. It is interesting to study the structure stability relationship between 2D and 3D systems by using more simple models that enable us to derive the vibrational frequencies at some symmetry points analytically.

In this paper, we focus on the Lennard-Jones (LJ) crystals and investigate how the dynamical stabilities in the 2D structures (HX, bHC, and bSQ) and 3D structures (FCC, BCC, and HCP) are correlated with each other. We categorize the stability property into four groups (Fig. 2 below) and discuss the stability relationship between the LJ crystals and the elemental metals in the periodic table. We also provide analytical expressions for studying the dynamical stability of the HX, bSQ, and BCC LJ crystals.

Among many models, the LJ potential is one of the simplest model for describing the interatomic forces between atoms in materials. The LJ potential is a central potential that depends on the interatomic distance $R$ only and can describe the dynamics of rare gas solids mermin . For positive integer pairs $(m,n)$, the generalized LJ potential is given by

with $C_{m}$ and $C_{n}$ being the positive values. The first and second terms in Eq. (1) describe the repulsive and attractive potential between atoms. The inequality $m>n$ must hold since an atom cannot penetrate into the volume of other atoms. It has been known that the FCC lattice is dynamically stable for all $m>n$, while the BCC lattice is unstable except for small $m$ and $n$ born ; wallace . Modified LJ potentials taking the Friedel oscillations into account have been used to study the stability of alloys mihalkovic . The Morse potential that is also a central potential has been used to model the chemical bonds in alloys alsalmi . For an accurate description beyond the central potential approximation, many-body corrections are required, as implemented in such as the embedded atom method finnis ; foiles and DFT KS . Therefore, we expect that the investigation regarding the LJ potentials would provide a simple understanding for the dynamical stability of solids.

The main results are summarized as follows: For any sets of $(m,n)$, the FCC and HCP structures are dynamically stable, which is consistent with the previous studies born ; wallace . However, the HX structure is unstable, which is due to the zero thickness of this system, as discussed below. The stability properties for the BCC, bHC, and bSQ structures are presented in Fig. 2. It is categorized into four groups: (A) The bHC and bSQ structures are dynamically stable (blue); (B) The bHC is dynamically stable only (green); (C) The BCC is dynamically stable only (red); and (D) No stable structures are found (orange).

As an example for each group, we study the $(12,6)$, $(8,4)$, $(6,3)$, and $(9,3)$-LJ crystals in the HX, bHC, bSQ, FCC, BCC, and HCP structures. The optimized parameters for these structures are listed in Table 1 and 2. The $R_{0}$s satisfying $\phi^{\prime}(R_{0})=0$ are also listed in Table 1. Figure 3 shows the $R$-dependence of the (12,6), (8,4), (6,3), and (9,3)-LJ potentials. For the case $(m,n)=(12,6)$, the LJ potential takes the minimum at $R=R_{0}\simeq a_{0}$ and decays to zero monotonically as $R$ increases, which will create the short range interaction forces between atoms. As the values of $(m,n)$ decrease, the value of $R_{0}$ shifts to large $R$, as listed in Table 1, and the minimum value of the LJ potential becomes shallow. The curvature around $R$ being the lattice parameter is important to understand the instability of the BCC and 2D structures, which will be discussed below.

The dispersion curves for the $(12,6)$, $(8,4)$, $(6,3)$, and $(9,3)$-LJ crystals are shown in Figs. 4, 5, 6 and 7, respectively. The vibrational frequency is expressed in units of $\omega_{0}=\sqrt{E_{0}/(Ma_{0}^{2})}$. The imaginary frequency is represented as negative value of $\omega$. As shown in Figs. 4(a)-7(a), the instability of the HX structure is attributed to the appearance of the imaginary frequencies of the out-of-plane (ZA) mode over the Brillouin zone. With Eq. (8), the frequency of the ZA mode for the HX structure is expressed as

where the summation is taken over the lattice sites having the magnitude of $L=a_{0},\sqrt{3}a_{0},2a_{0},\cdots$. The dominant contribution is from the sites with $L=a_{0}$, but $\phi^{\prime}(a_{0})$ is always negative because the minimum of $\phi(R)$ is located at $R=R_{0}>a_{0}$, as listed in Table 1. The second largest contribution is from the sites with $L=\sqrt{3}a_{0}$, but the size of $\phi^{\prime}(\sqrt{3}a_{0})(>0)$ is not large enough to stabilize the ZA mode.

When the buckling is assumed as in bHC and bSQ, the lattice parameter is shortened compared to the HX structure, resulting in lower total energy, as listed in Table 1. No imaginary frequencies are observed within the Brillouin zone for the (12,6)-LJ crystals in the bHC and bSQ structure (see Fig. 4(b) and (c)) and the (8,4)-LJ crystal in the bHC structure (see Fig. 5(b)). For the (8,4)-LJ crystals, the bSQ structure is unstable against the ZA modes along the $\Gamma$-X direction. The (6,3) and (9,3)-LJ crystals in both the bHC and bSQ are unstable.

To understand this, we derive analytical expressions for the vibrational frequencies at the point X for the bSQ structure. By considering up to the second nearest neighbor (NN) sites, one obtains an expression for the lowest frequency (doubly degenerate). The derivation is provided in Appendix A. The stability condition is given by

with $\tilde{A}=A(a_{\rm bSQ})+A(d_{12})$, $\alpha=a^{2}/(4d_{12}^{2})$, $\beta=\delta^{2}/d_{12}^{2}$, and $d_{12}^{2}=a^{2}/2+\delta^{2}$. For the case of $(m,n)=(12,6)$, the optimized lattice parameter ($a_{\rm bSQ}$) and the interatomic distance between the particle 1 and 2 ($d_{12}$) are almost equal to the location of the potential minimum, that is, $a_{\rm bSQ}\simeq R_{0}$ and $d_{12}=1.009a_{0}\simeq R_{0}$, as listed in Table 1. This leads to $\tilde{A}\simeq 0$ and $0<B(a_{\rm bSQ})\simeq B(d_{12})$, so that the inequality in Eq. (21) is easily satisfied. For the case of $(m,n)=(6,3)$, the size of $a_{\rm bSQ}=0.885a_{0}$ is smaller than $d_{12}=1.171a_{0}\simeq R_{0}$, leading to $A(a_{\rm bSQ})<0$ (see Fig. 3). In addition, since the LJ potential is shallow, the magnitude of $B(d_{12})(>0)$ is not so large enough to cancel negative $\tilde{A}$ in the same brackets in Eq. (21). This gives rise to the instability of the bSQ structure. This is due to the long-range nature of the (6,3)-LJ potential: The size of $a_{\rm bSQ}$ becomes small in order to increase the energy gain, whereas $a_{\rm bSQ}$ is shifted from $R_{0}$. Unfortunately, it was difficult to derive the stability condition analytically for the bHC structure because the off-diagonal elements in Eq. (7) are not zero.

The dynamical stability of the BCC structure shows opposite tendency: When the long and short range LJ potentials are used, the BCC structure is dynamically stable and unstable, respectively. As shown in Figs. 4(e)-7(e), the BCC is dynamically stable for the (6,3)-LJ crystal only, otherwise it is unstable against the vibrational modes along $\Gamma$-N direction. The lowest frequency at the point N is expressed by ono2019_jap ; ono2020_jap

where in Eq. (8) the force constants up to the second NN sites, $a_{1}=\sqrt{3}a_{\rm BCC}/2$ and $a_{2}=a_{\rm BCC}$, are considered. Since $A(a_{1})<0$, the positive value of $B(a_{2})$ must be large enough to overcome the negative contribution from $A(a_{1})$. When the set $(m,n)=(6,3)$ is used, the size of $a_{\rm BCC}$ is about $0.8a_{0}$, yielding $a_{1}\simeq 0.7a_{0}$ and $a_{2}=0.8a_{0}$, respectively. Since these are less than $R_{0}=1.187a_{0}$, a relation $B(a_{2})\gg|A(a_{1})|$ holds, so that the BCC becomes dynamically stable.

Finally, we discuss the stability relationship between the LJ crystals and the elemental metals in the periodic table. In the previous study ono2020_2 , it has been shown that most transition metals in the bHC and bSQ structures are dynamically stable, while those in the HX structure are unstable. In this respect, most transition metals are classified into group A in Fig. 2. Interestingly, the overall shape of the dispersion curves in the (12,6)-LJ crystals in the HX, bHC, and bSQ structures is quite similar to that in 2D Ni, Pd, and Pt (see Supplemental Material (SM) SM including Fig. 20 in Ref. ono2020_2 ). The HCP Co may be classified into group B because the bHC structure is dynamically stable only (see SM SM including Figs. 8 and 19 in Ref. ono2020_2 ). Recent DFT and empirical potential calculations have shown that alkali metals (Li, Na, and K) in the FCC, BCC, and HCP structures are dynamically stable takahashi ; nichol , while BCC Li is stabilized by anharmonic effects hellman . For 2D Li, Na, and K, the bHC structure is dynamically stable only ono2020_2 . This indicates that there are no counterparts for the LJ crystals in group C, showing the limitation of the central potential for describing lattice dynamics of simple metals correctly. Although the overall shape of the phonon dispersions in the other alkali metals (Rb and Cs), alkali earth metals, and noble metals in the bHC and bSQ structures (see SM SM including Figs. 11, 12, and 21 in Ref. ono2020_2 ) are similar to that in group B, the stability properties for either the HX, BCC, or HCP structures are different between the realistic metals and the LJ crystals: For Rb and Cs, BCC (HCP) structure is stable (unstable) takahashi , while for the alkali earth and noble metals, the HX structure is dynamically stable. The boundary between the groups B, C, and D in Fig. 2 would be corrected by considering the many-body interaction beyond the central potential approximation finnis ; foiles or by considering the Friedel-like oscillation effect simply nichol .

We have studied the lattice dynamics of 2D and 3D LJ crystals by diagonalizing the dynamical matrix numerically and analytically. For all $(m,n)$, the HX structure is unstable due to the zero thickness. The buckling is important for obtaining dynamically stable 2D crystals: The bHC and/or bSQ structures are dynamically stable if the interaction between the LJ particles is described by the short-range potential. We have categorized the stability property into four groups shown in Fig. 2 and shown that the stability property of the LJ crystals with short-range potential ($m+n\geq 17$ and $n\geq 6$) is similar to that of transition metals in the periodic table. More parameters for the potential function are required to model the dynamical stability of metals in a unified manner.