Theory for anisotropic local ferroelectric switching

Theoretical modeling of polarization switching around a biased tip contact is important for fundamental understanding and advanced applications of ferroelectrics. Here we propose a simple in-plane two-dimensional model that considers surface charge transport and the associated evolution of the electric field driving domain growth. The model reproduces peculiar domain shapes ranging from round to faceted in KTiOPO4 (C2v symmetry) and LiNbO3 (C3v symmetry). This is done through modulation of dielectric permittivity, which mimics domain wall pinning on the lattice. In contrast to previous works, which attempted to justify domain anisotropy by means of point symmetry invariants, here we illustrate the necessity of taking translational symmetry into account. The results are pertinent to ferroelectric racetrack memories and other applications requiring domain tailoring.


Introduction
Ferroelectric domains are regions that differ by their direction of spontaneous polarization.Under electrical, mechanical and optical stimuli, domains respond differently and this property is used in numerous applications including memories [1,2], optics [3,4] and electro-mechanics [5,6].Piezoresponse force microscopy (PFM) allows domain imaging with nanoscale spatial resolution [7][8][9][10][11].Moreover, when a large enough voltage is applied to the conductive tip of the piezoresponse force microscope, spontaneous polarization switches locally under the tip, giving the opportunity for precise domain writing.With the help of PFM desirable domain patterns may be tailored with nanoscale precision for biological and chemical sensors [12], data storage [13] and nonlinear optical applications [4].Physical grounds for local domain switching around a piezoresponse force microscope tip are of both applied and fundamental interest.Here the following two subproblems are important: (i) dynamics of the domain wall (DW) in the bulk of the material and (ii) screening of the bound charge on the surface.Dynamics of the DW may be simply described by so-called 'dry' and 'viscous' friction mechanisms [14], which originate from dissipative polarization dynamics, DW pinning on Pierls barriers and defects, etc, and result in a threshold electric field for DW motion and DW mobility [15].Screening on the surface is described by the dynamics of the free charges that, in the case of interest, originate mainly from the piezoresponse force microscope tip and are transported over the ferroelectric surface to the DW [14,16,17].Research on surface screening is of importance because it often limits the speed of ferroelectric devices.To give an example, ferroelectric memories could be better than magnetic ones because ferroelectric domains may be much smaller (10 nm compared to 1 mm) as are their DWs (1-5 nm compared to 10-100 nm), if there were fast and energy efficient ways for their writing and reading [16].During recent years, the important role of the adsorbed water layer has been recognized for the transport of screening charge.The influence of ambient air relative humidity (RH) on the switching kinetics in ferroelectrics has been studied in a number of works [18][19][20].It has been shown that RH strongly influences the domain kinetics in ferroelectrics [18,21].The process of charge transport by the adsorbed surface layer is complex and involves surface electrochemistry [22,23].Current models for surface conductivity exist in the form of hypotheses, one being based on a twodimensional (2D) electron gas occupying surface states [24].Whatever the mechanism, it should not be responsible for the anisotropy of the domain shape.Indeed, surface phenomena are not expected to depend much on lattice orientation.Because in the present work we focus on anisotropy, we restrict the description of this subproblem to a simple model.In our simulations, Ohm's law is assumed to hold on the surface where RH is taken into account by the proxy parameter of surface conductivity.On balance, the model for DW dynamics in the crystal should contain terms responsible for the anisotropy of domain shape.Previous works in this direction assumed the angular anisotropy of DW energy to be responsible for the specific domain shape [25].However, such an interpretation conflicts with experimental observations [26].Experimentally observed domain faceting clearly indicates the need for translational symmetry operators to be of importance for the model, whereas DW energy anisotropy [25] takes into account the point group of the symmetry only.
In the present work, we illustrate that pinning of DWs on the lattice is responsible for the anisotropic domain growth.Modeling results are compared with experiments for domain structure evolution in KTiOPO 4 (KTP) crystals with C 2v symmetry, where two symmetrically non-equivalent DW types, designated as slow and fast walls in [15], coexist.Also, modeling is performed for the LiNbO 3 (LN) crystal family with C 3v symmetry and compared with the results in [27,28].These are the main materials for nonlinear optical applications.Their simple domain structure with 180°DWs allows electro-mechanical coupling to be neglected and discarding of the mechanical degrees of freedom.

The theoretical model
The aim of the model is to simulate domain growth in thin films or thin crystal slices with a continuous bottom electrode by voltage pulses applied via a piezoresponse force microscope tip in contact with the top surface, as schematically shown in figure 1(a).The high-symmetry axis of the crystal is directed along the z-axis of the film.The propagation of the DW is studied during simulations by the finiteelement method using Comsol Multiphysics software.In comparison with existing three-dimensional models [29][30][31] and 2D models for the XZ plane [28,32] we introduce following simplifications.(i) The nucleation process is not of interest here, we focus on lateral domain growth which allows us to reduce the problem to a 2D statement in the XY plane.
(ii) Electro-mechanical coupling is not considered, in accordance with the non-ferroelastic nature of 180°DWs.Namely, both domains have the same value of spontaneous strain; the contrast of piezoelectric response across the DW may lead, however, to some interesting phenomena which are beyond the scope of the present article.(iii) The domain growth is limited with the transport of screening charge from the tip to DWs.This is expected to mainly occur through a nanometersthick layer of moisture on the surface, whose thickness and conductivity increase with increase in ambient humidity [28].
In the frame of our model, 2D conductivity thus carries information on the RH.The transition from three to two dimensions is done under the assumption that electric polarization P z and electric field E do not depend on z, which is a reasonable approximation for poorly conductive materials where DWs tend to preserve electrically neutral orientations parallel to the z-axis.The electrical potential V thus linearly drops through the thickness of the film, the value of the potential at the surface being assigned to a point in the 2D computational domain.
The following set of equations is solved: Here ε b is the background dielectric permittivity of the ferroelectric [33], which describes contributions of modes other than the ferroelectric soft mode, e 0 is the universal dielectric constant and s is the (effective) conductivity of the surface.The ferroelectric polarization dynamics obeys the Landau-Khalatnikov differential equation (5), where a, b, g are Landau coefficients, > g 0 is the gradient term, responsible for the energy and thickness of the DW, and t 0 is a characteristic time.
The function ( ) A x y , , hereafter called the pinning function, has translational symmetry in accordance with the lattice of the material under study.The corresponding term introduces modulation of dielectric permittivity, which results in pinning of DWs on effective Pierls barriers.In the frame of our model, domain growth anisotropy may be entirely described with this term.
The initial state of the system corresponds to a monodomain state with spontaneous polarization -P, s to within random fluctuations of magnitude 0.01%.The voltage at the tip with respect to the bottom electrode is set to > V 0 t , which is applied as a boundary condition to a circle with radius r .0 Zero polarization and zero potential boundary conditions at  ¥ r were implemented using an auxiliary virtually infinite domain (for more details see the Supplementary materials).

Domain growth in KTP
A series of simulations was carried out for the parameters of KTP, symmetry group C 2v (see table 1 for material properties).´-5 10 15  ´-1. 5  ´-1 10 10 ´- shows the domain growth process.The initial shape of the domain is close to round, then it elongates and attains faceting.This corresponds to the idea that a slowly moving DW is subject to stronger pinning effects.The pinning was described using a two-parameter function plotted and compared with the KTP lattice in supplementary figure S2.Note that, for computational reasons, the scale of modulation is p times larger than the lattice constant.The term = Á 4 10 1 8 m F −1 was set to fit the domain size in accordance with experimental observations [26].With increase in A 2 , DW propagation becomes more restricted in the x-direction compared with the y-direction, and thus the height-to-width aspect ratio may be controlled.An aspect ratio of approximately 2.5 in accordance with experiment [26] was achieved with = A 2. 2 The discrete nature of the lattice, modeled here with equation (6), manifests itself in figure 1(b).In addition, we plotted it as a guide for the eye on top of the potential plot in figure 1(d).On domain growth the DW makes a breakthrough, for example in the center of the vertical facet, as shown with the arrows in figure 1(d).Then, two steps are propagating in both directions (see supplementary video 1).Because pinning is weaker at inclined walls, a lower voltage is needed for them to propagate (2.20 V) compared with vertical walls (2.43 V) (see figure 1(d)).This is in accordance with the concept of 'fast' and 'slow' walls in [15].Despite the main goal of our simulations being to reproduce the shape of the domain and uncover the mechanisms leading to domain anisotropy, another important aspect is the time dependence of domain size.Here our model is clearly premature, as it assumes ohmic conductivity, which is a rough approximation; nevertheless, here we devote considerable attention to the interpretation of domain sizes in view of the great interest in this issue.The results are compared with experimental data and analytical formula, previously obtained for circular domains [14] where t is the time of the voltage pulse and R is the domain size.Figure 1(c) shows the dependence of the length and width of the domain versus time of the voltage pulse.The curves are nearly linear in logarithmic scale, which is consistent with experiments [26].Also, as noted in the Introduction, important work in the field is devoted to the dependence of domain shape and size on the RH of the environment during PFM switching experiments [18-21, 28, 34, 35].In the frame of our model, this is described through variation of the surface conductivity, which is a proxy parameter.The effect of conductivity on domain growth is summarized in figure 2. Figure 2   ´-3.22 10 S. 15 Material parameters are listed in table 1.
parameters: it gradually changes from a circle for short times and/or high conductivity to a faceted shape.Figure 2(d) shows a comparison between our calculations for KTP and formula (7).While the two curves are similar, the best fit was achieved for a slightly different value for the surface conductivity s = ´-3.22 10 S, 15 as compared with s = ´-5 10 S 15 used in the simulations.This is because formula (7) neglects pinning of the DW on the lattice and the anisotropy of the material.The simulations in this sense refine the results with respect to (7).
Further, we explore the details of charge transport.Figure 3 illustrates currents of the screening charge.The configuration essentially has a fixed charge source, which is the atomic force microscope tip in the center.The sinks of the charge are the points of domain growth, and these are multiple and moving with time.Therefore, there is a complex arrangement of currents, which, however, obey the following logic.Conduction currents are almost isotropic and pump the charge to and behind the domain with small regard for the switching process (figure 3(b)).Displacement currents (due to the in-plane dielectric permittivity of the system) adapt quickly and deliver charge to the sinks from the surrounding areas (figure 3(c)).Finally, figure 3(d) shows the z-component of displacement currents, which essentially corresponds to polarization dynamics.Parenthetically, a test simulation was run with zero in-plane permittivity to show that this value has only minor importance for polarization dynamics.In that case, the in-plane displacement current is null and the conduction current roughly becomes a superposition of the current densities shown in figures 3(b) and (c).

Domain growth in LN
Further, we simulated domain growth in trigonal systems such as lithium niobate (LN).The pinning function was used in the form of ( ) where = a 0.51 nm is the scale of modulation.Note that the scale of modulation is p 2 times larger than the lattice constant of the real LN material for computational reasons.The first term here has hexagonal symmetry, whereas the second one breaks it down to trigonal.Figure 4 shows domain growth for the case where A 3 = 4 10 8 m F -1 and A 4 = 0.The shape of the domain gradually changes from round to hexagonal.DWs also propagate by a breakthrough in the center of facets and then by steps flowing from the center to the corners.A few such steps are visible in the figure, and more details are give in supplementary video 2. A hexagonal domain corresponds to experimental observations for LiNbO 3 [36].The experimental data for domain size in KTP, adapted with permission from [26].Copyright © 2021 American Chemical Society.Solid lines are the theoretical fits using formula (9).(b) The experimental data for domain size in LN, reprinted from [36], with permission of AIP Publishing, and their theoretical fits using formula (9).Fitting parameters listed in supplementary tables s1 and s2.
experimental results [34,36].As for the charge transport, it is analogous to that discussed for KTP but now with trigonal symmetry (figures 4(c) and (d)).Notably, equipotentials are concave near the points of growth: this is consistent with sinks of conduction current at these points.Furthermore, our results for A 4 = 1/6 presented in the supplementary materials illustrate that the method may describe the more complex domain shapes observed, for example, for lithium tantalate [25].

Comparison with experimental data
As noted above, Ohm's law is only a rough approximation for the surface conductivity as applied to polarization reversal around a piezoresponse force microscope tip.Such a law would correspond to a uniform carrier density and mobility on the surface.In real systems, it is expected that more carriers accumulate near the tip where the electric field is maximal.Correspondingly, in real systems domain growth is expected to be even faster at the initial stage and then slow down far from the tip compared with the idealized case of ohmic conductivity.In order to better compare the results of our theory with the experimental data we have updated formula (6) to a more general dependence, where the surface conductivity follows a power law on radial coordinates Then, along the lines of the theory [14] (for derivations see supplementary note 3) one obtains where r 0 is the radius of contact with the piezoresponse force microscope tip, set for all cases to 2 nm, and l and a are the fitting parameters.Fits of the experimental results for KTP [26] and LN [36] are plotted in figures 5(a) and (b), respectively.This analysis is still rough, as it does not consider the kinetics of chemical reactions on the surface, but it does give an idea of the surface conductivity of the ferroelectric samples during local polarization reversal around a piezoresponse force microscope tip and how it may depend on the radial coordinate (see supplementary note 3 and figure s4).

Summary and applications
In summary, we theoretically describe anisotropic domain growth in z-cuts of ferroelectric materials, a process extensively studied experimentally in pursuit of the development of memories and agile electronic devices.In the frame of our model, the shape of the growing domain is determined by a combination of effects related to transport of screening charge and pinning of DWs on the crystal lattice.As a result, the domain shape ranges from round to faceted depending on time and surface conductivity in good agreement with experimental observations.This justifies the concept of DW pinning on the lattice being responsible for domain shapes, in contrast to earlier works on DW energy anisotropy which are in poor agreement with experiment.Speaking about potential applications of the results for memories, the model is more pertinent to race-track memories [16], where DW displacement is traced.Storage density for such memories will rely on reproducibility of domain shapes and sizes, and their speed will be related to surface conductivity, to which we make a theoretical contribution.Since we do not focus on the nucleation process, and thus do not compute minimal domain size, applications to bubble memories are limited.The model readily applies to materials with other symmetry provided that only two out-of-plane domain states are enabled by the switching process.This way, in cubic-to-tetragonal ferroelectrics such as BaTiO 3 or Pb(Zr,Ti)O 3 we expect octagonal faceting.We highlight that the faceted domain shapes are expected in pure crystals, while in materials with a large concentration of defects, or in materials with an irregular lattice such as polymer-based ferroelectrics (e.g.PVDF) isotropic round domains are expected or domains with irregular shapes, as determined by structural fluctuations.These may be readily modeled using the same approach through introduction of irregular pinning functions, point lattice defects or random fields in the computational domain.To conclude, we believe that the simple 2D model we developed may become the basis for more powerful resource-efficient computational tools.In particular it may be used to describe local domain growth in more complex multi-physical or even physicalbiological processes which are now in trend.

Figure 1 .
Figure 1.Domain shape simulated for KTP.(a) Schematic of the process showing a biased atomic force microscope tip in contact with the ferroelectric film and a growing domain.(b) Domain evolution for three consequent time moments: color shows polarization.(c) The width and length of domains during the numerical experiment.(d) Distribution of electric potential with equipotential lines for domain growth.Arrows indicate 'unit cell' breakthroughs of the domain wall.Material parameters are listed in table 1.
(a) shows the dependence of domain length and width on the surface conductivity for a 50-s long voltage pulse.Time dependence of domain length for a few values of conductivity is plotted in figure 2(b).
Figure 2(c) shows the domain shape for selected

Figure 2 .
Figure 2. Impact of surface conductivity on domain growth in KTP.(a) Dependence of domain size on conductivity.(b) Domain length evolution for different conductivities.(c) Domain shapes for various values of conductivity and time.(d) Comparison between simulations and formula (7) with s = ´-3.22 10 S.

Figure 3 .
Figure 3. 2D distribution and directions of the currents by streamlines in KTP.(a) Distribution of the absolute value for the current density.(b) Distribution of the conductive current density in plane: color, absolute value; arrow, direction.(c) Distribution of the displacement current density in plane: color, absolute value; arrow, direction.(d) Distribution of the z-component of displacement current density.Values of inplane current densities are converted from surface currents to A m -2 through division by 100 nm.Material parameters are listed in table 1.
Figure 4(b) shows a comparison between our calculations and formula (7) [14].The two curves are similar, and the best fit was achieved for only a slightly different value for the surface conductivity (s = ´-2.61 10 S 16 ) compared with s = ´-1.5 10 S 16 used in the simulations.The size of the hexagon plotted in figure 4(b) corresponds to the distance from the center to a vertex.The shape of the curve qualitatively reproduces

Figure 4 .
Figure 4. Domain growth simulated for LN.(a) Evolution of domain growth during the numerical experiment.(b) Comparison between our calculations of the domain size during the numerical experiment and formula (7) with s = ´-2.61 10 S. 16 (c) Distribution of electric potential with equipotential lines.(d) Distribution of the normal current density.Streamlines show directions of the displacement current.Values of in-plane current densities are converted to A m -2 through division by 100 nm.Material parameters are listed in table 1.

Figure 5 .
Figure 5.Comparison of the theory with experimental results.(a)The experimental data for domain size in KTP, adapted with permission from[26].Copyright © 2021 American Chemical Society.Solid lines are the theoretical fits using formula(9).(b) The experimental data for domain size in LN, reprinted from[36], with permission of AIP Publishing, and their theoretical fits using formula(9).Fitting parameters listed in supplementary tables s1 and s2.