Prediction of intrinsic two-dimensional ferroelectrics in In2Se3 and other III2-VI3 van der Waals materials

Interest in two-dimensional (2D) van der Waals materials has grown rapidly across multiple scientific and engineering disciplines in recent years. However, ferroelectricity, the presence of a spontaneous electric polarization, which is important in many practical applications, has rarely been reported in such materials so far. Here we employ first-principles calculations to discover a branch of the 2D materials family, based on In2Se3 and other III2-VI3 van der Waals materials, that exhibits room-temperature ferroelectricity with reversible spontaneous electric polarization in both out-of-plane and in-plane orientations. The device potential of these 2D ferroelectric materials is further demonstrated using the examples of van der Waals heterostructures of In2Se3/graphene, exhibiting a tunable Schottky barrier, and In2Se3/WSe2, showing a significant band gap reduction in the combined system. These findings promise to substantially broaden the tunability of van der Waals heterostructures for a wide range of applications.

F erroelectricity, a property of materials associated with the emergence of spontaneous electric polarization, has a wide range of technological applications, such as non-volatile memories, field effect transistors and sensors 1,2 . Previous studies of ferroelectric materials have mainly focused on complex oxides, such as ABO 3 perovskite compounds 3 . Driven by technological demand for device miniaturization, exploration of the ferroelectric properties of perovskite thin films has been made more intensively 1,4,5 . Separately, a rapidly increasing number of two-dimensional (2D) van der Waals materials have been discovered, exhibiting a rich variety of emergent physical properties [6][7][8][9][10] . These developments in principle may offer new and alternative opportunities for realizing ferroelectricity in the ultimate single-layer regime 11,12 , especially with regard to the most technologically relevant polarizability perpendicular to the film direction.
The existence of 2D ferroelectricity was predicted long time ago based on an idealized Ising model 13 , but realistic materials normally suffer from the fundamental constraint that ferroelectricity would disappear when the film thickness is below a critical value, due to the effects of surface energy, depolarizing electrostatic field and electron screening 4,5,[14][15][16] . In general, the emergence of electric polarization demands breaking of the structural centrosymmetry in the polarization direction. Yet in the pristine structures of all known 2D materials including the well-known graphene and transition-metal dichalcogenides, the projections of their atomic positions on the out-of-plane axis are exclusively centrosymmetric, seemingly excluding any possible out-of-plane polarization.
Here we present the discovery of a class of stable single-layer 2D ferroelectric materials based on III-VI compounds in the form of III 2 -VI 3 . Using first-principles density-function theory (DFT) calculations, we reveal that the ground state structures of an intrinsic prototypical In 2 Se 3 quintuple layer (QL) possess both spontaneous out-of-plane and in-plane electric polarization, which can be reversed via laterally shifting the central Se layer through readily accessible kinetic pathways with the assistance of a modest out-of-plane or in-plane electric field. Furthermore, we demonstrate tunability in the Schottky barrier height within an In 2 Se 3 /graphene junction and a significant band gap reduction in the van der Waals heterostructure of In 2 Se 3 /WSe 2 , in each case achieved by reversing the out-of-plane polarity of the In 2 Se 3 layer. These findings effectively classify the well-known III 2 -VI 3 compounds actually as the long-sought 2D ferroelectric materials.

Results
Structures of In 2 Se 3 layered phases. Before exploring the ferroelectric properties, we first systematically examine the detailed structures of In 2 Se 3 . Bulk In 2 Se 3 has been shown to exist in two layered crystalline phases named a and b (refs 17,18), formed by vertical stacking of two different types of covalently bonded 2D In 2 Se 3 QLs via weak van der Waals interactions (Fig. 1a,b). The van der Waals nature of the inter-QL force is supported by earlier experimental observations that few-layer In 2 Se 3 samples can be obtained by exfoliation 19,20 or chemical vapour deposition 21 . Each atomic layer in a QL contains only one elemental species, with the atoms in a given layer arranged in a triangular lattice. The five atomic layers in a QL then stack in the sequence of Se-In-Se-In-Se atomic layers. Despite extensive experimental studies of the bulk structures, the precise alignments of the atomic layers within the a and b phases are still controversial [17][18][19][21][22][23] .
In this study, as prerequisites we employed DFT calculations to explore all possible atomic configurations ( Fig. 1c-h and Supplementary Fig. 1) within a QL, including the ones derived from the most-common crystal structures, such as zincblende, wurtzite and face-centred cubic (fcc). From the calculated total energy versus lattice constant for each structure ( Supplementary  Fig. 2), we find that none of the zincblende, wurtzite or fcc structures is stable. In either the zincblende (ABBCC) or wurtzite (ABBAA) structure shown in Fig. 1c  FE-ZB 0 (ABBCA) and FE-WZ 0 (ABBAC) (Fig. 1f,g), with the total heights of one QL around 6.8 Å. The dynamical and thermal stability of each of the two structures is further examined by calculating the phonon band structures ( Supplementary  Fig. 4a,b), confirming the absence of imaginary phonon modes, and ab initio molecular dynamic simulations ( Supplementary  Fig. 5). In addition, we find that the highly symmetric fcc (ABCAB) structure ( Fig. 1e) is unstable, and a metastable structure, fcc 0 , can be derived by shifting the central Se layer slightly away from the ideal fcc positions (Fig. 1h). The total energy of this metastable fcc 0 structure is 0.057 eV per unit cell higher than the two degenerate ground state structures. The most stable FE-ZB 0 and FE-WZ 0 and the metastable fcc 0 structures are all semiconductors. Their calculated band structures are provided in Supplementary Fig. 6. Based on the structural results presented above, we also obtain that the two degenerate ground states of FE-ZB 0 and FE-WZ 0 have an in-plane lattice constant of 4.106 and 4.108 Å, respectively, while the metastable state of fcc 0 has an in-plane lattice constant of 4.048 Å. When compared with the experimentally measured lattice constants, we identify the ground states to be the a phase, while the fcc 0 state to be the b phase 19,22 . Furthermore, our detailed calculations confirm that the experimentally observed Raman active A1 mode undergoes a blue shift when the In 2 Se 3 structure transforms from the a phase to the b phase ( Supplementary Fig. 4) 19,24 . Additionally, we note that the metastable nature of the b phase is consistent with the experimental observation that it is reached at higher temperature from the a phase through a structural phase transition 17,25 .
In contrast, in a related recent DFT study, only one ground state was considered for the a phase, while the energetically higher fcc phase was identified to be the b phase 26 .
Ferroelectric nature of In 2 Se 3 . Next, we turn to the key prediction that each of the degenerate ground-state structures of the In 2 Se 3 QL is an intrinsic 2D ferroelectric material with both out-of-plane and in-plane electric polarization. As shown in Fig. 1f,g, the Se atoms in both the top and the bottom surface layers reside on the hollow sites of the respective second-layer In atoms, while each atom in the central Se layer is tetrahedrally coordinated by the two neighbouring In layers, with one Se-In bond connecting to one side vertically and three Se-In bonds to the other side. As a result, the interlayer spacing between the central Se layer and the two In layers is dramatically different, effectively breaking the centrosymmetry and providing the very underlying basis for the emergence of the spontaneous out-of-plane electric polarization. The calculated magnitudes of the electric dipoles for one QL In 2 Se 3 in the degenerate ground states are both around 0.11 eÅ per unit cell (calculated by HSE06, 0.094 eÅ per unit cell by generalized gradient approximation-Perdew-Burke-Ernzerhof (GGA-PBE)). In addition, each groundstate structure hosts two equivalent states with opposite electric polarizations, which only differ by the energetically degenerate positions of the central Se-layer atoms. Specifically, in the FE-ZB 0 structure illustrated in Fig. 2a, the atoms in the central Se layer are at the B sites vertically aligned with the lower In layer (left in Fig. 2a), and the resultant electric dipole points downwards; by moving the central Se layer to the neighbouring C sites aligned with the upper In layer (right in Fig. 2a), the resultant electric dipole points upwards. In addition, the a-phase FE-ZB 0 and FE-WZ 0 structures also have in-plane electric polarization due to the in-plane centrosymmetry breaking. The in-plane electric polarization is along the [110] direction defined by the lattice vectors as illustrated in the insets of Fig. 3b. The magnitude of the in-plane electric polarization is calculated to be 2.36 and 7.13 eÅ per unit cell for the FE-ZB 0 and FE-WZ 0 phases, respectively, using the Berry phase approach. This difference can be attributed to the ions that make the two structures deviate from the non-polar reference structure possessing different numbers of charge, as illustrated in Supplementary Fig. 10a. A critical issue for ultrathin ferroelectric materials is the depolarization effect. It is known that the ferroelectricity of conventional ferroelectric thin films is usually suppressed, as the films are thinner than a critical thickness due to the effects of a depolarizing field induced by uncompensated charges in the presence of metal electrodes 4,5,14,15 . To examine the influence of the depolarizing field on the stability of the ferroelectric phase of In 2 Se 3 , we performed calculations with supercells containing one QL of In 2 Se 3 sandwiched between two graphite electrodes in short-circuit, as illustrated in Supplementary Fig. 7a,b for the In 2 Se 3 layer in the ferroelectric FE-ZB 0 phase and non-polar fcc 0 phase, respectively. A detailed description of the calculations is provided in Supplementary Note 1. The calculated results confirm that the depolarization effects only slightly reduce the energy difference between the two phases of In 2 Se 3 , and the FE-ZB 0 phase is still more stable than the non-polar fcc 0 phase. The electrostatic potential plot, shown in Supplementary Fig. 7c, also indicates that the built-in electric field within the ferroelectric In 2 Se 3 layer still exists in the presence of the graphite electrodes. All these results predict a realistic material system that is energetically stable and possesses 2D ferroelectricity with out-of-plane electric polarization at the single-layer limit.
For multilayer In 2 Se 3 films, their net polarization increases as a function of the film thickness and saturates as the thickness is above two QLs, in the cases that the polarization of all the ferroelectric QLs is aligned along the same orientation at each thickness. Detailed calculation results and discussions are provided in Supplementary Fig. 8 and Supplementary Note 2.
It is known that all ferroelectric materials are also piezoelectric and pyroelectric. We have investigated the piezoelectric property of a single QL of ferroelectric In 2 Se 3 by calculating the variation of the electric dipole as a function of the in-plane lattice deformation (Supplementary Fig. 9a) and the height variation as a function of an external electric displacement field applied in the vertical direction (Supplementary Fig. 9b). The results indicate that the compressive strain has a more significant effect on the electric polarization than the tensile strain. In particular, a tensile strain can slightly enhance the dipole moment of the ferroelectric layer. Furthermore, the electric-field-induced height variation in one QL ferroelectric In 2 Se 3 is very small, on the order of 10 À 3 Å at the electric displacement field as large as 0.2 V Å À 1 , which would be challenging to be resolved by piezoresponse force microscopy.
Kinetics of polarization reversal. To further demonstrate the ferroelectric nature of the systems, we must show that the direction of the electric polarization can be readily reversed by the application of a physically realistic electric field. We address this issue in two stages: first, we identify the most effective kinetic pathway connecting the two degenerates states with different polarities in the absence of an electric field; secondly, we investigate the effect of an external electric field on further reduction of the activation barrier along the pathway.
As a brute force check in the first stage, we find that the activation barrier against direct shifting of the central Se layer is 0.85 eV per unit cell, as shown in Fig. 2a. More importantly, an alternative process with a significantly lower activation barrier to reverse the electric polarization is revealed via a three-step concerted motion of the upper three Se-In-Se layers, as illustrated in Fig. 2b. In the first step, the a-phase FE-ZB 0 structure transforms into the metastable b-phase fcc 0 structure by laterally shifting the top three atomic layers together along the same direction to neighbouring sites. In the second step, the central Se atoms rotates around the C sites by 60°to a degenerate fcc 0 structure. In the third step, only the top two layers laterally shift along a direction that is rotated away from the original shifting direction by 60°, finally reaching an equivalent FE-ZB 0 structure, but now with the electric polarization reversed. The overall activation barrier of this concerted process is much lower, with the highest barrier to be only 0.066 eV per unit cell along the first step, comparable to that of the popular ferroelectric material PbTiO 3 (ref. 27).
In the second stage, we examine how the application of a perpendicular electric field reduces the kinetic barrier by lifting the degeneracy of the two polarized states. Our detailed calculations show that the activation barrier associated with the three-step concerted mechanism decreases linearly with the electric displacement field in the range of the field strength less than 0.3 V Å À 1 (Fig. 3a). It is worthwhile to note that the electric a b  field induces much more dramatic variation in the energy difference between the two oppositely polarized states than in the activation barrier. As shown in Fig. 3a, an electric displacement field of 0.3 V Å À 1 gives rise to an energy difference as large as 0.056 eV per unit cell, which is expected to result in a nearly ten times population difference of the two oppositely polarized states at room temperature. Moreover, for a ferroelectric domain of a practical device, the energy difference between the two polarized states is proportional to the domain size. Therefore, the population difference increases exponentially with the domain size. Although the activation barrier is also proportional to the domain size, given the relatively small activation barrier of 0.066 eV per unit cell, it is still possible to have an optimal domain size that gives rise to not only a sufficiently large energy difference to drive the reversal of the electric polarization but also a moderate activation barrier to make the kinetic process accessible at room temperature. Experimentally, an electric displacement field as large as 0.3 V Å À 1 had been demonstrated previously, for example, to open a sizable band gap in bilayer graphene 28 , which is expected to have a lower electric breakdown voltage than the present systems. In addition, it is important to point out that the reversal of the out-of-plane electric polarization accompanies with the reversal of the in-plane electric polarization for the FE-ZB 0 phase, as illustrated in Fig. 3b. We also examine the effects of the application of an in-plane electric field on the kinetics of the electric polarization reversal process. The calculated results, as summarized in Fig. 3b, indicate that an in-plane electric field of 0.03 V Å À 1 applied in the [110] direction gives rise to an energy difference as large as 0.142 eV per unit cell between the two oppositely polarized states, which is expected to result in a more than 200 times population difference of the two oppositely polarized states at room temperature. These features may offer an alternative approach to switch the orientation of the out-of-plane polarization by the application of an in-plane electric field. So far we have limited ourselves to ideal 2D systems of infinite size. In physically realistic growth conditions, the systems are more likely to contain different types of defects, especially domain walls 3,5 . Here we show that the electric polarization reversal process can be further facilitated by the presence of a domain wall between two oppositely polarized domains. In doing so, we construct four possible domain wall structures, as shown in the initial states of Fig. 4a,b, by moving half of the Se atoms in the central layer initially aligned to one In layer to the neighbouring sites aligned to the other In layer. The total formation energy of the two domain walls in the configuration as shown in Fig. 4a (0.22 eV per unit cell) is much lower than that in Fig. 4b  cell against the direct shifting mechanism of the central Se layer discussed earlier.
In 2 Se 3 -based van der Waals heterostructures. Next, we demonstrate the device potential of the discovered 2D ferroelectric materials in van der Waals heterostructures, focusing on the electrical transport properties. As a reference, a single ferroelectric In 2 Se 3 QL is a semiconductor with an indirect band gap of 1.46 eV (calculated by HSE06, 0.78 eV by GGA-PBE) (Fig. 5a). Owing to the presence of the out-of-plane electric polarization of the ferroelectric layer, there is a built-in electric field within the material, leading to different alignments of the energy bands with respect to the vacuum level on different sides of a given ferroelectric QL. For a ferroelectric In 2 Se 3 QL, such a difference is as large as 1.37 eV (calculated by HSE06). As a van der Waals 2D material is stacked with a ferroelectric In 2 Se 3 layer, the energy bands of the two components are approximately aligned with respect to the vacuum level, due to their weak van der Waals interaction. Therefore, as different sides of the ferroelectric layer are in contact with the other 2D material, different band alignments result in different global electronic structures. As the first specific system, we consider a bilayer heterostructure by stacking a QL of ferroelectric In 2 Se 3 onto a single-layer graphene, which is a non-ferroelectric semimetal. As shown in Fig. 5b-e, the Schottky barrier across the interface can be altered by switching the electric dipole orientation of the In 2 Se 3 layer. The magnitude of the electric dipoles of the system is 0.11 and 0.03 eÅ per In 2 Se 3 unit cell for the two oppositely polarized configurations as shown in Fig. 5b,d, respectively. The next bilayer heterostructure system considered is formed by stacking a QL of ferroelectric In 2 Se 3 on a monolayer of WSe 2 , which is a non-ferroelectric semiconductor. As shown in Fig. 5f-i, the band shift leads to a significant band gap reduction when switching the electric dipole orientation of the In 2 Se 3 layer. The magnitude of the electric dipoles of the system is 0.10 and 0.06 eÅ per In 2 Se 3 unit cell for the two oppositely polarized configurations as shown in Fig. 5f,h, respectively. For both heterostructures, the reduction of the electric dipoles in one of the polarized configurations can be attributed to the screening effects due to the charge transfer between the two layers as indicated in Fig. 5e,i. The influence of the graphene and WSe 2 layers on the energetics and kinetics of the polarization reversal processes of the ferroelectric In 2 Se 3 layer is discussed in Supplementary Note 3. The observed tunable band alignments with the ferroelectric layer can be exploited for different technological applications, such as for non-volatile memory devices or in graphene-based electronics. It is particularly worthwhile to note that the tunability in the properties can be achieved by the application of an external field, but the desired functionalities can be preserved even after the external field is removed. To provide a generic guideline for the design of desirable heterostructures, a schematic diagram of the band alignments of a single ferroelectric In 2 Se 3 QL is provided in Supplementary Fig. 13.
Family of 2D ferroelectric III 2 -VI 3 compounds. So far we have limited our discussions on the intrinsic ferroelectric properties of In 2 Se 3 . Next, we show that ferroelectric 2D van der Waals materials can be harboured in a wider family of the III 2 -VI 3 materials. Our DFT calculations suggest that the ferroelectric phases, FE-ZB 0 and FE-WZ 0 , are also the ground states of Al 2 S 3 , Al 2 Se 3 , Al 2 Te 3 , Ga 2 S 3 , Ga 2 Se 3 , Ga 2 Te 3 , In 2 S 3 and In 2 Te 3 when such materials are prepared in the QL form. Their semiconducting electronic band structures and optimal lattice constants are shown in Supplementary Fig. 14, and their dynamic stability is confirmed by the lack of imaginary phonon modes in the calculated phonon band structures (Supplementary Fig. 15). We further note that all the In-containing compounds have both the stable ferroelectric (FE-ZB 0 and FE-WZ 0 ) and metastable fcc 0 structures, while all the Ga-containing compounds only possess the stable ferroelectric structures, with the fcc-derived structure being unstable.

Conclusions
In this work, we have discovered a class of stable single-layer van der Waals 2D ferroelectric materials rooted in III 2 -VI 3 compounds that possess both intrinsic out-of-plane and in-plane electric polarization, which can be reversed through readily accessible kinetic pathways with the assistance of a modest out-of-plane or in-plane electric field. In a broader prospective, these discoveries add an important branch to the family tree of 2D materials. Proper integration of these materials with other classes of 2D systems is expected to substantially broaden the tunability and device potential of van der Waals heterostructures.

Methods
Computational methods. The first-principles DFT calculations were performed using the Vienna Ab Initio Simulation Package 29 . Valence electrons were described using the projector-augmented wave 30,31 method. The plane wave expansions were determined by the default energy cutoffs given by the Vienna Ab Initio Simulation Package projector-augmented wave potentials. The exchange and correlation functional was treated using the PBE 32 parametrization of GGA for structural relaxations and total energy calculations. For the band structure calculations of pristine III 2 -VI 3 compounds, we also used the hybrid functional of Hyed-Scuseria-Ernzerhof (HSE06) (ref. 33). To model the 2D films, the supercells contain a unit cell of single QL structures with a vacuum region of more than 15 Å. A saw-like self-consistent dipole layer was placed in the middle of the vacuum region to adjust the misalignment between the vacuum levels on the different sides of the film due to the intrinsic electric polarization. A G-centred 12 Â 12 Â 1 Monkhorst-Pack 34 k-mesh was used for k-point sampling. Optimized atomic structures were achieved when forces on all the atoms were o0.005 eV A À 1 . The in-plane electric polarization was evaluated by using the Berry phase method 35 . The climbing image nudged elastic band method 36 is used to determine the energy barriers of the various kinetic processes. In the heterostructure calculations, we included the van der Waals corrections as parameterized in the semiempirical DFT-D3 method 37 . More details are provided in Supplementary Note 1.
Data availability. All relevant data are available from the authors.