Mechanical switching in ferroelectrics by shear stress and its implications on charged domain wall generation and vortex memory devices

Recently, the resurgence of interest in flexoelectricity in solids has promoted research on switching of ferroelectric domains via mechanical loads. In this work, combining thermodynamic calculation and phase-field simulation, we revealed that mechanical switching of polarization can be achieved in ferroelectrics by shear stress via a simple mechanism where the presence of flexoelectricity is not necessary. The switching is a consequence of the trilinear coupling between shear stress and the two polarization components that lie in the shear strain plane. Specifically, when the direction of one polarization component is fixed, switching of the other component can be induced by applying a shear stress. Moreover, when the shear stress is fixed, switching of one of the polarization components can lead to the switching of the other component. Phase diagrams of the stable polarization state as a function of shear stress, temperature and electric field are calculated. Furthermore, as demonstrated by phase-field simulation on a ferroelectric thin film, domain switching can indeed be realized by a local shear stress. Importantly, the combining effect of shear stress and electric field can lead to a deterministic writing of charged domain walls, which should be very useful for the fabrication of domain wall devices. The implication of such a shear-stress-modified polarization switching on the design of vortex memory device is also discussed.


Introduction
'Ferroelectrics' refers to a class of materials that exhibit a spontaneous polarization P below a temperature point and P is switchable under an electric eld E. This attribute makes them promising in developing nonvolatile memory devices. 1 Ferroelectrics also have wonderful applications such as in sensors, 2 actuators, 3 energy conversion/storage devices, 4-6 optical switches, 7 microwave devices, 8 etc., due to a wide spectrum of other functional properties including piezoelectricity, pyroelectricity, electrooptic effects, and nonlinear dielectric/optic behavior. Importantly, the properties of ferroelectrics generally exhibit dependence on the polarization state and the domain structure, and it is wellknown that the polarization state and the domain structure of ferroelectrics are strongly coupled to external thermal, electrical and mechanical elds. [9][10][11] These features indicate abundant controllability of ferroelectric devices via domain engineering.
Due to the direct coupling between polarization P and electric eld E, it is natural to manipulate (e.g., switch, reorient or erase) ferroelectric domains electrically. As the applied eld is required to be larger than the coercive eld, which is usually in order of 10 7 to 10 8 V m À1 , electric-eld-driven processes such as charge injection (extraction) from (into) the electrode, electrochemical process near the electrode/ferroelectric interface, and charged defect migration can be quite signicant. While it is still under debate, the injected charges or mobile defects are believed to be captured by the interface (and form a space charge layer which decreases the effective eld applied and depresses domain inversion) or by the charged domain walls that appear during domain switching (and consequently pin the domain walls and depresses domain inversion). As a consequence, fatigue problem is quite difficult to avoid. 12,13 While internal electric eld might also play an important role, fatigue problem, at least the contribution from charge injection, is expected to be largely eased if external electric eld is absent. Electrical switching also causes complex memristive dynamics, 14,15 leakage and even dielectric break-down problem. It is therefore desirable to identify possible effects that can facilitate domain switching and alternative switching strategies using nonelectrical elds.
Most ferroelectric phase transition behaviors are caused by the collective lattice distortions. Consequently, ferroelectric polarization is intrinsically coupled with stress and strain elds.
However, it is not a common idea to apply mechanical elds to switch ferroelectric domains. This is because stress s or strain 3 is coupled with ferroelectric polarization P in an even symmetry, i.e., ÀsP 2 or À3P 2 . The switching of P to ÀP would not change the system's energy. Nevertheless, it should be pointed out that this argument is only applicable to perfect 180 switching and to homogeneous system. On the one hand, it is well-known that ferroelectric domains can be reoriented by strain or stress, such as strain-induced transformation from a domain to c domain in a ferroelastic domain pattern. On the other hand, when the system is inhomogeneous, e.g., subjected to inhomogeneous stress/strain state, or with surfaces/interfaces, other effects can enter into the role in facilitating mechanical switching of ferroelectric domains. For example, it has been realized that strain gradient can lead to 180 switching of polarization in a similar way of electric eld, namely exoelectricity. This is because strain gradient is an odd-parity tensor without inversion symmetry, similar to an electric eld.
Recently, notable progresses have been achieved in issue of mechanical switching of ferroelectric domains. [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] On the experimental front, Lu et al. 16 for the rst time demonstrated that the stress gradient generated by the tip of an atomic force microscope (AFM) can switch the polarization locally in BaTiO 3 (BTO) thin lm. Aerwards, a number of works have been performed, demonstrating the feasibility of mechanical switching in ferroelectric polymer lms, 17,18 thick lms, 19 multiferroic thin lms, 20,21 semiconducting lms, 22 etc. The effects of substrate mist strain 23,24 and loading mode (i.e., perpendicular mechanical load vs. sliding load) 25 on the mechanical switching behavior, and the asymmetric feature of mechanical switching in the vicinity of an existing domain wall, 26 have been explored. On the theoretical aspect, a number of works based on methods such as nite element simulations, 16,25 molecular dynamics simulation, 26 phase-eld method, 27-29 rst-principles calculation, 30 etc., have also been performed to gain a deep insight into the mechanism of mechanical switching. It should be pointed out that previous works mainly focused on mechanical switching in thin lms under perpendicular tiploading mode or bending mode, where exoelectric eld due to the transverse and longitudinal strain/stress gradients is believed to be important. In such cases, the shear strain/stress gradient is at least an order smaller in magnitude and its inuence is negligible. As two exceptional cases, Očenášek et al. 25 have shown that shear strain gradient can play a dominant role over compressive strain gradient in the exoelectric domain switching process under sliding loading mode; Li et al. 30 have studied the shear strain gradient effect on the polarization reversal of BTO thin lms. Also noteworthy is that, for systems with surfaces or interfaces, other effects beyond exoelectricity such as surface piezoelectricity, 31 interfacial bonding and screening effect, 32,33 contact electrication, 34 and surface electrochemical process 35 are possible to be modied by mechanical elds to facilitate mechanical switching of ferroelectric domain. For example, a recent work has shown that mechanical switching driven by depolarization eld can be achieved in ferroelectric thin lms, with both up-to-down and down-to-up switching modes and a hidden nonlocal switching mode. 33 Considering the controversy in the literature on the magnitude of exoelectric coefficients and many other possible sources that may affect domain switching, it is of signicance to reveal other possible mechanical switching mechanisms.
In this work, based on thermodynamic calculation and phase-eld simulation, we showed an overlooked mechanism of mechanical switching of polarization in BTO via shear stress engineering. In this mechanism, the trilinear coupling between shear stress and the two polarization components that lie in the shear strain plane is exploited, and the presence of exoelectricity is not necessary. Specically, when the direction of one polarization component is xed, switching of the other component can be induced by applying a shear stress. Moreover, when the shear stress is xed, switching of one of the polarization components can lead to the switching of the other component. The stable polarization state as a function of shear stress, temperature and electric eld is calculated. Phase-eld simulation is then performed to verify the shear-stressinduced domain switching in ferroelectric thin lm. Besides the normal switching where neutral domain walls are involved, a deterministic writing of charged domain walls can be also achieved via the combining application of shear stress and electric eld. This should be very useful for the fabrication of domain wall devices. At the end, we discuss that such a shearstress-modied polarization switching can provide us a novel design of vortex memory device.

Thermodynamic model
We model the thermodynamics of ferroelectrics based on the phenomenological Landau-Devonshire theory. The quantities such as polarization and stress are described in a Cartesian coordinate system (Fig. 1a). To consider the response of a thermodynamic system in response to external stress and electric eld, here we adopt the Gibbs free energy. Taking the state of the paraelectric crystal at zero stress and electric eld as the reference, the free energy density of a ferroelectric crystal as a function of polarization, stress and electric eld is written as, 36 a ijkl P i P j P k P l þ 1 6 a ijklmn P i P j P k P l P m P n where ellipsis '.' means higher order terms, a ij , a ijkl and a ijklmn are the phenomenological Landau-Devonshire coefficients, s ijkl is the elastic compliance tensor, Q ijkl is the electrostrictive constant tensor, s ij is the stress tensor, and E i is the electric eld vector. Note that all the indexes are in range 1-3, and the number 1, 2 and 3 correspond to the coordinate x, y and z, respectively. For BTO, an eight-order Landau-Devonshire potential is usually employed. 37 For simplicity, we consider that the BTO crystal is subjected to a uniform shear stress s zx , and an external electric eld along the z-axis E z . The explicit expression of the free energy density is then given by where a 1 , a 11 , a 12 , a 111 , a 112 , a 123 , a 1111 , a 1112 , a 1122 , and a 1123 are the nontrivial Landau-Devonshire coefficients; s 44 and Q 44 are the relevant elastic compliance constants and the electrostrictive constants written in Voigt notation. It can be seen from eqn (2) that the two polarization components, P x and P z , and the shear stress s zx are coupled with each other in a trilinear form À 1 2 Q 44 s zx P x P z . This coupling term, together with the last term, ÀE z P z , enforces a condition of the sign of these three quantities. The other terms of the free energy density do not depend on the signs of P x , P z , and s zx . That is to say, once the signs of two of them are xed, the sign of the third one cannot change arbitrarily, as required by the energy minimum principle. This condition provides us the mechanism of polarization switching by shear stress as schematically shown in Fig. 1b. Specically, two possible switching behaviors of a polarization vector (P x , P z ) can be expected, depending on the relative stability of the two polarization components. In the rst case, the polarization component P z is more stable, and the application of shear stress s zx leads to the switching of P x . In the second case, the polarization component P x is more stable, and the application of shear stress s zx leads to the switching of P z . Note that it is the trilinear term that determines the relative signs of P x , P z , and s zx , as it is the only term that links the three quantities, meanwhile the relative stability of the two polarization components can be modied by external electric eld via the last term. Other factors, such as temperature, mist strain, screening condition, etc., can also modify the relative stability of the two polarization components. The equilibrium polarization state can be obtained by minimizing the free energy density G in eqn (2) The application of shear stress s zx and electric eld E z lowers the degeneracy of the ferroelectric phases of BTO. To classify the different phases, we introduce the following notations: (i) the cphase, where P x ¼ P y ¼ 0 and P z s 0; (ii) the O-phase, where |P x | ¼ |P z | s 0 and P y ¼ 0; (iii) the r 1 -phase, where |P x | ¼ |P y | ¼ |P z | s 0; (iv) the r 2 -phase, where P x ¼ ÀP z s 0, and P y s 0; (v) the r 3phase, where P x ¼ P z s 0, and P y s 0; (vi) the ac 1 -phase, where |P x | s |P z |, P x P z < 0 and P y ¼ 0; (vii) the ac 2 -phase, where P x ¼ ÀP z s 0 and P y ¼ 0; (viii) the ac 3 -phase, where |P x | s |P z |, P x P z > 0 and P y ¼ 0; and (ix) the ac 4 -phase, where P x ¼ P z s 0 and P y ¼ 0. Values of expansion coefficients of the eight-order approximated potential for BTO are from ref. 37, and those of the electrostrictive coefficients and the elastic compliance coefficients are from ref. 38.

Phase-eld simulation
We construct a phase-eld model to investigate the domain switching behavior of BTO thin lm subjected to a local shear stress. In this model, the temporal evolution of the spontaneous polarization eld is described by the time-dependent Ginzburg-Landau (TDGL) equations, where F is the total free energy of the system, M is a kinetic coefficient and t is time.
The total free energy is a sum of the Landau free energy F Land , elastic energy F elas , gradient energy F grad , electrostatic energy F elec and surface energy F surf , that is, where f Land , f grad , f elas , f elec and f surf are the corresponding energy densities, and V and S are the volume and the surface of the nanolm, respectively. The Landau energy density f Land is taken to be at zero stress and zero electric eld, that is, the one given in eqn (2) with the elastic and the electric terms excluded. Under the condition of applied stress, the mechanical stress eld (applied and internal) and its coupling with the polarization eld contribute to the elastic energy density, which is given by, In the absence of body forces, the inhomogeneous stress eld s ij (r) is calculated by solving the mechanical equilibrium equation, i.e., s ij,j (r) ¼ 0, where the comma in the subscript denotes spatial differentiation. Here the thin lm is assumed to be epitaxially grown on an elastic substrate. The lm/substrate system is subjected to the contact force of an AFM tip with a trend of relative motion along the lm surface. The mechanical boundary at the top surface is given by n i s ik ¼ s k , where s k is the local surface traction caused by the tip force and n i is the ith component of the unit vector normal to the surface. Specically, the local surface traction consists of a normal component s z and a lateral component s x , corresponding to a compression force F p and a friction force F f , respectively. In calculating the stress eld, a thick-enough layer of substrate is included and the bottom boundary of the lm/substrate system is xed in displacement. For simplicity, the elastic property of the substrate is taken to be the same with that of the ferroelectric thin lm.
The gradient energy, or called domain wall energy, is caused by the inhomogeneity of the polarization eld. For a perovskite ferroelectric such as BTO, to the lowest order of the Taylor expansion, the gradient energy density takes the form as where G 11 , G 12 , G 44 , and G 0 44 are gradient energy coefficients, and P i,j denotes the derivative of the ith component of the polarization vector, P i , with respect to the jth coordinate.
The electric energy density of a given polarization distribution is given by where 3 b of about 503 0 is the dielectric constant of the background material. 39,40 For a nanolm subjected to an external electric eld E a , the total electric eld E is equal to sum of the external eld E a and the depolarization eld E d . E can be calculated by the electrostatic equilibrium equation for a freecharge-absent body as Due to truncation at the surface of the nanolm, the spontaneous polarizations are inhomogeneous across the out-ofplane direction. Thus an additional surface energy is necessary to describe this intrinsic effect. Using the effective extrapolation length d eff i , the surface energy density can be approximately given by where D 11 , D 22 and D 44 are the material coefficients related to the gradient energy coefficients.
Phase-eld simulation of domain switching of BTO thin lm is performed by numerically solving the TDGL equations (eqn (4)) together with the mechanical and electrostatic equilibrium equations. For simplicity, a two-dimensional (2D) mesh along x and z direction is adopted. Thickness of the BTO lm is taken to be 12 nm, and that of the substrate layer is taken to be 36 nm. The lateral size of the simulation cell is set to be 200 nm. A periodic boundary condition is applied along the x direction. The TDGL equations are solved by using an explicit Euler iterative scheme, whereas the mechanical and electrostatic equilibrium equations are solved by a nite element method (FEM). Values of the expansion coefficients of the Landau free energy, the electrostrictive coefficients and the elastic compliance coefficients are the same as those adopted in the thermodynamic calculation. Values of the gradient energy coefficients are taken to be:

Results and discussion
Polarization state of BTO crystal at zero external electric eld We start our investigation on the temperature evolution of the polarization state of BTO crystal subjected to different values of shear stress s zx but at zero external electric eld. For comparison, in Fig. 2a Fig. 2b depicts the phase diagram of polarization state of BTO as a function of temperature and shear stress in the absence of external electric eld. It can be seen that the application of shear stress lowers the degeneracy of the ferroelectric phases of BTO. In the negative stress regime, phases with the polarization components P x and P z in opposite signs are stable. They are ac 1 -phase, ac 2phase and r 2 -phase. In contrast, in the positive stress regime, the stable phases are those with the polarization components P x and P z in the same signs. They are ac 3 -phase, ac 4 -phase and r 3phase. The increase of shear stress magnitude |s zx | from zero to 4 GPa increases the paraelectric-ferroelectric transition temperature from about 395 K to above 520 K, meanwhile it decreases the ac 2 /r 2 and the ac 4 /r 3 -phase transition temperatures from about 200 K to below 40 K. Note also that ac 1 and ac 3phases are stable at relatively high temperature and small stress. With the increase of stress, the stable temperature window of the two phases decreases, forming two inverted triangular regions in the phase diagram. To see more clearly the change of polarization, in Fig. 3 we depict the temperature evolution of polarization state of BTO under different values of shear stress s zx and at zero external electric eld, with Fig. 3a-d corresponding to s zx ¼ À0.4 GPa, 0.4 GPa, À2 GPa, and 2 GPa, respectively. Paraelectric-ferroelectric phase transition occurs at a temperature of about 400 K for the cases of |s zx | ¼ 0.4 GPa, whereas for |s zx | ¼ 2 GPa the system exhibits paraelectric-ferroelectric phase transition at about 455 K. Moreover, the paraelectric-ferroelectric phase transition is rst-order at |s zx | ¼ 0.4 GPa and it turns to be second-order at |s zx | ¼ 2 GPa. The ferroelectric phases are quite different due to the shear stress. Specically, at s zx ¼ À0.4 GPa, the system adopts ac 1 -phase at T ¼ [325 K, 400 K], ac 2 -phase at T ¼ [325 K, 280 K], and r 2 -phase at temperature below 170 K. At s zx ¼ 0.4 GPa, the system adopts ac 3 -phase at T ¼ [325 K, 400 K], ac 4 -phase at T ¼ [325 K, 280 K], and r 3 -phase at temperature below 170 K. At s zx ¼ À2 GPa, the system evolves into ac 2 -phase at about 455 K, then into r 2 -phase at about 95 K. At s zx ¼ 2 GPa, the system evolves into ac 4 -phase at about 455 K, then into r 2phase at about 95 K. For negative stress s zx < 0, the two components P x and P z are always in opposite signs, i.e., P x P z < 0, whereas for positive stress s zx > 0, they are always in positive signs, i.e., P x P z > 0. The change of direction of the two polarization components due to the change of sign of shear stress is clearly seen.

Polarization state of BTO crystal at nonzero external electric eld
Note that while the above result indicates the feasibility of changing the relative sign of the two polarization components by shear stress, it is not deterministic which polarization component would change its sign when shear stress applies. This is due to the fact that, phases with P x > 0, P z < 0 or P x < 0, P z > 0 are possible at negative shear stress, meanwhile, phases with P x > 0, P z > 0 or P x < 0, P z < 0 are possible at positive shear stress. In practice, the direction of which polarization component will be switched depends on the relative stability of the two polarization components (see, Fig. 1b). Such relative stability can be modied by many factors, such as temperature, depolarization eld, external electric eld, and mist strain. Take ferroelectric thin lm as an example. The more stable polarization direction can be tuned from out-of-plane to in-plane when the mist strain changes from compression to tension. Another possible effect that can lead to a deterministic switching is that the degeneracy of the polarization states with opposite polarization component is broken. For example, due to the application of an external eld along positive z direction, states with P z > 0 should be more stable than states with P z < 0. As a result, a switching from states with P z > 0 to states with P z < 0 is not likely to trigger by the shear stress, hence P x rather than P z would change sign. To see it more clearly, Fig. 4 depicts the stable polarization vector of BTO as a function of shear stress s zx and external electric eld E z at room temperature. Result shows that due to the external electric eld, the direction of P z component is controlled by the external eld. Consequently, the direction of P x component depends on the sign of the shear stress. A deterministic switching of the P x component is thus achieved. Note that, besides the application of external eld, the penalty of domain wall energy can also cause a broken degeneracy of the polarization states. A deterministic switching of the polarization can occur in a natural way without the need of an external eld. This point will be claried later by our phase-eld simulation.
Moreover, Fig. 4 also indicates that when the shear stress s zx is xed, the change of external eld E z can lead to switching of P z component together with the switching of P x component. Therefore, shear stress can lead to a switching of in-plane polarization by an out-of-plane electric eld. To see more clearly, in Fig. 5 we depict the polarization state of BTO as a function of external electric eld under different values of shear stress s zx and at room temperature, with Fig. 5a-d corresponding to s zx ¼ À0.4 GPa, 0.4 GPa, À2 GPa, and 2 GPa, respectively. It can be seen that while P z component changes its sign with the inversion of the external eld, P x also changes its sign with the inversion of the external eld. Moreover, the signs of P x and P z are either the same or opposite depending on the shear stress. Note also that while the magnitude of P z increases as the external eld increases, the magnitude of P x decreases with the external eld. Such a change of the magnitude of P x can exhibit either a linear feature (see, Fig. 5c and d for |s zx | ¼ 2 GPa) or a nonlinear feature (see, Fig. 5a and b for |s zx | ¼ 0.4 GPa) depending on the magnitude of the shear stress.

Shear-stress-induced domain switching in ferroelectric thin lm
The effects of shear stress/strain on polarization stability have been discussed by several works. 25,30,[43][44][45] Particularly, using a thermodynamic model, Zembilgotov et al. 43 have calculated the phase diagrams of the stable polarization state in epitaxial lms subjected to an in-plane mist shear strain. Gruverman et al. 44 observed that polarization switching could occur far beyond the tip contact area in polycrystalline lms. Shear stress deformation of the grain underneath the tip was believed to cause this abnormal switching behavior. More recently, it has been shown that shear strain gradient can play a dominant role over compressive strain gradient in the exoelectric domain switching process. 25 First-principle calculation has been also applied to reveal the shear-strain gradient induced polarization reversal in BTO thin lms. 30 These works indicate that shear stress/strain (gradient) can be effective to control the polarization stability of ferroelectrics. Here we emphasize that polarization switching can be induced by shear stress, due to the trilinear coupling effect between shear stress and the two polarization components that lie in the shear strain plane. This shear stress effect has been overlooked in recent discussion on First of all, this effect can be utilized to design generator or sensor devices based on ultrathin ferroelectric lms. The merit of such devices lies in the fact that, in ultrathin ferroelectric lms, depolarization eld tends to force the polarization along the inplane direction, which leads to a depression of piezoelectric response that utilizes the coupling of out-of-plane strain and polarization. In contrast, piezoelectric response based on the shear mode can enter into the role. This should has potential applications such as power generation or state detection of uids where shear ow is important. Another potential application of the shear-stress-modied polarization switching effect is memory device. One possible working mode of the memory device is that exploits mechanical switching of the domain information without the need of external electric eld. Domain switching can be induced in a thin lm with orthorhombic or rhombohedral single-domain state, by exerting a local shear stress to the lm. The local shear stress may be generated by the friction force via an AFM tip working at the sliding loading mode, or even under perpendicular loading mode. This is due to the fact that two shear stress regions with signicant values can be induced inside a lm when it is subjected to a compression tip force. Note that the contribution of shear stress to domain switching would be much less important for the switching of tetragonal single- Fig. 4 The stable polarization vector of BTO as a function of shear stress s zx and external electric field E z at room temperature (left panel). Different colors correspond to different directions of the polarization vector as defined in the right panel. domain state, as it is based on the coupling between two polarization components. Another possible working mode is that exploits the complicate domain switching behavior under the combining application of local shear stress and external electric eld, which may be applied by a biased AFM tip working at the sliding mode. Note that for such a working mode, head-to-head or tail-to-tail charged domain walls can be generated in a deterministic way according to the result shown in Fig. 4. Such charged domain walls can be stabilized in samples with enough free charge carriers. This should be very useful for the design of domain-wall-based memory devices.
To verify the feasibility of domain switching via shear stress, we further perform phase-eld simulation on a 12 nm-thick BTO thin lm under the contact of an AFM tip (with a radius of 50 nm) at the lm surface. The AFM tip is subjected to a compression force F p along the lm thickness direction (i.e., negative z direction), which imposes a nonuniform surface traction s z within the contact region. Moreover, the tip has a trend of relative motion along the lm surface (i.e., positive or negative x direction). This leads to a lateral surface traction s x due to the static friction. Assuming that the static friction coefficient is m, we have s x ¼ ms z , and the static friction force F f ¼ mF p . Fig. 6 depicts the distributions of stress components s zz and s zx in the lm under a tip force F p ¼ 1.8 mN with a trend of relative motion along the negative x direction. It can be seen that the tip force induces a compressive stress state s zz along the lm direction right below the tip (Fig. 6a). Moreover, it also causes two shear stress regions with s zx taking opposite signs in the lm at the two sides of the tip (Fig. 6b). Due to the trend of relative motion along the negative x direction, the magnitude of shear stress at the le side (in front of the tip) is much larger than that at the right side (behind the tip).
In the following, we explore if domain switching can occur at the enhanced shear stress region. The electric boundary condition at the top and bottom surfaces of the lm is set to be shortcircuited. The simulation temperature is set to be 100 K and the mist strain is À0.002 so that the lm can be stabilized in a rhombohedral phase. The thin lm is initially set into a singledomain state with P x > 0 and P z > 0. The evolution of this singledomain state in response to the tip force shown above is simulated and the result is shown in Fig. 7a. In response to the local shear stress s zx < 0, the polarization eld nearby the stressapplied region is expected to be switched into P x < 0 and P z > 0, or P x > 0 and P z < 0. However, in fact, the polarization eld is denitely switched into P x > 0 and P z < 0; otherwise, head-to-head and tail-to-tail charged domain walls will appear, which violates the energy minimization principle. An up-to-down domain switching thus can occur at the shear stress region with s zx < 0. This is exactly the case we found in Fig. 7a. One can see that, as the tip force is applied, the polarization eld below the tip becomes mainly along the in-plane direction due to the large compressive stress. Moreover, the polarization eld near the le side of the tip (i.e., at the enhanced shear stress region with s zx < 0) is switched to negative z direction. Aer the tip force is switched off, a switched domain with P x > 0 and P z < 0 is formed slightly ahead of the tip. Similarly, we also simulate the case where the thin lm is initially set into a single-domain state with P x > 0 and P z < 0, and the trend of relative motion is changed to be along the positive x direction. As shown in Fig. 7b, a down-to-up domain switching has been triggered by a local tip force. For this case, the switched domain with P x > 0 and P z > 0 is formed, and its location is also slightly ahead of the tip. Note that our result shows that shear stress can lead to both up-to-down and down-to-up domain switching, which is determined by the sliding direction and the polarization direction of domain. This is different from exoelectric domain switching in tip-lm architecture, where the directions of strain gradients and exoelectric eld are xed, thus one can only expect either up-to-down or down-to-up switching. Note also that, for a lm under sliding load, previous work showed that the domain switching due to exoelectric eld includes a region behind the sliding tip, which is determined by the distribution of the shear strain gradient. 25 This is in contrast to our case, where the switched domain is slightly ahead of the tip, as determined by the  distribution of the shear stress and the direction of the polarization eld (Fig. 6b). In practice, the feasibility of bi-directional switching and the location of switching domain could be evidences of mechanical switching beyond exoelectricity.
As pointed out previously, it is also possible to generate charged domain walls in thin lm under the combining effect of shear stress and external electric eld. To see this, we perform another simulation, with all of the setting similar to the case shown in Fig. 7a, except that now an external electric eld E z ¼ 3 Â 10 8 V m À1 is also applied to the thin lm (at the middle region of the lm with a width 80 nm), mimicking a biased tip working at the sliding mode. The simulation result is shown in Fig. 8. It can be seen that the polarization eld nearby the stress and eld applied region is switched into P x < 0 and P z > 0, with the appearance of tail-to-tail and head-to-head alike domain walls between the switched domain and the unswitched domains, in contrast to the case shown in Fig. 7a where charge neutral domain walls are formed instead. The difference is due to the fact that, the direction of polarization component P z is xed here due to the external electric eld, therefore the polarization component P x has to be switched by the shear stress according to the trilinear coupling. It should be pointed out that the shear-stressinduced charged domain walls we found here is not strictly tailto-tail and head-to-head to avoid the depolarization energy (as no free charge carriers is assumed), and the domain walls would not maintain stable aer tip force and eld is switched off. Nevertheless, such charged domain walls in practice could be stabilized by free charges and mobile charge defects such as oxygen vacancies in the ferroelectric lm, similar to those reported in previous experimental works. [46][47][48][49][50] Our result thus indicates a novel method to write charged domain walls.

Shear-stress-induced vortex switching in ferroelectric nanoring
Finally, it is noteworthy that the trilinear coupling between the shear stress and the polarization components also indicates a potential strategy of switching the chirality of ferroelectric vortices. Although it is very attractive to utilize the chirality of ferroelectric vortices to develop high-density information storage, switching of the vortex chirality is challenged so far. This is due to the lack of a bilinear coupling between the electric eld without curl and the toroidization. Here we show that shear stress might help solve the problem. The basic idea is schematically illustrated in Fig. 9, where switching of the vortex formed in a ferroelectric nanoring with a radius R is taken as an example. The system is placed in a cylindrical coordinate system (r, 4, z). For simplicity, we assume that the nanoring has a uniform polarization P 4 along the angular direction, and is subjected to a shear stress s 4z , which can be exerted by a torque to the nanoring. The toroidization of the nanoring is given by G z ¼ V À1 ∭ V r Â PdV xRP 4 . Then we apply an external eld along the z direction so that the polarization component P z is induced. The coupling energy density between the shear stress and the two polarization components is f c ¼ ÀQs 4z P 4 P z , with Q being the electrostrictive constant. Aer integrating f c over the volume of the nanoring, one can obtain the coupling energy F c ¼ ∭ V f c dV ¼ ÀAQs 4z G z P z , with A being a constant related to the geometry parameters of the ring. Therefore, the toroidization G z , the polarization P z and the shear stress s 4z are coupled in a trilinear form. For a xed shear stress s 4z , a change sign of P z , which is controlled by the external electric eld, can lead to a change sign of the toroidization G z . Moreover, if the external electric eld is xed so that P z is xed in sign, a change sign of shear stress s 4z can also lead to a reversal of the toroidization G z . Note that a polar-toroidal state (G z s 0, P z s 0) is necessary for the appearance of the trilinear coupling. Usually one can apply an external eld to induce such a polartoroidal state. In some systems, polar-toroidal state can be spontaneously formed, e.g., in ferroelectric/paraelectric nanowires. 51 In a future simulation work, we would like to show that these systems are promising to realize electric and mechanical switching of the vortex chirality. The feasibility of a simple switching strategy of the chirality of ferroelectric vortices should be very useful for vortex memory devices.

Conclusions
Thermodynamic calculation and phase-eld simulation have been performed to reveal the effect of shear stress on the polarization switching behavior in ferroelectric with orthorhombic or rhombohedral phase. Using BTO as an example, we show that switching of polarization can be achieved by shear Fig. 8 Phase-field simulation result of domain switching in a 12 nmthick BTO thin film with the formation of charged domain walls under the combining application of a tip force F p ¼ 1.8 mN and an external electric field E z ¼ À3 Â 10 8 V m À1 . The tip has a trend of relative motion along the negative x direction. stress via a simple mechanism. In this mechanism, the trilinear coupling between shear stress and the two polarization components that lie in the strain plane is exploited. Phase diagrams of polarization state as a function of shear stress, temperature and electric eld are calculated. The shear-stressinduced domain switching in ferroelectric thin lm is veried by phase-eld simulation. In addition to the normal domain switching where neutral domain walls are involved, abnormal domain switching with the formation of charged domain walls can be achieved via the combining application of shear stress and electric eld. The role of shear stress in the feasibility of chirality switching of ferroelectric vortices is also discussed. Our work sheds light on the potential use of shear stress or strain in domain engineering of ferroelectrics.

Conflicts of interest
There are no conicts to declare.