Purely one-dimensional ferroelectricity and antiferroelectricity from van der Waals niobium oxide trihalides

Intrinsic one-dimensional (1D) ferroelectric materials are rarely reported but are highly sought to break the size limit of nanostructured conventional ferroelectrics. Herein, we report a class of inborn 1D ferroelectric nanowires, namely 1D NbOX3 (X = Cl, Br, and I), that can be directly obtained from experimentally realized van der Waals crystals. In addition to the sizable spontaneous polarization, 1D NbOX3 exhibits low ferroelectric switching barriers, small coercive electric fields, and high critical temperature, governed by the hybridization of the Nb empty d orbitals and the O p orbitals (d0 rule). Moreover, the double-channel structure of 1D NbOX3 also enables the emergence of 1D antiferroelectric metastable states. Our findings not only propose a class of 1D ferroelectric materials toward the development of miniaturized and high-density electronic devices, but also pave an avenue of obtaining intrinsic 1D ferroelectrics from van der Waals crystals.


INTRODUCTION
Ferroelectric (FE) materials sustaining a spontaneous electrical polarization that can be switched by an external electric field, have broad technological applications, such as non-volatile randomaccess memories 1-3 , field-effect transistors 4,5 , sensors 6,7 , and photovoltaics 8,9 . Driven by the increasing demand for device miniaturization, significant progress has been made in two/onedimensional (2D/1D) ferroelectrics. The exploration of 2D FE materials started with the direct downscaling of conventional bulk FE materials, a strategy found numerous obstacles such as enhanced depolarizing electrostatic fields, severe surface reconstruction, dangling bonds at the surface, etc., thus imposing great limits on their use for practical applications [10][11][12][13] . The recent emergence of van der Waals (vdW) layered 2D FE materials 14,15 has stimulated extensive research, both experimentally and theoretically, by providing a possible solution to overcome these challenges thanks to their inborn atomic-thin structures [16][17][18][19][20][21][22][23][24][25][26][27] . Compared to 2D ferroelectrics, 1D FE materials are expected to be superior for building high-density FE devices. Several 1D FE nanostructures, including nanowires, nanotubes, nanoribbons, and belts, have been proposed. These are all based on manipulating conventional FE compounds [28][29][30][31][32][33] , meaning that such structures could suffer from the same drawbacks of their 2D counterparts. Moreover, the diameters of these 1D structures usually range from several tens to hundreds of nanometers, namely they are far from a nearly atom-scale size 28,34 . Ideally, 1D ferroelectrics should not only possess robust ferroelectricity but also simultaneously harbor small diameters, strong chemical bonding within the chain, and high stability. In order to overcome these challenges, it is highly important to explore new routes for the discovery of 1D ferroelectrics.
Experimentally a number of vdW crystals made of weakly bonded 1D building blocks (chains or molecular wires) have been successfully fabricated [35][36][37][38] . Inspired by the exfoliation of 2D functional materials from bulk layered structures, such type of vdW crystals could offer a great platform to explore the rich physics within the pure 1D limit. To date, research on 1D materials obtained from vdW bulk compounds is still limited, but some interesting properties have been already discussed, such as charge density waves in NbSe 3 39 , nontrivial band topology in TiCl 3 40 , and power-law dependent tunneling conductance in MoSe 41 . Importantly, the isolation of pure 1D structures from vdW crystals has been demonstrated possible in molybdenum based polyoxometalates 42 . If the 1D building block in a vdW crystal possesses broken centrosymmetry along the chain direction, then intrinsic 1D ferroelectrics at ultrasmall dimensions could be obtained via exfoliation from the corresponding vdW crystals.
In this work we report a family of intrinsic 1D ferroelectrics, namely 1D NbOX 3 (X = Cl, Br, and I), which are analyzed by mean of density functional theory (DFT) and model Hamiltonian Monte Carlo (MC) simulations. 1D NbOX 3 could be easily exfoliated from the experimentally synthesized vdW crystals due to the small binding energy. Moreover, 1D NbOX 3 exhibits sizable spontaneous polarization above room temperature, low FE transition barriers, and low coercive electric fields. The two-channel 1D geometry of NbOX 3 also enables the emergence of AFE metastable states, which has never been reported in the purely 1D regime. Our findings highlight an interesting avenue to realize intrinsic 1D ferroelectricity and antiferroelectricity from the vdW crystals.

RESULTS AND DISCUSSION Structural exploration and stability
The structures of the experimental bulk phases 43,44 of NbOX 3 (X = Cl, Br and I) are presented in Fig. 1a-c. Bulk NbOCl 3 and NbOBr 3 crystalize in a tetragonal lattice with the P42 1 m space group, while NbOI 3 displays a monoclinic structure with a lowered C2 symmetry. The PBE-calculated lattice parameters are summarized in Table 1 and they are in great agreement with the experimental values 43,44 . Interestingly, bulk NbOX 3 possesses a peculiar unit cell consisting of two parallel 1D NbOX 3 nanowires, where the Nb ions show an off-center displacement, d Nb , along the O-Nb-O atomic chain direction. In the 3D crystals, d Nb in neighboring 1D nanowires aligns antiparalleley in NbOCl 3 /NbOBr 3 , and parallely in NbOI 3 , leading to anti-polarized and polarized bulk structures for NbOCl 3 /NbOBr 3 and NbOI 3 , respectively.
Despite the different overall polar orders of bulk NbOX 3 , each 1D NbOX 3 nanowire within the bulk phase is perfectly polarized and weakly vdW-bound. These features provide the opportunity to exfoliate the long-sought 1D polarized structures with perfect radical boundaries from experimentally accessible bulk phases. The binding energy E b is then evaluated to determine the strength of the vdW interactions in the bulk. The relevant quantity is E b = E 1D −E bulk /2, where E 1D and E bulk are the energies of 1D and bulk NbOX 3 , respectively. The binding energies of NbOCl 3 , NbOBr 3 , and NbOI 3 are calculated to be 62.2, 78.0, and 97.4 meV per atom, respectively, and they are comparable to those of graphene (52 ± 5 meV per atom), MoS 2 (77 meV per atom) 45 , and phosphorene (81 ± 5 meV per atom) 46 . This suggests that the extraction of 1D NbOX 3 from bulk is highly possible. The cohesive energy E c is calculated by and E X represent the energy of 1D NbOX 3 , single Nb, O, and X atoms, respectively. The E c of 1D NbOCl 3 , NbOBr 3 , and NbOI 3 is calculated to be −4.82, −4.46, and −4.09 eV/atom, respectively, indicating their strong bonding characteristics.
1D NbOX 3 possesses two parallel NbOX 4 channels with the shared halogen atoms in the center as shown in Fig. 1d. In order to extensively explore the polar configurations possible in these double-channel 1D structures, the centrosymmetric phase (Pmmm symmetry) in which there is no Nb displacement along the b axis  Table 1. The lattice parameters (a, b, and c), the Nb polar displacement (d Nb ) relative to the center of the NbO 2 X 4 octahedra, and the spontaneous polarization (P) for bulk and 1D NbOX 3 . For bulk NbOX 3 , the differences between calculated lattice parameters and the experimental values are shown in parentheses (in percent). The P (in eÅ per unit cell) is also given in parentheses for clarity. Note that bulk NbOI 3 contains two 1D nanowires in its unit cell thus its P (in eÅ per unit cell) is around twice than that of 1D FE NbOI 3 .
is chosen as paraelectric (PE) configuration and the associated phonon spectra are examined first. Figure 1e presents the calculated phonon spectrum of 1D PE NbOBr 3 as an example. Clearly, two optical phonon modes show imaginary frequency (Γ 1 = −3.4 THz and Γ 2 = −3.1 THz) at the Γ point, indicating a dynamical instability of 1D PE NbOBr 3 . Further analysis of the phonon eigenvectors confirms that the Γ 1 and Γ 2 soft modes correspond to the off-center displacement, d Nb , with either identical or opposite sign in the two NbOBr 4 channels. Such soft modes spontaneously drive the system into either a FE (space group of Pmm2) or an AFE (space group of P2/m) state, whose dynamical stabilities are confirmed by the phonon spectra (see Fig. 1f, g) showing no negative frequencies. Similarly, we can identify the stable FE and AFE phases of 1D NbOCl 3 and NbOI 3 , as shown in Supplementary Fig. 1. We also consider head-to-head antipolar configurations as shown in Supplementary Fig. 2, which are energetically not favorable compared with the FE and AFE phases. The elastic constants are calculated to be 76.93, 73.75, and 54.67 GPa for 1D FE NbOCl 3 , NbOBr 3 , and NbOI 3 respectively by using density-functional perturbation theory (DFPT) method, suggesting their flexible 1D structures. The thermal stability of the ferroic states is confirmed by the small energy fluctuations and structural variations observed during 5ps-long AIMD simulations at 300 K (see Supplementary Fig. 3). It is worth noting that the 1D AFE states are only slightly higher in energy (<3 meV) than the FE ones, suggesting that metastable 1D AFE NbOX 3 might also be experimentally accessible.

Ferroelectric polarization and switching barrier
The central symmetry of 1D FE NbOX 3 is broken by the cooperative displacements of the Nb ions in both NbOX 4 channels, leading to a significant spontaneous polarization, P, along the chain direction. The magnitude of P is calculated by Berry phase through the modern theory of polarization 47,48 3 54 , and some lead-zirconate-titanate (PZT) 55,56 . The piezoelectric coefficients along the polar direction are calculated to be 2.69, 2.71, and 3.01 C per m 2 for 1D NbOCl 3 , NbOBr 3 , and NbOI 3 respectively by using DFPT method. In the 1D AFE phases, the Nb opposite polar displacements in the two NbOX 4 channels cancel out, giving an overall vanishing P. Their polar strength can be estimated from the magnitudes of d Nb in the NbOX 4 channels, which are close to what found for the FE phases (see Table 1).
In a FE material, the direction of the electric polarization can be switched by an external electric field, hence it is of paramount importance to understand the kinetics of the polarization reversal process in 1D NbOX 3 . In order to explore the most favorable transition pathway connecting the two energetically degenerate states with opposite polarization directions, we systematically scanned the energy surface of 1D NbOX 3 with respect to two independent reaction coordinates, d Nb1 and d Nb2 , namely the polar displacements of Nb in the individual channels. As shown in Fig. 2a-c, two low-energy FE and AFE states with opposite polarized directions (namely, S FE /S′ FE and S AFE /S′ AFE ) are found along the diagonal directions, and they are separated by a highenergy PE state. This is consistent with the predicted two soft modes of 1D PE NbOX 3 , which spontaneously drive the PE state to a lower-energy state (either FE or AFE). FE switching would then proceed by crossing the AFE metastable state, instead of being a direct transition going through the PE configuration, due to the large energy difference between the FE and PE states. The NEB method is further exploited to refine the transition pathways from S FE to S′ FE and to determine accurately the transition barriers. As shown in Fig. 2d, d Nb in one of the S FE channels is first reversed, while that in another channel remains unchanged, leading the system to an intermediate AFE state. Then, d Nb in the unchanged channel gets reversed as well, realizing the transition from S FE Fig. 2 Energy map and NEB results. a-c The energy contour plot (in meV) of the 1D NbOX 3 unit cell as a function of the polar displacements of two Nb ions for NbOCl 3 , NbOBr 3 , and NbOI 3 , respectively. The energy of the PE phases is set to zero. d Energy profiles along the NEB pathway for the polarization switching of 1D FE NbOX 3 . The inset shows the structure evolution from FE to AFE states through the transition state (TS).
to S′ FE . The overall transition barriers are calculated to be only 47, 34, and 16 meV for 1D FE NbOCl 3 , NbOBr 3 , and NbOI 3 respectively. These values are much smaller than those of many widely studied FE materials, such as PbTiO 3 , BaTiO 3 , and 2D In 2 Se 3 57-59 .

Electronic properties
Having explored the polar nature of 1D NbOX 3 , we further investigate the origin of the 1D ferroelectricity and antiferroelectricity. 1D NbOX 3 nanowires are found to be semiconductors with band gaps ranging between 1.4 eV and 4.3 eV and display similar band compositions, as shown in Fig. 3a-f and Supplementary Figs. 4-6. By taking 1D FE NbOBr 3 as an example, we plot the orbitalresolved band structures for Br, O, and Nb atoms in 1D NbOBr 3 (see Fig. 3d-f). The Nb-d x 2 y 2 /Br-p x/y , Nb-d xz/yz /O-p x/y , and Nb-d z 2 / O-p z hybridization can be clearly seen below the Fermi level and it is expected from the Nb-Br and Nb-O bonding geometry. The mixing of Nb 5+ empty d orbitals and O 2p orbitals along the polar direction thus dominates the emergence of the FE/AFE states in 1D NbOBr 3 , namely the d 0 principle is found for this class of compounds, in analogy with the well-known FE perovskite oxides 60,61 . Figure 3g and Supplementary Fig. 7 present the calculated -COHP (crystal orbital Hamilton population) integral for 1D NbOBr 3 and NbOCl 3 /NbOI 3 , which further confirms that the FE/ AFE states are stabilized by the enhanced orbital hybridization of Nb-d xz/yz /O-p x/y and Nb-d z 2 /O-p z . Furthermore, we find the pdπ interactions of Nb-d xz/yz /O-p x/y to show more prominent enhancements than the pdσ interaction of Nb-d z 2 /O-p z . This is a signature also found for the O-Ti interactions in BaTiO 3 and PbTiO 3 62 .

Monte Carlo simulations
For practical applications, the critical temperature, T C , of 1D FE NbOX 3 needs to be high enough so that their polarization can persist above room temperature. The effect of finite temperature is then investigated via Monte Carlo (MC) simulations based on a Landau-Ginzburg model 20 , in which the energy, E, of 1D NbOX 3 can be expressed by an expansion of the order parameter d i (d Nb of ith NbO 2 Br 4 octahedron), where the first term describes the on-site potential energy and the parameters A, B, and C can be obtained by fitting the energy-d Nb double-well curve (see Fig. 4a). The other two terms represent the harmonic interactions between two neighboring dipole moments ( Supplementary Fig. 8) along the x and y directions, respectively. The coefficients D x and D y can be estimated based on the nearestneighbor approximation as shown in Supplementary Fig. 9. All the parameters (A, B, C, D x , and D y ) for 1D FE NbOX 3 are summarized in Supplementary Table 1. The limitation of this model is discussed in Supplementary Note 1.
With the above effective Hamiltonian, Eq. (1), the phase transition at finite temperature can be investigated by MC simulations. Figure 4b plots the averaged polar displacements of Nb ions as a function of temperature. Clearly, this quantity drops abruptly at a temperature close to the T C , suggesting the occurrence of a phase transition. The T C are then calculated from the singular point of the specific heat (see Supplementary Fig. 10) to be 630, 520, and 320 K for NbOCl 3 , NbOBr 3 , and NbOI 3 , respectively. The robust room-temperature ferroelectricity for 1D NbOX 3 can also be validated by the averaged d Nb during the last 2 ps in our AIMD simulations at 300 K ( Supplementary Fig. 11), where no significant reduction of d Nb is observed under thermal perturbation. Notably, the Tc of 1D FE NbOX 3 follows the NbOCl 3 > NbOBr 3 > NbOI 3 order, mirroring the order of their double-well potential depths.
Finally, we investigate the electric field induced FE transition by MC simulations. When an electric field, E, is applied along the polar direction, an additional energy term -E(∑d i Z i * ) must be incorporated into Eq. (1), where Z i * is the Born effective charge. Figure 4c presents the simulated FE hysteresis loops. The critical electric fields that trigger the depolarization are determined to be around 0.61, 0.33, and 0.04 MV per cm for NbOCl 3 , NbOBr 3 , and NbOI 3 , respectively, which are comparable or even smaller than those of reported FE materials including HfO 2 (~1 MV per cm) 63 , HfZrO 4 (~1.2 MV per cm) 64 , CuInP 2 S 6 (~10 MV per cm) 65 , and 2D In 2 Se 3 (6−10 MV per cm) 66 . The small coercive fields of 1D FE NbOX 3 will substantially facilitate the FE switching under low electric voltage, promising great potential for low energy-cost FE devices.
In summary, we have predicted the long-sought intrinsic 1D ferroelectricity in 1D NbOX 3 nanowires. 1D NbOX 3 is highly likely to be exfoliated from the experimentally synthesized vdW bulk phases due to the small binding energies. Notably, 1D NbOX 3 possesses great dynamical and thermal stabilities, sizable spontaneous polarizations, low switching barriers and coercive fields, and above room temperature Tc, holding great potentials for applications in nanoscale FE devices such as high-density nonvolatile memories. In addition, the double-channel 1D structure also enables AFE metastable states in 1D NbOX 3 , offering a great platform to explore complex ferroic orders down to the 1D limit. The polarized nature of 1D NbOX 3 originates from the d 0 rule, namely from the hybridization of the Nb-d xz/yz /O-p x/y and Nb-d z 2 / O-p z orbitals, a mechanism similar to that found in conventional FE materials, like PbTiO 3 and BaTiO 3 . Our findings highlight a class of intrinsic 1D FE compounds with extraordinary ferroelectricity and point to a feasible route for the exploration of exotic 1D physics from the vdW crystals containing 1D building blocks.

Geometry optimization and electronic structure calculations
Our DFT calculations were carried out using the Vienna ab initio simulation package (VASP) [67][68][69] . The electron-core interaction was described by using the projector augmented wave (PAW) method 70 and the exchange and correlation energy was treated with the generalized gradient approximation (GGA) parameterized by Perdew, Burke, and Ernzerhof (PBE) 71,72 . The vdW interactions were included throught the DFT-D3 scheme 73 . The cutoff energy for the plane-wave expansion was set to 520 eV and the Brillouin zone was sampled by a 1 × 12 × 1 mesh. The geometries were fully optimized until the residual forces and the energy were converged to 0.005 eV per Å and 1 × 10 −6 eV, respectively. A vacuum region always >15 Å was introduced to avoid the interaction between the neighboring periodic images. The Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) was also utilized to obtain an accurate description of the electronic properties 74 .

Phonon, NEB, and AIMD calculations
Phonon spectra were calculated based on 1 × 4 × 1 supercells by using the finite displacement method as implemented in the PHONOPY code 75 . The threshold for energy is tightened to 10 −9 eV to get the accurate forces. The nudged elastic band (NEB) method 76 was applied to study the FE phase transition, with 17 images used in total. Ab initio molecular dynamics (AIMD) simulations were carried out using 1 × 8 × 1 supercells for a total of 5 ps with a time step of 1.0 fs. The NVT ensemble is used in the simulations with the temperature controled by the Langevin thermostat 77,78 .

Monte Carlo simulations
Metropolis-algorithm MC simulations were performed for a periodic 1D supercell containing 20,000 unit cells. The first 2 × 10 8 MC steps were used for equilibration, followed by additional 2 × 10 8 MC steps to obtain the thermal averages. Each simulation was repeated 100 times and then the results were averaged to eliminate numerical errors.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. Fig. 4 Double-well potential and MC simulation results. a Energy-d Nb double-well plots connecting the PE and the two degenerate FE phases (-P and P states) in 1D NbOX 3 . The DFT data is fitted using the Landau-Ginzburg model. b The averaged Nb polar displacements of 1D FE NbOX 3 as a function of temperature as obtained from MC simulations. c The ferroelectric hysteresis loops of 1D FE NbOX 3 simulated by MC calculations at room temperature.