Current vortices and magnetic fields driven by moving polar twin boundaries in ferroelastic materials

Ferroelastic twin boundaries often have properties that do not exist in bulk, such as superconductivity, polarity etc. Designing and optimizing domain walls can hence functionalize ferroelastic materials. Using atomistic simulations, we report that moving domain walls have magnetic properties even when there is no magnetic element in the material. The origin of a robust magnetic signal lies in polar vortex structures induced by moving domain walls, e.g., near the tips of needle domains and near domain wall kinks. These vortices generate displacement currents, which are the origin of magnetic moments perpendicular to the vortex plane. This phenomenon is universal for ionic crystals and holds for all ferroelastic domain boundaries containing dipolar moments. The magnetic moment depends on the speed of the domain boundary, which can reach the speed of sound under strong mechanical forcing. We estimate that the magnetic moment can reach several tens of Bohr magnetons for a collective thin film of 1000 lattice planes and movements of the vortex by the speed of sound. The predicted magnetic fields in thin slabs are much larger than those observed experimentally in SrTiO3/LaAlO3 heterostructures, which may be due to weak (accidental) forcing and slow changes of the domain patterns during their experiments. The dynamical multiferroic properties of ferroelastic domain walls may have the potential to be used to construct localized magnetic memory devices in future.


INTRODUCTION
Interface physics has attracted much interest in the field of ferroelastic materials. Known as domain boundary engineering, ferroelastic materials can be functionalized by designing and optimizing domain walls 1,2 . Many studies have observed that interface possess properties that do not exist in bulk, such as the superconductivity 3-5 , ferroelectricity 6,7 . For example, twin walls are known to be superconductive in WO 3-x 8,9 . Another common effect is that polarity is generated near domain walls via the flexoelectric effect 10,11 , e.g., in CaTiO 3 12-14 . Besides the isolated domain walls, emerging functional properties have also been found near the interfaces between different oxides [15][16][17] or different orientations of the same oxide 18 . For instance, the two-dimensional electron gas at the interface between SrTiO 3 /LaAlO 3 based heterostructures has been verified to be responsible for the occurrence of the interfacial conductivity [19][20][21] and superconductivity 22,23 . Among these functional properties, magnetism has rarely been reported to emerge from twin walls. Although some ideas have been proposed to explore the possible magnetisms inside ferroelastic heterostructures 24,25 , the corresponding magnetic signals were too weak for practical applications or even to be measured experimentally. Recently, within the framework of dynamical multiferroicity 26,27 , magnetic signals were shown to emerge from inside moving ferroelectric domain walls. In a similar approach, we expand this idea to the moving kinks in twin walls and needle tips in ferroelastic materials. We argue, using atomistic simulations, that topological nano-structures (stripe twins or needles domains [28][29][30][31][32][33][34][35][36] etc.) can, thus, generate noticeable magnetism. This behavior relies on the nucleation and movement of kinks in twin walls 11,[37][38][39] . These kinks were shown before to be highly mobile and can be accelerated by external strain. They form polar vortices near the kinks, which, in turn, are the origin of displacement currents and hence of magnetism.
Using atomistic simulations, we report a robust magnetic signal induced by moving ferroelastic domain walls based on a twodimensional theoretical model (the sample is charge neutral, no free electrons exist) 11 . The displacement current and magnetic moment are comparable with mobile polar structures and the formation of polar vortices proposed in refs 26,27 , where magnetic moments amounted to~10 −6 μ B in the domain walls. Here we describe much stronger magnetism by moving ferroelastic needle domains although the magnetic moment density is smaller than that in ref. 26 . The primary magnetic field is not generated at the twin walls, but near twin walls, where the displacement current vortices are induced by kinks in the domain walls. Our findings are potentially important for designing nano-devices based on the magnetic field produced by the moving polar ferroelastic structures in thin slabs instead of interfaces in heterostructures.

RESULTS
Displacement current and magnetic field driven by a moving single kink First we analysed toy models based on Landau springs to describe the interatomic potential in ferroelastic materials (see "Methods" section). The ferroelastic shear angles at the equilibrium state are chosen to be 2°and 4°, corresponding to materials like SrTiO 3 and LaAlO 3 (see "Methods" section). In this paper, we show the results using the model with shear angle of 2°(see the simulation results for 4°model in Supplementary materials). In order to illustrate the basic mechanism how a magnetic moment is generated, we initially focus on a single kink motion inside a ferroelastic twin wall under free boundary conditions ( Supplementary Fig. S1). When the system relaxes, the kink moves from left to right. The moving kink induces relative displacements of nearby cations and anions at a time interval Δt, as shown in Fig. 1a. The relative 1 displacements of the cations and anions on upper and lower part of the domain are symmetrical. The moving kink generates displacement currents (see "Methods" section and Supplementary  Fig. S2), where the displacement currents (per time interval) form vortices on either side of the twin wall. The important observation is that these vortices have the same sense of rotation (Fig. 1b). As the two displacement currents rotate in the same direction, their induced magnetic moments are additive. Using the "Biot-Savart" law, the magnetic moment is calculated to be~2.7 × 10 −4 μ B (see "Methods" section and Supplementary Fig. S3), where μ B = 9.274 × 10 −24 JT −1 is the magnetic moment of an electron. This fact is the key point of our argument. It is different from the idea in refs 26,27 , where the displacement current is considered inside domain walls. In our approach we consider displacement currents outside domain walls where displacement currents are induced by kinks and other topological defects. Our walls are ferroelastic and need not to be ferroelectric. This opens applications to many more materials outside ferroelectrics.
Displacement current and magnetic field generated by moving needle domains We now concentrate on topological structures where kinks are predominant. One such configuration is a needle domain where kinks occur inside the needle tip-but not in the shaft 28,29,37,[40][41][42][43][44] . The simulation box has the size 5.6 nm in x direction and 15.0 nm in y direction. We first explore the behavior of a single needle domain. Moving needle domains are extremely common in ferroelastic materials. Their dynamics have been described for several decades 28,29,37,[40][41][42][43][44][45] . The needle domains contain kinks near the tip 46 , where the domain walls appear mesoscopically bent. We first consider three kink pairs inside the domain structure ( Supplementary Fig. S4). The upper and lower sample surfaces were fixed. The local strain distributions near the middle kink pair exhibit typical Eshelby strain patterns 46 (Supplementary Fig. S4). Electric dipoles occur inside twin walls ( Supplementary Fig. S5) due to the flexoelectric effect in the twin walls and at the location of large local inhomogeneous strains ( Supplementary Fig. S4). The out-of-twin-plane dipoles are activated by the flexoelectric effect of the strain gradients ( Supplementary Fig. S5). When needle domains advance and retreat, they carry these dipoles with them and generate dipolar fields that change coherently. The needle movement, and hence the kinks and dipolar fields, oscillate under oscillatory forcing.
We then applied a simple shear to drive the movement of the needles and the kinks in the domain boundaries with a strain rate of 2.7 × 10 8 s −1 (Supplementary Fig. S4). The kinks have the same orientation and produce repulsive force between them. Under weak stress, kink pair 1 is firstly activated to move, followed by the other two kink pairs. Figure 2a shows the atomic configuration of the moving kink pair 3. The relative displacement of cations and anions are generated around the kink pair, as shown in Fig. 2b. Figure 2c shows the corresponding displacement current and the magnetic moment perpendicular to the atomic plane. The displacement current is of the order of~10 −18 A. The current density map for the kink pair is composed of one displacement current (current vortex) inside the needle and two half-displacement currents on each side of the kink (Fig. 2c). The magnetic moment corresponding to the middle full current vortex is 1.2 × 10 −3 μ B . We roughly estimated the magnetic moment density and found a value of 60% of the moving ferroelectric twin walls 26 (see "Methods" section). The magnetic moment depends on the distance between kinks and the needle thicknesses ( Supplementary Fig. S6). The increase of the magnetic moment with increasing kink-kink distance is attributed to the increase of the displacement current while the flexoelectric effect between kink pairs decreases so that the simulated magnetic moment is a compromise between these two factors. In our simulations, kink pairs reach a constant speed of 0.72 km s −1 (pair 1), 0.80 km s −1 (pair 2), and 1.03 km s −1 (pair 3), as shown in Fig. 2d. The small differences manifest that moving needle domains slightly change their shape under high speeds, known as "sharpening the needle". Needles slow down when approaching the lower fixed surface. Two snapshots (Δt = 0.5 ps) of the moving trajectory of every kink pair were then chosen to calculate the current density near every kink pair 26,27,47 . In comparison with the 2°model, the 4°degree model produces stronger displacement currents and stronger magnetic fields ( Supplementary Fig. S6). The ferroelastic switching is achieved by moving kinks in the 4°model, which induces a larger relative displacement of cations and anions.
Another typical microstructure in ferroelastic materials is the comb pattern 28,29,37,[40][41][42][43][44][45]48 , which consists of needle arrays with various distances between the needles. In our model, the stressdriven propagation of comb pattern consists of three vertical thin needles with 4 atomic layers for each ( Supplementary Fig. S7). By controlling the external stress, vertical needles are stabilized at the same distance from the surface with equal needle distances of 1.5 nm. The local strain fields show typical Eshelby patterns ( Supplementary Fig. S7), similar to that in the kink pair system. The profile of each needle consists of two parts: the needle tip and double kinks inside the walls near the tip. The planar dipole configurations now show antiparallel dipole configurations with out-of-twin-plane dipoles ( Supplementary Fig. S8).
A shear deformation drives the propagation of the comb pattern with a strain rate 2.7 × 10 8 s −1 . All needle tips are energetically equivalent. The activation barriers for driving needle tips movement are therefore the same. The tips are activated at the same time with a gradual increase of speed before reaching a steady state at a velocity of 0.24 km s −1 . They finally slow down when moving close to the lower surface. Two atomic configurations with a time interval of Δt = 3.0 ps during the steady-state motion were chosen to calculate the local displacement current density and the corresponding magnetic moment ( Supplementary  Fig. S9). Displacement currents are observed between the vertical needles, as shown in Supplementary Fig. S9. The strength of the local displacement current density is weaker than that in the kink pair system, with magnetic moments of the same magnitude when the geometry effect of the displacement current (larger vortex diameters) are taken into account. The weaker displacement current can be explained by the slower speed of the needle tips and a smaller flexoelectric effect for needle tips with larger tip-tip distances. By comparing the magnetic moment in the 2°model and the 4°model ( Supplementary Fig. S9), we confirm that the larger flexoelectric effect increases the resulting magnetic field.
To generate stronger magnetic fields, we shortened the intertip distances to 0.7 nm (Fig. 3a). The comb pattern now consists of 5 vertical needle domains with a thickness of the needle of only 3 atomic layers. The symmetry breaking of every single needle is enhanced by a strong strain interaction between nearby vertical needles. The application of shear strain moves the comb pattern downwards. Two snapshots on the moving trajectory were used to calculate the displacement current and magnetic moment. The local displacement current increases relative to the sparse pattern because of stronger strain interactions between vertical needles (Fig. 3b). Simultaneously, the diameter of each single current vortex becomes smaller compared with the previous model. The calculated magnetic moment is 7.3 × 10 −4 μ B for each displacement current vortex (Fig. 3c), which is smaller than the value of 1.1 × 10 −3 μ B for individual needles in Supplementary Fig. S9. However, the total magnetic moment of the entire comb pattern becomes larger due to a higher density of needle domains (Fig. 3c). The speed of the comb pattern in Fig. 3d is 0.24 km s −1 . The system with larger shear angle of 4°generates an enhanced displacement current and magnetic moment (Supplementary Fig. S10).

DISCUSSION
The displacement current and the corresponding magnetic field depend on the distances between kink pairs or the needle density inside the comb patterns. Such tunability stems from the generation of different dipolar moments between needle tips that vary the diameter of the displacement current vortices. The generated magnetic field does not only depend on the moving speed of the polar twin boundaries, but also on the flexoelectricity and geometric configurations of the twin structures. The flexoelectricity characterizes the coupling between strain gradients and polarization. Structural defects carrying strong strain gradients, like moving kinks, needles, and junctions, induce large displacement currents and hence a large magnetic moment. For example, enhancing the flexoelectric polarization by increasing the ferroelastic shear angle produces stronger magnetic fields. Magnetic fields produced by the fast movement of our ferroelastic twins are much larger than the recent investigations on the heterogeneous structures near LaAlO 3 and SrTiO 3 (LAO/STO) interfaces 24,25 . In these observations, the occurrences of the magnetic signals were explained in terms of the couplings between lattices and spins or free electrons near the LAO/STO interfaces or the domain walls below T_c = 105 K for SrTiO 3 . Such couplings are normally very weak in ferroelastic materials, and the observed magnetic signals are therefore very small. The strong increase of the magnetic field in our simulations stems from the fast movement of kinks in ferroelastic walls, which produces significant displacement currents and hence potentially bigger magnetic signals than those observed in the LAO/STO interface.
Ferroelastic thin slabs often consist of more than 1000 atomic layers. In addition, the movement of the polar ferroelastic twin walls is often much faster than those shown in our simulations (even supersonic speed has been found 46 ). In a rough estimation we use the speed of sound of aluminium (6.42 km s −1 ) as the speed of the patterns and extend the slab thickness to 1000 layers. Under these circumstances we find a very large magnetic moment (m = 78.11 μ B ) during the motion of the ferroelastic comb structures (Fig. 3). Such large magnetic moments are sufficient for applications as magnetic memory devices. Nevertheless, the We finally proposed an experimental arrangement to measure the magnetic field. Recently developed experimental technique could prepare ferroelastic materials with a high density of domain walls. With such a sample, the first step is to apply a local stress in an atomic force microscope (AFM) to move atomic kinks and needle domains inside the polar twin walls, as shown in Fig. 4a, b. Simultaneously the magnetic field is measured using a superconducting quantum interference device (SQUID) 24,25 . This experimental arrangement is similar to that proposed to measure the magnetization in a moving ferroelectric domain wall 26 , while has a major difference. In order to measure the magnetic signal generated by the ferroelectric walls, a pick-up loop is positioned on top of the ferroelectric wall. Such high-resolution SQUID pickup loops have typically one micro-meter diameter and are expected to be miniaturized further in future. The magnetic fields emerge from the displacement current vortices near the moving twin walls where the diameter of the current loop can be approximated by the dimension of the pick-up coil. A measurable magnetic signal is then expected to be observable directly in SQUID experiments.
In ferroelectric/ferroelastic material, such as BaTiO 3 , the required thin stripes with periodic twins (similar to Fig. 3, but with larger twin spaces) have been successfully fabricated and carefully controlled 49 . Other materials like LaAlO 3 and Pb 3 (PO 4 ) 2 were equally used to control the ferroelastic domain patterns 50,51 . The design of nano-devices based on the magnetic field produced by the moving polar ferroelastic twin walls then might be possible based on our proposed coupling mechanism between ferroelastic strain, polarity, and magnetic fields.

Interatomic potential and model construction
Our simulations are based on a two-dimensional ionic spring model 11 with two base atoms (A and B), carrying negative and positive charges, respectively. The needle domain structures (A sublattice) were constructed by using nonlinear Landau springs while the interactions between B atoms and A-B atoms are purely harmonic to exclude any additional polar instabilities in bulk. We constructed two models with shear angles of 2°and 4°, respectively. The construction of the 2°model were inspired by SrTiO 3 with T_c = 105 K. The construction of the 4°model is inspired by the fact that most typical ferroelastic shear angels are below 4°, such as LaAlO 3 43 , CaTiO 3 [12][13][14] . The 4°model can be used to explore the possible maximum value of magnetic signal we may detect. All parameters for the two models are shown in Table I in the Supplementary Materials. The polarity occurs along the twin walls due to strong flexoelectricity. Additional long-range Coulomb forces were considered between A-A, A-B, and B-B atoms. Free boundary conditions were used in both x and y directions. Both directions are charge free (no excess charges) by changing charges of the surfaces layers to be one half of that in the inner layers 33 . The temperature of the sample was controlled by the Nosé-Hoover thermostat 52,53 . All the calculations were carried out with the LAMMPS code 54 .

Calculation of displacement current density
The displacement current density in each unit cell is calculated as Δt þQ À sanion Δt (ref. 47 ), where j cation and j anion are the current density contributed by cation and anion, Q + and Q − are the amount of charge carried by cation and anion, v cation and v anion are the moving velocity of cation and anion, s cation and s anion are the moving distances of cation and anion during a time interval Δt. More details are shown in Supplementary Fig. S2 and refs 26,27,47 . Calculation of magnetic moment produced by the displacement current The magnetic moment perpendicular to the current vortex plane can be calculated as m ¼ 1 47   the atomic mass and L i is the orbital angular momentum for the ith moving particle. The direction of L i is proportional to the vector product r i × v i . The magnetic moment can be further calculated as m ¼ 1 where r i is the position vector, v i is the moving velocity and θ i the angle between the position vector and moving direction for the ith particle. More details are shown in Supplementary Fig. S3 and refs 26,27,47 .

Estimation of magnetic moment density
We chose a kink pair to estimate the magnetic moment density (Fig. 2). For a moving kink pair with velocity 1.03 km s −1 , the induced magnetic moment is m = 1.2 × 10 −3 μ B . The area of displacement current vortex contains ca. 81 unit cells. The relevant volume V is then one atomic layer times 81 unit cells and the moment density can be approximated as m/V = 1.2 × 10 −3 μ B /(1 atomic layer × 81 unit cells). Comparing with the results in ref. 26 , their magnetic moment is m = 2.47 × 10 −5 μ B per unit cell right at the moving ferroelectric domain wall where the volume is one unit cell, i.e., m/V = 2.47 × 10 −5 μ B /(1 unit cell). The ratio between the kink model can then be roughly estimated to be 60% of the domain walls 26 .

DATA AVAILABILITY
The data that support the findings within this paper are available from the corresponding authors upon reasonable request. Fig. 4 Experimental arrangements for measuring displacement currents and the corresponding out-of-twin-plane magnetic field. The signals are generated by (a) moving kinks and (b) moving needles in a ferroelastic material. The local stress is applied by the AFM tip to drive the motion of the polar ferroelastic twin walls.