Ultrahigh electromechanical response from competing ferroic orders

Thin film growth

The conventional solid-state reaction method was used to fabricate the ceramic targets of NaNbO₃ using Na₂CO₃ (99.0%) and Nb₂O₅ (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₃ (111) substrates by sputtering at 790 °C with Ar:O₂ = 2:1, power of 100 W and a deposition time of 4 h. The Nb-doped SrTiO₃ substrate (Shinkosha) was cleaned in acetone before putting it into the sputtering chamber. The NNO film has a stoichiometry of Na_1.16NbO_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 1s, O 1s and Nb 3d were collected. The background type used is Shirley. The spectra were processed using CasaXPS v.2.3.23 software (ref. ⁴⁶).

Rutherford backscattering spectroscopy

RBS experiments were performed using a nuclear microprobe at the Centre for Ion Beam Applications, National University of Singapore⁴⁷. 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. ⁴⁸). 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⁴².

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⁻¹ 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² area with a +10 V tip bias, and then followed by a magnification of scanning (2.5 × 2.5 µm²) with a −10 V tip bias.

First-principles calculations

We studied the structural phase transitions and electronic structures of NaNbO₃ by using the first-principles calculations within GGA-PBESol approximation⁴⁹ as implemented in the WIEN2K software package, the fully localized limit using full-potential (linearized) augmented plane wave (FP-LAPW) + local orbital method⁵⁰. 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_minK_max = 8.0. Fourier expansion of the charge density was performed up to R_minG_max = 20.0. Here, R_min is the minimum sphere radius, that is, 1.61 a.u. The Na 2s²2p⁶3s¹, Nb 4s²4p⁶4d⁴5s¹ and O 2s²2p⁴ 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.⁻¹. 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.⁻¹. The polarization of R3c structure is obtained by the Berry phase simulation, which is implemented in the WIEN2K code⁵¹, P = 6.24 μC cm⁻², at a = 3.889 Å. The schematics of atomic structures demonstrated in this paper were drawn by VESTA software⁵².

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

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⁵³. 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,

The energy contributions are given by

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 ε⁰(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. ⁵³.

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₁ = −6.5 × 10⁷ J m C⁻², a₁₁ = 0.9 × 10⁸ J m⁵ C⁻⁴, a₁₂ = 8.0 × 10⁸ J m⁵ C⁻⁴, a₁₁₁ = 3.3 × 10⁹ J m⁹ C⁻⁶, a₁₁₂ = −3.5 × 10⁹ J m⁹ C⁻⁶, a₁₂₃ = −1.0 × 10⁹ J m⁹ C⁻⁶, a₁₁₁₁ = −3.1 × 10¹⁰ J m¹³ C⁻⁸, a₁₁₁₂ = 0.2 × 10¹⁰ J m¹³ C⁻⁸, a₁₁₂₂ = 4.2 × 10¹⁰ J m¹³ C⁻⁸, a₁₁₂₃ = −5.0 × 10¹⁰ J m¹³ C⁻⁸, b₁ = −6.5 × 10⁵ J m C⁻², b₁₁ = 0.9 × 10⁸ J m⁵ C⁻⁴, b₁₂ = 7.3 × 10⁸ J m⁵ C⁻⁴, b₁₁₁ = 3.3 × 10⁹ J m⁹ C⁻⁶, b₁₁₂ = −3.5 × 10⁹ J m⁹ C⁻⁶, b₁₂₃ = 25.0 × 10⁹ J m⁹ C⁻⁶, b₁₁₁₁ = 3.1 × 10¹⁰ J m¹³ C⁻⁸, b₁₁₁₂ = 0.2 × 10¹⁰ J m¹³ C⁻⁸, b₁₁₂₂ = −4.0 × 10¹⁰ J m¹³ C⁻⁸, b₁₁₂₃ = 15.0 × 10¹⁰ J m¹³ C⁻⁸ and d₁₁ = d₁₂ = 1.0 × 10⁹ J m⁵ C⁻⁴. The gradient energy coefficients are obtained from the domain wall width, taken as g₁₁ = 3.2 × 10⁻¹¹ J m³ C⁻², g₁₂ = 0 and g₄₄ = 1.6 × 10⁻¹¹ J m³ C⁻². The background dielectric constant takes a common value across ferroelectric perovskites⁶², \({{\kappa }}_{11}^{{\rm{b}}}=40\). The elastic stiffness is c₁₁ = 2.30 × 10¹¹ J m⁻³, c₁₂ = 0.90 × 10¹¹ J m⁻³ and c₄₄ = 0.76 × 10¹¹ J m⁻³ (ref. ⁵⁹). 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₁₁ = 0.160 m⁴ C⁻², Q₁₂ = −0.072 m⁴ C⁻² and Q₄₄ = 0.084 m⁴ C⁻².

The phase-field simulations are performed in a two-dimensional simulation system with a total size taken as l₁ × l₃ = 256 nm × 128 nm, which is discretized into an array of 256 × 128 grids. The vertical dimension of the system (with l₃ = 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⁶³. 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.