Two-dimensional ionic crystals: The cases of IA-VII alkali halides and IA-IB CsAu

The alkali halides, known as ionic crystals, have the NaCl-type or CsCl-type structure as the ground state. We study the structural, vibrational, and electronic properties of two-dimensional (2D) ionic crystals from first-principles. Two potential structures that are hexagonal and tetragonal are investigated as structural templates. Through phonon dispersion calculations, 8 and 16 out of 20 alkali halides in the hexagonal and tetragonal structures are dynamically stable, respectively. The electron energy gaps range from 6.8 eV for LiF to 3.9 eV for RbI and CsI in the tetragonal structure within the generalized gradient approximation. By considering the Madelung energy and the core-core repulsion, we propose a hard sphere model that accounts for the nearest-neighbor bond length and the cohesive energy of 2D alkali halides. The 2D CsAu in the tetragonal structure is also predicted to be stable as an ionic crystal including only metallic elements, showing a band gap of 2.6 eV that is higher than that of the 3D counterparts.


I. INTRODUCTION
The solids are classified into metals and insulators, and the latter is further classified into covalent, ionic, and molecular crystals. The covalent bonding has played an important role in creating two-dimensional (2D) crystals, such as graphene, hexagonal boron-nitride, and transition metal dichalcogenides [1]. The honeycomb lattice with zero or finite thickness is a structural template in 2D covalent crystals, while the interlayer coupling is described by a van der Waals interaction. As the 3D metals have the fcc and hcp structures, the 2D metals can form close packed structure of hexagonal or buckled honeycomb lattices [2,3], while most of them are left to be synthesized. The 2D ionic crystals have also attracted significant interest recently. Tikhomirova et al. predicted the stable structures of 2D NaCl by using density functional theory (DFT) and evolutionary approach, and synthesized 2D NaCl experimentally, where the stability of the 2D NaCl relies on the chemical interaction with the diamond substrate [4]. Shi et al. synthesized Na-rich phases of 2D Na-Cl crystals that are stabilized on graphene surface [5]. Zhao et al. reported the square lattice formation of 2D NaCl using molecular dynamics (MD) simulation [6]. Chen et al. found 2D EuS as an ionic crystal having ferromagnetism within DFT [7].
Several calculations within DFT have been performed to predict stable structures and interesting properties of 2D compounds including group-IV [8], III-V [8][9][10], III-VI [11], II-VI [12], IV-VI [13], V-IV [14,15], and V-IV-III-VI [14] elements. While the group-IV elements have perfect covalent bonds, any other compounds can have ionic bonds as well due to the electronegativity difference [12]. In this sense, the 2D IA-VII compounds have strong ionic bonds. More recently, Kumar et al. investigated the structural and electronic properties of 2D alkali halides except for the Cs-based compounds, whereas only the planar honeycomb structure was assumed [16]. * shota o@gifu-u.ac.jp The alkali halides as 3D ionic crystals can have the NaCl-type (B1) or CsCl-type (B2) structures in the ground state. Whereas the relative phase stability has been studied in detail using first-principles approach [17] and Landau free energy [18], such a preference can be simply determined by the interplay between the Madelung energy and the ratio of the ionic radii R = r > /r < , where r > and r < are the radii of the larger and smaller ions, respectively [19,20]. The Madelung constant of the CsCl-type structure (α B2 ≃ 1.7627) is larger than that of the NaCl-type structure (α B1 ≃ 1.7476), indicating that the former is more stable if the nearestneighbor interatomic distance is the same in both structures. On the other hand, when the value of R is increased, the stable structure changes from the CsCl-type to NaCl-type structure across the value of R = ( √ 3 + 1)/2 ≃ 1.37. This is speculated from the 3D packing of hard (i.e., impenetrable) spheres: the nearest-neighbor interatomic distance is √ 3r > (between the atoms with the same charges) rather than r > + r < (between the atoms with different charges) [21]. For the 2D ionic crystals, however, no simple models have been established for understanding the stable structures.
In this paper, we explore the stable structures of 2D alkali halides AB (IA-VII compounds). Two crystal structures with the hexagonal and tetragonal lattices are assumed, as shown in Fig. 1. These structures have been identified as the 2D CuI-type and SnSe-type structures, respectively, within a computational search based on a topology-scaling algorithm [22], and the latter structure has also been assumed in studying 2D ionic EuS crystals [7]. Based on the Madelung constant calculations and the 2D sphere packing, we predict that the tetragonal structure is more stable than the hexagonal structure for R ≤ √ 2 + 1. By using DFT, we confirm that the 2D alkali halides AB (A = Li, Na, K, Rb, Cs; B = F, Cl, Br, I) follow this trend, that is, the 2D alkali halides have the tetragonal phase except for the Li-based compounds. We demonstrate that 8 and 16 out of 20 AB in the hexagonal and tetragonal structures, respectively, are dynamically stable by performing phonon dispersion calculations. The 2D alkali halides have an electron en-TABLE I. The 2D AB identified to be dynamically stable as the formula of A1B1 or A2B2. These ABs exhibit semiconducting property except for III-V GaSb [9], IV-IV C [23], IV-V CBi, GeBi, SnBi, PbN, PbSb [15], and bimetallic systems M-M [24,25]. For the IV-IV compounds, the cases of A = B indicate graphene, silicene, germanene, and stanene for A = C, Si, Ge, and Sn, respectively.  [24,25] ergy gap that ranges from 6.8 eV for LiF to 3.9 eV for RbI and CsI in the tetragonal structure, and exhibit a relatively flat band near the Fermi level. We also predict that the 2D tetragonal CsAu, an ionic crystal including only metallic elements, have an energy gap of 2.6 eV that is higher than that of the 3D counterpart by a factor of more than two.
In Table I, we list the 2D ABs that have been identified to be dynamically stable within phonon calculations. We also note that the Computational 2D Materials Database (C2DB) [26], created from high-throughput DFT calculations, includes 225 AB compounds that are dynamically stable. The C2DB includes 8 alkali halides having the formula of A 2 B 2 (A = K, Rb, Cs; B = Cl, Br, I) except for Cs 2 I 2 . The stable structure is the same as the 2D CuBr-type structure [22]: the Cu and Br atoms form the square structure, where the Br atoms are displaced alternately above and below the plane of Cu atoms. The present work is expected to stimulate more investigations for understanding the stabilities of 2D ionic crystals.
The paper is organized as follows: in Sec. II we first discuss the stable structures based on a hard sphere model (Sec. II A), next provide the structural, vibrational, and electronic properties of 2D alkali halides through firstprinciples calculations (Sec. II B), and compare the cohesive energies within DFT to those within a hard sphere model (Sec. II C). The physical properties of 2D CsAu are provided in Sec. II D, and our conclusion is presented in Sec. III. The computational details are included in Appendix.

II. RESULTS AND DISCUSSION
A. Hard sphere model or square structures; and those monolayers are stacked along the surface normal to form the AB bonds. Such an ionic bonding along the z direction is expected to stabilize the 2D crystals against the out-of-plane vibrations. The basis vectors for the A1 and B1 atoms are τ A1 = (0, 0, −z A ) and τ B1 = (0, 0, z B ), respectively, while those for the A2 and B2 atoms are given by τ A2 = (0, a/ √ 3, z A ) and τ B2 = (0, a/ √ 3, −z B ) for the hexagonal structure and τ A1 = (0, 0, −z A ) τ B1 = (0, 0, z B ), τ A2 = (a/2, a/2, z A ), and τ B2 = (a/2, a/2, −z B ) for the tetragonal structure, respectively, with the lattice constant a. We define the thickness of the 2D materials as h = z A + z B . Figure 2 shows the h-dependence of the Madelung constant α j for the 2D structures j including the hexagonal ( Fig. 1(a)) and tetragonal ( Fig. 1(b)) structures, where the equality of z A = z B is assumed. For comparison, two other structures were also assumed: the buckled honeycomb structure with τ A = (0, 0, −z A ) and τ B = (0, a/ √ 3, z B ) and the buckled square structure with τ A = (0, 0, −z A ) and τ B = (a/2, a/2, z B ) [24]. The values of α j were calculated by using pymatgen code [27]. The hexagonal and tetragonal structures have the maxi- mum value of α j ≃ 1.621 at h/a = 1/ √ 3 and α j ≃ 1.682 at h/a = 1/ √ 2, respectively. The optimal values of h/a are obtained when the interlayer distance is equal to the interatomic distance between atoms A and B in each monolayer. The buckled honeycomb and square structures have the maximum values of α j ≃ 1.542 and 1.615 at h/a = 0, respectively, so that the buckled structures with finite h/a will be energetically unstable, which is in contrast to the structure-stability property of the 2D metals [3]. Since the Madelung energy is given by E = −α j e 2 /d with the electron charge e and the nearestneighbor distance d, the 2D ionic crystals are predicted to have the tetragonal structure when the values of d are the same for all js.
We next investigate a potential transformation from the tetragonal to hexagonal structures by considering the difference of the ionic radii between atoms A and B within a 2D sphere packing. For the tetragonal structure, the interatomic distance is equal to d = r > + r < . However, when r > = a/ √ 2, the smaller ions may not touch the larger ones, yielding d = √ 2r > . Such a situation occurs at a critical ratio given by In a similar manner, the critical ratio R hex of the hexagonal structure is given by R hex = √ 3(2 + √ 3) ≃ 6.464. For the alkali halides AB, the values of R are listed in Table II. The Li atom has a small ionic radius of 0.6Å, so that the R is larger than R tet except for LiF. Other alkali halides satisfies R < R tet . From the calculations of the α j and R, we predict that the Na-, K-, Rb-, and Csbased compounds (except for the Li-based compounds) prefer the tetragonal structure.

B. DFT and DFPT
To show the validity of the hard sphere model for the structure-stability relationship of 2D alkali halides, we performed DFT calculations to calculate the energy difference ∆E = E hex − E tet , where E hex and E tet indicate the total energy (per unit cell) of AB in the hexagonal and tetragonal structures, respectively. The DFT calculations were performed by using Quantum ESPRESSO (QE) [28]. The computational details are provided in Appendix. The values of ∆E for the alkali halides AB are plotted in Fig. 3. The negative values of ∆E are observed for the Li-based compounds, which is consistent with our model calculations.
We next studied the dynamical stability of 2D alkali halides by performing the phonon dispersion calculations within density-functional perturbation theory (DFPT) [29]. The phonon dispersions for 20 AB in the hexagonal and tetragonal structures are shown in Fig. 4. For the hexagonal and tetragonal structures, 8 and 16 out of 20 AB were found to be dynamically stable, respectively. Although E hex is larger than E tet for the Li-based compounds, the hexagonal structure was unstable because imaginary phonon frequencies were observed along Γ-M and Γ-K lines in the Brillouin zone (BZ). On the other   hand, the Li-based compounds in the tetragonal structure is dynamically stable except for the LiI. Figure 5 shows the electron band structure of 2D alkali halides in the hexagonal and tetragonal structures. Basically, the magnitude of the energy gap decreases when the atoms A and B are changed from Li to Na to K to Rb to Cs and F to Cl to Br to I, respectively, except for a few relationships such as between RbF (CsF) and RbCl (CsCl). A similar trend of the energy gap variation has also been reported in the 2D alkali halides having the planar honeycomb structure, where such a variation was correlated with an increase in the ionic radii [16].
The conduction band minimum (CBM) of the 2D alkali halides is located at the Γ point, while the valence band maximum (VBM) is located at along the Γ-K line for the hexagonal structure and at the Γ or M point or the Γ-X line for the tetragonal structure. Note that the dispersion curve around the Fermi level is relatively flat, that is, small k dependence of the single-particle energy. As suggested in the 2D II-V compounds such as SrSe [12], ferromagnetism may be induced by the strong Coulomb interaction when the amount of carriers is tuned. Table III lists the optimized lattice parameters used in the DFT and DFPT calculations. The values of a, z A , and z B for 2D AB increase with the sum of the ionic radii listed in Table II, while the increase in z A is more moderate than that in z B when the atom A is fixed. For the cases of NaF, KF, RbCl, and CsBr, the difference be-tween z A and z B is small or zero, so that the RbF, CsCl, and CsF show the relationship of z A > z B . Table III also lists the 2D bulk modulus B 2D (see Appendix for the definition) and the in-plane elastic constants c 11 and c 12 [2,25] by assuming the tetragonal structure. The positive value of ∆ indicates the elastic stability of 2D AB. The B 2D , c ij , and ∆, in turn, decrease with the sum of the ionic radii. The LiF has the largest values of B 2D and c 11 and the small c 12 , resulting in the largest value of ∆. Table IV lists the values of d and the cohesive energies of E within the DFT and the hard sphere model for the tetragonal structure, where we assumed d 2 DFT = a 2 /2 + (z A − z B ) 2 because the optimized structure may satisfy z A = z B . The agreement of d between the DFT and model calculations is good, where an averaged relative error is 1.3 %. However, the Madelung energy of E coul sp (r) = −αe 2 /r at r = d sp overestimates the cohesive energy of E DFT . To remedy it, we include the core-core repulsive forces due to the Pauli principle into the total energy per ion pair as [20] E tot sp (r) = E coul sp (r) + C r m . By assuming that E tot sp takes a minimum value at r = d sp , the parameter C can be expressed as a function of d sp , and the total energy is formally expressed as

C. DFT versus hard sphere model
hence the cohesive energy is smaller by a factor of (m − 1)/m. The value of m is related to the B 2D = S(∂ 2 E/∂S 2 ) 0 , where S is the total area of the surface, and the subscript "0" indicates the derivative at the equilibrium. For the tetragonal structure, the m is written as (see Appendix for the derivation) In the present work, we calculated the B 2D by using DFT, and estimated the value of m using Eq. (4). As listed in Table IV, the corrected total energies given in Eq. (3) agree with the DFT results, where an averaged relative error was reduced from 28 % for E coul sp to 10 % for E tot sp . Note that the magnitude of m estimated in 2D AB is smaller than that estimated in 3D AB; for example, m = 5.88, 6.73, 7.24, and 7.06 for LiF, LiCl, LiBr, and LiI in the NaCl-type structure, respectively [20]. This implies that the core-core repulsive forces are long-ranged in the 2D systems. The time-evolution of the total energy (per atom) at T = 300 K and 500 K for 2D CsAu in the tetragonal structure. Right: top and side views for the atomic distribution after the first-principles MD simulation of 5 ps for T = 500 K.

D. CsAu
As an application of the structural templates in Fig. 1, we investigate the stability of 2D CsAu (IA-IB compound). It is known that the Cs and Au atoms form an ionic crystal in the CsCl-type structure. The CsAu shows an optical bandgap of 2.6 eV that was extracted from the photoemission spectroscopy experiments [30] and an indirect bandgap of 0.9 eV (from the point X to R in the Brillouin zone) that was predicted by DFT calcula-tions [31,32]. The MD simulations have predicted that the liquid CsAu has a frequency gap between the longitudinal and transverse optical branches in the limit of long wavelength, also showing a non-metallic character [33]. The 2D CsAu in the buckled honeycomb structure has been predicted to be a wide-gap semiconductor, although the free-standing structure is unstable against the out-of-plane deformations [25]. Such an instability will be avoided by adapting the structural templates for 2D ionic crystals.
The CsAu compound has the ratio of R = 1.23 < R tet by assuming r = 1.37Å for the Au atom [20], so that the tetragonal structure is preferred. In fact, the value of ∆E was estimated to be 0.15 eV/cell within the DFT. The optimized lattice parameters in units ofÅ are as follows: (a, z A , z B ) = (6.15, 1.79, 1.83) and (5.16, 1.82, 1.84) for the hexagonal and tetragonal structures, respectively. For the tetragonal structure, we obtain d DFT = 3.65Å that is larger than d sp = 3.06Å. Figure 6 shows the phonon dispersions of 2D CsAu in the hexagonal and tetragonal structures, showing that the tetragonal structure is dynamically stable. For the tetragonal phase, we performed first-principles MD simulations by assuming the temperature at T = 300 K and 500 K, and confirmed that 2D CsAu is thermodynamically stable as well (see Fig. 7). The electron band structure for the 2D CsAu is also shown in Fig. 6 (bottom panels). The indirect band gap from around K to M point is observed for the hexagonal structure, while the direct band gap at X point is observed for the tetragonal structure. For the hexagonal and tetragonal structures, the magnitude of the energy gap is estimated to be 3.1 eV and 2.6 eV within the generalized gradient approximation [34], respectively, which are higher than that of the 3D CsAu by a factor of more than two. When the spin-orbit coupling (SOC) is included, the lattice parameters of z A and z B were almost the same as those without the SOC, while the flat bands around ε ≃ −1 eV split into two bands that are separated by 1.5 eV.

III. CONCLUSION
In conclusion, we investigated the dynamical stability of 2D alkali halides in the hexagonal and tetragonal structures by performing first-principles phonon calculations, and showed that 8 and 16 out of 20 AB were dynamically stable, respectively, indicating that the tetragonal structure can serve as a structural template for the free-standing 2D ionic crystals. The electron energy gaps range from 6.8 eV for LiF to 3.9 eV for RbI and CsI in the tetragonal structure, and a relatively flat band can be observed near the Fermi level. The cohesive energies using a phenomenological model based on the hard spheres with positive or negative charges are in agreement with those using DFT within an average error of 10 %. The tetragonal CsAu, an ionic compound including only metallic elements, was also predicted to be dynamically and ther- modynamically stable and show a band gap of 2.6 eV that is higher than that of the 3D counterparts. The present work is expected to stimulate more investigations for interesting combination of A and B atoms, yielding 2D materials that are partially ionic and partially covalent. Exploring suitable substrates [10,35,36] for the 2D ionic crystals and a phase transition between 2D structures, as an analog of the B1-B2 transition in 3D systems [17,18], is left for future work.

ACKNOWLEDGMENTS
This work was supported by JSPS KAKENHI (Grant No. 21K04628). A part of numerical calculations has been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo and the supercomputer "Flow" at Information Technology Center, Nagoya University. calculations were performed within the generalized gradient approximation of Perdew-Burke-Ernzerhof [34]. In the self-consistent field calculations, the ultrasoft pseudopotentials provided in pslibrary.1.0.0 [37] were used; 60 Ry and 600 Ry were used for the cutoff energies for the wavefunction and the charge density, respectively; and a 20×20×1 k grid was used [38]. The structures were optimized within an accuracy of 10 −5 Ry for the total energy and 10 −4 a.u. for the total force. The thickness of the vacuum layer was set to be 15Å to avoid the artificial interactions between the 2D crystals along the z direction. Spin-polarized calculations were performed in the structure optimization, whereas no magnetic moments were observed for the 2D crystals.
We defined the cohesive energy per ion pair as E DFT (AB) = ε a (A) + ε a (B) − ε tet (AB)/2, where ε a (X) and ε tet (AB) are the total energy of the atom X and the 2D crystal AB in the tetragonal structure, respectively. To calculate the ε a (X), we performed atom-in-a-box calculations of a single atom X in a unit cell with a volume of 15×15×15Å 3 within spin-unpolarized approximation.
The phonon dispersion calculations were done within DFPT [29], where a 4×4×1 q grid was used (4 and 6 q points in the BZ for the hexagonal and tetragonal structures, respectively), and the convergence parameter of tr2 ph was set to be 10 −14 . Due to the insulating characteristics of each system, the nonanalytic term accounting for the long-range Coulomb forces is included in the force constant matrix. For the electron band structure calculations, the improved tetrahedron method [39] was adapted for creating the charge density. To study the relativistic effect on the electron dispersions of the 2D CsAu, the SOC was included within the scheme described in Refs. [40,41]. The first-principles MD simulations were done at 300 K and 500 K for 5 ps, where the velocity scaling method and a time step of 1 fs were used. A 3×3×1 supercell was used for 2D CsAu in the tetragonal structure. Calculations were performed using QE [28].
We defined the 2D bulk modulus as where S is the total area of 2D ionic crystals. When there are N c cells in the crystal, S = N c s with s being the area per cell. By using the total energy per cell, e = E/N c , the bulk modulus is expressed by B 2D = s(∂ 2 e/∂s 2 ) 0 . For the tetragonal unit cell with the lattice constant a, s = a 2 and thus δs ≃ 2aδa. Using the finite difference method, we obtain B 2D ≃ e(s + δs) + e(s − δs) − 2e(s) 4(δa) 2 .
(A2) This is evaluated by the total energy calculations within DFT. We next express the B 2D by using the interatomic distance r. Assuming r = a/ √ 2, S = 2N c r 2 , and the total number of ions N = 4N c , we obtain where E tot sp is given by Eq. (2), and the derivative is evaluated at the equilibrium r = d sp . We then obtain the expression of m given by Eq. (4).
The in-plane elastic constants of c 11 (= c 22 ) and c 12 were calculated by using the formulation described in Refs. [2,25]. The elastic stability of the 2D AB was studied by using the determinant of the Hessian matrix of the strain energy, i.e., ∆ = c 2 11 − c 2 12 . We also investigated the dynamical stability of the 2D alkali halides in the planar honeycomb and square structures (i.e., the case of h = 0) using DFPT [29], where a 6×6×1 q grid (7 q points) and a 4×4×1 q grid (6 q points) were used, respectively. The optimized lattice parameters, cohesive energies, and phonon dispersions of 20 AB in the planar honeycomb and square structures are provided in Table V and Fig. 8. The 2D alkali halides were unstable except for the planar honeycomb LiF, KCl, KBr, RbCl, RbBr and planar square KF. Such an instabilities may be due to the absence of the ionic bonds along the surface normal, in contrast to the hexagonal and tetragonal structures depicted in Fig. 1.