Niobium oxide dihalides NbOX 2: a new family of two-dimensional van der Waals layered materials with intrinsic ferroelectricity and antiferroelectricity

Two-dimensional (2D) ferroelectric (FE) materials displaying spontaneous polarizations are promising candidates for miniaturized electronic and memory devices. However, stable FE orderings are only found in a small number of 2D materials by experiment so far. In the current work, based on high-throughput screening of a 2D van der Waals layered materials database and first-principles calculations, we demonstrate niobium oxide dihalides NbOX 2 (X = Cl, Br and I), a group of experimentally synthesized yet under-explored van der Waals layered compounds, as a new family of 2D materials that simultaneously exhibit intrinsic in-plane ferroelectricity and antiferroelectricity. Similar to FE perovskite oxides, polar displacement of Nb cations relative to the center of the anion octahedral cage can lead to experimentally measurable FE polarizations up to 27 l C cm (cid:2) 2 in layered NbOX 2 . The presence of low-lying antiferroelectric (AFE) phases can eﬀectively reduce the energy barrier associated with polarization switching, suggesting switchable ferroelectricity is experimentally achievable. In addition, the mechanism driving FE phase transitions in NbOX 2 monolayers around Curie temperature T C is clearly revealed by our finite-temperature simulations. NbOCl 2 monolayer is predicted to be a stable ferroelectric with T C above room temperature. Moreover, application of NbOBr 2 and NbOI 2 monolayers as 2D dielectric capacitors is further developed, where electrostatic energy storage of nearly 100% eﬃciency can be achieved in the 2D single-layer regime. Two-dimensional (2D) ferroelectric (FE) materials with stable spontaneous polarizations above room temperature hold great promise in miniaturized electronic and memory devices. Here, we demonstrate a new strategy to uncover a wide spectrum of experimentally synthesizable yet unreported in PbTiO 3 . Based on our material screening scheme, we not only obtain the known 2D FE materials, such as group-IV monochalcogenides, but also identify a new class of layered FE materials – NbOX 2 .

Introduction ABO 3 perovskite oxides with broken spatial inversion symmetry are commonly known ferroelectric (FE) materials, which exhibit spontaneous polarizations below Curie temperature T C . The continuous demands for miniaturized FE devices have motivated extensive explorations of FE perovskite oxide thin films with reduced thickness. 1 In the past decade, the emergence of various two-dimensional (2D) materials with diverse functionalities has offered a new platform to investigate ferroelectricity and cooperative FE properties at the nanoscale. In particular, newly discovered 2D FE materials can have stable out-of-plane (vertical) FE orderings even in atomic-thin layers, 2-4 with interlayer spacing below the critical thickness where ferroelectricity in most FE perovskite oxides would disappear. 5,6 In addition, strong coupling of 2D ferroelectricity with ultra high elastic strain 7,8 or superior optoelectronics reponses 9 make 2D FE materials display rich coupled physical properties, enabling promising applications such as miniaturized electronic and memory devices. 10 Guided by first-principles calculations, a number of 2D FE materials have been successfully synthesized and demonstrated with stable FE orderings by experiment. Typically, in-plane ferroelectricity and FE phase transitions were predicted in group-IV monochalcogenide (GeS, SnS and SnSe, etc.) monolayers, [11][12][13][14] and then experimentally demonstrated in atomic-thick SnTe. 15,16 Computational discovery of 2D FE a-In 2 Se 3 layers 17 initiated the extensive experimental studies on intercorrelated in-plane and out-of-plane ferroelectricity within this system. [18][19][20][21][22] Spontaneous symmetry breaking in the 1T phase of MoS 2 monolayers can lead to robust in-plane ferroelectricity, 23 while its structural polarity was detected by second harmonic generation measurement very recently. 24 The material discoveries mentioned above largely relied on researchers' experience or chemical intuition, which was usually restricted to identifying 2D FE materials with structure or composition analogous to the known ones. For example, group-IV monochalcogenide monolayers were investigated as structural analogues to phosphorene. 11 Therefore, providing rational strategies for a wide spectrum of new 2D FE material discovery will largely accelerate research progress within this field, which is also beneficial for a complete and insightful understanding regarding the unique properties of 2D ferroelectricity at the nanoscale.
In this work, instead of using a traditional intuition guided material design approach, we demonstrate a new strategy for discovering experimentally synthesizable 2D layered FE materials. Recent data-mining studies have provided the complete database covering almost all experimentally reported 2D van der Waals layered materials (including 2D layered FE materials). 25,26 Further screening of these 2D materials databases for acentric layered semiconductors or insulators with switchable polarizations can lead to the discovery of new classes of 2D layered FE materials which are not simply composition or structural analogues to any experimentally reported ones. After searching the 2D weakly bonded layered materials database provided by Cheon et al. 25 using basic structural and energetic criteria ( Fig. 1), we identified niobium oxide dihalides NbOX 2 (X = Cl, Br and I) as long-sought-after 2D van der Waals layered FE materials.
Bulk NbOX 2 have been synthesized since the 1960s, 27,28 and they can be potentially exfoliated into 2D layered forms. Despite their structural polarities at room temperature, 28-31 NbOX 2 have not been considered or investigated as 2D layered ferroelectrics. Except for earlier tight-binding band structure calculations, 32 accurate predictions regarding electronic and FE properties for NbOX 2 bulk and 2D layers have never been reported. In our work, based on first-principles calculations, model Hamiltonian and molecular dynamics simulations, we identify the presence of stable FE and metastable antiferroelectric (AFE) phases in both 3D bulk and 2D layered forms of NbOX 2 . The coexistence of FE and AFE phases can lead to unique polarization switching and structural phase transition features in NbOX 2 , distinct from the currently known 2D FE materials. Our findings demonstrate NbOX 2 as a practical system to study the intrinsic ferroelectricity and antiferroelectricity down to the monolayer limit, and to explore the novel functionalities associated with 2D FE/AFE phase transitions.

Results and discussion
Ground state structures and ferroelectric properties Niobium oxide dihalides NbOCl 2 , NbOBr 2 and NbOI 2 belong to a group of transition metal oxide halide compounds with chemical formula MOX 2 (M = V, Nb, Mo and Ta, and X = Cl, Br and I). 27,28 Bulk NbOX 2 is a typical van der Waals layered material, 25 which is formed by stacking of NbOX 2 monolayers along the out-of-plane direction. Due to the weak van der Waals interlayer interactions (the interlayer binding energy E b = 12-13 meV Å À2 in Table 1, while the E b of graphite is predicted to be 20.3 meV Å À2 33 ), it should be eminently possible to obtain 2D NbOX 2 multi-layers or monolayers via mechanical exfoliation techniques. Similar to the perovskite crystal structure, the NbOX 2 monolayer is composed of NbO 2 X 4 octahedra with mixed edge-and corner-sharing connectivity. Typically, as shown in Fig. 2(a), NbO 2 X 4 octahedra build up a 2D structural network through extensive interconnecting X-X edges along one planar direction (crystallographic a axis) and cornered O atoms along the other (b axis).
As determined by experimental structural characterizations, 29-31 bulk NbOX 2 can crystalize in a ferroelectric C2 phase with monoclinic symmetry, where Nb cations exhibit non-zero off-center displacement along the Nb-O-Nb atomic chain direction (b axis). In our work, in order to explore all possible NbOX 2 polar configurations, we choose the Nb polar-distortion free, centrosymmetric C2/m phase (ground state structures for MoOCl 2 34 and TaOI 2 35 ) as the paraelectric (PE) reference for bulk NbOX 2 . Fig. 2(a) displays the crystal structure for the PE bulk phase (crystallographic parameters shown in Table S1 of the ESI †), in which there exits one-dimensional Peierls distortion 36 of Nb atoms in each NbO 2 X 4 octahedra, leading to an alternation of Fig. 1 High-throughput screening scheme for 2D FE materials. We start by looking for non-metallic 2D monolayers without inversion symmetry, which restrict the search to about 260 types of 2D insulating or semiconducting monolayers with C 1 , C 2 , C 2v , C 3v , C 4v and C 6v polar point groups. We then analyze polarization for each candidate, and identify two potential FE states with opposite polarizations (AEP states). The potential PE phase is obtained after ''averaging'' the atomic structures of the AEP states. After structural optimization, those 2D materials with unstable FE or PE structures are excluded. More importantly, the energy difference between the FE and PE phases (DE) should be small enough so that the polarization can be switchable. The energy criteria for D E is chosen to be 40 meV per atom, corresponding to the PE-FE energy difference in PbTiO 3 . Based on our material screening scheme, we not only obtain the known 2D FE materials, such as group-IV monochalcogenides, but also identify a new class of layered FE materials -NbOX 2 .
two unequal Nb-Nb distances along the Nb-X-Nb direction (a axis). As we will discuss later, Nb-Nb Peierls distortion is crucial for determining the electronic structures and band gaps of NbOX 2 . After examining the phonon spectra of PE phases, the soft phonon modes and corresponding polar structural instabilities for bulk NbOX 2 can be identified.
We use bulk NbOI 2 as an example and display the calculated phonon spectrum of its PE phase in Fig. 2(b). Three types of soft optical phonon modes with imaginary frequencies at the G point (o = i 83, 82 and 71 cm À1 ) are identified, corresponding to FE (all Nb atoms have exactly the same polar displacement), antipolar (two neighboring monolayers have anti-parallel Nb polar displacement) and AFE (Nb atoms from neighboring octahedra within each monolayer displace oppositely) Nb polar displacement patterns relative to the center of the NbO 2 I 4 octahedra. It is noted that Nb polar displacement in all three modes is restricted along the Nb-O-Nb direction. We next inject each soft mode into PE NbOI 2 , followed by structural optimization to obtain lower symmetry polar phases. FE phases with C2 symmetry (Table S2 in the ESI †) are obtained after relaxation of FE-mode related structure. Cooperative displacement of all Nb atoms (polar displacement amplitude d Nb E 0.14 Å) in the FE phase can break inversion symmetry and generate non-zero spontaneous polarizations along the b axis, which can be qualified using Berry phase calculations. The typical polarization-energy double well curves are obtained for bulk NbOX 2 FE phases ( Fig. S1 in the ESI †), after recording the variation of polarization and total energy of the system with respect to the amplitude of the FE mode connecting the FE C2 and PE C2/m phases. Based on our calculations, bulk NbOCl 2 exhibits the largest spontaneous polarization, highest FE potential depth and largest Nb polar displacement amplitude (Table 1), due to the significant electronegativity of the Cl anion. Besides the FE phase, antipolar phases with P2/c symmetry and AFE P% 1 phases (crystallographic parameters shown in Tables S3 and S4 of the ESI †) with anti-parallel interlayer or intralayer Nb polar displacement are also identified. Each Nb cation in the antipolar or AFE phases has non-zero polar displacement, but their contributions to the overall polarization of the system cancel out.
After investigation of bulk NbOX 2 , we now turn to 2D NbOX 2 multi-layers and monolayers, where we will demonstrate these NbOX 2 layers as intrinsic 2D FE (AFE) materials. The soft-mode adopted structural optimization scheme is also used to obtain the polar ground state structures for 2D NbOX 2 layers. FE, antipolar and AFE soft modes are found in the PE phase of NbOX 2 multi-layers, while the monolayers only have FE and AFE modes ( Fig. S2 and Table S5 in the ESI †). Without a vertical stacking sequence and out-of-plane periodicity, all 2D NbOX 2 layers are optimized into orthorhombic symmetry. Fig. 3(a) summarizes our calculated energetic results for all polar phases in NbOI 2 (results for NbOCl 2 and NbOBr 2 are provided in Fig. S3 of the ESI †). The FE phase is the most stable structure for all NbOX 2 layered systems, independent of the layer numbers. Except for the monolayer, the antipolar phase is the next most stable phase, and is very close to the FE phase in energy. Table 1 Calculated interlayer binding energy (E b ), Nb polar displacement (d Nb ) relative to the center of the NbO 2 X 4 octahedra, FE potential depth D E FE , spontaneous polarization P and nominal energy band gaps E g for bulk and monolayer FE NbOX 2 . We assume the thickness of the NbOX 2 monolayer is same as the interlayer spacing of bulk NbOX 2 , when simulating its P. The optical transition allowed E g are given in parentheses Even the less stable AFE phase has a small energy difference (o2 meV f.u. À1 ) relative to the FE phase. FE phases of NbOX 2 also show weak layer-dependent in-plane ferroelectricity due to the van der Waals interlayer interactions. 37 Shown in Fig. 3(b) are the calculated double-well potential curves for FE NbOI 2 , where NbOI 2 monolayer, bilayer and bulk exhibit almost identical Nb polar displacement amplitude and FE potential depth in their FE phases (similar results are also found in NbOCl 2 and NbOBr 2 ). Based on layerindependent energetic stability and ferroelectricity, it is expected that in-plane FE polarization of almost the same magnitude (Table 1) can be obtained in all structural forms of NbOX 2 , including 3D bulk and 2D layers. This unique feature makes 2D NbOX 2 distinct from group-IV monochalcogenides (e.g. GeS, SnS and SnTe), a well studied 2D FE material class with in-plane ferroelectricity. As it is demonstrated by experiment, 2D orthorhombic SnTe layers (g phase) exhibit strong antipolar interlayer coupling, leading to a vanishing of ferroelectricity in the bulk phase and those multi-layers containing even numbers of layers. 16 NbOX 2 monolayer can be stabilized into two polar phases, either a stable FE or metastable AFE phase, which contain two NbO 2 X 4 octahedra in their unit cell. To evaluate the possibility of FE-to-AFE transition, we calculate 2D energy contour plot of NbOI 2 monolayer as a function of polar displacement of two Nb atoms within a single unit cell. As shown in Fig. 3(c), two degenerate FE (AFE) phases locate along the diagonal directions, separated by the PE phase. Transition between two FE phases with opposite polarizations does not need to go through the high-energy PE phase. Instead, it will cross a low-lying AFE phase, where the polar displacement of one Nb atom is reversed, while keeping the other almost unchanged.

Electronic structure and optical absorption properties
Layered NbOX 2 , including 3D bulk and 2D layers, have FE ground states with spontaneous polarizations. The coupling between FE ordering and the electronic and optical properties of 2D materials can lead to interesting phenomena, such as the bulk photovoltaic 9,38 and ferrovalley effects. 12, 39 We will investigate the electronic structures and optical absorption properties for the FE phases in bulk and monolayer NbOX 2 . Fig. 4 displays our calculated energy band structures and projected density of states (PDOS) for FE NbOX 2 , obtained from HSE hybrid density functional calculations. All NbOX 2 systems, including 3D bulk and 2D monolayer, are predicted to be indirect semiconductors with nominal E g (the energy difference between the filled and empty band edges, shown in Table 1)   40. For NbOX 2 monolayer, the contribution of Nb-d z 2 orbital to energy bands are highlighted by green circles. The variation of the energy separation between the occupied and empty Nb-d z 2 orbitals at G point (DE Nb-d ) as a function of Nb-Nb paring distance (d Nb-Nb ) in FE NbOX 2 monolayers is also presented. NbOCl 2 monolayer with the shortest d Nb-Nb has the strongest Peierls distortion intensity and therefore the largest DE Nb-d . around 1.8 eV. Due to the hybridization of Nb-4d with anion p orbitals, largely dispersed energy bands around the conduction band minimum (CBM) in NbOX 2 are mainly contributed to by empty Nb-4d states. Below the Fermi level, the highest occupied valence band is almost dispersionless through the entire Brillouin zones of the monoclinic bulk and orthorhombic monolayer. For NbOCl 2 and NbOBr 2 , such a flat band is isolated from other valence bands of hybridized anion p orbitals. Further orbital component analysis indicates that the dispersionless valence bands (marked as green circles in the energy bands of the monolayers) originate from localized Nb-d z 2 states. Meanwhile, similar weakly dispersed energy bands also appear in conduction states (above CBM), corresponding to empty Nb-d z 2 orbitals.
Formation of the energy gap between the filled and empty Nb-d z 2 states (DE Nb-d ) is a direct consequence of 1D Peierls distortion. 36 Each Nb 4+ cation in NbOX 2 has one unpaired 4d electron. Paring of two 4d electrons simultaneously occurs in the NbOX 2 single unit cell, after periodic dimerization of neighboring Nb cations created by the 1D Peierls distortion. The energy gap is then formed as two paired electrons are favored to fully occupy the low-energy Nb-d z 2 orbital, leaving the high-energy d z 2 orbital empty. We double checked the electronic structures for NbOX 2 using spin polarized calculations, which also predict nonmagnetic ground states with paired electrons on the Nb-d z 2 orbitals. In fact, the intensity of the 1D Peierls distortion in NbOX 2 can be well qualified by energy separation between the occupied and empty Nb-d z 2 orbitals (DE Nb-d ). As shown in the last plot of Fig. 4, the strongest Peierls distortion with the shortest Nb-Nb distance in NbOCl 2 monolayer leads to the largest DE Nb-d . The overall extraordinarily large DE Nb-d in all three systems indicates that 1D Peierls distortion in NbOX 2 is ultra-stable, able to survive over thermal fluctuations at high temperature.
In layered NbOX 2 , Nb-d z 2 orbitals can hardly hybridize with any planar anion p orbitals. As a result, Nb-d z 2 states correspond to the dispersionless energy bands that are quite localized in energy. More importantly, optical transitions from the localized Nb-d z 2 valence band to other empty Nb-d orbitals from conduction bands are completely forbidden. The effective optical absorption in NbOX 2 comes from the transition between hybridized anion-p orbitals of the valence band and hybridized Nb-4d orbitals (except d z 2) of the conduction bands. Therefore, as far as optical absorption properties are considered, NbOCl 2 and NbOBr 2 are insulators with optical transition allowed E g 4 3.0 eV (Table 1). Due to the weak electronegativity of the I anion, bulk and monolayer NbOI 2 can exhibit semiconducting electronic and optical properties. Typically, as a ferroelectric semiconductor with an indirect E g , bulk NbOI 2 can effectively absorb visible light with a photon energy above 2.0 eV (Fig. S4 in the ESI †), making it suitable for FE photovoltaic applications. 41,42 2D NbOX 2 layers have orthorhombic symmetry with anisotropic planar crystal structures, where Nb atoms bond with X-X edges along one planar direction and O atoms along the other. Similarly to many other 2D FE materials, 12,43 strong anisotropy in crystal structures of 2D NbOX 2 will be reflected in their energy bands and optical absorption properties as well. Using semiconducting monolayer NbOI 2 as an example, we will illustrate its highly anisotropic and strongly coupled electronic and optical absorption properties. Fig. 4(a) shows the band structure for FE NbOI 2 monolayer, where all the energy bands around the Fermi energy level are marked by colored circles whose radii are proportional to the contribution of the corresponding atomic orbitals. I anion has a weaker electronegativity than O, therefore the valence bands around the Fermi level mainly come from I-p orbitals. Except for the localized d z 2 bands, most valence and conduction bands show anisotropic dispersion along two planar directions: valence I-p y bands and low-lying conduction Nb-d yz bands are highly dispersed along G-Y (reciprocal lattice equivalent of the crystallographic b axis), but weakly dispersed along the G-X direction.
Besides electronic structures, we have also investigated the anisotropic optical properties of NbOI 2 monolayer by simulating its optical absorption coefficient a, excited by incident light polarized along the x and y directions (corresponding to the a and b axes, respectively). As shown in Fig. 5(b), both a x and a y have optical absorption edges around 2.0 eV. Above the absorption edge, a strong optical absorption peak appears at 3.15 eV, when incident light is polarized along the x direction. The strength of this absorption peak in a x is larger than that of a y by almost 10 times, indicating a strong direction-dependent photo-absorption in NbOI 2 monolayer. As energy bands around the Fermi level are mainly contributed to by hybridized I-p and Nb-d orbitals, the low-energy optical absorption in NbOI 2 can be assigned to optical transitions from I-p to Nb-d orbitals. When NbOI 2 monolayer is excited by incident light polarized along the x direction, optical transitions between all weakly dispersed Nb-d yz and I-p y energy bands along the whole G-X direction (indicated by dashed red arrows in Fig. 5(a)) can contribute to optical absorption at 3.15 eV in a x . While along the y direction, the weaker absorption at 3.15 eV comes from the optical transition (blue dashed arrow) of Nb-d yz and I-p y bands around the G point only. As a result, when monochromatic light with photon energy of 3.15 eV is incident on FE NbOI 2 monolayer, optical absorption along the nonpolar a axis is significantly larger than that along the polar axis, leading to a nearly linearly polarized optical absorption selectivity in NbOI 2 monolayer. 39

Ferroelectric polarization switching and domain-wall properties
It is clearly shown that layered NbOX 2 systems have in-plane ferroelectricity at the monolayer limit. As a fundamental nature of ferroelectricity, FE polarizations in 2D NbOX 2 should be switchable under an applied electric field. Unlike other 2D FE materials with planar ferroelectricity (e.g. SnS and SnSe), FE ordering in 2D NbOX 2 monolayers indeed has a ''1D'' nature, as the ferroelectricity or antiferroelectricity is restricted along the polar b axis, but forbidden along the other nonpolar axis. Based on its 1D collinear nature, FE polarization switching in 2D NbOX 2 monolayers can be achieved through the direct reversal of the polarization direction by 1801, rather than rotation of polarization across the nonpolar axis.
During the polarization reversal (switching) process, FE domain-walls between oppositely orientated FE mono-domains will be formed in NbOX 2 monolayers. Meanwhile, the polarization reversal process can be carried out through the motion of FE domain-walls. 44,45 Fig. 6(a) displays the supercell configuration containing two FE mono-domains with opposite polarization directions, separated by a 1801 FE domain-wall located along the I-I lattice plane in NbOI 2 monolayer. Each mono-domain is composed of four NbOI 2 unit cells, stacking along the a axis. To preserve the periodic boundary condition required by first-principles plane wave calculations, we must include two coherent FE domain-walls (twin walls) in one supercell, so that the whole supercell can be repeated periodically along both axes.
In order to evaluate how polarization is reversed across the FE domain-wall, we analyze the evolution of Nb polar displacement amplitude (d Nb ) as a function of Nb distance (r) away from the domain-wall. As shown in Fig. 6(a), polarization is reversed as d Nb changes its sign abruptly across the domain-wall. d Nb of each Nb cation almost fully recovers to its mono-domain value (indicated by dashed lines). As a result, the 1801 FE domainwall supercell configuration is geometrically similar to a NbOI 2 superlattice, where AFE phases are sandwiched between two FE phases with opposite polarizations. As the energy difference between AFE and FE phases of NbOI 2 monolayer is as low as 1.24 meV f.u. À1 , the 1801 FE domain-wall has a very small domain wall energy. The d Nb -r curve shown in Fig. 6(a) can be quantitatively described using a function corresponding to the soliton solution of one-dimensional fourth-order Landau-Ginzburg domain wall theory 44 as:  44,45 the absence of PE phase in the NbOI 2 domain-wall configuration is responsible for its ultra narrow domain-wall width.
We then explore the kinetic process for domain-wall motion by simulating the corresponding minimum energy pathway (MEP) trajectory and the associated energy barrier height, using a generalized solid-state NEB method. 46 Such a method allows both atomic positions and lattice parameters to relax along the pathway. As shown in Fig. 6(b), transition from the initial to the final domain-wall configurations resembles the motion of the domain wall from one I-I lattice plane to the next, accompanied by the reversal of polarization direction of one group of Nb cations by 1801. At the saddle point configuration along the MEP, the domain-wall moves to the lattice plane where Nb cation is located. This very Nb cation is then constrained to have zero d Nb , while d Nb of other Nb cations remain unchanged. As a result, the MEP associated with domain-wall motion is similar to the energy pathway for transition from single FE to AFE phase, as we obtained earlier from the energy contour plot in Fig. 3(c). The overall barrier height for polarization reversal of the Nb cation in the NbOI 2 monolayer is predicted to be 14.9 meV, which is much smaller than the polarization rotation barrier (B40 meV) in FE perovskite PbTiO 3 . 47 Therefore, the presence of the low-lying AFE phase can lead to a small energy barrier, which is beneficial for easy domain-wall motion and polarization switching in NbOI 2 monolayer.
Using the simulated energy barrier, we can estimate how large the electric field needs to reverse FE polarization in NbOI 2 monolayer. The critical electric field E C is given by: E C C E barrier /PÁV, where P is the spontaneous polarization and V is the normalized cell volume for NbOI 2 monolayer. The estimated critical electric field E C C 0.63 MV cm À1 , which is readily accessible by experiment. 48 Therefore, electric field induced switching of FE polarization can be achievable under laboratory conditions.

Structural phase transition at finite temperature
Up to now, we have simulated the structural, electronic and polarization switching properties for FE NbOX 2 layers using zero-temperature DFT calculations. For practical applications of ferroelectrics, 2D NbOX 2 should have relatively high Curie temperature T C , so that their polarizations can still persist above room temperature. In the following, we will investigate the ferroelectricity and structural phase-transition of NbOX 2 monolayers at finite temperature using both Monte Carlo (MC) and ab initio molecular dynamics (MD) simulations.
The geometry of 2D NbOX 2 monolayers can be represented by a 2D lattice grid containing a number of Nb cations, where each Nb cation at 2D grid point i [i is a collapsed index of its 2D position (m, n)] has a unique polar displacement d i . Due to the 1D nature of polarization, d i of Nb cations are collinearly arranged into columns along the polar axis. Using polar displacement d i as the order parameter, the configuration of NbOX 2 monolayer at any condition can be unambiguously specified. Therefore, we can express the free energy of NbOX 2 monolayer as Landau-Ginzburg expansion of order parameter d i as: 13 GðdÞ where the first two terms describe the FE energy for each Nb cation and parameters A and B can be obtained after fitting a double-well potential for FE NbOX 2 monolayers. The last two terms correspond to the interactions of d i and d j between two nearest neighboring Nb cations. Here, we include the parameters C x and C y , which measure the interactions along the non-polar and polar axes respectively (Fig. S5 in the ESI †). Due to the high energy cost to form a head-to-head/tail-to-tail polar configuration (Nb cations within the same column have antiparallel d i ) in NbOX 2 , parameter C y c C x . As shown in Table 2, the absolute values for most of the parameters A-C follow the trend: NbOCl 2 4 NbOBr 2 4 NbOI 2 , consistent with the fact that FE NbOCl 2 has the strongest ferroelectricity. Based on the effective Hamiltonian we developed for NbOX 2 monolayers, we can investigate temperature induced structural transition using MC simulations (computational details in the ESI †). For comparison, we also perform parameter-free ab initio MD simulations, where the time evolution of the instantaneous temperature, total energies and cation polar displacement of NbOX 2 can be obtained (Fig. S6 in the ESI †). Due to the large computational cost, we restrict our MD simulations of NbOX 2 monolayers to selected target temperature. Fig. 7 shows our simulated temperature dependent macroscopic average polar displacement hd(T)i in NbOX 2 monolayers. The structural transition from FE to PE phase is identified, as the simulated hd(T)i drops abruptly to zero around T C . Moreover, at the selected temperature, ab initio MD predicted ensemble average hd(T)i are in overall good agreement with MC results (the largest derivation occurs in NbOI 2 , with DT B 60 K), which validates the effective Hamiltonian and parameters we used for MC simulations. Curie temperature T C can be quantitatively determined after fitting MC simulated hd(T)i as follows: 13 where d is critical exponent for the FE-PE phase transition. As shown in Table 2, the tendency of T C for FE NbOX 2 monolayers also follows: NbOCl 2 4 NbOBr 2 4 NbOI 2 . NbOCl 2 monolayer is predicted to have the highest T C of 396 K, and can exhibit considerable polarization at room temperature. In order to understand the microscopic mechanism governing FE-PE phase transition, we further examine temperature dependent cation displacement evolution. The lower panel of Fig. 7 plots the real space mapping of d i from the whole NbOI 2 lattice grid during MD simulations after the system reaches thermal equilibrium at target temperature below, approaching and above T C , respectively. Below T C , the magnitude of d i starts to decrease from its zero temperature value as temperature increases. The polarization reversal event (changing sign of d i ) Table 2 Parameters used for MC simulation in eqn (2). Except the fourthordered parameter B (in meV Å À4 per f.u.), parameters A, C x and C y are all in units of meV Å À2 per f.u. Curie temperatures T C are obtained via fitting to MC data based on eqn (3)  does not occur until T E T C , where |d i | of some Nb cations become small enough so that thermal energy is comparable with the energy cost to form a head-to-head/tail-to-tail polar configuration. The polarization reversal process around T C is initiated by switching d i of a single Nb cation from one cation column. Then the whole column simultaneously switches its polarization to the opposite direction, to avoid the unfavorable head-to-head/tail-to-tail configuration. The FE-PE phase transition is achieved through dynamic reversal of column polarization as mentioned above. As AFE phase in NbOX 2 is more stable than the distortion-free PE phase. Even when T 4 T C , the system stays in a thermal-equilibrium PE state, where the overall polarization is almost zero, but individual columns can still have non-zero and randomly oriented polarization along the polar axis. As a result, PE state above T C is more like a disordered AFE phase. Moreover, the energy cost associated with ''collective'' reversal of whole column polarization is negligible compared to the ''isolated'' cation dipole switching energy (C y c C x , Table 2), thus T C of NbOX 2 monolayer is mainly determined by ''isolated'' rather than ''collective'' cation dipole switching events. 49 Therefore, T C is not simply proportional to the energy difference associated with ''collective'' polarization reversal (FE-to-AFE phase transition). Even though AFE phase is close to FE phase in energy, NbOX 2 monolayers can still have T C around or even above room temperature.

Experimental outlook and FE device by design
We have provided strong evidence that 2D NbOX 2 possess intrinsic in-plane ferroelectricity and antiferroelectricity. Experimentally, bulk NbOX 2 have already been prepared using solid-state sintering methods, while their structure polarity at room temperature was verified by XRD measurement. [28][29][30][31] Direct identification of ferroelectricity or antiferroelectricity in bulk NbOX 2 requires further measurement of their polarization-electric field (P-E) hysteresis loops. 2D NbOX 2 layers, including multilayers and monolayers can be obtained by direct exfoliation of bulk samples or chemical vapor deposition. As long as the temperature is below T C , 2D NbOX 2 layers can crystalize in a ground state FE phase with polarization comparable with that of bulk NbOX 2 , after poling by electric field (above critical field E C ) applied along the polar axis. Based on our simulation, 2D NbOX 2 monolayers exhibit variable FE polarizations and T C . Chemical alloys made of individual NbOX 2 end members can offer a wide range of material choice to achieve optimal T C and FE properties. For example, a 2D NbOCl x I 2Àx monolayer is expected to possess T C above room temperature and semiconducting properties comparable with NbOI 2 . More specifically, increasing concentration x can lead to a higher T C and larger optical transition allowed E g in the NbOCl x I 2Àx monolayer. Alternatively, one can also use experimentally accessible tensile strain to engineer the FE properties of NbOX 2 monolayers. As shown in Fig. S7 of the ESI, † more than 6 times enhancement of the FE potential depth and twice of the Nb polar displacement can be achieved with 3% tensile strain applied along the polar axis of NbOX 2 monolayers.
As a result, tensile strained NbOBr 2 and NbOI 2 monolayers are expected to have elevated T C above room temperature. Due to linearly polarized optical selectivity of the semiconducting NbOI 2 monolayer, its in-plane polar and non-polar axes can be easily distinguished using incident polarized light. As NbOX 2 have the interconnected octahedral structural framework similar to perovskite oxides, after imaging of Nb off-center polar displacement (0.13-0.14 Å) using annular bright-field scanning transmission electron microscopy (ABF-STEM), 50,51 one can directly observe and quantify FE polarization appearing in 2D NbOX 2 layers. Moreover, measuring the shear piezoelectric effect arising from in-plane polarization by lateral piezoelectric force microscopy (PFM) 37 can also provide a direct experimental proof of FE polarization in 2D NbOX 2 . The 1D polarization nature and coexistence of FE/AFE phases make 2D NbOX 2 exhibit the unique polarization reversal and structural phase-transition features. In experiment, AFE perovskite oxides have been widely used as dielectric capacitors for electrostatic energy storage through electric charging and discharging processes. 52 In particular, high energy densities are achieved in perovskite solid solutions near the phase boundaries where FE and AFE phases coexist. 53 We therefore propose the application of NbOX 2 monolayers as 2D dielectric capacitors. Fig. 8(a) displays a schematic diagram for a 2D single-layered capacitor obtained by lateral growth of a heterostructure between a NbOX 2 monolayer and another 2D metallic material, where an electric field can be applied along the polar axis through a metal electrode. In order to evaluate the performance of a 2D NbOX 2 based capacitor, we choose NbOI 2 monolayer as an example, and investigate the electric-field induced phase transition by simulating its P-E loop around room temperature using MC simulations. When electric field E is applied along the polar axis, the additional energy term À(EÁP)V will be incorporated into the effective Hamiltonian in eqn (2). The effective polarization P can be computed from Nb polar displacement d i and Born effective charge Z i * by: Therefore, the electric field related energy term can be simplified as: ÀE(S i d i ÁZ i *). Based on our predictions, NbOI 2 monolayer adopts a disordered AFE (PE) phase at room temperature. The simulated P-E loop in Fig. 8(b) shows that under the charging process (increasing E) P increases abruptly around 0.15 MV cm À1 , indicating that a small critical field E C can trigger the AFE-to-FE transition in NbOI 2 monolayer at room temperature. Above E C , P increases linearly within FE phase upon further charging. Remarkably, an almost identical P-E curve associated with the FE-to-AFE transition is obtained under the discharging process (decreasing E). As a result, under the applied electric field, NbOI 2 monolayer possesses P-E loops with almost no hysteresis, which is consistent with its small FE/AFE energy difference and low polarization switching energy barrier at room temperature. 54 Similarly to NbOI 2 , NbOBr 2 monolayer also crystalizes in a disordered AFE phase at room temperature and exhibits nearly hysteresis-free P-E loops. While NbOCl 2 monolayer has a non-zero FE polarization at room temperature, it features P-E double hysteresis loops of typical ferroelectrics. Due to the energy loss associated with the hysteresis, NbOCl 2 monolayer is less suitable for electrostatic energy storage than NbOBr 2 or NbOI 2 . A 2D capacitor based on NbOBr 2 or NbOI 2 monolayers can exhibit nearly zero energy loss during the charging and discharging processes, yet has the advantage of operating at a low electric field, which can potentially enable nearly 100% efficiency for electrostatic energy storage in the 2D monolayer.

Conclusions
In summary, we have identified niobium oxide dihalides NbOX 2 as a new family of 2D van der Waals layered materials, exhibiting intrinsic in-plane ferroelectricity and antiferroelectricity. Comprehensive investigations regarding structural, ferroelectric, electronic, optical absorption and domain-wall properties of 2D NbOX 2 have been carried out using firstprinciples calculations. NbOX 2 adopts a 2D structural framework composed of interconnected NbO 2 X 4 octahedra, where considerable cation off-center polar displacement along one planar axis can lead to 1D ferroelectricity or antiferroelectricity. Along the other planar axis, the ultra-stable Peierls distortion between Nb 4+ cations gives rise to the semiconducting or insulating electronic properties, which are preferable for stable polarizations. In-plane ferroelectricity or antiferroelectricity in 2D NbOX 2 can be sustained down to the monolayer limit, and the energetic stabilities and polar magnitudes are almost independent of the number of layers. Coexistence of FE and AFE phases makes 2D NbOX 2 exhibit easily switchable polarizations. Moreover, the strong coupling between FE order and optical absorption leads to a nearly linearly polarized optical selectivity in the semiconducting NbOI 2 monolayer. Besides zero temperature DFT calculations, we also develop an effective Hamiltonian to simulate the structural phase transitions of NbOX 2 monolayers under finite temperature and an applied electric field. The NbOCl 2 monolayer is predicted to be a robust FE material with room temperature stable ferroelectricity, while NbOBr 2 and NbOI 2 monolayers adopt a PE state composed of disordered AFE phases at room temperature, but will undergo a structural transition into a stable FE phase under a small applied electric field. Based on our simulations, P-E loops characterizing AFE-FE transitions have nearly zero hysteresis, therefore a capacitor composed of NbOBr 2 or NbOI 2 has a great advantage for highly efficient electrostatic energy storage. We hope our work will motivate future experimental exploration of 2D ferroelectricity and antiferroelectricity in layered NbOX 2 , and realization of superior energy storage performance that our simulations predicted.

Computational methods
Our first-principles calculations are performed based on density functional theory (DFT) as implemented in the Vienna ab initio Simulation Package (VASP-5.4.1). 55,56 A plane-wave basis set within the projector augmented-wave method 57 is employed, using a 600 eV plane-wave energy cutoff. We simulate NbOX 2 bulk and 2D layered structures using the strongly constrained and appropriately normed (SCAN) meta-GGA functional. 58 After inclusion of the revised Vydrov-van Voorhis nonlocal correlation (rVV10), 33,59 SCAN + rVV10 functional can accurately predict the interlayer spacings as well as intralayer lattice parameters for bulk NbOX 2 ( Table S2 in the ESI †). 2D NbOX 2 layers are modeled as slabs with a vacuum region of more than 20 Å along the out-of-plane direction. A G-centered Monkhorst-Pack k-point grid of about 40 k-points per Å À1 spacing is used for k-point sampling. NbOX 2 bulk and layers are fully optimized until the residual Hellmann-Feynman forces are smaller than 1 meV Å À1 and the stresses less than 0.1 kbar. Phonon frequencies and eigenvectors are calculated based on the finite difference method implemented in the Phonopy package. 60 The electronic structures and optical absorption spectra of NbOX 2 systems are further calculated using the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional. 61,62 As band gap changes due to spin-orbital coupling (SOC) effects are quite small (DE g o 0.02 eV for bulk NbOI 2 ), the SOC effect is excluded during HSE calculations. The electronic contribution to the polarization is calculated following the Berry phase formalism. 63 The generalized solid-state nudged elastic band (ss-NEB) method 46,64,65 is used to simulate the minimum energy path (MEP) and associated energy barriers for various kinetic processes.

Conflicts of interest
There are no conflicts to declare.