Emergence and transformation of polar skyrmion lattices via flexoelectricity

Topological matter beyond the Landau paradigm hosts rich phases in real and momentum spaces, with magnetic skyrmions serving as a prominent example due to their novel physics and device potential. While magnetic skyrmion lattices (SkLs) are well established, polar skyrmions in ferroelectrics, though observed, have not yet been experimentally demonstrated to form long-range ordered SkLs despite theoretical predictions of hexagonal polar SkLs. Building long-range ordered polar SkLs is crucial for collective properties and applications (e.g., metamaterials, optoelectronics). The flexoelectric effect, which couples polarization and strain gradients, is significant in thin films and can induce orientation anisotropy in domain walls. This suggests flexoelectricity may stabilize and transform different SkL types and enable deterministic control. The study addresses whether SkLs beyond hexagonal can be stabilized, how flexoelectricity governs their stabilization and transformation, and whether one can deterministically control SkL orientation and lattice type via flexoelectricity and related anisotropies.

Magnetic SkLs exhibit distinct responses and emergent phenomena (e.g., topological magnon bands, Nernst effect) and can undergo rotation, distortion, and structural transformations between hexagonal and tetragonal lattices. Polar skyrmions in ferroelectric systems have been reported with properties such as negative capacitance and strong nonlinear optical responses, and theoretical works predict hexagonal polar SkLs and reversible transitions among vortex and skyrmion crystals. Collective dynamics of polar vortex lattices indicate potential device applications. Flexoelectricity, often negligible in bulk, becomes significant in thin films due to large strain gradients and is known to impose orientation anisotropy on non-Ising domain walls (Ising–Néel-like walls). Prior studies on flexoelectric effects and domain architectures motivate examining how flexocoupling anisotropy could shape long-range ordering and lattice types of polar skyrmions.

The authors employ three-dimensional phase-field simulations for PbTiO3 (PTO) thin films on substrates. The order parameter is the spontaneous polarization vector P = (P1, P2, P3), evolving via time-dependent Ginzburg–Landau (TDGL) dynamics: ∂P/∂t = −M δF/δP. The total free energy F follows Landau–Ginzburg–Devonshire theory and includes Landau, elastic, electrostatic, gradient, and flexoelectric energy densities, plus surface energy. Elastic energy uses stiffness tensor Cijkl and electrostriction eij = Qijkl Pk Pl; electrostatic energy includes depolarization with partial screening (Ed = βESC + (1−β)EOC) and external electric fields; gradient energy penalizes polarization inhomogeneity. Flexoelectric energy is expressed in Lifshitz-invariant form with three independent bulk flexocoupling coefficients for cubic symmetry: longitudinal f11, transverse f12, and shear f44; the associated flexoelectric field is Eflexo = fijk ∂εkl/∂xl. Mechanical equilibrium is solved with the top surface stress-free and the bottom clamped; total strain εij includes displacement gradients and any imposed bending strain distributions. Electrostatic equilibrium is solved with screening boundary conditions set by β. Numerical implementation: TDGL solved by Euler iteration with time step Δt = 0.001; simulation supercell 128 × 128 × 12 grid points with spacing Δh = 0.4 nm; in-plane periodic boundary conditions; biaxial misfit strain εm = −0.01; screening factor β = 0.6; random initial polarization noise amplitude 0.001 C/m^2 at room temperature. Bending operations are applied by adding fixed, precomputed strain distributions (pure bending along [010] or [100]). Out-of-plane (OP) electric fields E[001] are varied (e.g., 0, −0.5, −1.25 MV/cm). In-plane (IP) electric fields are applied in some simulations to isolate field-induced anisotropy; in those tests the flexoelectric effect is switched off to assess pure IP-field effects. Characterization includes: structural factor via 2D FFT of (Pz − ⟨Pz⟩) to identify ordering symmetry; Pontryagin density and layer topological charge Q(z) to quantify skyrmion topology. Phase diagrams are constructed in f11–E[001] space and f12–f11 space (with specified f44 and E[001]).

• Diverse long-range ordered polar skyrmion lattices (SkLs) emerge in PTO films under combined OP fields and flexoelectricity: tetragonal SkL (T-SkL) and hexagonal SkLs (H-SkLs) with multiple in-plane orientations—d1H-SkL and d2H-SkL (one basis vector ∥ ⟨110⟩), and v1H-SkL and v2H-SkL (one basis vector ∥ [100] or [010]). • Without flexoelectricity: at E[001] = 0 MV/cm the film forms a labyrinthine state; at −0.5 MV/cm it transforms to a hexagonal SkL (d1H-SkL); at −1.25 MV/cm skyrmions shrink and defects appear (imperfect H-SkL, iH-SkL), with loss of clear six-fold symmetry in the structural factor. • With finite f11 (e.g., f11 = 3 V; f12 = f44 = 0): at E[001] = 0 MV/cm stripes with nearly 90° angles form; at −0.5 MV/cm elongated skyrmions align along [100]/[010]; at −1.25 MV/cm a double-q tetragonal SkL (T-SkL) emerges. Increasing f11 strengthens localized IP flexoelectric fields around skyrmion edges, enhancing skyrmion–skyrmion interactions and introducing orientation anisotropy that drives a hexagonal-to-tetragonal lattice transformation. • Phase diagrams (f11 vs E[001]; f12 vs f11 at E[001] = −1.25 MV/cm) identify eight polar states: C (uniform down polarization), H-L (labyrinthine), T-L (labyrinthine with ~90° stripes), TH-L (intermediate labyrinthine), iH-SkL, H-SkL variants (d1, d2, v1, v2), TH-SkL (mixed hexagonal–tetragonal), and T-SkL. TH-L and TH-SkL appear as transitional states with increasing flexoelectricity, evidencing a gradual structural transformation. • Flexocoupling tensor symmetry governs lattice type: along f11 = f12 with f44 = 0 (isotropic tensor), d1H-SkL is stabilized. Deviations (cubic symmetry with anisotropy along [100]/[010]) transform d1H-SkL into T-SkL. With only f44 ≠ 0, T-SkL cannot be stabilized; combined nonzero f11/f12 with f44 modulates lattice symmetry and thresholds, yielding H-SkL, TH-SkL, or T-SkL depending on values (e.g., f11 = f44 = 3 V, f12 = 0 → H-SkL; f11 = 3 V, f12 = f44 = 0 → T-SkL; f12 = f44 = 3 V, f11 = 0 → T-SkL; f12 = 3 V, f11 = f44 = 0 → TH-SkL). Flexocoupling also raises OP field thresholds between labyrinthine, SkL, and uniform C states. • H-SkL in-plane orientation is reconfigurable by inducing IP anisotropy. Bending along [010] (or [100]) introduces an IP easy axis via electrostriction (Q11 > −Q12 > 0), rotating H-SkL from d1H-SkL to v1H-SkL (or v2H-SkL) upon annealing under E[001] = −1.25 MV/cm with f11 = f12 = 4 V, f44 = 0. Stripe-like intermediate states precede final H-SkL formation. Applying IP electric fields alone (with flexoelectricity off) similarly stabilizes v1H-SkL or v2H-SkL for fields along [100] or [010], and d2H-SkL or d1H-SkL for fields along [110] or [1̄10]. • Nonreciprocal bending response of T-SkL: starting from T-SkL (E[001] = −1.25 MV/cm, f11 = 3 V, f12 = f44 = 0), downward bending B1+ transforms T-SkL to v1H-SkL and the lattice reverts to T-SkL after removing bending; upward bending B1− does not trigger the transformation. The reciprocity reverses when Q × f11 < 0. Mechanism: flexoelectricity stabilizes Ising–Néel-like walls and, due to top–bottom asymmetry of skyrmion bubble polarity (central divergence vs convergence for Q = ±1), creates asymmetric Eflex along [001], yielding an anisotropic energy landscape that depends on bending direction and topological charge. • Overall, flexoelectricity-modified anisotropy stabilizes multiple SkL lattice types and governs continuous transformations (hexagonal → tetragonal), while external strain gradients or IP electric fields enable deterministic rotation of hexagonal SkLs.

The study directly addresses whether polar skyrmion lattices beyond the theoretically predicted hexagonal type can be stabilized and controllably transformed. It shows that flexoelectricity, prominent in thin films with strong strain gradients, modifies in-plane anisotropy to stabilize a spectrum of skyrmion lattices including tetragonal lattices and multiple hexagonal orientations. Phase diagrams elucidate how flexocoupling coefficients and OP electric fields determine phase stability and transformation pathways. Furthermore, macroscopic in-plane anisotropy, induced by bending (via electrostriction) or by direct IP electric fields, rotates hexagonal SkLs among crystallographic orientations, enabling reconfigurable skyrmion crystals. The discovery of a nonreciprocal bending response unique to T-SkL, arising from the coupling between flexoelectricity and skyrmion topology, highlights new electromechanical functionalities. These findings are significant for engineering polar topological textures as tunable metamaterials or photonic/optomechanical platforms, where lattice symmetry and orientation critically influence optical dispersion, nonlinear responses, and periodic strain fields.

Phase-field simulations reveal that flexoelectricity plays a critical role in stabilizing and transforming polar skyrmion lattices in PTO thin films. Beyond hexagonal SkLs, tetragonal SkLs and mixed transitional lattices (TH-SkL) are stabilized as flexocoupling-induced anisotropy increases, with comprehensive phase diagrams mapping stability across f11, f12, f44 and E[001]. Hexagonal SkLs can be deterministically rotated via strain gradients (bending) or in-plane electric fields by tuning in-plane anisotropy, and T-SkLs exhibit a flexoelectricity-driven nonreciprocal bending response linked to skyrmion topology. These results provide design guidelines for creating and manipulating ordered polar skyrmion crystals, with implications for reconfigurable optical and electromechanical devices. Future work should experimentally validate these predictions, quantify material-specific flexocoupling tensors, and explore reciprocal impacts of polar topologies on macroscopic flexoelectric responses and device-level performance.

Findings are based on simulations with specific boundary conditions and parameters (e.g., partial screening β = 0.6, biaxial misfit strain εm = −0.01, supercell size 128×128×12, periodic in-plane boundaries). The in-plane model size affects skyrmion size though not the main qualitative conclusions. Accurate values and symmetry of flexocoupling coefficients (f11, f12, f44) and electrostrictive constants are critical and may vary among materials or growth conditions. Bending-induced rotations of H-SkLs required annealing in simulations; experimental protocols may differ. Results focus on PTO films and may not directly generalize without material-specific calibration. Experimental demonstration and parameter extraction are needed to confirm phase boundaries, thresholds, and nonreciprocal responses.