Ferroelectricity in strained Hf₂CF₂ monolayer

In general, MXenes are exfoliated from the precursor MAX phases by selective chemical etching. During the exfoliation process, the surfaces of MXenes are usually passivated by fluorine, oxygen, or other chemical ligands [9, 27]. Here, Hf₂C monolayer can be obtained from its parent bulk Hf₂AlC, which has been synthesized in experiment [28]. For fully fluorinated Hf₂C monolayer, i.e., Hf₂CF₂, there are four most possible configurations for chemical terminations, as proposed in Ref. [10].

According to our structural optimization, the model 1, as shown in Fig. S1 in Supplemental Materials (SM) [29], is dynamic unstable, which will spontaneously transform to the model 2. Therefore, only the models 2, 3, and 4 will be discussed in the following. The differences among these three configurations can be clearly distinguished from the side view, as shown in Fig. 1(a-c), which mainly lie in the positions of surface fluorine ions relative to the central carbon. The top view of model 2 is also illustrated in Fig. 1(d).

Basic physical properties of Hf₂CF₂ monolayer obtained from DFT calculations are summarized in Table 1. For comparison, the parent Hf₂AlC bulk is also calculated, whose lattice constants match the experimental values ($a$=$3.271$ & $c$=$14.363$ Å [28]) very well. Also our DFT results of Hf₂CF₂ monolayer are highly consistent with previous theoretical values (e.g., $a=3.264$, $3.232$, and $3.179$ Å for the models 2, 3, and 4 in Ref. [30]).

For the free standing Hf₂CF₂ monolayer, the model 2 owns the lowest energy among all considered configurations, while the model 3 is slightly higher in energy. Both models 2 and 4 are metallic, while the model 3 is a semiconductor. All these results agree well with previous first-principles studies [10, 30], which indicates the reliability of our calculations. It is noted that the space group of model 3 is a sub-group of model 2, by breaking the inversion symmetry along $c$-axis. Thus, the model 3 is polar, while the model 2 is nonpolar.

The proximate energies of the models 2 and 3 provide the possibility to pursuit the switching functionality of their contrastive physical properties (e.g., nonpolar vs polar, metallic vs semiconducting). Figure 2(a) shows the evolutions of models 2-4 as a function of in-plane lattice constant $a$, i.e., to simulate the biaxial strain. Obviously, the compressive strain is more advantageous for the models 3 and 4, since their optimized in-plane lattice constants are smaller than that of model 2. Indeed, the energy of model 3 becomes lower than that of model 2 when the biaxial strain is over $-1.3\%$. In contrast, the energy of model 4 is always higher than other two configurations within the considered range of biaxial strain. The phase transition between the models 2 and 3 is accompanied by a metal-semiconductor transition, with a sudden change of band gap from $0$ to $0.39$ eV. With further compression, the band gap of model 3 increases first and then decreases to $0$ gradually, as shown in Fig. 2(a).

Besides the factor of energy, the dynamic stability is also essential for structural transition. To determine the stability of Hf₂CF₂ monolayer, the phonon spectra are calculated, as shown in Fig. 2(b). No imaginary phonon frequency emerges in all three cases: unstrained model 2, strained models 2 and 3 (at $\varepsilon=-3\%$), implying their dynamic stability. In other words, the model 2 remains (meta)-stable in the strained condition, which is an evidence of the first-order transition.

To explore the strain effect on electronic properties of Hf₂CF₂ monolayer, the density of states (DOS) and band structures under different compressive strains are shown in Fig. 3. Initially, the model 2 exhibits metallic behavior, and the bands around the Fermi level is mostly contributed by Hf’s $5d$ orbitals [Fig. 3(a)]. Since the spin-orbit coupling (SOC) of $5d$ orbitals is prominent, the band structure is recalculated with SOC enabled, presented in Fig. 3(b) for comparison. However, the SOC effect is almost negligible and the system remains metallic. Under a moderate biaxial strain, the model 3 becomes the ground state, which opens an indirect band gap, as shown in Fig. 3(c). For comparison, the HSE06 band structure is also shown in Fig. 3(c), which is similar to the GGA one with a slightly enlarged band gap. The HSE06 band structure for unstrained model 2 can be found in Fig. S2 of SM [29], which does not alter its electronic structure too much. Thus, in the following, only the GGA calculations are presented.

For transition metal elements, partially filled $d$ orbitals may lead to magnetism. However, our DFT calculations do not find any local magnetic moment on Hf³⁺ ion, even if a moderate $U$ is applied. The absence of local magnetic moment is due to the large bandwidth of Hf’s $5d$ orbitals (see Fig. 3 for example).

In practice, strain is usually achieved by proper substrates. Here the MoS₂ is chosen as the substrate to mimic the strain effect in real heterostructures, because of their approximate lattice constant ($3.181$ Å for MoS₂, and $3.263$ Å for Hf₂CF₂) and similar structure (hexagonal in-plane symmetry). The slight smaller lattice constant of MoS₂ can provide a compressive strain to the attached Hf₂CF₂ monolayer, as requested to stabilize the model 3. Our DFT calculation shows that the model 3 configuration indeed has a lower energy than model 2 for Hf₂CF₂ monolayer on MoS₂ substrate. More details regarding the MoS₂ substrate can be found in SM (Fig. S3 and Table S1) [29].

According to the symmetry analysis, the structure of model 3 is polar but model 2 is nonpolar. For the model 3, the upper Hf is just above the lower F, while the lower Hf is not below the upper F. In contrast, the F and Hf columns are always synchronized in the model 2. Such in-plane shift of unilateral F layer breaks the inversion symmetry along the $c$-axis: the distances between Hf and C layers become different for the upper ($d_{1}$) and lower ($d_{2}$) ones. For example, at $\varepsilon$=$-1.9\%$, $d_{1}$=$1.333$ Å and $d_{2}$=$1.156$ Å, while the original $d_{1}$=$d_{2}$=$1.219$ Å, as visualized in Fig. 4(a). Then an out-of-plane polarization is induced. The atomic positions at $\varepsilon$=$-1.9\%$ for the models 2 and 3 are summarized in Table 2.

In this sense, such origin of ferroelectricity in the model 3 is similar to the so-called interlayer-sliding ferroelectricity [31], although here the sliding occurs at the surface layers, not between the vdW layers in their proposal. Due to the stronger interactions of F-Hf-C bonding, here a moderate polarization $7.8$ $p$C/m (equals to $1.48$ $\mu$C/cm² in the 3D case by considering its thickness) is obtained at $\varepsilon$=$-1.9\%$, larger than the vdW interlayer-sliding ferroelectric BN ($\sim 0.68$ $\mu$C/cm² [31]) and thus easier to be detected in future experiments. This ferroelectric polarization can be slightly enhanced with increasing compressive strain, as shown in Fig. 4(b). This tendency is also reasonable: stronger compressive strain, looser structure along the $c$-axis, which is advantageous for ferroelectric displacements.

The ferroelectric switching process is also simulated using the NEB method, as shown in Fig. 4(c). Since the paraelectric state is dynamically stable, the energy profile is a triple well one, instead of conventional double well. At $\varepsilon$=$-1.9\%$, the energy barrier between the ferroelectric state and paraelectric state is relative higher $\sim 410$ meV/f.u., although the energy difference between these two phases is only $17$ meV/f.u.. Noting that this value should be considered as a theoretical upper limit, since the real switching process probably finds an easier path than our NEB one. Even though, this value is already lower than that of Sc₂CO₂ ($\sim 520$ meV [15]), and comparable to the most extensively studied multiferroic BiFeO₃ ($427$ meV) [32]. Furthermore, with increasing compressive strain, the middle well becomes slightly shallower. The ferroelectric properties, including the polarization and switching process, are not affected by the MoS₂ substrate [29].

Since the paraelectric state (as the intermediate state of ferroelectric switching path) is metallic, it is necessary to clarify the feasibility of switching by electric field. Usually, the ferroelectric switching needs a full insulating path. However, for some improper ferroelectric systems, metallicity can be allowed to coexist with polarity [33]. In addition, the concept of “2D hyperferroelectric metals” was also proposed, with a polarization pointing out-of-plane [34]. Furthermore, Fei et al. experimentally realized the ferroelectric switching of a 2D metal WTe₂ using gate electrodes [35]. Thus, for 2D ferroelectrics, the switching path can be compatible with the metallicity, at least in the nanoscale limit.

Finally, it should be noted that the in-plane shifting of F layers will not induce any in-plane polarization, which is beyond the intuitional expectation and different from previous claim [15]. Both the paraelectric model 2 and ferroelectric model 3 have the $C_{3}$ rotation symmetry, which forbids any in-plane polarization. In fact, the displacement of F ion starts from a high-symmetric point to another high-symmetric point, which can create a polarization quantum and will not induce a real polarization [36].

Considering the $C_{3}$ (i.e., the in-plane triple rotation) symmetry of paraelectric phase, there are three possible shifting directions for each F plane, denoted as A/B/C, as shown in Fig. 5(a). These triple directions plus two F layers provide the $3\times 2$ degrees of freedom, similar to the hexagonal improper ferroelectric $R$MnO₃ and $R$FeO₃ [37]. In those hexagonal oxides, the triple choices of Mn or Fe trimerizations (denoted as $\alpha$, $\beta$, $\gamma$) and ferroelectric polarization ($+$ or $-$) generate six possible domains, which form the $\mathbb{Z}_{3}\times\mathbb{Z}_{2}$ topological antiphase domain walls (i.e. the $\alpha^{+}$-$\beta^{-}$-$\gamma^{+}$-$\alpha^{-}$-$\beta^{+}$-$\gamma^{-}$ vortices). The origin of such topological domain vortices in those hexagonal oxides is their Mexican-hat like energy landscape [38].

Then it is interesting to ask whether similar topological domain structure can appear in Hf₂CF₂ monolayer? To preliminarily answer this question, here the energy landscape as a function of F layer shifting vectors are calculated using the NEB method, as shown in Fig. 5(b-c). Unexpectedly, this energy landscape is totally different from the Mexican-hat one [38]. The energy barriers of direct paths between adjacent energy minima, i.e., those broken lines in Fig. 5(b-c), are rather high. Then all preferred domain walls will be between the ferroelectric model 3 and paraelectric model 2, i.e., those solid lines in in Fig. 5(b-c), instead of walls between antiphase domains of the model 3. Thus, no topological domain structure is expected here.

To further investigate the domain wall profiles, a ferroelectric domain wall between the A^- and B⁺ domains is studied, as shown in Fig. 6(a). Starting from a sharp domain boundary, the relaxed structure shows a continuous transition between polarization up and down, as shown in Fig. 6(b). Such continuous modulation of local dipoles is a result of the energy landscape [Fig. 5(c)], different from the very sharp domain wall predicted in WO₂Cl₂ which owns the quadruple energy wells [39]. Since the domain wall is close to the model 2 state, it is expected to be metallic, as shown in Fig. 6(c). The electronic state around the Fermi level is visualized in Fig. 6(d), which further confirms the metallicity of domain wall.

The metallicity of ferroelectric domain wall in Hf₂CF₂ is intrinsic and free from defects. For comparison, the previously studied conductive domain wall in BiFeO₃ [40, 41] is mostly due to the accumulation of charged defects and remains semiconducting (but with smaller local band gaps). Thus, the underlying physics involved in BiFeO₃ domain wall and Hf₂CF₂ domain wall is conceptually different. Such metallic wall provides the possibility to pursuit domain wall nanoelectronics [42].