# Tunable vertical ferroelectricity and domain walls by interlayer sliding in
$\beta$-ZrI2

Xiaonan Ma  Physics Department, State Key Laboratory of Advanced Special
Steel, and International Center of Quantum and Molecular Structures, Shanghai
University, Shanghai 200444, China  Materials Genome Institute, and Shanghai
Key Laboratory of High Temperature Superconductors, Shanghai University,
Shanghai 200444, China  Chang Liu  Physics Department, State Key Laboratory of
Advanced Special Steel, and International Center of Quantum and Molecular
Structures, Shanghai University, Shanghai 200444, China  Materials Genome
Institute, and Shanghai Key Laboratory of High Temperature Superconductors,
Shanghai University, Shanghai 200444, China  Wei Ren<EMAIL_ADDRESS>Physics Department, State Key Laboratory of Advanced Special Steel, and
International Center of Quantum and Molecular Structures, Shanghai University,
Shanghai 200444, China  Materials Genome Institute, and Shanghai Key
Laboratory of High Temperature Superconductors, Shanghai University, Shanghai
200444, China  Sergey A. Nikolaev<EMAIL_ADDRESS>Tokyo Tech
World Research Hub Initiative (WRHI), Institute of Innovative Research, Tokyo
Institute of Technology, 4259 Nagatsuta, Midori-ku, Yokohama, Kanagawa
226-8503, Japan  Laboratory for Materials and Structures, Tokyo Institute of
Technology, 4259 Nagatsuta, Midori-ku, Yokohama, Kanagawa 226-8503, Japan
National Institute for Materials Science, MANA, 1-1 Namiki, Tsukuba,Ibaraki
305-0044, Japan

###### Abstract

Vertical ferroelectricity where a net dipole moment appears as a result of in-
plane ionic displacements has gained enormous attention following its
discovery in transition metal dichalcogenides. Based on first-principles
calculations, we report on the evidence of robust vertical ferroelectricity
upon interlayer sliding in layered semiconducting $\beta$-ZrI2, a sister
material of polar semimetals MoTe2 and WTe2. The microscopic origin of
ferroelectricity in ZrI2 is attributed to asymmetric shifts of electronic
charges within a trilayer, revealing a subtle interplay of rigid sliding
displacements and charge redistribution down to ultrathin thicknesses. We
further investigate the variety of ferroelectric domain boundaries and predict
a stable charged domain wall with a quasi-two-dimensional electron gas and a
high built-in electric field that can increase electron mobility and
electromechanical response in multifunctional devices. Semiconducting
behaviour and a small switching barrier of ZrI2 hold promise for novel
ferroelectric applications, and our results provide important insights for
further development of slidetronics ferroelectricity.

##  I Introduction

Ferroelectric materials with spontaneous electric dipole moments switchable by
an external electric field offer a broad range of technological applications,
such as non-volatile memories, field-effect transistors, and active elements
in electromechanical and electrooptical devices.ferro1 ; ferro2  The ability
to switch electric polarization is a key ingredient in modern nanotechnology,
where the need for further reduction of individually polarized domains towards
the atomic scale, as well as for the ease of their switchability has been been
constantly growing.

Thinning down ferroelectrics is one promising direction to miniaturize
electronic devices, and several ferroelectric materials have been found to
maintain macroscopic polarization in ultrathin films.thin1 ; thin2 ; thin3 ;
thin4 ; thin5  Furthermore, the development of van der Waals assembly has
enabled heterostructure engineering,assembly1 ; assembly2  where physical
properties are tuned to a desired functionality by combining different
individual layers, which can be used to design two-dimensional ferroelectrics
from non-ferroelectric parent compounds.boron1 ; boron2  Other venues to
overcome the challenges can be found in van der Waals layered materials that
possess properties favourable for tailoring ferroelectricity, such as
durability against strain and surface functionalization.vdwreview  Robust
polarization in these systems was shown to sustain down to atomic thicknesses
and can provide new building blocks for functional heterostructures.thin33 ;
thin44 ; wte2switch ; polarmote2

Following the immense growth of activity in two-dimensional systems, layered
transition metal dichalcogenides (TMDs) have recently drawn great attention
due to their diverse physical properties, ranging from extremely large
magnetoresistance magneto1  and superconductivitysupermote2  to the
topological electronic states.zhijun ; higherorder  Recent studies have
demonstrated that the out-of-plane switchable polarization originating from
interlayer sliding exists in polar Weyl semimetals MoTe2 and WTe2.polarmote ;
wte1  Although the value of spontaneous polarization is small and can be
partially screened by metallic states, its rigidity upon interlayer sliding
and importance for potential applications as a ferroelectric memory have
prompted the urge to search for novel vertical ferroelectrics. Despite rising
research activity in this field dubbed slidetronics, only a few materials have
been discovered so far, such as MoTe2, WTe2, CuInP2S6,cuinps  In2Se3,inse1 ;
inse2 , and VS2.renwei  Exploration of new layered ferroelectrics elucidating
the microscopic origin of electric polarization and their intrinsic
properties, such as domain walls, are in high demand for further development
of slidetronics.

Figure 1: Relation of the T0, Td and 1T′ phases of ZrI2. a Crystal structures:
1T′, designed T0, and Td phases. Side and top views presented to the right
show the interlayer zigzag I–I bonds (black dotted lines) and the intralayer
zigzag Zr–Zr chains (blue solid lines). Trilayers with the clockwise and
counterclockwise rotated ZrI6 octahedra are denoted as M (Minus) and P (Plus),
respectively. The $+/-$ signs stand for positive and negative displacements of
the interlayer zigzag I–I bonds (with respect to the lower trilayer). The
$\Gamma_{2}^{-}$ and $\Gamma_{4}^{+}$ modes relate the $Pnma$ structure to the
$Pmn2_{1}$ (Td) and $P2_{1}/m$ (1T′) phases, respectively. b Band structure of
the designed T0 phase and partial charge densities. The Zr $d$ and I $p$
states are highlighted with blue and pink colours, respectively. c Energy
profile obtained from the devised T0 phase as a function of interlayer
displacement and the monoclinic $a-c$ angle. d Enlarged phonon spectrum of the
T0 phase of ZrI2.

In this work, we report on the first-principles evidence of robust vertical
ferroelectricity in layered ZrI2, a semiconducting counterpart of
isostructural semimetal TMDs MoTe2 and WTe2. Being several orders of magnitude
smaller than in conventional ferroelectrics, the out-of-plane polarization in
ZrI2 is found to be rigid upon interlayer sliding, and the low energy barrier
for its ferroelectric switching combined with a small band gap can hold out
the prospect for novel slidetronics applications. Our theoretical study
reveals that the ferroelectric activity in ZrI2 stems from a subtle interplay
of charge redistribution and ionic displacements, providing important insights
on the origin of vertical ferroelectricity in layered materials and justifying
its persistence in the ultrathin limit. We investigate the complexity of
domain boundaries in multidomain structures of ZrI2 that arise due to the
breaking of stacking sequences. Our results demonstrate the formation of
stable charged domain walls with a quasi-two-dimensional electron gas and a
high built-in electric field that can be put to good use in multifunctional
devices.

##  II RESULTS

###  II.1 Ferroelectric activity in ZrI2

Historically, there have been identified three polymorph forms of ZrI2 studied
by Guthrie and Corbett: $\alpha$ ($P2_{1}/m$)struc1 , $\beta$
($Pmn2_{1}$)struc2 , and $\gamma$ ($R\bar{3}$)struc3  phases. The
$\gamma$-phase is comprised of the antiprism clusters formed by six Zr atoms
which are surrounded by twelve iodine atoms, whereas the $\alpha$\- and
$\beta$-phases have layered structures and are reported to be isostructural,
respectively, to the 1T′ and Td phases of MoTe2 and WTe2. The latter phases
consist of the buckled I–Zr–I trilayers coupled by weak van der Waals
interactions, where the Zr atoms form the zigzag chains in each layer, as
shown in Fig. 1a. Among the three polymorphs, only the $\beta$-phase has a
polar structure and can be expected to reveal ferroelectric properties.

Similar to its sister compounds MoTe2 and WTe2, the $\alpha$\- and
$\beta$-phases of ZrI2 (hereafter referred to as 1T′ and Td, respectively) can
be derived from the devised parent T0 phase with the $Pnma$ structure by
distorting the orthorhombic $a-c$ angle or by sliding adjacent trilayers,
respectively (Fig. 1a). All three phases are found to be semiconducting with
an indirect energy band gap $E_{g}\sim 0.20$ eV (see Supplementary Figure 1).
As shown in Fig. 1b, the states near the Fermi level correspond to the Zr $d$
states that form strong metal-metal bonding along the zigzag chains, and the
bands below represent the I $p$ states. In contrast to semimetals MoTe2 and
WTe2 which demonstrate a sizeable mixing of the anionic $p$ and metal $d$
states at the Fermi level that alongside with spin-orbit coupling stabilizes
the Weyl points,zhijun  the I $p$ states in ZrI2 are pushed way below the
Fermi level due to much longer intralayer I–I bonds, without affecting a small
gap opened in the Zr $d$ states.

In Fig. 1c, we computed the energy landscape as a function of the monoclinic
distortion and interlayer sliding starting from the devised T0 structure and
fixing the unit cell volume and atomic arrangement in each trilayer. For both
types of displacement, a double-well structure is clearly seen in the
calculated energy profile, where the center with higher energy corresponds to
the parent T0 phase. By this means, interlayer sliding leads to two local
minima at $\Delta\sim\pm 0.71$ Å representing two polar Td structures, whereas
rotating the monoclinic $a-c$ angle stabilizes two ferroelastic 1T′-I and
1T′-II phases at 84.5∘ and 95.5∘, respectively. The corresponding energy
barrier connecting the twin structures via direct pathways through the T0
phase is $\sim 5.3$ meV/u.c for both the Td and 1T′ phases, and the energy
difference between the 1T′ and Td phases is found to be small, less then 0.05
meV/u.c. ($E_{T_{\mathrm{d}}}-E_{1T^{\prime}}=-0.17$ meV/u.c. for the fully
optimized structures).

The calculated phonon spectra of the T0 phase presented in Fig. 1d clearly
indicate two types of lattice instabilities: an optical zone-centered mode
$\Gamma_{2}^{-}$ and a linear phonon mode $\Gamma_{4}^{+}$ along the
$\Gamma$-X direction carrying an elastic instability. From group-theoretical
analysis (see Supplementary Table 3), it follows that the two modes with
irreducible representations $\Gamma_{2}^{-}$ and $\Gamma_{4}^{+}$ transform
the $Pnma$ phase to the $Pmn2_{1}$ and $P2_{1}/m$ structures, respectively.
Indeed, the atomic displacement corresponding to the $\Gamma_{2}^{-}$
instability include an in-plane sliding displacement of the alternating
trilayers, while the $\Gamma_{4}^{+}$ mode amounts to a shear distortion of
the unit cell, resulting in a ferroelastic structure.

Figure 2: Ferroelectric polarization in the Td phase of ZrI2. a Td↑ and Td↓
structures obtained from the devised T0 phase by interlayer sliding. The M
trilayers are denoted with blue rectangles, and the P trilayers with positive
and negative sliding displacements are highlighted with orange and yellow
rectangles, respectively. b Double-well energy profile and electric
polarization along the $c$ axis calculated from the Berry phase theory as a
function of interlayer sliding in the devised T0 phase.

The T0, 1T′ and Td phases can be characterized by the way the trilayer
sequences are stacked in a layered structure. In the T0 phase with the
$D_{2h}$ symmetry, the I–Zr–I trilayers are invariant under the mirror
reflection with respect to the $b$ axis accompanied by a half-lattice vector
translation, $\\{M_{b}|\frac{\bm{b}}{2}\\}$, and the adjacent trilayers can be
transformed by $\\{M_{a}|\frac{\bm{a}}{2}+\frac{\bm{b}}{2}\\}$ or
$\\{M_{c}|\frac{\bm{a}}{2}\\}$. Consequently, the trilayers are stacked with
two alternating orientations where the ZrI6 octahedra twist either clockwise
or counterclockwise, denoted by taking notations of Ref. polarmote  as M
(Minus) and P (Plus), respectively. In the Td phase, the interlayer sliding
along the $a$ axis breaks inversion symmetry reducing the point symmetry to
$C_{2v}$, so that the adjacent trilayers are invariant under
$\\{M_{b}|\frac{\bm{b}}{2}\\}$ and can be connected by
$\\{M_{a}|\Delta+\frac{\bm{a}}{2}+\frac{\bm{b}}{2}\\}$. On the other hand, the
monoclinic distortion along the $a$ axis in the 1T′ phase preserves inversion
symmetry but lowers the symmetry to $C_{2h}$ with the mirror plane along the
$b$ axis. The twin structures in both phases driven by the $\Gamma_{2}^{-}$
and $\Gamma_{4}^{+}$ modes can be further classified by the way the I ions in
adjacent trilayers move relative to each other: in a positive ($+$) or
negative ($-$) direction, if seen with respect to the lower layer, as shown in
Fig. 1a. Then, if going up from the bottom trilayer, the patterns M$+$P$+$ and
M$-$P$-$ will correspond to the ferroelastic 1T′-I and 1T′-II structures,
while the patterns M$-$P$+$ and M$+$P$-$ will stand for the polar Td↑ and Td↓
phases, respectively. With that said, the 1T′ phase has only one type of the
interlayer displacement, whereas the Td phase contains both $+$/$-$ types that
break inversion symmetry.

Figure 3: Schematics of electronic charge displacements in ZrI2. a Wannier
functions (top and side views) constructed for the occupied I $p$ (left) and
Zr $d$ (middle and right) states. The Wannier functions are presented for one
trilayer, and only one $p$ Wannier function at each I site is shown for
clarity. Centers of the $d$ orbital Wannier functions are denoted with red
circles. b Electronic charge displacements in the centrosymmetric and non-
centrosymmetric phases due to interlayer sliding. Vertical shifts of the
Wannier centers off their specified atomic positions are schematically shown
with arrows. Centers of the $d$ orbital Wannier functions are denoted with
green circles. The alternating longer $d_{\mathrm{L}}$ and shorter
$d_{\mathrm{S}}$ zigzag I–I bonds are shown with blue and pink colours,
respectively. Inversion center is marked with the black cross.

To investigate ferroelectric activity of the Td phase, we considered the
devised T0 phase modulated along the $\Gamma_{2}^{-}$ mode while keeping fixed
the unit cell volume and atomic arrangement in each trilayer. Electric
polarization calculated from the Berry phase theory as a function of
interlayer sliding is presented in Fig. 2 and shows that a net dipole moment
appears strictly along the vertical axis and changes its direction within the
double-well profile when going between positive and negative sliding
displacements. Since interlayer sliding amounts to rigid ionic translations,
electric polarization is purely of electronic origin. The absolute value of
$P_{\mathrm{el}}$ corresponding to the energy minimum of the constrained Td
structure is 0.22 $\mu$C$\cdot$cm-2 (0.24 $\mu$C$\cdot$cm-2 for the fully
optimized Td phase). Despite being a few orders of magnitude smaller compared
to conventional ferroelectrics, electric polarization in ZrI2 is caused solely
by interlayer sliding and is switchable without vertical ionic displacements.
Importantly, the switching barrier between two ferroelectric Td structures of
ZrI2 is much lower than in conventional ferroelectrics (34 meV/atom for BaTiO3
and 67 meV/atom for PbTiO3).rabe

Given the rigidity of a single trilayer, intralayer antiferroelectric
structures are found to be highly unstable relaxing to the non-polar T0 phase
and confirming the robustness of ferroelectric properties in the Td phase.
From molecular dynamics simulations, the lower bound for the Curie temperature
of the Td phase can be estimated as $\sim 400$ K.

It is worth noting that interlayer sliding alone is not sufficient to explain
the appearance of the net dipole moment in the Td phase. Several previous
studies attributed the ferroelectric activity in WTe2 to the vertical charge
transfer between adjacent trilayers that occurs upon interlayer sliding. liju
; ren1  While the Berry phase analysis presented above shows that the
ferroelectric activity in ZrI2 is purely of electronic origin, the calculated
change in the charge density demonstrates that the P and M trilayers are rigid
in the sense that a sheer ionic displacement is accompanied by the same
intralayer shift of the electronic charges, and the corresponding change in
the interlayer charge density representing interlayer bonding indicates that
interlayer sliding results solely in the charge density redistribution between
the trilayers (see Supplementary Figure 2).

The microscopic origin of the uncompensated dipole moment in the Td phase of
ZrI2 can be elaborated by examining displacements of the electronic charges
driven by interlayer sliding. The electronic contribution to the net dipole
moment can be expressed as a sum of the centers of the Wannier functions,
$\mathbf{r}_{n}$, as $\mathbf{d}_{el}=-2e\sum_{n}\mathbf{r}_{n}$ ($e$ is a
positive electron charge). As shown in Fig. 3a, the Wannier functions of the I
$p$ states are well localized in the vicinity of the I ions, whereas the
Wannier functions corresponding to the Zr $d$ states are not centered at the
Zr sites and are, instead, located at the Zr triangles, reflecting strong
metallic bonding in the zigzag chains. Due to the hybridization effects, the
Wannier functions move off the atomic sites, and their vertical displacements
can be described by the shift of the Wannier centers relative to their
symmetry specified positions, $\delta_{d}$ and $\delta_{\mathrm{I}}$, as
depicted in Fig. 3b. As a result of the intralayer hybridization, the
$p$-orbital Wannier functions at the top and bottom I sites tend to shift
inwards the trilayer, which is accompanied by the shift of the $d$ orbital
Wannier functions at the corresponding Zr triangles. The intralayer
hybridization is strong and determines an overall displacement of electronic
charge centers. However, the charge density distribution is also affected by
interlayer bonding, in particular by steric I–I interactions. In the
centrosymmetric phase, the top and bottom I1 and I2 sites are equivalent in
each trilayer resulting in symmetric shifts of the Wannier centers, so that
all corresponding $\delta$’s are equal giving zero $\mathbf{d}_{el}$. The
calculated charge centers shifts for the centrosymmetric T0 phase are
$\delta_{d}=0.0712$ Å, $\delta_{\mathrm{I1}}=0.1660$ Å, and
$\delta_{\mathrm{I2}}=0.2197$ Å (for $\delta_{\mathrm{I}}$ is an average over
three $p$ orbitals). When sliding adjacent trilayers, inversion symmetry is
broken, and the top and bottom I sites in each trilayer become non-equivalent.
Notably, interlayer sliding modifies steric bonding in the interlayer zigzag
I–I chains leading to the alternating shorter and longer I–I bonds that change
depending on the sliding direction (Fig. 3b). These bonds reverse their order
between adjacent trilayers, so that the top and bottom I ions in each trilayer
experience asymmetric bonding environment giving rise to non-equivalent shifts
of their charge centers. The shorter the bond length, the more the electronic
charge density is pulled outwards the trilayer reducing its vertical shift.
Since interlayer bonding is weak, the resulting shifts and the net dipole
moment are small. The calculated shifts in the Td↑ phase are
$\delta_{d1}=-0.0756$ Å, $\delta_{d2}=0.0687$ Å, $\delta_{\mathrm{I1}}=0.1661$
Å, $\delta_{\mathrm{I2}}=0.2176$ Å, $\delta_{\mathrm{I3}}=-0.2201$ Å,
$\delta_{\mathrm{I4}}=-0.1662$ Å, and the net dipole moment is $0.0585$
$e\cdot$Å, in perfect agreement with the value of $0.0579$ $e\cdot$Å ($0.243$
$\mu$C$\cdot$cm-2) obtained from the Berry phase theory. It should be noted
that in the 1T′ phase, the interlayer zigzag I–I chains also alternate between
the trilayers, but the corresponding top and bottom I ions are always
connected by the same type of the I–I bonds and have symmetric shifts of their
charge centers (see Supplementary Figure 3).

The hybridization also causes the in-plane shifts of the charge centers within
a trilayer. Since the interlayer zigzag I–I chains alternate periodically in
the bulk structure, the in-plane shifts will sum up to zero. Nevertheless, in
the case of only two trilayers, such cancellation will not be perfect giving
rise to the in-plane electric polarization (see Supplementary Figure
9).prbzri2

###  II.2 Domain walls

When cooled down below the transition temperature, ferroelectric materials
tend to form complex multidomain structures with different orientations of
macroscopic polarization, and nanometric-scale domain walls (DW) develop at
the domain boundaries. Structural and functional properties of DWs, such as
width and energy barrier across the DW, substantially impact polarization
switching process and the DW mobility.dwreview

Figure 4: First-principles calculations of the charged DW in $\beta$-ZrI2. a a
$1\times 1\times 8$ supercell with the Td/Td arrangement containing the head-
to-head (HH) and tail-to-tail (TT) DWs at the center and boundaries,
respectively. The M trilayers are denoted with blue rectangles, and the P
trilayers with positive and negative sliding displacements are highlighted
with orange and yellow rectangles, respectively. b Electronic band structure
and charge densities corresponding to the valence band maximum (VBM) and
conduction band minimum (CBM1 and CBM2). The isosurfaces correspond to 0.01
$e\cdot$Å-3. c Planar and macroscopic average electrostatic potentials across
the supercell and planar average charge densities of VBM and CBM.

As follows from the results of first-principles calculations, the energy
difference between the Td and 1T′ phases of ZrI2 is small, and both forms can
coexist in a sample producing different types of DW structures. Following the
symmetry analysis reported in Ref. polarmote , a stacking sequence
$($T${}_{\mathrm{d}})_{n}/($1T${}^{\prime})_{m}$ has the $Pm$ space symmetry
($C_{1h}$ point group) that can be derived from both the $P2_{1}/m$ and
$Pmn2_{1}$ space groups by taking the corresponding $\Gamma_{2}^{-}$ and
$\Gamma_{3}$ modes, respectively, which represent the breaking of the $+/-$
patterns. Owing to the rigidity of a single trilayer, there exist several
possible ways of vertical stacking at the boundaries, which adopt the patterns
of the Td and 1T′ stacking sequences. For example, when two ferroelectric Td↑
and Td↓ domains are stacked along the $c$ axis ($...$M$-$P$+$/M$+$P$-$$...$),
the 1T′-I pattern will be formed at the boundary. Or, similarly, for a
1T′-I/1T′-II stacking ($...$M$+$P$+$/M$-$P$-$$...$) the DW will have the Td↑
pattern. Because the $Pm$ space group is polar, there are two possibilities
for vertical DWs in ZrI2: a (polar) DW with the Td pattern at the 1T′/1T′
stacking will exhibit ferroelectric properties due to the local inversion
symmetry breaking at the boundary, whereas the Td/Td and Td/1T′ arrangements
will form a (non-polar) charged DW with the 1T′ pattern arising from
polarization discontinuity.

When a normal component of electric polarization changes across the domain
boundary as in the latter case of the Td/Td and Td/1T′ arrangements, the
interface will accommodate a high density of bound charges that give rise to
large depolarizing fields suppressing ferroelectricity. If these charges are
compensated by free carriers or electron-hole transfer across the gap, a
stable charged DW can form as a movable ultrathin conductive layer. Taking a
nominal value of $P_{\mathrm{el}}$ in the Td phase of ZrI2 and assuming a DW
width $w$ of several nanometers, the bound charge concentration can be
estimated as $2P_{\mathrm{el}}/ew\sim 3\cdot 10^{19}$ cm-3, which is an
appreciable value for a small band gap semiconductor.cdwreview

The results obtained for the Td/Td arrangement are summarized in Fig. 4, where
the $\mathrm{P}+\mathrm{M}+$ and $\mathrm{P}-\mathrm{M}-$ patterns at the
boundaries form, respectively, the 180∘ head-to-head (Td↑/Td↓) and tail-to-
tail (Td↓/Td↑) charged DWs. Due the presence of bound charges at the
interfaces, the bottom of the conduction states and the top of the valence
states approach the Fermi level, thus providing electrons and holes for
screening at the head-to-head and tail-to-tail DWs, respectively (Fig. 4b).
This band bending driven by electron-hole transfer across the gap further
leads to a potential difference. The calculated potential profile across the
heterostructure reveals a small barrier of $\sim 0.24$ eV with the minimum and
maximum corresponding to the head-to-head and tail-to-tail DWs (Fig. 4c).
Assuming that a domain width can be varied up to tens of nanometers, the
induced potential will give rise to a large built-in electric field of
hundreds kV$\cdot$cm-1 across the domain. Combined with a strong dielectric
anisotropy of ZrI2 (see Supplementary Table 4), the dielectric and
piezoelectric responses in ZrI2 can be greatly enhanced and controlled by
changing the density of charged DWs.enhanced  The calculated formation energy
of the charged DWs well converges to the estimation
$E_{\mathrm{DW}}=2P_{\mathrm{el}}E_{g}/e\sim 1.0$ mJ$\cdot$m-2 as a function
of the supercell size (see Supplementary Figure 11). This value is
significantly smaller than in conventional ferroelectrics (35 mJ$\cdot$m-2 for
the 90∘ DWs and 132 mJ$\cdot$m-2 for the 180∘ DWs in PbTiO3,dw1  71
mJ$\cdot$m-2 for the 180∘ DW in BiFeO3dw2 , $\sim 70$ mJ$\cdot$m-2for anti-
phase DWs in SmFeO3dw3 ), allowing for the ease to control the DW motion (see
Supplementary Figure 12).

As a direct consequence of strong charge compensation, the charged DWs exhibit
metallic behaviour and can feature noticeable inwall conductivity.quasigas
The tail-to-tail DW at the Td↓/Td↑ interface is found to be less pronounced
slowly decaying over several layers beyond its core, because the band bending
by the valence states is rather weak (Fig. 4c). On the contrary, the head-to-
head DW at the Td↑/Td↓ stacking demonstrates a much higher concentration of
bound charges that are strongly localised at the potential cavity of the
interfacial trilayers forming a quasi-two-dimensional electron gas. Since the
bound charges at the interface are proportional to the normal component of
electric polarization, varying conductivities are expected for different types
of DWs at the Td/Td and Td/1T′ boundaries. Thus, the electron confinement can
be modulated by interlayer sliding that allows to controlling both the
conductivity and carrier mobility at the interfaces.gas1

Taking the estimated charge density at the DW boundary, one can calculate the
Debye length to be $\sim 0.5$ nm at room temperature. The value is rather
small, but this estimate should be taken with care. As shown in Fig. 4, the
charged DWs of ZrI2 reveal a generic quantum structure where the electron
motion perpendicular to the DW is quantized. Since the number of subbands
corresponding to the electron’s inwall propagation is small, the semiclassical
approximations for screening effects can be violated. This feature is
different from DWs in typical oxide ferroelectrics, where the number of
subbands is much larger and the Thomas-Fermi approximation holds
true.debyelength

Figure 5: Polar DW in ZrI2. Electronic band structure and planar average
electrostatic potential for a $1\times 1\times 8$ supercell with the
1T′-II/1T′-I arrangement. The M and P trilayers are denoted with blue and grey
rectangles. The yellow and orange rectangles denote the Td↑ and Td↓
boundaries, respectively.

An example of the polar DW structure is presented in Fig. 5 for the 1T′/1T′
arrangement. The Td pattern formed at the boundary breaks inversion symmetry,
and the $+/-$ stacking discontinuity gives rise to the uncompensated dipole
moment at the interfacial trilayer. Since the 1T′ phase in the bulk is non-
polar, there are no boundary charges, and the 1T′/1T′ stacking remains
semiconducting with the ferroelectric interface.

We have considered the variety of DW structures in ZrI2 that can be promising
for novel slidetronics applications. The presented analysis is limited to the
vertical arrangements and based on the breaking of stacking sequences.
However, it is worth noting that the DW structures have also been reported for
the in-plane boundaries in MoTe2.polarmote  While our calculations show that
the antiferroelectric arrangements in a single trilayer are highly unstable
(see Supplementary Figure 14), the formation of neutral DWs at the Td/Td
structures in ZrI2 is also possible by introducing a minimal strain mismatch
between adjacent domains.

##  III Discussions

Based on extensive first-principles calculations, we have established that
vertical ferroelectricity in $\beta$-ZrI2 is a robust consequence of lattice
instability resolved by interlayer sliding between the adjacent trilayers. The
microscopic origin of electric polarization in ZrI2 is attributed to a subtle
interplay of ionic displacements and charge redistribution, leading to
asymmetric shifts of the electronic charge centers within each trilayer and
thus refuting the previously proposed scenario of interlayer charge transfer.
While the out-of-plane polarization in ZrI2 is found to be several orders of
magnitude smaller than in conventional ferroelectrics, the low energy barrier
for its ferroelectric switching combined with a small band gap can provide
opportunities for further developments in the field of ferroelectric
semiconductors, such as ferroelectric tunnel junctions.junction

Our study can shed light on some aspects of ferroelectric activity in
isostructural TMDs MoTe2 and WTe2, whose properties are summarized in Table
I.expmote2 ; expwte2 ; expwte22  The essential difference between these polar
semimetals and semiconducting ZrI2 is that electric polarization of the latter
is not screened by metallic states, and $\beta$-ZrI2 itself can be regarded as
a new alternative to the layered TMDs, thus enriching the family of
slidetronics ferroelectrics. The presented analysis of ferroelectric activity
applies to semimetallic TMDs, explaining the persistence of electric
polarization down to the ultrathin limit, and, in principle, can be extended
to other vertical ferroelectrics.

We have established the complexity and variability of domain boundaries in
multidomain structures of ZrI2. Importantly, the same types of DWs and
superlattice-like arrangements were previously reported experimentally in
MoTe2polarmote  between the Weyl semimetal Td and higher-order topological 1T′
phases,zhijun ; higherorder  which can provide a proving ground for exploring
interfacial topological states and related quantum phenomena. Given their
semimetallic behaviour, no static charge accumulation is expected at the
boundaries in TMDs. In contrast, the DWs in semiconducting ZrI2 offer another
diversity related to the breaking of stacking sequences and polarization
discontinuity, that greatly adds up to multifunctional aspects of vertical
ferroelectrics. In particular, we have predicted a stable charged DW at the
phase boundaries with a pronounced metallic behavior and a high built-in
electric field that can be used to enhance dielectric and piezoelectric
properties. A head-to-head charged DW was shown to form a high density quasi-
two-dimensional electron gas confined at the interfacial trilayers. Such DWs
can be manipulated in a controllable way and allow to creating ultrathin
conductive layers embedded in a semiconducting matter and to increase electron
mobility in electronic devices, that can bring new functionalities for future
slidetronics applications.

Table 1: Experimental lattice parameters and theoretically predicted electric polarization for the Td phase of ZrI2, MoTe2, and WTe2.  | ZrI2 | MoTe2 | WTe2
---|---|---|---
$a$ (Å) | 3.744struc2  | 3.468expmote2 , 3.458zhijun  | 3.496expwte2 , 3.477expwte22
$b$ (Å) | 6.831struc2  | 6.310expmote2 , 6.304zhijun  | 6.282expwte2 , 6.249expwte22
$c$ (Å) | 14.886struc2  | 13.861expmote2 , 13.859zhijun  | 14.070expwte2 , 14.018expwte22
$P_{\mathrm{el}}$ ($\mu$C$\cdot$cm-2) | 0.243 | 0.058polarmote2  | 0.190wte1
Band gap (eV) | 0.20 | – | –

##  IV Methods

Structural parameters of the $Pmn2_{1}$ and $P2_{1}/m$ phases of ZrI2 were
adopted from previous experimental studies,struc1 ; struc2  and the $Pnma$
phase was devised by symmetrising the $Pmn2_{1}$ structure. One should bear in
mind that notations of the $Pnma$ and $Pmn2_{1}$ structures have different
orders of the lattice vectors, which are related as $(a,b,c)$ and $(b,a,c)$,
respectively.

Group-theoretical analysis was performed using the Bilbao Crystallographic
Server (see Supplementary Table 3).bilbao1 ; bilbao2  The crystal structures
were visualized with VESTA.vesta

###  IV.1 Electronic structure

First-principles calculations were performed using the Vienna ab-initio
simulation package (VASP) vasp  within the framework of projected augmented
waves paw  and Quantum-ESPRESSO (QE) realized in the basis of plane waves.qe
The calculations were carried out using local density approximation (LDA) lda
and generalized gradient approximation (Perdew-Burke-Ernzerhof, PBE,pbe  and
Perdew-Burke-Ernzerhof revised for solids, PBEsolpbesol ) for the exchange-
correlation potential. The calculations in QE were performed with ultrasoft
and norm-conserving pseudopotentials.corso  For all calculations, the valence
state configurations were taken as $4s^{2}4p^{6}5s^{2}4d^{2}$ for Zr and
$5s^{2}5p^{5}$ for I. The plane wave cutoff was set to 500 eV and 900 eV for
VASP and QE, respectively. The Brillouin zone was sampled by a Monkhorst-Pack
$k$-point mesh,mp  $6\times 10\times 3$ for the T0 and 1T′ phases and
$10\times 6\times 3$ for the Td phase. The convergence criteria for the total
energy calculations was set to 10-8 eV, and all structures were optimized with
the total force convergence criteria of 10-8 eV$\cdot$Å-1. The optimized
crystal structures used in the main text are given in Supplementary Table 1.
Electronic band structures of the T0, Td, and 1T′ phases calculated with LDA,
PBE and PBEsol are shown in Supplementary Figure 1. The effect of spin-orbit
coupling was shown to give minor changes (see Supplementary Figure 4a).

To investigate the effect of van der Waals interactions on the structural and
ferroelectric properties of ZrI2, we performed a comparative analysis using
the DFT-D2 correction method of Grimme,grimme1  the DFT-D3 method with Becke-
Jonson damping,grimme2 ; grimme3  the vdW-DF and vdW-DF2 functionals,dion1 ;
dion2 ; dion3 , and the optimized vdW functionals (optPBE, optB88, and
optB86b).optdion1 ; optdion2  The optimized structure parameters for the Td
phase are shown in Supplementary Table 2, and the corresponding band
structures are presented in Supplementary Figure 5. One can see that most of
the functionals show a reasonable agreement with the experimental crystal
structure parametersstruc2  and with the band gap of $\sim 0.1$ eVstruc1
(available for the 1T′ phase only). Significant deviations were obtained for
the PBE and vdW-DF functionals that give larger values of the $c$ axis and
band gaps with smaller values of electric polarization compared to the other
functionals.

The effect of electronic correlations in the Zr $d$ shell was checked within
the LDA$+U$ methodldau  and shown to give minor changes (see Supplementary
Figure 4b). We have also considered the effect of long-range Coulomb
interactions using the $GW$ methodgw1  as implemented in VASP:gw2  a single
shot $G0W0$ approach with 150-200 unoccupied bands and 100-150 frequencies.
The results are presented in Supplementary Figure 5. While the $GW$ method may
seem applicable owing to the large spatial extension of the Zr $4d$ orbitals,
one can see that the band gap is largely overestimated compared to the
experimental estimates, which could be related to the implementation or other
known issues of the $GW$ method in polar materials.gw3

_Given its computational flexibility, PBEsol without spin-orbit coupling was
used for all the results presented in the main text, including the analysis of
ferroelectric properties and DW structures._

High-symmetry $k$-points used for the band structure and phonon dispersion
calculations are $\mathrm{X}=(\frac{1}{2},0,0)$,
$\mathrm{S}=(\frac{1}{2},\frac{1}{2},0)$, $\mathrm{Y}=(0,\frac{1}{2},0)$,
$\mathrm{Z}=(0,0,\frac{1}{2})$, $\mathrm{U}=(\frac{1}{2},0,\frac{1}{2})$,
$\mathrm{R}=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$,
$\mathrm{T}=(0,\frac{1}{2},\frac{1}{2})$ for the $Pnma$ (#62) and $Pmn2_{1}$
(#31) structures; $\mathrm{Z}=(0,\frac{1}{2},0)$,
$\mathrm{B}=(0,0,\frac{1}{2})$, $\mathrm{Y}=(\frac{1}{2},0,0)$,
$\mathrm{C}=(\frac{1}{2},\frac{1}{2},0)$,
$\mathrm{D}=(0,\frac{1}{2},\frac{1}{2})$,
$\mathrm{A}=(-\frac{1}{2},0,\frac{1}{2})$,
$\mathrm{E}=(-\frac{1}{2},\frac{1}{2},\frac{1}{2})$ for the $P2_{1}/m$ (#11)
structure.

###  IV.2 Phonon spectra

The phonon spectra for the T0 phase of ZrI2 were calculated using the method
of frozen phonons as implemented in VASP and Phonopy phonopy  and density
functional perturbation theorydfpt  (DFPT) as implemented in QE. The
calculations within frozen phonons were carried out for a $2\times 4\times 1$
supercell on a $3\times 3\times 3$ $k$-point mesh. The calculations within
DFPT were performed on a $3\times 6\times 2$ $k$-point mesh. The calculated
phonon spectra are summarized in Supplementary Figure 6.

The Born effective charges and the dielectric matrix with and without local
field effects were calculated using DFPT as implemented in VASP.bornvasp

###  IV.3 Electric polarization

Electric polarization was calculated within the Berry phase theory. berry1 ;
berry2  For the results shown in Fig. 2, we have considered the T0 phase
modulated along the $\Gamma_{2}^{-}$ mode while keeping fixed the unit cell
volume and atomic arrangement in each trilayer. Electric polarization of the
fully optimized Td phase obtained within different approximations for the
exchange-correlation potential is summarized in Supplementary Table 2 and
Supplementary Figure 7.

The analysis of electronic charge centers was carried out using maximally
localized Wannier functions as implemented in the Wannier90 package.wan90  The
resulting Wannier functions are obtained by projecting the states below the
Fermi level onto the Zr $d$ and I $p$ atomic orbitals. Contributions from all
sites to the net dipole moment can also be seen in the non-diagonal components
of the Born effective charge tensors (see Supplementary Figure 8).

###  IV.4 Molecular dynamics

Ab-initio molecular dynamics simulations were performed on a $2\times 2\times
1$ supercell for the structure optimized with PBEsol. The NVT ensemble with
the Nosé-Hoover thermostatnose  was chosen to simulate the effects of
temperature, and the calculations were run up to 2000 fs with a time step of 2
fs.

In order to estimate the Curie temperature of ZrI2 while taking into account
the energetic proximity of the Td and 1T′ phases, one has to consider large
supercells and allow for elastic changes in the unit cell’s shape and volume.
In this study, molecular dynamics was performed with the fixed volume of the
supercell, and a possible structural instability of the Td phase was
associated with the deviation of the interlayer I-I distances from their
values at 0 K. The results are summarized in Supplementary Figure 10. At low
temperatures, the interlayer I-I bonds between the middle and upper/lower
trilayers evolve on average around their equilibrium 0 K values. At
temperatures higher than 400 K, the middle trilayer starts to float
alternately between the upper and lower trilayers (the values of the
interlayer I-I bonds with the upper and lower trilayers decrease and increase
by turns). While this behaviour does not unambiguously imply the transition to
the paraelectric T0 phase, it can be associated with the onset of structural
instability. We believe that such approach allows us to consider 400 K as an
effective lower bound for the Curie temperature in the Td phase of ZrI2.

###  IV.5 Domain walls

Calculations for the Td↑/Td↓ heterostructure were performed on a $1\times
1\times 8$, $1\times 1\times 10$, and $1\times 1\times 12$ supercells in the
$Pmn2_{1}$ notation. Calculations for the 1T′/1T′ heterostructure were carried
out for a $1\times 1\times 8$ supercell in the $Pnma$ notation. The Brillouin
zone was sampled by a $10\times 5\times 1$ Monkhorst-Pack $k$-point mesh, and
the heterostructures were optimized with the total force convergence criteria
of 10-5 eV/Å. The DW energy is calculated as
$E_{\mathrm{DW}}=(E_{\mathrm{MD}}-E_{\mathrm{SD}})/2A$, where
$E_{\mathrm{MD}}$ is the energy of a supercell containing the DW,
$E_{\mathrm{SD}}$ is the corresponding single domain energy, and $A$ is the DW
cross-sectional area. To reduce systematic errors, the single domain
structures were optimized for the same supercell. Since the
$($T${}_{\mathrm{d}})_{n}($1T${}^{\prime})_{m}$ heterostructure has the $Pm$
symmetry, there may be a tendency towards monoclinic distortions during
structural optimization. We have compared our results for the fully optimized
structures and the constrained structures with a fixed shape and found only
minor changes in the calculated DW energies and electrostatic potentials. The
results obtained with PBEsol are summarized in Supplementary Figure 11,
showing the convergence of the DW energy and electrostatic potentials with the
respect to the supercell size.

The formation of charged DWs largely depends on the values of the band gap and
electric polarization, which specify the degree of the band bending.
Electronic spectra of the Td↑/Td↓ heterostructure obtained for a variety of
exchange-correlation functionals considered in our study are compared in
Supplementary Figure 13. One can see that the band bending by the conduction
states and the formation of the head-to-head DWs are well defined and
pronounced for most of the functionals. In contrast, the band bending from the
valence states is rather weak, and the formation of the tail- to-tail DWs is
less definite and may be unfavourable. The optPBE and vdW-D2 functionals only
reveal the head-to-head DW, which can be attributed to the reduced value of
electric polarization compared to the other functionals. The calculations with
PBE and vdW-DF do not show the formation of any charged DWs, which can be
related to the fact that these functionals largely overestimate the $c$ axis
and give larger band gaps with smaller values of electric polarization
compared to the other functionals.

The Debye length was estimated using the Debye-Hückel model:debye

$\lambda=\sqrt{\frac{\epsilon_{r}\epsilon_{0}k_{B}T}{ne^{2}}},$ (1)

where $\epsilon_{r}$ is the static dielectric constant along the $c$ axis,
$\epsilon_{0}$ is the vacuum permittivity, $k_{B}$ is the Boltzmann constant,
$T$ is temperature, $n$ is the carrier charge density at the DW, and $e$ is
the elementary charge.

## References

  * (1) Dawber, M., Rabe, K. M., Scott, J. F. Physics of thin-film ferroelectric oxides. _Rev. Mod. Phys._ 77, 1083–1130 (2005).
  * (2) Scott, J. F. Applications of Modern Ferroelectrics. _Science_ 315, 954–959 (2007).
  * (3) Bune, A. V. et al. Two-dimensional ferroelectric films. _Nature_ 391, 874–877 (1998).
  * (4) Fong, D. D. et al. Ferroelectricity in ultrathin perovskite films. _Science_ 304, 1650–1653 (2004).
  * (5) Böscke, T. S., Müller, J., Bräuhaus, D., Schröder, U., Böttger, U. Ferroelectricity in hafnium oxide thin films. _Appl. Phys. Lett._ 99, 102903 (2011).
  * (6) Martin, L. W. & Rappe, A. M. Thin-film ferroelectric materials and their applications. _Nat. Rev. Mater._ 2, 16087 (2017).
  * (7) Wang, H. et al. Direct observation of room-temperature out-of-plane ferroelectricity and tunneling electroresistance at the two-dimensional limit. _Nat. Commun._ 9, 3319 (2018).
  * (8) Geim, A. K. & Grigorieva, I. V. Van der Waals heterostructures. _Nature_ 499, 419–425 (2013).
  * (9) Frisenda, R. et al. Recent progress in the assembly of nanodevices and van der Waals heterostructures by deterministic placement of 2D materials. _Chem. Soc. Rev._ 47, 53–68 (2018).
  * (10) Stern, M. V. et al. Interfacial Ferroelectricity by van der Waals Sliding. _Science_ 372, 1462–1466 (2021).
  * (11) Yasuda, K., Wang, X., Watanabe, K., Taniguchi, T., Jarillo-Herrero, P. Stacking-engineered ferroelectricity in bilayer boron nitride. _Science_ 372, 1458–1462 (2021).
  * (12) Duong, D. L., Yun, S. J., Lee, Y. H. Van der Waals Layered Materials: Opportunities and Challenges. _ACS Nano_ 11, 11803–11830 (2017).
  * (13) Shirodkar, S. N. & Waghmare, U. V. Emergence of ferroelectricity at a metal-semiconductor transition in a 1T monolayer of MoS2. _Phys. Rev. Lett._ 112, 157601 (2014).
  * (14) Fei, R. X., Kang, W., Yang, L. Ferroelectricity and phase transitions in monolayer group-IV monochalcogenides. _Phys. Rev. Lett._ 117, 097601 (2016).
  * (15) Fei, Z. et al. Ferroelectric switching of a two-dimensional metal. _Nature_ 560, 336–339(2018).
  * (16) Yuan, S. et al. Room-temperature ferroelectricity in MoTe2 down to the atomic monolayer limit. _Nat. Commun._ 10, 1775 (2019).
  * (17) Ali, M. N. et al. Large, non-saturating magnetoresistance in WTe2. _Nature_ 514, 205–208 (2014).
  * (18) Qi, Ya. et al. Superconductivity in Weyl semimetal candidate MoTe2. _Nat. Commun._ 7, 11038 (2016)
  * (19) Wang, Z. et al. MoTe2: A Type-II Weyl topological Metal. _Phys. Rev. Lett._ 117, 056805 (2016).
  * (20) Wang, Z., Wieder, B. J., Li, J., Yan, B., Bernevig, B. A. Higher-Order Topology, Monopole Nodal Lines, and the Origin of Large Fermi Arcs in Transition Metal Dichalcogenides $X{\mathrm{Te}}_{2}$ ($X=\mathrm{Mo},\mathrm{W}$). _Phys. Rev. Lett._ 123, 186401 (2019).
  * (21) Huang, F.-T. et al. Polar and phase domain walls with conducting interfacial states in a Weyl semimetal MoTe2. _Nat. Commun._ 10, 4211 (2019).
  * (22) Sharma, P. et al. A room-temperature ferroelectric semimetal. _Sci. Adv._ 5, eaax5080 (2019).
  * (23) Liu, F. et al. Room-temperature ferroelectricity in CuInP2S6 ultrathin flakes. _Nat. Commun._ 7, 12357 (2016).
  * (24) Zhou, Y. et al. Out-of-plane piezoelectricity and ferroelectricity in layered $\alpha$-In2Se3 nanoflakes. _Nano Lett._ 17, 5508–5513 (2017).
  * (25) Cui, C. et al. Intercorrelated In-Plane and Out-of-Plane Ferroelectricity in Ultrathin Two-Dimensional Layered Semiconductor In2Se3. _Nano Lett._ 18, 1253–1258 (2018).
  * (26) Liu, X., Pyatakov, A. P., Ren, W. Magnetoelectric Coupling in Multiferroic Bilayer VS2. _Phys. Rev. Lett._ 125, 247601 (2020).
  * (27) Guthrie, D. H. & Corbett, J. D. Synthesis and structure of an infinite-chain form of ZrI2 ($\alpha$), _J. Solid State Chem._ 37, 256–263 (1981).
  * (28) Corbett, J. D. & Guthrie, D. H. A second infinite-chain form of zirconium diiodide ($\beta$) and its coherent intergrowth with $\alpha$-zirconium diiodide. _Inorg. Chem._ 21, 1747–1751 (1982).
  * (29) Guthrie, D. H. & Corbett, J. D. Two zirconium iodide clusters. Hexazirconium dodecaiodide (Zr6I12) and cesium hexazirconium tetradecaiodide (CsZr6I14). _Inorg. Chem._ 21, 3290–3295 (1982).
  * (30) Rabe K. M., Ahn, C. H., Triscone, J.‐M. Physics of Ferroelectrics: A Modern Perspective. Springer‐Verlag Berlin Heidelberg, Vol. 105 (2007).
  * (31) Yang, Q., Wu, M., Li, J. Origin of Two-Dimensional Vertical Ferroelectricity in WTe2 Bilayer and Multilayer. _J. Phys. Chem. Lett._ 9, 7160–7164 (2018).
  * (32) Liu, X. et al. Vertical ferroelectric switching by in-plane sliding of two-dimensional bilayer WTe2. _Nanoscale_ 11, 18575–18581 (2019).
  * (33) Zhang, T. et al. Ferroelastic-ferroelectric multiferroics in a bilayer lattice. _Phys. Rev. B_ 103, 165420 (2021).
  * (34) Catalan, G., Seidel, J., Ramesh, R., Scott, J. F. Domain wall nanoelectronics. _Rev. Mod. Phys._ 84, 119–156 (2012).
  * (35) Bednyakov, P. S., Sturman, B. I., Sluka, T. et al. Physics and applications of charged domain walls. _npj Comput. Mater._ 4, 65 (2018).
  * (36) Sluka, T., Tagantsev, A. K., Damjanovic, D., Gureev, M., Setter, N. Enhanced electromechanical response of ferroelectrics due to charged domain walls. _Nat. Commun._ 3, 748 (2012).
  * (37) Meyer, B. & Vanderbilt, D. Ab initio study of ferroelectric domain walls in PbTiO3. _Phys. Rev. B_ 65, 104111 (2002).
  * (38) Ren, W. et al. Ferroelectric Domains in Multiferroic ${\mathrm{BiFeO}}_{3}$ Films under Epitaxial Strains. _Phys. Rev. Lett._ 110, 187601 (2013).
  * (39) Yang, Ya. et al. Improper ferroelectricity at antiferromagnetic domain walls of perovskite oxides. _Phys. Rev. B_ 96, 104431 (2017).
  * (40) Sluka, T., Tagantsev, A. K., Bednyakov, P., Setter, N. Free-electron gas at charged domain walls in insulating BaTiO3. _Nat. Commun._ 4, 1808 (2013).
  * (41) Ohtomo, A. & Hwang, H. A high-mobility electron gas at the LaAlO3/SrTiO3 heterointerface. _Nature_ 427, 423–426 (2004).
  * (42) Sturman, B., Podivilov, E., Stepanov, M., Tagantsev, A., Setter, N. Quantum properties of charged ferroelectric domain walls. _Phys. Rev. B_ 92, 214112 (2015).
  * (43) Garcia, V. & Bibes, M. Ferroelectric tunnel junctions for information storage and processing. _Nat. Commun._ 5, 4289 (2014).
  * (44) Tamai, A. et al. Fermi Arcs and Their Topological Character in the Candidate Type-II Weyl Semimetal ${\mathrm{MoTe}}_{2}$. _Phys. Rev. X_ 6, 031021 (2016).
  * (45) Brown, B. E. The crystal structures of WTe2 and high-temperature MoTe2. _Acta Cryst._ 20, 268–274 (1966).
  * (46) Mar, A., Jobic, S., Ibers, J. A. Metal-metal vs tellurium-tellurium bonding in WTe2 and its ternary variants TaIrTe4 and NbIrTe4. _J. Am. Chem. Soc._ 114, 8963–8971 (1992).
  * (47) Aroyo, M. I., Perez-Mato, J. M., Orobengoa, D., Tasci, E., de la Flor, G., Kirov, A. Crystallography online: Bilbao Crystallographic Server. _Bulg. Chem. Commun._ 43, 183–197 (2011)
  * (48) Capillas, C. et al. SYMMODES: a software package for group-theoretical analysis of structural phase transitions. _J. Appl. Cryst._ 36, 953–954 (2003).
  * (49) Momma, K. & Izumi, F. VESTA3 for three‐dimensional visualization of crystal, volumetric and morphology data. _J. Appl. Crystallogr._ 44, 1272–1276 (2011).
  * (50) Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals. _Phys. Rev. B_ 47, 558–561 (1993).
  * (51) Blochl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 50, 17953 (1994).
  * (52) Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _J. Phys.: Condens.Matter_ 21, 395502 (2009).
  * (53) Ceperley, D. M. & Alder, B. J. Ground State of the Electron Gas by a Stochastic Method. _Phys. Rev. Lett._ 45, 566–569 (1980).
  * (54) Perdew J. P., Burke, K., Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 77, 3865–3868 (1996).
  * (55) Csonka, G. I. et al. Assessing the performance of recent density functionals for bulk solids. _Phys. Rev. B_ 79, 155107 (2009).
  * (56) Dal Corso, A. Pseudopotentials periodic table: From H to Pu. _Computational Material Science_ 95, 337–350 (2014).
  * (57) Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations. _Phys. Rev. B_ 13, 5188–5192 (1976).
  * (58) Grimme, S. Semiempirical GGA-type density functional constructed with a long‐range dispersion correction. _J. Comp. Chem._ 27, 1787–1799 (2006).
  * (59) Grimme, S., Antony, J., Ehrlich, S., Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 132, 154104 (2010).
  * (60) Grimme, S., Ehrlich, S., Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. _J. Comp. Chem._ 32, 1456–1465 (2011).
  * (61) Dion, M., Rydberg, H., Schröder, E., Langreth, D. C., Lundqvist, B. I. Van der Waals Density Functional for General Geometries. _Phys. Rev. Lett._ 92, 246401 (2004).
  * (62) Román-Pérez, G. & Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. _Phys. Rev. Lett._ 103, 096102 (2009).
  * (63) Lee, K., Murray, É. D., Kong, L., Lundqvist, B. I., Langreth, D. C. Higher-accuracy van der Waals density functional. _Phys. Rev. B_ 82, 081101(R) (2010).
  * (64) Klimeš, J., Bowler, D. R., Michaelides, A. Chemical accuracy for the van der Waals density functional. _J. Phys.: Condens. Matter_ 22, 022201 (2009).
  * (65) Klimeš, J., Bowler, D. R., Michaelides, A. Van der Waals density functionals applied to solids. _Phys. Rev. B_ 83, 195131 (2011).
  * (66) Liechtenstein, A. I., Anisimov, V. I., Zaanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. _Phys. Rev. B_ 52, R5467–R5470 (1995).
  * (67) Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. _Phys. Rev._ 139, A796–A823 (1965).
  * (68) Shishkin, M. & Kresse, G. Implementation and performance of the frequency-dependent $GW$ method within the PAW framework. _Phys. Rev. B_ 74, 035101 (2006).
  * (69) Botti, S. & Marques, M. A. L. Strong Renormalization of the Electronic Band Gap due to Lattice Polarization in the $GW$ Formalism. _Phys. Rev. Lett._ 110, 226404 (2013).
  * (70) Togo, A. & Tanaka, I. First principles phonon calculations in materials science. _Scr. Mater._ 108, 1–5 (2015).
  * (71) Baroni, S., de Gironcoli, S., Dal Corso, A., Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. _Rev. Mod. Phys._ 73, 515–562 (2001).
  * (72) Gajdoš, M., Hummer, K., Kresse, G., Furthmüller, J., Bechstedt, F. Linear optical properties in the projector-augmented wave methodology. _Phys. Rev. B_ 73, 045112 (2006).
  * (73) King-Smith, R. D. & Vanderbilt, D. Theory of polarization of crystalline solids. _Phys. Rev. B_ 47, 1651–1654 (1993).
  * (74) Resta, R. Macroscopic polarization in crystalline dielectrics: The geometric phase approach. _Rev. Mod. Phys._ 66, 899–915 (1994).
  * (75) Mostofi, A. A. et al. An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions. _Comput. Phys. Commun._ 185, 2309–2310 (2014).
  * (76) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. _J. Chem. Phys._ 81, 511–519 (1984).
  * (77) Debye P. & Hückel E. Zur Theorie Der Elektrolyte. I. Gefrierpunktserniedrigung Und Verwandte Erscheinungen. _Phys. Z._ 24, 185–206 (1923).

##  V ACKNOWLEDGEMENTS

The authors thank Xingen Liu, Minglang Hu, Wei Wu, and Yue-Wen Fang for
stimulating discussions.

This work was supported by the Tokyo Tech World Research Hub Initiative (WRHI)
Program of the Institute of Innovative Research, Tokyo Institute of
Technology, the National Natural Science Foundation of China (Grants No.
51861145315, No. 11929401, No. 12074241, No. 52130204), the Science and
Technology Commission of Shanghai Municipality (Grant No. 19010500500, No.
20501130600, No. 19DZ2270200), the Independent Research and Development
Project of State Key Laboratory of Advanced Special Steel, Shanghai Key
Laboratory of Advanced Ferro Metallurgy, Shanghai University (Grant No. SKLASS
2020-Z07), Austrian Research Promotion Agency (FFG, Grant No. 870024, project
acronym MagnifiSens), and High Performance Computing Center, Shanghai
University.

##  VI AUTHOR CONTRIBUTIONS

W.R. and S.A.N conceived the study. X.M., S.A.N., and C.L. performed first-
principles calculations of electronic structures, lattice dynamics, and
ferroelectric properties. S.A.N. performed the analysis of structural
stability and calculations of the Wannier functions. S.A.N. and X.M.
investigated the properties of domain walls. S.A.N. analyzed the results and
wrote the manuscript with the input from all authors.

##  VII DATA AVAILABILITY

The authors declare that all source data supporting the findings of this study
are available within the article and the Supplementary information file.

##  VIII COMPETING INTERESTS

The authors declare no competing interests.