### Thin film growth

The conventional solid-state reaction method was used to fabricate the ceramic targets of NaNbO_{3} using Na_{2}CO_{3} (99.0%) and Nb_{2}O_{5} (99.9%) as the raw materials. The 3-inch target with 10 mol% excess of Na was used in this study. Epitaxial NNO films with 200 nm thickness (Fig. 2e, top) were grown on Nb-doped SrTiO_{3} (111) substrates by sputtering at 790 °C with Ar:O_{2} = 2:1, power of 100 W and a deposition time of 4 h. The Nb-doped SrTiO_{3} substrate (Shinkosha) was cleaned in acetone before putting it into the sputtering chamber. The NNO film has a stoichiometry of Na_{1.16}NbO_{3.08} derived from the X-ray photoemission spectroscopy (XPS) and Rutherford backscattering spectroscopy (RBS) measurements (Extended Data Fig. 3).

### Characterization

#### X-ray photoelectron spectroscopy

Surface analysis was performed using XPS Kratos AXIS Supra monochromatic Al-Kα (1486.6 eV) X-ray beam. Charge neutralization was used to correct the charge shift by irradiating low-energy electrons and ion beams onto the sample. XPS narrow scan spectra for Na 1*s*, O 1*s* and Nb 3*d* were collected. The background type used is Shirley. The spectra were processed using CasaXPS v.2.3.23 software (ref. ^{46}).

#### Rutherford backscattering spectroscopy

RBS experiments were performed using a nuclear microprobe at the Centre for Ion Beam Applications, National University of Singapore^{47}. A 1.5 MeV He^{+} beam was focused to micrometre dimensions by a quadrupole lens system. The count frequency on a silicon surface barrier detector was about 1,500 Hz by controlling the incident beam current. Single measurement continuously proceeds for about 30 min to collect enough signals. The experimental RBS spectrum was fitted by a program named SIMNRA (v.7.02) (ref. ^{48}). The experimental settings were validated using a standard WSi sample.

#### Synchrotron X-ray diffraction

The synchrotron-based X-ray diffraction experiments were conducted at room temperature in the X-ray development and demonstration beamline of the Singapore Synchrotron Light Source and the surface diffraction beamline (BL02U2) at Shanghai Synchrotron Radiation Facility. The lattice parameters of the orthorhombic antiferroelectric phase and rhombohedral ferroelectric phase in the NNO film were calculated based on the reciprocal space vector method^{42}.

#### Electron microscopy characterizations

A dual beam focused ion beam (FIB) (FEI DB Helios 450 S Nanolab) was used to prepare an ultrathin sample (about 50 nm) for TEM and scanning transmission electron microscopy (STEM) measurements. The specimen cutting was conducted using 30 kV Ga ions, and the specimen cleaning was conducted using 1 kV Ga ions. Cross-sectional and atomic structural images were taken using TEM with a 200-keV field-emission gun and a probe current of 0.5 nA at 1 nm. The instrument used was an FEI Tecnai G2 TF20 TEM. Atomic resolution images were taken using a 200-kV cold field-emission probe-corrected transmission electron microscope (JEM-ARM200F) with a STEM resolution of 0.078 nm. Images were monitored and taken using a Gatan 4,000 × 4,000 one-view CMOS camera. DigitalMicrograph software was used to analyse the TEM and STEM data.

#### Electrical properties

A manual probe station is used for electrical tests. Top electrodes Au/Pd with diameters of about 200 μm were deposited by an e-beam evaporator. Dynamic ferroelectric polarization-electric field loops were measured at an interval of 0.3 ms with top–top geometry using a ferroelectric tester (Precision Multiferroic II, Radiant Technologies). A laser Doppler vibrometer was used to obtain the effective piezoelectric coefficient by measuring the surface displacement on the electrode of the tested thin film under an applied electric field. The electric field consists of a small a.c. voltage with an amplitude of 1 V (*V*_{a.c.} = sin *ωt*) and a d.c. bias used for inducing the phase transition. Only the first harmonic strain was measured during the test.

#### PFM measurements

PFM measurements were conducted on a commercial scanning probe microscope system (MFP-3D, Asylum Research) at room temperature under ambient conditions. A Pt-coated conductive probe (240AC-PP, OPUS) with a force constant of around 2 N m^{−1} and free resonance frequency of about 70 kHz was used. During PFM scanning, the piezoresponse signal was amplified by using contact resonance, and the Dual AC Resonance Tracking mode was enabled to keep the resonance tracked. An a.c. drive of 2–8 *V*_{pp} was typically used for PFM mapping. The FE domain writing was performed by firstly scanning a 5 × 5 µm^{2} area with a +10 V tip bias, and then followed by a magnification of scanning (2.5 × 2.5 µm^{2}) with a −10 V tip bias.

### First-principles calculations

We studied the structural phase transitions and electronic structures of NaNbO_{3} by using the first-principles calculations within GGA-PBESol approximation^{49} as implemented in the WIEN2K software package, the fully localized limit using full-potential (linearized) augmented plane wave (FP-LAPW) + local orbital method^{50}. The simulations have been calculated with the augmented plane wave + local orbital (APW + lo) basis set and the muffin-tin radii (*R*_{MT}) for the atomic spheres were chosen as 2.12 a.u., 1.78 a.u. and 1.61 a.u. for Na, Nb and O, respectively. The partial waves were expanded up to *l*_{max} = 10 in the muffin-tin spheres and outside a plane-wave cutoff *K*_{max} defined by *R*_{min}*K*_{max} = 8.0. Fourier expansion of the charge density was performed up to *R*_{min}*G*_{max} = 20.0. Here, *R*_{min} is the minimum sphere radius, that is, 1.61 a.u. The Na 2*s*^{2}2*p*^{6}3*s*^{1}, Nb 4*s*^{2}4*p*^{6}4*d*^{4}5*s*^{1} and O 2*s*^{2}2*p*^{4} were treated as semi-core states and were included in the valence. The Brillouin zone sampling was done using a Kmesh of 10 × 10 × 10, and 10 × 10 × 3 *k*-points in the full Brillouin zone for *R3c* and *Pbcm* structures, respectively. All structures are fully relaxed by using energy minimization with Hellmann–Feynman force of 1 mRy a.u.^{−1}. For each strain state, that is, different in-plane lattice constant *a*, lattice constant *c* and all atomic positions were fully relaxed until the final Hellmann–Feynman forces were less than 1 mRy a.u.^{−1}. The polarization of *R3c* structure is obtained by the Berry phase simulation, which is implemented in the WIEN2K code^{51}, *P* = 6.24 μC cm^{−2}, at *a* = 3.889 Å. The schematics of atomic structures demonstrated in this paper were drawn by VESTA software^{52}.

### Phase-field simulations

In the phase-field model, the microstructure of the NNO film is described by the ferroelectric polarization field **P**(**x**) and the antiferroelectric order parameter **A**(**x**), where **x** is the spatial position vector. The stable state of the system under given conditions as well as the response of the system to given external fields is simulated by evolving the fields, **P** and **A**, to equilibria as governed by the time-dependent Ginzburg–Landau equation^{53,54}, written as

$${{\boldsymbol{\gamma }}}_{{\rm{P}}}\frac{\partial {\bf{P}}}{\partial t}=-\frac{\partial F}{\partial {\bf{P}}}+{{\bf{E}}}^{{\rm{P}},{\rm{therm}}}.$$

(1)

$${{\boldsymbol{\gamma }}}_{{\rm{A}}}\frac{\partial {\bf{A}}}{\partial t}=-\frac{\partial F}{\partial {\bf{A}}}+{{\bf{E}}}^{{\rm{A}},{\rm{therm}}}.$$

(2)

Equations (1) and (2) describe a relaxation process of the system driven by a minimization of the free energy of the system under given conditions. Here **γ**_{P} and **γ**_{A} are the damping coefficients for ferroelectric polarization and antiferroelectric order, respectively, *t* is the time, and *F* is the free energy of the system. **E**^{P,therm} and **E**^{A,therm} are random fields arising from thermal fluctuations, which obey a normal distribution with the strength given by the fluctuation–dissipation theorem^{54,55}.

The free energy *F* is formulated as a functional of the ferroelectric polarization field, antiferroelectric order parameter and external conditions such as the external electric field **E**^{ext}, that is, *F* = *F*[**P**(**x**), **A**(**x**), **E**^{ext}], following existing phase-field theories^{53}. It takes a sum of the Landau free energy for ferroelectric polarization, the Landau free energy for antiferroelectric order, the coupling energy between ferroelectric polarization and antiferroelectric order, the gradient energy, the electrostatic energy and the elastic energy, that is,

$$F={F}_{{\rm{Landau}}}^{{\rm{P}}}+{F}_{{\rm{Landau}}}^{{\rm{A}}}+{F}_{{\rm{coupling}}}+{F}_{{\rm{gradient}}}+{F}_{{\rm{electrostatic}}}+{F}_{{\rm{elastic}}}.$$

(3)

The energy contributions are given by

$${F}_{{\rm{Landau}}}^{{\rm{P}}}=\int \left({a}_{i}{{P}_{i}}^{2}+{a}_{ij}{{P}_{i}}^{2}{{P}_{j}}^{2}+{a}_{ijk}{{P}_{i}}^{2}{{P}_{j}}^{2}{{P}_{k}}^{2}+{a}_{ijkl}{{P}_{i}}^{2}{{P}_{j}}^{2}{{P}_{k}}^{2}{{P}_{l}}^{2}\right){\rm{d}}{x}^{3},$$

(4)

$${F}_{{\rm{Landau}}}^{{\rm{A}}}=\int \left({b}_{i}{{A}_{i}}^{2}+{b}_{ij}{{A}_{i}}^{2}{{A}_{j}}^{2}+{b}_{ijk}{{A}_{i}}^{2}{{A}_{j}}^{2}{{A}_{k}}^{2}+{b}_{ijkl}{{A}_{i}}^{2}{{A}_{j}}^{2}{{A}_{k}}^{2}{{A}_{l}}^{2}\right){\rm{d}}{x}^{3},$$

(5)

$${F}_{{\rm{coupling}}}=\int {d}_{ij}{{P}_{i}}^{2}{{A}_{j}}^{2}{\rm{d}}{x}^{3},$$

(6)

$${F}_{{\rm{gradient}}}=\int {g}_{ijkl}^{{\rm{P}}}\frac{\partial {P}_{i}}{\partial {x}_{j}}\frac{\partial {P}_{k}}{\partial {x}_{l}}\,{\rm{d}}{x}^{3}+\int {g}_{ijkl}^{{\rm{A}}}\frac{\partial {A}_{i}}{\partial {x}_{j}}\frac{\partial {A}_{k}}{\partial {x}_{l}}\,{\rm{d}}{x}^{3},$$

(7)

$${F}_{{\rm{electrostatic}}}=\int \left(-\frac{1}{2}{E}_{i}^{{\rm{d}}}{P}_{i}-{E}_{i}^{{\rm{ext}}}{P}_{i}\right){\rm{d}}{x}^{3},$$

(8)

$${F}_{{\rm{elastic}}}=\int \frac{1}{2}{c}_{ijkl}\left({\varepsilon }_{ij}-{\varepsilon }_{ij}^{0}\right)\left({\varepsilon }_{kl}-{\varepsilon }_{kl}^{0}\right){\rm{d}}{x}^{3},$$

(9)

Indices *i*, *j*, *k*, l = 1, 2, 3, 4 indicate components of a vector or a tensor in a Cartesian coordinate. An Einstein summation convention over repeated indices *i*, *j*, *k* and *l* is used. *a*_{i}, *a*_{ij}, *a*_{ijk} and *a*_{ijkl} are the Landau coefficients for ferroelectric polarization. *b*_{i}, *b*_{ij}, *b*_{ijk} and *b*_{ijkl} are the Landau coefficients for antiferroelectric order, **d** is the ferroelectric–antiferroelectric coupling coefficient and **g**^{P} and **g**^{A} are the gradient energy coefficients for ferroelectric polarization and antiferroelectric order, respectively. The Landau free energy function is constructed with local or global minima corresponding to the AFE orthorhombic *Pbcm* phase (that is, **P** = **0** and **A** along ⟨110⟩) and the FE rhombohedral *R3c* phase (that is, **P** along ⟨111⟩ and **A** = **0**) so that only these two phases can be stabilized, thus avoiding the formation of other phases. **E**^{d}(**x**) is the depolarization field obtained by solving the electrostatic equilibrium equation at every evolution time step. **c** is the elastic stiffness tensor, **ε**(**x**) is the strain field obtained by solving the mechanical equilibrium equation, and **ε**^{0}(**x**) is the eigenstrain field given by \({\varepsilon }_{ij}^{0}={Q}_{ijkl}{P}_{k}{P}_{l}\), with ** Q** being the electrostrictive coefficient. For details of the mechanical and electrostatic equilibrium equations, please refer to ref.

^{53}.

The material constants for NNO used in the phase-field simulations are obtained from refs. ^{56,57,58,59,60,61,62}. The Landau coefficients for ferroelectric polarization and antiferromagnetic order as well as the ferroelectric–antiferroelectric coupling coefficients are constructed to reproduce the experimental and DFT data of spontaneous polarization, P–E loop, dielectric constant, lattice constant and free energy of the FE and AFE phases^{56,57,58,60,61}, taken as *a*_{1} = −6.5 × 10^{7} J m C^{−2}, *a*_{11} = 0.9 × 10^{8} J m^{5} C^{−4}, *a*_{12} = 8.0 × 10^{8} J m^{5} C^{−4}, *a*_{111} = 3.3 × 10^{9} J m^{9} C^{−6}, *a*_{112} = −3.5 × 10^{9} J m^{9} C^{−6}, *a*_{123} = −1.0 × 10^{9} J m^{9} C^{−6}, *a*_{1111} = −3.1 × 10^{10} J m^{13} C^{−8}, *a*_{1112} = 0.2 × 10^{10} J m^{13} C^{−8}, *a*_{1122} = 4.2 × 10^{10} J m^{13} C^{−8}, *a*_{1123} = −5.0 × 10^{10} J m^{13} C^{−8}, *b*_{1} = −6.5 × 10^{5} J m C^{−2}, *b*_{11} = 0.9 × 10^{8} J m^{5} C^{−4}, *b*_{12} = 7.3 × 10^{8} J m^{5} C^{−4}, *b*_{111} = 3.3 × 10^{9} J m^{9} C^{−6}, *b*_{112} = −3.5 × 10^{9} J m^{9} C^{−6}, *b*_{123} = 25.0 × 10^{9} J m^{9} C^{−6}, *b*_{1111} = 3.1 × 10^{10} J m^{13} C^{−8}, *b*_{1112} = 0.2 × 10^{10} J m^{13} C^{−8}, *b*_{1122} = −4.0 × 10^{10} J m^{13} C^{−8}, *b*_{1123} = 15.0 × 10^{10} J m^{13} C^{−8} and *d*_{11} = *d*_{12} = 1.0 × 10^{9} J m^{5} C^{−4}. The gradient energy coefficients are obtained from the domain wall width, taken as *g*_{11} = 3.2 × 10^{−11} J m^{3} C^{−2}, *g*_{12} = 0 and *g*_{44} = 1.6 × 10^{−11} J m^{3} C^{−2}. The background dielectric constant takes a common value across ferroelectric perovskites^{62}, \({{\kappa }}_{11}^{{\rm{b}}}=40\). The elastic stiffness is *c*_{11} = 2.30 × 10^{11} J m^{−3}, *c*_{12} = 0.90 × 10^{11} J m^{−3} and *c*_{44} = 0.76 × 10^{11} J m^{−3} (ref. ^{59}). The electrostrictive coefficients are taken to be isotropic for simplicity with values consistent with our DFT data of the lattice constant differences between *R3c* FE and *Pbcm* AFE phases (Fig. 1c), with *Q*_{11} = 0.160 m^{4} C^{−2}, *Q*_{12} = −0.072 m^{4} C^{−2} and *Q*_{44} = 0.084 m^{4} C^{−2}.

The phase-field simulations are performed in a two-dimensional simulation system with a total size taken as *l*_{1} × *l*_{3} = 256 nm × 128 nm, which is discretized into an array of 256 × 128 grids. The vertical dimension of the system (with *l*_{3} = 128 nm) along the thickness direction of the film consists of a bottom substrate region of 24 nm, a middle film region of 100 nm and a top vacuum region of 4 nm. An electric potential boundary condition is used at the top and bottom surfaces of the film for the electrostatic equilibrium equation^{63}. The piezoelectric performance of the film is modelled by simulating the strain response of the system on applying a given voltage to the top surface of the film while keeping the bottom surface of the film electrically grounded.