Fluorite-type materials in the monolayer limit

By using Materials Project database materalsproject and pymatgen code pymatgen , we first extracted 94 fluorite-type materials (spacegroup of $Fm\bar{3}m$). These materials have been found in inorganic crystal structure database (ICSD). We excluded $f$-electron systems (lanthanide and actinide compounds), some of hydrogen compounds (Li₂NH, K₂PtH₄, Ca₂RhH_5.4, and Sr₂RhH₅), high-temperature phase materials (Cu₂S, Cu₂Se, and Al₂O), and oxides (PbO₂ and BiO₂). H₂S and H₂Se were also excluded because of a large discrepancy between experimental and calculated lattice constants. Zirconia (ZrO₂) and hafnia (HfO₂) in the fluorite-type structure are known to be stable above 2650 K and 2870 K, respectively shin2018zirconia , but these are studied in the present work.

We study 53 fluorite-type materials that consist of $AB_{2}$ (19), $A_{2}B$ (14), and others (20), where $A$ is a metallic element and $B$ is H, F, Cl, O, S, Se, and Te. The $AB_{2}$-type system includes hydrides (8), fluorides (7), chlorides (2), and oxides (2). The $A_{2}B$ materials consist of alkali metals $A=$ Li (4), Na (4), K (3), and Rb (3) and chalcogen atoms $B=$ O, S, Se, and Te. The other materials consist of alkali earth metals Be (2) and Mg (4) and other metallic elements Rh (2), Ir (2), Pt (4), Au (3), Ni (1), Co (1), and Pd (1).

We consider four crystal structures in the monolayer limit: 1T, 2H, 1T^′, and PCK structures. These are illustrated in Fig. 1(b). The 1T structure corresponds to a monolayer that is truncated from the (111) surfaces of 3D phase. The 1T and 2H structures have a hexagonal unit cell with lattice constant $a$. The atoms are located at $A(0,0,0)$, $B(0,a/\sqrt{3},h)$, and $B(a/2,a/(2\sqrt{3}),-h)$ for 1T structure, and $A(0,0,0)$, $B(0,a/\sqrt{3},\pm h)$ for 2H structure, where $2h$ is the thickness of the monolayer. For the case of $A_{2}B$ systems, $A$ and $B$ atoms should be interchanged. The 1T^′ structure has a rectangular unit cell with lattice constants of $a$ and $b$. This is regarded as a distorted 1T structure, as $A$ atoms are dimerized along $y$ direction and $B$ atoms follow the movement of $A$ atoms. The 1T^′ structure may be buckled along the $z$ direction. The PCK structure is regarded as a monolayer truncated from the (110) surface. $A$ atom in the center of rectangular cell is displaced to $z$ direction, i.e., $A(a/2,b/2,\delta)$. Each $A$ atom is surrounded by four $B$ atoms in the same planes, $z=0$ and $\delta$.

We used Quantum ESPRESSO (QE) code qe to perform density-functional theory (DFT) calculations. We used exchange-correlation energy functional within generalized-gradient approximation (GGA) parametrized by Perdew, Burke, and Ernzerhof (PBE) pbe . The electron-ion interactions were treated by using ultrasoft pseudopotentials in pslibrary1.0.0 dalcorso . The energy cutoff for wavefunction was set to be $E_{\rm cut}=\max(E_{A},E_{B})$ plus 20 Ry, where $E_{j}(j=A,B)$ is the suggested value for atom $j$, and the energy cutoff for charge density was set to be $10E_{\rm cut}$. Spin-polarized calculations were performed. A vacuum layer was set to be 15 Å for 2D systems. Convergence thresholds for the total energy in the self-consistent field calculations was set to be $10^{-12}$ Ry, and those for the total energy and forces for structure optimization were set to be $10^{-5}$ Ry and $10^{-4}$ a.u., respectively.

We first optimized the crystal structure by using smearing parameter of 0.015 Ry smearingMP and $k$-point distance $\Delta k$ smaller than 0.1 Å^-1. The lattice constant of $0.95a_{3D}/\sqrt{2}$ was assumed as an initial guess for 1T and 2H structures, where $a_{3D}$ is the lattice constant of the cubic phase and the factor of $0.95$ accounts for the in-plane contraction due to the lack of atoms along the out-of-plane direction. For 1T^′ structure, an initial structure was prepared by referring to the crystal structure of WTe₂ monolayer. For PCK structure, $b/a\simeq\sqrt{2}$ and $\delta=0.4a$ were assumed as an initial guess. Ferromagnetic phase was assumed in the initial spin configurations. For semiconducting systems, we optimized the structure again by using no smearing parameters and assuming $\Delta k\leq 0.15$ Å^-1.

The stability of $A_{n}B_{m}$ system is studied by calculating the formation energy

where $n$ and $m$ are integers, $\varepsilon_{\alpha}(A_{n}B_{m})$ is the total energy of $A_{n}B_{m}$ in the structure $\alpha$, and $\varepsilon(X)$ with $X=A,B$ is the total energy of element $X$. To obtain $\varepsilon(X)$, we extracted the ground state structure of $X$ from the Open Quantum Materials Database (OQMD) oqmd , as done in the construction of Computational 2D Materials Database (C2DB) C2DB2021 . We optimized the structure and calculated $\varepsilon(X)$.

For the bulk and the 2D system having negative formation energy and the lowest energy among four structures (1T, 2H, 1T^′, and PCK), we calculate the phonon dispersions within density-functional perturbation theory dfpt . The long-range Coulomb correction is included in the force constant matrix for semiconductors. The Coulomb interaction in the $z$ direction is truncated for 2D systems iso2D . The $8\times 8\times 1$ $q$ grid is used for the 1T and 2H structures, $6\times 4\times 1$ for the 1T^′ and PCK structures, while $4\times 4\times 4$ for the 3D structure. The variable for the acoustic sum rule is set to be “crystal” qe .

For the systems having NPR, first-principles molecular dynamics (MD) simulations are performed by using QE qe . A $4\times 4\times 1$ supercell is assumed, the volume of the unit cell is fixed, and the ionic temperature is kept to 300 K by adapting the velocity scaling. The Newton’s equation is integrated by using the Verlet algorithm with a time step of 1 fs.

To investigate the ionic character of 2D $AB_{2}$ and $A_{2}B$ systems, we estimate the Madelung constant $M$ of 1T and 2H structures. The $M$ is then characterized by $h/a$ only. The electrostatic energy is calculated by using pymatgen code pymatgen .

We calculate the surface energy $\gamma$ for the (111) and (110) planes of fluorite-type materials. This is defined as

where $E_{\rm slab}$ is the total energy of the slab including $n$ unit cells, $E_{\rm bulk}$ is the total energy of 3D bulk, and $A$ is the area of the surface. The factor of $1/2$ accounts for the presence of two surfaces on top and bottom sides of the slab. For the (111) and (110) surfaces, 7 and 9 layer-thick slabs (21 and 27 atoms) including a vacuum layer of 15 Å are assumed and the atomic position of the middle (fourth and fifth) layer is fixed in the geometry optimization. The slab models are constructed by using atomic simulation environment (ASE) package ase .

Before studying the 2D systems, we calculated phonon dispersions of 53 compounds in the fluorite structure to check whether the methodology used in the present work correctly predicts the dynamical stability of 3D counterparts. We have found that HgF₂, PIr₂, ZrO₂, and HfO₂ are unstable within PBE. When Perdew-Zunger pz and PBEsol pbesol functionals are used, HgF₂ and PIr₂ are dynamically stable.

Table 1, 2, and 3 lists the $E_{\rm form}$ of $AB_{2}$, $A_{2}B$ systems, and other ordered alloys, respectively, in the 1T, 2H, 1T^′, and PCK structures for 53 compounds. The $E_{\rm form}$ of bulk is also tabulated in these Tables. For $AB_{2}$ and $A_{2}B$ systems except for ScH₂, YH₂, TiF₂, PbF₂, and Li₂Te, 1T phase is the most stable 2D structure because 2H and PCK phases are higher in energy and 1T^′ phase is relaxed to 1T structure, as listed in Table 1 and 2. The $E_{\rm form}$ of the 1T phase is larger than of bulk by less than 0.2 eV, while hydrides tend to be less stable. These 1T phases are dynamically stable. However, some of them show imaginary frequencies around $\Gamma$ point in the Brillouin zone. This is attributed to an instability against the flexural (out-of-plane) vibrations. Phonon dispersions for the lowest energy phase are provided in Supplemental Materials SM .

The stability preference of ordered alloys is different from $AB_{2}$ and $A_{2}B$ systems. As listed in Table 3, the 1T^′ or PCK phases are preferred rather than 1T phase. However, many of them have positive $E_{\rm form}$, while they have negative $E_{\rm form}$ in the bulk. The systems having negative $E_{\rm form}$ are $X_{2}$Pt ($X=$ Al, Ga, In, and Sn) in the 1T^′ structure, PRh₂, Ga₂Au, and CoSi₂ in the PCK structure, and Al₂Pd in the 2H structure. These structures are all dynamically stable except for Sn₂Pt. The phonon dispersions are provided in Supplemental Materials SM .

We next calculate the Madelung constant $M$ to study the ionic character of $AB_{2}$ and $A_{2}B$ systems. The hydrides are excluded due to the metallic band structure. We assume the oxidation states of $B^{-}$ for halogens (F and Cl) and and $B^{2-}$ for chalcogens (O, S, Se, and Te). Figure 2 plots $M$ as a function of $h/a$ for $AB_{2}$ and $A_{2}B$ systems. The 1T structures have $h/a\in[0.15,0.32]$ and $M>2.3$, whereas the 2H structures have $M<2.3$. Therefore, the ionic character of 1T structure is more significant than of 2H structure.

Note that TiF₂ monolayer has a large value of $h/a=0.48$ and prefers 2H structure. It is interesting that 2H TiF₂ is more stable than the 3D bulk morita (see Table 1). The 2H structure has an indirect band gap of 1.3 eV (PBE) from K to $\Gamma$ in the Brillouin zone, while the 1T structure is metallic. This relationship between the crystal structure and electronic property is similar to that in MoS₂ monolayer magneMoS22018 . Such a structure-property relationship reflects the electronic configuration of transition metals. For Ti atom [Ar]$(3d)^{2}(4s)^{2}$, and a possible oxidation state is Ti²⁺. For Mo atom [Kr]$(4d)^{5}(5s)^{1}$, and a possible oxidation state is Mo⁴⁺ by assuming S^2-. The remaining two electrons in the Ti and Mo atoms should be occupied into $d$ energy levels. As the 2H structure has the D_3h symmetry, $d_{z^{2}}$ state becomes the lowest energy state. By considering the spin degeneracy, the $d_{z^{2}}$ state is completely filled, while the other $d$ states are empty, resulting in a finite band gap. On the other hand, the 1T structure has the O_h symmetry, and $d_{xy}$, $d_{yz}$, and $d_{zx}$ states become stable in energy. These states are partially filled with two electrons, resulting in a metallic phase.

As we have obtained stable structures of fluorite-type materials in the monolayer limit, we next explore physical properties of PCK structures below.

We study the elastic property of the PCK structures: ScH₂, YH₂, PbF₂, Li₂Te, PRh₂, Ga₂Au, and CoSi₂. The PCK structure is stretched along the $y$ direction. We define the out-of-plane strain as $\varepsilon_{z}=(d-d_{0})/d_{0}$, where $d$ and $d_{0}$ are the thickness of the monolayer with and without the strain, respectively. Figure 3(a) shows $\varepsilon_{z}$ as a function of $\varepsilon_{y}$. The atom position and the lattice constant $a$ are optimized, while the lattice constant $b$ is fixed. Interestingly, when $\varepsilon_{y}$ is increased, $\varepsilon_{z}$ also increases for Ga₂Au, PbF₂, and PRh₂ monolayers. This indicates an NPR in the out-of-plane direction, $\nu=-\partial\varepsilon_{z}/\partial\varepsilon_{y}<0$. As shown in Fig. 3(b), $\nu$ takes a minimum value of $-0.8$ at $\varepsilon_{y}=0.02$ in Ga₂Au monolayer. First-principles MD simulations show that they are thermodynamically stable at 300 K except for PRh₂, as shown in Fig. 3(c). For the case of 2D semiconductors (PbF₂ and Li₂Te), we plotted the potential energy and $\varepsilon_{z}$ as functions of $(a-a_{0})/a_{0}$ and $(b-b_{0})/b_{0}$, where $a_{0}$ and $b_{0}$ are the lattice constants in equilibrium. The NPR is observed only when the PbF₂ monolayer is stretched along the $y$ direction. This is provided in the Supplemental Material SM .

We have investigated whether the NPR also appears for thicker PbF₂ films. The PCK structure is constructed by stacking two PbF₂ layers [see Fig. 1(a)]. Therefore, we have constructed three (9 atoms) and four (12 atom) PbF₂ layers. The strain is increased from $\varepsilon_{y}=0$ to 0.06. However, no NPR is observed in these systems because the thickness monotonically decreases from 4.017 and 6.134 Å to 3.947 and 6.056 Å for three and four layers, respectively.