Influence of kinks on the interaction energy between ferroelastic domain walls in membranes and thin films

INTRODUCTION

Ferroelectric and ferroelastic materials spontaneously split into domains where the order parameter (spontaneous polarisation or spontaneous strain, respectively) is uniform. The boundaries between these domains are called domain walls. They move in response to an applied external field (electric field or stress field)^[1-6]. They also exhibit emergent properties that do not exist in bulk, such as a spontaneous polarisation in ferroelastics^[7] or anomalous electrical conductivity in ferroelectrics^[8-11].

Landau theory predicts that domain walls are smooth with a hyperbolic tangent profile. However, experimental observations demonstrate a different behaviour, showing that domain walls can exhibit complex profiles with meanders and atomic steps called kinks^[12-15]. This phenomenon is particularly evident in various materials such as membranes^[16], ceramics^[17,18] and thin films^[19]. This internal structure has direct consequences on the emergent properties of domain walls. In lithium niobate, scanning transmission electron microscopy images on a lamella cut from a bulk single crystal revealed kinks and antikinks (i.e., atomic steps in the opposite direction compared to kinks) at domain walls. It was proposed that these kinks and antikinks lead to localised electric charges that influence the dielectric response of the material^[20].

Kinks are also essential elements that form during switching under an applied external field^[21-23]. Clear experimental evidence supports the notion that kink formation and local bending constitute the first stage of domain wall motion^[24]. As such, the collective motion of domain walls in avalanche-like processes is triggered by kinks^[25]. Kinks seem to play an even more important role in switching along non-polar directions or in non-polar materials. In lithium niobate, switching on non-polar cuts is governed by the generation and propagation of charged kinks^[13]. Simulations also indicate that kink movements dominate the switching mechanism of polar domain walls in non-polar ferroelastic^[26]. In addition, the movement of kinks has been predicted to be supersonic, opening possibilities for materials applications at GHz frequencies^[27].

In samples with high densities of domain walls, interactions between walls are governed by their junctions^[28,29]. In simpler but not uncommon configurations of parallel domain walls separated by distances larger than the boundary thickness, the interactions may arise from kinks^[30] that are known to lead to enhanced areas of strain^[18]. For nanoscopic sizes, the interaction energy between domain walls was found to be dipolar, i.e., to decay as 1/d² where d is the distance between the walls^[30], which agrees well with the energy decay of surface steps^[31,32]. For large samples, the interaction energy was found to decay as 1/d^[33].

In this work, using reasonable interatomic potentials, we show that a crossover between “dipolar” 1/d² interactions and “monopolar” 1/d interactions in free-standing samples occurs at lateral sizes in the order of 1,000 l.u. (l.u. = lattice unit), which corresponds to around 0.5 μm for a unit cell parameter of 0.5 nm. This behaviour changes greatly if the sample is clamped on one side, such as in thin films, since we find that the kink-kink interactions stay monopolar even for very thin films (~ 100 l.u.). Our findings are important for understanding the organisation of domain walls and their response to an applied external field in membranes, transmission electron microscopy lamellae, ceramics, and thin films.

METHODS

Ferroelastic domains and confined atomic kinks residing inside domain walls are described by a model for ferroelastic transitions based on a Landau-type double-well potential^[34], as schematically shown in Supplementary Figure 1. The potential energy U(r) contains three terms: the harmonic first nearest atomic interactions U(r) = 20(r - 1)² (black springs), the anharmonic second-nearest interactions U(r) = -25(r - √2)² + 20,000(r - √2)⁴ (yellow springs) along diagonals in the lattice unit and the fourth-order third-nearest interactions U(r) = 8(r - 2)⁴ (green springs), where r is the distance between atoms. The first- and third-nearest interactions are related to the elastic interactions and constitute the elastic background in ferroelastic materials. The model parameters are inspired by SrTiO₃ with the energy scale determined by T_c = 105 K and typical ferroelastic shear angle of 2°. The atomic mass is M = 50 amu. The relevant strain components ε_xx, ε_yy, and ε_xy are calculated from the appropriate interatomic distances relative to those of the monodomain. The sample was relaxed for ca. 10⁵ computational steps. All simulations were performed using the LAMMPS program^[35], which minimises the potential energy of the total system. Although our formulated interatomic potential is simple, the main elastic properties of ferroelastic materials have been considered such that the simulated results can successfully reveal the elastic interactions between kinks, which would also be predicted by advanced methods, such as analytical potentials^[36] and DFT calculations, where more parameters and physical processes are considered. Two different boundary conditions are used in our simulations. The first one is the open or Dirichlet boundary condition applied in both x and y directions, where sample relaxations, including shape changes and rotations, are allowed. The initial lattice parameters in x and y directions are set to a = 1 l.u. and relax to 1.0001297 l.u. and 0.9995027 l.u. in x and y directions, respectively. A single domain wall with a kink residing at its centre was constructed to investigate the size dependence of tilt angles and self-energies of the kink [Figure 1A and B]. Two parallel domain walls containing one kink in each wall, i.e., a kink-kink pair and a kink-antikink pair [Figure 1A] with various separations [Figure 1C] were constructed to investigate the effect of the sample size on the kink interactions.

Figure 1. Domain wall configurations for kink interactions. (A) Domain wall configurations containing one or two kinks. The structural details of the atomic kinks, kink pairs and strain fields near the kinks are shown in (B and C). The strain maps in (A-C) are colour-coded by the atomic-level strain: ɛ_xx. Green dashed lines in (A) indicate the position of domain walls. Red solid lines in (A) indicate the atomic steps of the kink structures.

The second boundary condition was constructed with the bottom layer fixed while all other surfaces were free to relax. This configuration represents an extreme case of hard interfacial bonding without any lattice defect between the sample and the substrate. This situation is encountered in thin films on “hard” substrates. All other parameters were identical to the first case of open boundary conditions. The kinks were initially created inside domain walls, and the system was then relaxed using a conjugate gradient method followed by 5 × 10⁶ (5 × 10³ ps) molecular dynamics (MD) simulation steps to obtain the full ferroelastic domain structure. Ferroelastic domain structures were obtained by averaging structural snapshots every 1,000 MD steps (1 ps). To avoid the movement of kinks in the domain wall, the temperature was kept very low at T = 0.001 K using a Nosé-Hoover thermostat^[37]. All simulations were performed using the LAMMPS code. The atomic configurations were displayed using the OVTIO software^[38].

RESULTS AND DISCUSSION

Free standing sample (membrane, lamella)

We first construct a single domain wall with a stable kink located at the centre within a cell with open boundary conditions with a constant lateral size (L_x) of 1,601 l.u. in the x direction and variable vertical sizes (L_y) between 101 l.u. and 1,601 l.u. in the y direction (sample thickness indicated by Δ in Figure 2A). The kink-induced distortion is observed in the strain ɛ_xx map, as shown in Figure 2A. The strain fields generated by the kinks manifest obvious compressive and tensile regimes on the top and bottom of the kink, similar to those of dislocations^[39]. A local bending near the kink, together with a macroscopic tilt of the domain wall on the left and right sides of this kink, is observed in Figure 2A and B. Such local bending can also be discussed under the framework of the Helfrich model^[40]. The tilt angle is defined as the macroscopic angle between the domain wall and the horizontal direction, with the angles θ₁ and θ₂ shown in Figure 2B. All atomic layers from the bottom (green lines in Figure 2A) to the top surface (blue lines in Figure 2A) show similar tilt angles. Figure 2C shows the tilt angles as a function of the sample thickness. Both θ₁ and θ₂ decrease as the system thickness increases, following θ ~ Δ^-1. All these results are in agreement with our previous results on kink interactions^[30].

Figure 2. Dependence of the kink-induced wall tilt and kink self-energy on the sample thickness Δ. (A) Strain distribution induced by the kink for a sample with a thickness Δ of 101 l.u. (B) Bending of the local lattice with tilt angles θ₁ and θ₂. The green and blue lines indicate the bottom and top surfaces, respectively. (C) The variation of macroscopic tilt angles as a function of sample thickness Δ. The data points in (C) are fitted by θ = a + b × Δ^c with a ~ 0 for θ₁ and θ₂, b = 2.246 for θ₁ and b = -2.246 for θ₂, and c = -1 for θ₁ and θ₂. The fitted line in (C) shows the scaling exponent of -1 between the sample thickness and the tilt. (D) Relationship between kink self-energy and sample thickness Δ. The data points in (D) are fitted by E_kink = E_core - A × Δ^B with E_core = 0.17 eV, A = 6 eV l.u. and B = -1. The fitted line in (D) shows a scaling of ~1/Δ while logarithmic scaling equally fits the data for thicknesses > 1,000 l.u.

We then calculate the self-energy of kinks in domain walls with different sample thicknesses. The self-energy of kinks is calculated as the energy difference between samples with and without kinks inside the wall. The self-energy is dependent on the sample thickness, as shown in Figure 2D. In analogy with dislocations, the self-energy of kinks consists of two components, namely the core energy of the kink and the elastic energy around the core. The core energy is found to be 0.17 eV, similar to values found for kinks in dislocations in silicon (~0.12 eV)^[41]. The evolution of the total energy with the sample thickness Δ follows a logarithmic size dependence for thicknesses > 1,000 l.u [Figure 2D]. At lower distances, the energy is lower and follows a power law E_kink = E_core - A/Δ with E_core = 0.17 eV and A = 6 eV l.u. The kink energy per unit cell is defined as E_kink divided by the total number of unit cells in the sample and decays with a power law dependence from 10^-7 eV/unit cell in the smallest sample with Δ = 150 l.u. to 2 × 10^-8 eV/unit cell for a sample with a thickness of Δ = 1,601 l.u.

The size dependence of kink-kink and kink-antikink pair interactions is then investigated. Two parallel domain walls with one kink or one antikink in each wall are created [Figure 1A and C]. The two walls were initialised at symmetric positions with respect to centre of the sample, with a kink (or an antikink) at the centre of each wall. The samples range from small (L_x = 201 l.u., L_y = 200 l.u.) to large sizes (L_x = 1,601 l.u., L_y = 1,600 l.u.). To calculate the interaction energy of the kink-kink configuration, the total potential energy is reduced by the potential energy of two noninteracting kinks E_kink-kink = E_total - 2E_kink. Figure 3 shows the sample size dependence of kink-kink interactions. For small sample sizes, the kink-kink interaction energies (black symbols and fitted lines in Supplementary Figure 2) are similar to Lu et al., and the wall-wall interaction shows an atypical “dipolar” character with a scaling exponent of -2 [Figure 3A]^[30]. As the sample size increases, the scaling exponents decrease and “monopolar” interactions with an exponent of -1 are reached when the sample size is larger than 1,001 l.u. This “monopolar” wall-wall interaction is consistent with previous theoretical predictions^[42]. Our energy scaling for the sample thickness reveals a wide crossover regimen, which had not been found before, that occurs near thicknesses of 1,000 l.u [Figure 3B]. To explore further the physical processes leading to this crossover, lattice profiles [Figure 4A-D] and strain maps [Figure 4E] are analysed. Strong strain deformations observed near the surfaces of small-sized sample [Figure 4E(a-c)] decay when the system size increases [Figure 4E(d-f)] and almost disappear when the size is over 1,400 l.u [Figure 4E(g and h)]. The simulated energies as a function of the wall-wall distances are power laws for all thicknesses. They do not represent superpositions of two power laws of more complex functions, which could have represented the intermediate range between the surface-dominated and bulk-dominated interactions. The evolution of this size-dependent surface strain is, thus, seen as the change of exponents of the wall-wall interaction in Figure 3B. The top layers bend nearly parabolically in opposite directions so that the negative displacement of the bottom layer [Figure 4A] is the same as the positive bending of the top layer [Figure 4B]. Small samples bend very strongly, while thick samples are more rigid. The lower and upper domain walls show more local deformations at the kink positions. The same thickness dependence is found for interactions between kinks and antikinks [Supplementary Figure 3].

Figure 3. Interaction energies of kink-kink configurations residing inside two parallel walls as a function of the wall-wall distance d. (A) Interaction energy on logarithmic scales with the fitted scaling exponents. (B) Scaling exponents as a function of sample sizes. The thickness scaling changes from d^-2 for thin samples to d^-1 for thick samples.

Figure 4. Lattice profiles and strain maps of kink-kink configurations with different system sizes and wall-wall distance d. (A-D) Lattice profiles for the bottom surfaces, top surfaces, lower walls and upper walls of samples with d = 10 l.u., 82 l.u. and 200 l.u. Strain maps in (E) are colour-coded according to the atomic-level normal strain: ɛ_xx with a wall-wall distance of d = 200 l.u. The sample sizes for (A-D) are 301 l.u. × 300 l.u., 501 l.u. × 500 l.u., 701 l.u. × 700 l.u., 901 l.u. × 900 l.u., 1,001 l.u. × 1,000 l.u., 1,201 l.u. × 1,200 l.u., 1,401 l.u. × 1,400 l.u., and 1,601 l.u. × 1,600 l.u. The black lines indicate the bottom surface, lower wall, upper wall, and top surface.

The change of the interaction energy exponent from -2 for small samples to -1 for thick samples extends over an interval of sample thicknesses between 700 l.u. and 1,400 l.u. In real materials, the interval extends from 350 nm to 700 nm for a lattice unit of 0.5 nm. The weaker wall-wall interactions could lead to a higher concentration of domain walls.

Sample on a substrate (thin film)

A thin film on a substrate cannot bend on the fixed surface but can relax the free surface. To explore the effect of a fixed boundary on the kink interactions, we simulate this scenario by clamping the lower surface in the limit of hard interfaces with no lattice misfit. The model has a sample size of 601 l.u. in the x direction and 100 l.u. in the y direction [Figure 5A]. Similar to the cases of samples with free boundaries, the two parallel domain walls are symmetric with respect to the centre of the sample, and the kink is stabilised at the centre of each domain wall. We then relaxed the entire sample. The resultant strain fields are shown in Figure 5A. The strain fields do not extend to the bottom surface, while strain fields deform the top surface. Figure 5B shows the local displacements of the fixed bottom surface, lower wall, upper wall and top surface of samples with wall-wall distances of d = 10 l.u., 22 l.u., and 42 l.u. The bottom surface remains flat, while the top surface shows local ridge-and-valley deformations. These local deformations are in the order of 0.02 l.u. or 2% strain for d = 42 l.u. The local deformations of the domain walls are of the same order of magnitude, with a sharp singularity at the kinks. The upper wall exhibits a larger deformation than the lower wall.

Figure 5. Interaction energies of kink-kink configurations with a clamped bottom surface. (A) Strain fields of the thin film with a thickness of Δ = 100 l.u. and wall-wall distances of d = 42 l.u. The strain map was colour-coded according to the atomic-level strain ɛ_xx. (B) Lattice displacements of fixed bottom surface, lower wall, upper wall, and top surface due to the kink-kink interactions. (C) The variation of interaction energies as a function of wall-wall distance. The data points in (C) are fitted by using the equation E_kink-kink = E₀ - A × d^B with E₀ ~ 0 (noninteracting kinks), A = 0.832 eV l.u. and B = -1. The scaling exponent of -1 is shown in (D).

The fundamental difference between the clamped and the free samples is that no macroscopic bending can occur under clamping conditions, while the local deformations are visible in both cases. The decay of the valley structure away from the centre is exponential and extends over some 50 l.u. The decay resembles the deformation caused by intersections between domain walls and surfaces^[43,44].

The interaction between the two kinks, and hence between the walls, is shown in Figure 5C. The scaling function is a power law with an exponent of -1 [Figure 5D]. In straight samples without any macroscopic bending, the scaling E ~ d^-1 holds even though only one surface is clamped and the opposite surface is free to deform. Similar results were found for the kink-antikink interactions, as shown in Supplementary Figure 4.

CONCLUSIONS

Kink-kink interactions in bulk samples interact as “monopoles” with a d^-1 dependence when they are separated by the distance d. As the sample size decreases, the interaction for thin samples decays following a characteristic d^-2 trend similar to that of dipoles. This behaviour of any singularity (dielectric, dislocations, interstitials, etc.) is commonly described analytically by the concept of “image force”. The construction is based on the calculation of the surface relaxation as having the same energy as if a fictitious image force was placed outside the sample. Such image forces have also been used to describe the dynamics of dislocation movements^[45]. Our results clarify the role of the crossover regime near d = 1,000 l.u., which is rather wide.

Detailed investigations of wall profiles are often attempted by transmission electron microscopy where the typical sample thickness is 50 l.u., well within the “dipolar” range. We demonstrate the role of image forces that are responsible for the dipolar relaxation by comparing the relaxation patterns of free-standing samples and thin films. They are closely related to the bending of the sample and much less to the bulging of the surface to form ridge-and-valley structures. This corresponds not only to substrate effects but also to thin wedges in transmission electron microscopy investigations where the imaged part of the sample is very thin but constrained by the thicker part of the specimen. Still, even in clamped samples, such as thin films on rigid substrates, rather large surface structures, such as ridges and valleys, are observed in our simulation and could be detected experimentally.

Our simulated results provide a comprehensive understanding of the elastic interactions between kinks in ferroelastic and ferroelectric domain walls. The existence of a wide crossover regime in free-standing samples indicate that their domain structures should organise differently and exhibit abnormal behaviours in response to external fields, leading to unusual functionalities (e.g., superelasticity). Future research on membranes with in-situ TEM or dynamic piezoresponse force microscopy should provide valuable insights and investigations on novel nanomaterials.