Introduction

Ferroelectrics have a switchable spontaneous polarization which is used, for example, in memory applications1, and their large piezoelectric and dielectric responses are essential for capacitors, sensors, and actuators2,3. Proper ferroelectrics, whose primary order parameter is the spontaneous polarization, such as Pb(Zr,Ti)O3 and BaTiO3, exhibit a qualitatively well-defined electric-field-driven polarization reversal process4 that is determined by the elastic and electrostatic boundary conditions5,6. At the local scale, the reversal is governed by domains, following material-dependent nucleation and growth processes that involve the movement of ferroelectric ___domain walls7. In general, for small electric fields reversible ___domain wall movements can occur (Rayleigh’s law), whereas larger electric fields lead to nonlinearity, manifesting as Barkhausen jumps8. Characteristic distances reported for reversible ferroelectric ___domain wall movements according to the Rayleigh’s law are usually in the sub-10-nm range9. To increase the regime in which reversible ___domain wall movements occur, different approaches have been tested, including point defect10,11,12 and ___domain engineering13, as well as cycling of ferroelectric ___domain walls between pinning centers14. Following these approaches, electric-field-driven reversible ___domain wall motions over distances up to several micrometers have been achieved. In contrast to proper ferroelectrics, the local dynamics and reversible displacements of ___domain walls in improper ferroelectrics, where the electric polarization arises as a secondary order parameter15,16, are relatively unexplored.

Here, we investigate the ___domain wall dynamics in the improper ferroelectric model system ErMnO3. By applying band-excitation piezoresponse force microscopy (BE-PFM)17, we map the ferroelectric domains and resolve their local response to applied electric fields. We find that away from the topologically protected structural vortex lines, which are characteristic for the ___domain pattern in hexagonal manganites18,19,20, ___domain walls readily move under the electric field. Most intriguingly, the ___domain walls consistently return to their original positions after bipolar electric field cycling, even after being displaced by distances exceeding 250 nm. Complementary phase field simulations exhibit the same switching behavior, indicating that this phenomenon is characteristic for hexagonal manganites and linked to their improper ferroelectric order.

Results

ErMnO3 is a geometrically driven improper ferroelectric, where the primary order parameter relates to a structural trimerization mode. Its fundamental ___domain19,20 and ___domain-wall physics21,22 are well-understood, making it an ideal model system for dynamical investigations at the nanoscale. The system is a uniaxial ferroelectric with a spontaneous polarization, P ≈ 6 µC/cm², parallel to the hexagonal c-axis (space group P63cm). The ___domain structure consists of ferroelectric 180° domains that meet at topologically protected six-fold structural vortex lines, which from at the improper ferroelectric phase transition as described elsewhere23. The family of hexagonal manganites is attracting broad attention due to its unusual physical properties, including multiferroicity24,25, functional ___domain walls22,26, Kibble-Zurek scaling27, and inverted ___domain-size / grain size scaling28,29, with potential applications in barrier layer capacitors30. In this work, we use polycrystalline ErMnO3 from the same batch as studied in ref. 28, where the interested reader can find relevant details concerning the processing conditions of the sample and a crystallographic and microstructural characterization.

We begin by measuring the local electric-field response using BE-PFM17,31 as presented in Fig. 1. In contrast to conventional PFM, which operates at a single frequency, the response is detected in BE-PFM within a specified frequency range (here 360 to 420 kHz), allowing to capture the broad mechanical resonance frequency of the tip-sample system32. Fig. 1a shows a schematic illustration of the measurement principle. To image the polarization-reversal process, a bipolar triangular voltage signal sequence (Fig. 1a) regularly used in SS-PFM33, is applied point by point to the surface via the probe tip. The BE-PFM response is measured after pulse application, as shown in the inset to Fig. 1a (the protocol applied to correct for electrostatic background contributions and electric-field screening effects34,35 is described in Supplementary Figs. S1 and S2). The initial vertical BE-PFM phase response of our polycrystalline ErMnO3 specimen is displayed in Fig. 1b (for corresponding amplitude data see Supplementary Fig. S3). We find that the BE-PFM response is maximized in the vertical channel, which leads us to the conclusion that the grain has primarily an out-of-plane polarization component, i.e., the c-axis is approximately perpendicular to the imaged surface. Ferroelectric \(\pm\)P domains can be distinguished based on their characteristic BE-PFM phase response, with -P domains showing a phase shift of 0 rad and +P domains showing a -π rad phase shift, consistent with their opposite polarization directions. The BE-PFM data shows the well-established ___domain structure of ErMnO318,19,22, consisting of ferroelectric 180° domains that come together by forming a characteristic sixfold line defect. The surface intersection point of the line is marked by the red arrow in Fig. 1b.

Fig. 1: BE-PFM data showing the voltage-dependent ___domain-structure evolution of ErMnO3.
figure 1

a Schematic time-dependence of the applied DC bipolar voltage sequence, corresponding to a square triangular function. The inset schematically visualizes the BE-PFM measurement applied to monitor the piezoresponse around the tip-sample resonance frequency. Details are provided in the methods section. Corrected off-field BE-PFM phase images of the ___domain structure at b maximum negative and c maximum positive voltage. White areas represent domains with upward polarization, +P, whereas black areas indicate downward polarization, −P. The red arrow indicates the position of the intersection of a vortex line with the surface. Green dots indicate the position of the ferroelectric ___domain wall prior to the application of the bipolar volage sequence (at t = 0 ms). A corresponding BE-PFM amplitude image is displayed in Fig. S1. The continuous change of the ferroelectric ___domain structure is displayed in a movie in the Supplementary Information (Supplementary Movie 1). A zoom-in into the region marked by the blue square, featuring its voltage-dependent ___domain-structure evolution, is displayed in d). The time-dependent BE-PFM phase (after the application of the DC pulses) is extracted along the dashed white line and displayed in Fig. 3d. The dimension of the square is 0.66 × 0.66 µm².

Full size image

Representative BE-PFM phase images recorded at maximum negative (−12.0 V) and positive (+12.0 V) voltages are displayed in Fig. 1b and c, respectively. We find pronounced spatial variations in the response of the ferroelectric domains, depending on the local ___domain structure. For example, the ___domain wall movements are largely suppressed in the vicinity of the vortex cores, corroborating that the vortices act as anchor points for the ferroelectric ___domain walls36,37. In contrast, substantial changes in the ___domain wall position are observed away from the vortex lines as presented in Fig. 1d. Figure 1d shows the voltage-dependence of the BE-PFM phase response for the area marked by the blue square in Fig. 1b, c. To demonstrate the effect of the electric field on the ferroelectric domains, the initial ___domain wall position is highlighted by green dots. Under the application of a negative voltage the −P domains contract, followed by an expansion under a positive voltage. Most interestingly for this work, we find that the initial and final ___domain-wall positions in Fig. 1d are almost indistinguishable, indicating that the ___domain walls revert to their starting position after applying the bipolar electric-field cycle. The behavior is consistent over the sample as confirmed by Fig. S4, which displays the analysis at a different position. Our results align with previous scanning electron microscopy studies, which demonstrate that the initial and final state after ion- or electron-beam irradiation are similar38. In contrast to those studies, we apply voltage via an electrical contact, enabling precise monitoring of ___domain wall positions under controlled conditions.

To understand the local polarization reversal, we examine the voltage-dependent BE-PFM phase response on the nanoscale in the region of interest (Fig. 1d). The BE-PFM phase response can be classified as a function of the distance from the ___domain wall. Representative responses are displayed in Fig. 2a–c. At distances exceeding 150 nm from the original ___domain wall position, the BE-PFM phase response is independent of the applied voltage pulse (Fig. 2a), showing that the local polarization does not switch up to the maximum applied voltage of \(\pm\)12.0 V. As we approach the ___domain wall, highly asymmetric hysteretic BE-PFM phase responses emerge as a function of the voltage, with the effect being position-dependent, varying between 40 nm and 150 nm from the initial ___domain wall position. The sign of the asymmetry in the hysteresis depends on the side from which the ___domain wall is approached (Fig. 2b). For both sides, a change of the BE-PFM phase by \({{\mathrm{\pi}}}\) rad under the maximum applied voltage can be measured, which is an indication of local polarization reversal. Depending on the side of the ___domain wall, the asymmetry in the hysteresis loop is expressed by the absence of a negative or positive coercive voltage. Mechanistically, the BE-PFM data indicates that polarization reversal occurs by attraction or repulsion of the ferroelectric ___domain wall towards or away from the AFM tip39,40, implying that the ___domain wall reverts back to its initial position after the bipolar electric-field cycle. Finally, at a position close to the ferroelectric ___domain wall (Fig. 2c), a hysteresis loop develops as expected for electric-field-driven polarization reversal. The behavior can be quantified by measuring the width of the hysteresis loop, \({V}_{c}^{+}-{V}_{c}^{-}\), and its shift, \(({V}_{c}^{+}-{V}_{c}^{-})/2\). Voltages \({V}_{c}^{+}\) and \({V}_{c}^{-}\) are displayed as dashed lines in Fig. 2b and c. As depicted in Fig. 2d, the width of the hysteresis loop reaches up to ~12 V at the ___domain wall and gradually decreases away from it. The hysteresis shift reaches maximum values of −6 V and + 6 V away from the ___domain wall and gradually decreases towards the position of the wall. As expected, the hysteresis shift is positive when approaching the ___domain wall from the +P ___domain and has a negative value when approaching the ___domain from the P ___domain. Figure 2e displays the average hysteresis width alongside the hysteresis shift plotted as a function of the distance from the ___domain wall position (the position is indicated by a dashed straight line). The shift of the hysteresis loop is symmetric with respect to the initial ___domain wall position. The asymmetry in the hysteresis width can be explained based on surface charges, which can lead to a substantial bias in ferroelectric switching41. The asymmetric switching behavior is consistent with previous spatially resolved measurements and plays an important role at the local scale38. The observed switching behavior (Fig. 2) is interesting as it shows that distinct switching behaviors can be achieved within the same electric field range at the nanoscale, depending on the relative positioning of the electrode with respect to the ___domain walls and vortex lines.

Fig. 2: Nanometric evaluation of the switching behavior in ErMnO3.
figure 2

The voltage-dependent BE-PFM phase response illustrates distinct switching types relative to the distance from the ferroelectric ___domain wall: a the BE-PFM phase response is constant for all DC offset voltages and no switching is observed at distances larger than 150 nm from the ___domain wall. The –P domains exhibit a voltage-independent phase of 0, whereas the +P domains exhibit a phase of \(-{{\mathrm{\pi}}}\). b Asymmetric hysteresis loop measured at two different sites of the ___domain wall at a distance of ~30 nm. The light green loop represents the response measured in a +P ___domain, and the dark green loop corresponds to the response measured in –P ___domain. c Symmetric hysteresis observed at the position of the ferroelectric ___domain wall. Dashed lines in b and c represent the positive and negative coercive voltage Vc+ and Vc−, respectively. d Illustration of  the width and horizontal shift of the hysteresis loop, showcasing the validity of our classification (a–c) for the area of interest (Fig. 1d). Green dots indicate the position of the ferroelectric ___domain wall prior to the application of the bipolar voltage sequence (t = 0 ms). e Hysteresis shift and width, assessed in the vicinity of the ___domain wall in the section shown in Fig. 2d, depicted as a function of the distance from the ___domain wall position (displayed as a black dashed line) averaged between t = 0 ms and t = 250 ms. The error bars indicate the standard deviation of all pixels grouped into the corresponding data point.

Full size image

To evaluate whether the response to the electric field is specific to the ErMnO3 polycrystal or representative for the family of hexagonal manganites in general, we compare the experimental results with phase-field simulations42. Following the established phase-field model for hexagonal manganites, the system is described by a Landau expansion of the trimerization tilt amplitude \(Q\), the azimuthal tilt phase \(\varPhi\), and the polar mode \({{{\mathcal{P}}}}\)20,43,44, as elaborated in the method section. The model reproduces the characteristic ___domain structure of ErMnO3 observed in our experiment (Fig. 1 and S4), including the topologically protected vortex lines, as displayed in Fig. 3a. Importantly, the data allows for investigating the impact of an electric field on different surfaces and within the bulk material, giving a three-dimensional (3D) model of the electric-field-induced ___domain structure. For this purpose, we consider the coupling between the energy density and the electric field as an additional term in the free energy, \(F={F}_{{{{\rm{Landau}}}}}-{PE}\), where \(P\) is the polarization and \(E\) is the external electric field, as explained in more detail in the method section.

Fig. 3: Phase-field simulations showing the relation between the electric field and the ___domain structure.
figure 3

The influence of an electric field on the ferroelectric domains is displayed for a the initial configuration without an applied electric field and b at the maximum negative electric field. The continuous change of the simulated ferroelectric ___domain structure is displayed in a movie in the Supplementary Material (Supplementary Movie 2). c The polarization is extracted along the dotted line in panels a and b as a function of an equivalent triangular-shaped voltage sequence, analogous to our experiment (Fig. 1a). d Voltage-dependent data of the BE-PFM phase, obtained from experiment, extracted along a dotted line of the area of interest in Fig. 1d is displayed. The green dots represent a guide to the eye for the electric-field-dependent positions of the ___domain wall, determined through thresholding. Jump-like changes in the position are highlighted by white arrows.

Full size image

The simulated ___domain structure for a triangular electric field sequence analogous to the experiment (Fig. 1a) is displayed in Fig. 3a and b, visualizing the electric-field driven evolution of the ferroelectric ___domain structure from the initial to the maximum electric field, respectively. The simulation reproduces the experimentally observed contraction and expansion of the ferroelectric domains in response to the applied electric field. Close to the vortex lines, variations in the ___domain structure are much less pronounced than elsewhere, corroborating that these lines act as anchor points for the ___domain walls (Fig. S6a)36,37. Interestingly, in agreement with our experimental observations, the simulation shows that the ferroelectric ___domain walls return to their initial position after the bipolar electric-field cycle. Given its 3D nature, the phase-field simulations extend the reversible long-range motion observed in the BE-PFM data to the bulk of the material. The importance of the improper ferroelectric nature and the topologically protected vortex lines for the reversibility of the ___domain structure is visualized in Fig. S6b and c. The observed behavior is fundamentally different from proper ferroelectric materials, which often show irreversible ___domain-structure changes during electric-field cycling, where the ___domain structure is either transformed into a mono-___domain state45 or spatially unrelated configurations46. The reversibility is displayed by the electric-field-dependent evolution of two representative walls in Fig. 3c. Green dots are used to indicate the positions of the ferroelectric ___domain walls, visualizing how they evolve as a function of the electric field. For comparison, experimental data extracted from Fig. 1d are presented in Fig. 3d. Different from the simulation, jump-like changes in position, exemplarily highlighted by white arrows, arise in the experiment, which we attribute to imperfections in the crystallographic structure of the material. For example, it has been calculated that oxygen interstitials have a lower formation energy at the ___domain walls and, hence, promote pinning26,47. Despite such jump-like changes, experiment and theory qualitatively exhibit the same behavior. The ___domain walls relax back to their initial positions after the bipolar electric-field cycle. This behavior suggests that the reversibility in the ___domain switching is linked to the improper ferroelectric order of ErMnO3, with the relaxation of the primary structural order parameter acting as an additional resorting force, overcoming the pinning potential of defects. The phase-field simulations incorporate restoring forces through the stiffness term, which is directly associated with the primary structural order parameter (Eq. 1). Remarkably, the ___domain walls revert to their original position after having traversed distances of up to 284 ± 25 nm (see Fig. 1d), moving between displacements under maximum positive (+12 V) and negative voltages (\(-\)12 V). Such reversible long-range ___domain wall motions are consistently observed at various positions across the sample. Another example is presented in Fig. S4 (Supplementary Movie 3), showing the case of neighboring stripe-like domains, where the ___domain walls are displaced by about 150 ± 25 nm under application of the electric field as evaluated in Fig. S5. In general, even larger traveling distances may be expected if higher voltages are applied with the upper limit being constrained by the proximity of neighboring ___domain walls. In addition, by optimizing the synthesis and reducing the amount of lattice imperfections, the emergence of defect-related jump-like ___domain wall motions may be suppressed, facilitating deterministic control with even larger traveling distances. Most importantly, ErMnO3 shows this reversibility in its pristine state, that is without the need for engineering by, e.g., defect dipoles10 or point defects12,14. To gain additional insight into the microscopic interactions that govern the ___domain wall motion in ErMnO3, extended atomic-level calculations are required, which could be performed as a next step. Most importantly for the present study, the agreement between the BE-PFM measurements and the phase-field simulations reflects that the observed behavior is common to improper ferroelectric hexagonal manganites and not specific to the ErMnO3 polycrystal under investigation. Furthermore, this finding highlights that chemical defects play a secondary role in this system.

In conclusion, we have visualized the electric-field-driven local polarization reversal in improper ferroelectric ErMnO3 using BE-PFM with nanometric precision. We attribute this phenomenon to the improper nature of the ferroelectric order, with the structural order parameter acting as a restoring force that promotes a relaxation of the ___domain walls to their initial position. Our observation of reversible long-range ___domain wall movements in ErMnO3 is consistent with reports on the electric-field-induced contraction of domains into meandering lines19,21 and extends electric-field-dependent transmission electron microscopy studies36 to the level of ferroelectric domains. We expect that the observed reversibility of the switching process is not limited to ErMnO3 and the family of hexagonal manganites. Other possible candidate materials include isostructural ferrites48, indates49, and gallates50, as well as hexagonal tungsten bronzes51 or 2H compounds52. To further increase the traveling distance of the ___domain walls or reduce the amplitude of the driving electric field, parameters such as temperature and mechanical pressure may be leveraged to reduce the coercive field53, representing a pathway for optimizing the dynamical responses towards potential device applications. The demonstrated controllable and reversible motion of ___domain walls holds promise for advanced devices such as sensors54 and tunable capacitors55. For sensors, the demonstrated reversibility enables dynamic adjustment of the figures of merit, whereas in capacitors, it facilitates precise fine-tuning of the dielectric constant. In both cases, the electric field stimulates a reversible dynamic shift in ___domain wall position. In addition, the dynamic responses may be utilized to translate electric input signals into complex ___domain wall displacements with well-defined nonlinear relaxation behavior, providing all key characteristics required for reservoir computing and giving new opportunities for in-materio computing56.

Methods

Synthesis and sample preparation

Synthesis of ErMnO3 powder was done by a solid-state reaction of dried Er2O3 (99.9% purity; Alfa Aesar, Haverhill, MA, USA) and Mn2O3 (99.0% purity; Sigma-Aldrich, St. Louis, MO, USA) powders. The powders were mixed and ball-milled (BML 5 witeg Labortechnik GmbH, Wertheim, Germany) for 12 hrs at 205 rpm using yttria-stabilized zirconia milling balls of 5 mm diameter and ethanol as dispersion medium. The reaction to ErMnO3 was done by stepwise heating to 1000 °C, 1050 °C, and 1100 °C with a dwell time of 12 hrs. More details on the powder processing can be found in ref. 28. The powder was isostatically pressed into samples of cylindrical shape at a pressure of 200 MPa (Autoclave Engineers, Parker-Hannifin, Cleveland, OH, USA). Sintering was carried out in a closed alumina crucible at a temperature of 1350 °C for 4 hrs with a heating and cooling rate of 5 K/min (Entech Energiteknik AB; Ängelholm, Sweden).

PFM measurements

Prior to PFM measurements, the samples were lapped with a 9-µm-grained Al2O3 water suspension (Logitech Ltd., Glasgow, UK) and polished using silica slurry (SF1 Polishing Fluid, Logitech AS, Glasgow, Scotland). BE-PFM measurements were conducted with the same Cr/Pt coated Multi75E-G AFM probe (force constant, k≈3 N/m, Budget Sensors, Sofia, Bulgaria), using a Cypher atomic force microscope (Oxford Instrument, Abingdon, UK). For BE-PFM imaging the electrical bias was applied to the surface via the probe tip, while the sample was grounded. More information on the measurement technique are provided in ref. 32. The tip nominal radius was ~25 nm and the pixel size was ~47 nm. For the analysis of the data in Figs. 2 and 3, we considered that the regular grid structure of the scan is non-collinear with the features of interest under investigation. BE-PFM was achieved by arbitrary-wave generator and data-acquisition electronics. A custom software program was employed to generate the probing signal and record local amplitude and phase hysteresis loops. BE-PFM measurements were performed at a frequency range of 360 to 420 kHz, a power spectral density of 0.57 V2/kHz, and an amplitude of 5 V, respectively. Simple harmonic oscillator fits were subsequently applied to the measured spectra to extract amplitude and phase information at the resonance frequency. The SS-PFM triangular-square function had a maximum of ±12 V, applied over 600 ms with a sequence of 64 steps of on-field and off-field measurements.

Phase-field simulations

Phase-field simulations were performed based on the Landau expansion of the free energy of ferroelectric hexagonal manganites as described in ref. 20 as

$$\tt {F}_{{Landau}}\left(Q,\varPhi {{,}}{{{\mathcal{P}}}}\right)=\frac{a}{2}{Q}^{2}+\frac{b}{4}{Q}^{4}+\frac{{Q}^{6}}{6}(c+{c}^{{\prime} }\cos 6\varPhi )-g{Q}^{3}{{{\mathcal{P}}}}\cos 3\varPhi -\frac{{g}^{{\prime} }}{2}{Q}^{2}{{{{\mathcal{P}}}}}^{2} \\ +\frac{{a}_{{{{\mathcal{P}}}}}}{2}{{{{\mathcal{P}}}}}^{2}+\frac{1}{2}{\sum}_{i=x,y,z}[{s}_{Q}^{i}({\partial }_{i}Q{\partial }_{i}Q+{Q}^{2}{\partial }_{i}\varPhi {\partial }_{i}\varPhi )+{s}_{{{{\mathcal{P}}}}}^{i}{\partial }_{i}{{{\mathcal{P}}}}{\partial }_{i}{{{\mathcal{P}}}}]{{,}}$$
(1)

where \(Q\) is the amplitude and \(\phi\) is the angle of the lattice-trimerizing bipyramidal MnO5 tilt, and \({{{\mathcal{P}}}}\) is a displacement field that corresponds to the polar mode \({\Gamma }_{2}^{-}\). This displacement field \({{{\mathcal{P}}}}\) is proportional to the polarization, P, of the system. To describe the coupling between the polarization, P, and the external electric field, \(E\), an additional term \(-{PE}\) is included. The term \({s}_{Q}^{i}({\partial }_{i}Q{\partial }_{i}Q+{Q}^{2}{\partial }_{i}\varPhi {\partial }_{i}\varPhi )\) constitutes the stiffness term representing the structural order parameter and accounting for inhomogeneities in the structural order. The structural order is coupled to the ferroelectric order through the terms (\(-g{Q}^{3}{{{\mathcal{P}}}}\cos 3\varPhi -\frac{{g}^{{\prime} }}{2}{Q}^{2}{{{{\mathcal{P}}}}}^{2}\)). The parameters for the Landau expansion were chosen as in ref. 20, as \(a=-2.626\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-2}\), \(b=3.375\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-4}\), \(c=0.117\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-6}\), \({c}^{{\prime} }=0.108\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-6}\), \({a}_{P}=0.866\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-2}\), \(g=1.945\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-4}\), \({g}^{{\prime} }=9.931\,{{{\rm{eV}}}}{\mathring{\rm A} }^{-4}\), \({s}_{Q}^{z}=15.40\,{{{\rm{eV}}}}\), \({s}_{Q}^{x}=5.14\,{{{\rm{eV}}}}\), and \({s}_{p}^{z}=52.70\,{{{\rm{eV}}}}\). To ensure stability, the gradient energy coefficient was set to \({s}_{p}^{x}=+ 8.88\,{{{\rm{eV}}}}\) as validated in ref. 44. Our system was simulated with a uniform Cartesian computational mesh with spacing \({d}_{x}={d}_{y}=0.2\,{{{\rm{nm}}}}\) and \({d}_{z}=0.3\,{{{\rm{nm}}}}\). Our simulations were performed on a mesh with size \({n}_{x}={n}_{y}={n}_{z}=64\). Periodic boundary conditions were chosen in all three principal directions to simulate a bulk system. The Ginzburg-Landau equations were integrated with a Runge-Kutta 4 integrator with time steps of \(\Delta t=5\cdot {10}^{-4}\). The system was initialized with random values of all three order parameters. The initial ___domain pattern was then generated by evolving the system for \(8\cdot {10}^{4}\) time steps, with no external electric field applied, i.e., \(E=0\). From this initial ___domain pattern, the system was integrated with an external field applied parallel to the polarization direction, i.e.,

$$E\left(n\right)=\left\{\begin{array}{cc}\frac{4n}{N}\cdot {E}_{\max } & \frac{n}{N}\le 1/4\\ \left(1-4\left(\frac{n}{N}-\frac{1}{4}\right)\right)\cdot {E}_{\max } & \frac{1}{4} < \frac{n}{N}\le \frac{3}{4}\\ \left(-1+4\left(\frac{n}{N}-\frac{3}{4}\right)\right)\cdot {E}_{\max } & \frac{3}{4} < \frac{n}{N}\end{array}\right.$$
(2)

which is similar to the experimentally applied bipolar voltage sequence displayed in Fig. 1a, yet with a flipped sign of \(E(n)\). Here, \(N=2\cdot {10}^{4}\) is the number of the time steps of a full cycle and \(n\) is a given time step of the cycle. The maximal electric field has been chosen to \({E}_{\max }=0.11\,{{{\rm{V}}}}/{{\text{\AA}}}\).