Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials

Converse flexoelectricity is a mechanical stress induced by an electric polarization gradient. It can appear in any material, irrespective of symmetry, whenever there is an inhomogeneous electric field distribution. This situation invariably happens in piezoresponse force microscopy (PFM), which is a technique whereby a voltage is delivered to the tip of an atomic force microscope in order to stimulate and probe piezoelectricity at the nanoscale. While PFM is the premier technique for studying ferroelectricity and piezoelectricity at the nanoscale, here we show, theoretically and experimentally, that large effective piezoelectric coefficients can be measured in non-piezoelectric dielectrics due to converse flexoelectricity.


Supplementary Note 1. Constitutive equations and numerical approximation
Following a linear continuum theory of piezoelectricity [1], augmented with flexoelectricity [2], the total electrical enthalpy of a dielectric solid possessing piezoelectricity and flexoelectricity is written as where Ω is the solid domain, C is the fourth-rank tensor of elastic moduli, e is the third-rank tensor of piezoelectricity, k is the second-rank dielectric tensor, and g is the sixth-rank strain gradient elasticity tensor. The last term in Supplementary Eq. (1) is usually discarded, but it is important to guarantee the thermodynamic stability of the model in the presence of flexoelectricity [3,4,5]. The flexoelectric coupling between the strain gradient ∇ε and the electric field E is through the third term, where the fourth-rank flexoelectric tensor µ describes both direct and converse flexoelectric effects [6,7]. The constitutive equations for the mechanical stress σ and the electric displacement D are derived from the electrical enthalpy density H as [8] σ ij =σ ij −σ ijk,k = C ijkl ε kl − e kij E k + µ lijk E l,k − g ijklmn ε lm,nk , where the third term in Supplementary Eq. (2) is the mechanical stress induced by converse flexoelectricity and the last term in Supplementary Eq. (3) is the polarization induced by direct flexoelectricity. The usual stressσ, the higher-order (hyper) stressσ and the electric displacement D are given bŷ These equations are subject to appropriate electrical boundary conditions as whereφ and ω are the prescribed electric potential and surface charge density and Γ φ ∪ Γ D = ∂Ω is the boundary of the domain Ω with a unit normal n i . The mechanical boundary conditions in the form of displacements or traction t can be specified as below: where u i and t k are the prescribed mechanical displacements and tractions, D j = ∂ j − n j D n is the surface gradient operator, D n = n k ∂ k is the normal gradient operator, Γ u ∪ Γ t = ∂Ω, and Γ u ∩ Γ t = ∅. It is clear that the traction boundary condition in Supplementary Eq. (7) is affected by the higher-order stresses. In addition to these, the strain gradients result in other types of boundary conditions as [5]: where υ is the prescribed normal derivative of displacement, r is the higher-order traction, Γ v ∪ Γ r = ∂Ω, and Γ v ∩ Γ r = ∅. Here, we assume natural boundary conditions on Γ v and Γ r , i.e. υ = r = 0, and neglect edge forces arise from a non-smooth boundary [9].
The constitutive equations (2) and (3) in the presence of flexoelectricity are fourth-order which demand at least C 1 continuous basis functions for a direct Galerkin method. Therefore, the weak forms of the extended flexoelectric model cannot be discretized with standard finite elements. Following a recent work [8], we deal with the fourth-order nature of the partial differential equations by approximating displacements and electric potential using a meshfree method with smooth basis functions [10]. Meshfree methods provide an approximation to continuum field equations based on basis functions that do not rely on a mesh or its connectivity. Essentially, these methods allows us to determine a set of smooth basis functions p a (x), each localized around its corresponding node of the grid. We expand the continuum displacement and electric potential fields as and their derivatives as where the first and second variations involve the gradient and Hessian of the LME basis functions, respectively. From now on, we omit the arguments of the basis functions for simplicity, i.e. u = N a=1 p a u a . Introducing these expansions in Supplementary Eq. (1) and considering a 2D axisymmetric problem dealing with rotationally symmetric structures under axisymmetric loading, we obtain its discrete representation as where r and z are the cylindrical coordinates. The stiffness tensor C, the dielectric tensor k, the piezoelectric tensor e, the flexoelectric tensor µ, and the strain-gradient tensor g are written in Voigt form as where The longitudinal, transversal, and shear flexoelectric coefficients are represented by µ 11 , µ 13 , and µ 44 respectively and k 11 is the dielectric constant for an isotropic medium. The gradient operators B u and B φ and the Hessian operators H u and H s for an axisymmetric case are written in Voigt form as The discrete algebraic equations for the equilibrium can be derived following the usual Galerkin procedure as where the local contribution of each quadrature point to the matrix of system has the structure where the derivatives of the basis functions are evaluated at the corresponding quadrature point.

Supplementary Note 2. Contact model
The Signorini-Hertz-Moreau model [12,13] is written for a contact point pair (x; x 0 ) as: where n stands for the outward unit normal vector to the contact boundary Γ c . The normal gap function g n is given by where u n = u·n, u is the mechanical displacement, x is a point on the flat surface and x 0 is the corresponding point of the rigid sphere which is coming into contact. The simplest method to enforce the non-penetration constraints in Supplementary Eq. (31) consists on transforming the constrained minimization problem into an unconstrained minimization problem by the penalty method, which adds an extra surface energy term over the contact boundary, penalizing deviation from the conditions in Supplementary Eq. (31), to the total energy in Supplementary Eq. (1) as [12] H t = H + where β is the penalty parameter and • is the Macaulay brackets which returns the argument itself if the argument is positive, otherwise it returns zero (no contact). The rotational symmetry of the spherical tip-sample contact allows us to employ a two-dimensional axisymmetric model, as depicted in Fig. 2(b) of the manuscript. The contact force F is then obtained as where a is the contact radius, i.e. the distance of the farthest point from the symmetry axis where g n = 0. The discrete form of the total energy H t in Supplementary Eq. (33) in then obtained as where H is given in Supplementary Eq. (15) and F a is the contact force given by where R c is the contact radius and t is the imposed mechanical traction due to contact, obtained as t = −β g n n. The right-hand side contribution in Supplementary Eq. (30) then changes to f U = tp a r and the total contact force F is obtained as F = a F a .
To verify this contact model, we perform simulations for different given vertical positions of the tip and compute the magnitude of the resulting contact force. These results are compared with the classical Hertzian contact model which provides an analytical solution for the relationship between the indentation depth d, the tip radius R, the contact radius R c , and the applied load F as [14] where E * is the apparent Young's modulus. Assuming that the Young's modulus of the tip is noticeably higher than that of the sample, E * is approximated as E * = E/(1 − ν 2 ), where E and ν are the Young's modulus and Poisson's ratio of the sample, respectively. Supplementary Figure 1 shows the results considering elastic parameters E = 230 GPa and ν = 0.3, a tip radius of R = 100 nm, a domain length of L = 2 µm, a domain height of H d = 1 µm and a value of 5 × 10 17 Nm −3 for the penalty parameter β. A good agreement is observed between the simulation and the Hertzian contact model. To perform the simulations of non-piezoelectric SrTiO 3 (STO), we consider the material parameters given in Ref. [11] and presented in Supplementary Table 1. The length-scale l 1 has been computed as 10 nm using a formulation in Ref. [5]. Simulations are performed considering a tip radius of R = 100 nm, a domain length of L = 2 µm, a domain height of H d = 1 µm and a tip voltage V = -5 V. We first progressively indent (increase the vertical position of the tip), until reaching the desired force F in the absence of the tip voltage (V = 0). Then, in the second simulation, we apply the tip voltage V (see Eqs. (2) and (3) in the manuscript) and change the position of the tip until we obtain the same calculated force as in the first simulation. This change corresponds to the value of ∆u z . The effective piezoelectric coefficient is then calculated as d eff 33 = −∆u z /V . We perform these simulations for different values of the applied force F .

Supplementary Note 3. Electrochemical strain
Another effect that can yield an electromechanical response is electrochemical strain [15]. The intense electric fields generated by the tip can sometimes attract ionic defects, such as oxygen vacancies, towards the contact region, thereby inducing a local swelling of the lattice. This effect is proportional to ionic conductivity and it should therefore be expected to increase with temperature [16]. In contrast, the converse flexoelectric coefficient is proportional to the dielectric constant, which in SrTiO 3 (STO) decreases with increasing temperature [17]. We have looked at the temperature dependence of the STO piezoresponse to distinguish both mechanisms. As can be seen in Supplementary Figure 2, the piezoresponse decreases as temperature increases, consistent with a flexoelectric origin. The ionic migration contribution to the electromechanical response was neglected due to the lack of hysteresis in the cycles of the remanent PFM response as a function of the applied V dc field (Supplementary Figure 3). The absence of oxygen vacancies that could lead to ionic contribution to electromechanical response [18,19] was also probed by XPS (Supplementary Figure 4).

Supplementary Note 4. XPS analysis of the SrTiO 3 Single Crystal
STO single crystals were analysed by XPS at room temperature and different water gas pressures. XPS experiments were performed at the Near Ambient Pressure XPS end station of the CIRCE beamline (BL-24) at the ALBA synchrotron radiation facility. The photon energy range at CIRCE is 100 eV-2000 eV. The end station is equipped with a Phoibos 150 NAP electron energy analyzer. The analyzer is provided with four differentially pumped stages connected by small apertures. A set of electrostatic lenses focuses the electrons though the apertures to maximize transmission. Such set-up allows keeping the detector in Ultra High Vacuum while the sample is at a maximum pressure of 20 mbar [20,21]. The total electron energy resolution in the experimental conditions used for the high-resolution spectra was better than 0.3 eV. Further details about the system can be found elsewhere [22]. Surface oxygen vacancies can be monitored by checking the presence of shoulders in the low BE side of the Ti 4+ 2p 3/2 peak (corresponding to Ti 4+ ) due to reduced Ti species, i.e. Ti 3+ (Supplementary Figure 4b).
In our experiments the RH conditions change from vacuum to 10% and surface oxygen vacancies cannot be observed in the Ti 2p 3/2 peak in any case. Since our samples were all previously exposed to atmospheric conditions and we did not perform any annealing treatment to create vacancies, we conclude that all oxygen vacancies close to the surface area within a penetration depth of several nm have been compensated before the experiments began. No differences in the shape of Ti 2p 3/2 neither the Sr 3d peaks are observed before and after water exposure.

Supplementary Note 5. Topographical damage after electromechanical measurements
In order to account for possible electrochemical effects in the shape of permanent damage in the material produced by the application of an ac voltage by an Atomic Force Microscope (AFM) conductive tip, the topography of the sample was monitored before and after the experiments. Supplementary Figure 5 shows that there is no relevant difference in sample topography before and after measuring the electromechanical response with a diamond coated tip.
Supplementary Note 6. Electromechanical response of TiO 2 The effective piezoelectric coefficient was measured for TiO 2 in the same conditions as for SrTiO 3 . TiO2 shows a rutile structure as compared to the cubic STO structure and is known to also show a reasonable flexoelectric response, even its flexoelectric coefficient µ eff 13 =1.7 nCm −1 is reduced by a factor of 0.68 with respect of that of SrTiO 3 . This lower flexoelectric performance leads to a reduction of the measured effective piezoelectric coefficient due to converse flexoelectric effects as expected. In this case, change of contact radius was achieved by using AFM tips with nominally different radius, and applying forces of the order of F ≈ 100 nN for each of the measurements. In every measurement, the tip force was calibrated, tip radius was determined from SEM images after the measurements, and effective contact radius was calculated a posteriori using the Hertzian contact model. Measurements were done for each tip alternatively on each sample (Supplementary Figure 6).    Figure 6. Study of electromechanical response of TiO 2 as compared to SrTiO 3 under the use of tips with different radius. The radius of the tips was determined after the nanoscale electromechanical measurements from scanning electron microscopy (SEM) images of the used tips. a SEM image of a Nanosensor EFM tip of medium stiffness (k), showing the sharpest radius. b SEM image of Nanosensor PtSi tip also with a medium stiffness k. c SEM image of a Nanosensor NCH CDT diamond coated tip, with high stiffness and big tip radius. d Effective electromechanical response of TiO 2 (full diamonds) and SrTiO 3 (empty circles) as a function of the contact radius a. In this case, the different experimental contact radius were obtained by varying the tip radius instead of varying the force. The experimental contact radius a is calculated from the tip radius and the applied force following the Hertz contact model as stated in Supplementary Eq. (37). Under the same applied force, the contact radius will scale with the tip radius. The black points corresponding to the smaller contact radius are obtained from measurements performed with the sharper tip as shown in a, red points are obtained from measurements with the tip shown in b and finally green points correspond to measurements with the tip shown in c. The error bars correspond to the error of the linear fitting of the experimental data, which correlates the measured electromechanical amplitude of oscillation ∆h with the V ac applied voltage.