Preisach modelling of nonlinear response in electrically biased lead zirconate titanate-based piezoceramics

Appl Phys A (2013) 112:1081–1088 DOI 10.1007/s00339-012-7492-3 Preisach modelling of nonlinear response in electrically biased lead zirconate titanate-based piezoceramics Diego A. Ochoa ·Rafel Pérez ·Jose E. García Received: 21 September 2012 / Accepted: 3 December 2012 / Published online: 11 December 2012 © Springer-Verlag Berlin Heidelberg 2012 Abstract The alteration of the high-field electrical permittivity (nonlinear response) of PZT-based ceramics when an electrical bias field is applied is reported in this work. Large differences are observed between soft and hard PZT behaviours. While in soft PZT a bias field does not modify the nonlinear behaviour, a notable dependence is verified in hard PZT. The Preisach model is satisfactorily used to describe experimental results. A distribution function containing the first terms of the Maclaurin development series of a function composed by two Gaussian-like functions of different amplitudes is proposed. The model gives a satisfactory explanation for the fact that the permittivity depends not only on the amplitude of the applied electric field, but also on the bias field, both for soft and hard ceramics and for poled or unpoled samples. 1 Introduction The response of piezoceramics to an electrical or mechanical excitation can be separated into two contributions, known as intrinsic and extrinsic. The intrinsic contribution is associated with the linear lattice distortion of the unit cell, whereas the extrinsic contribution is mostly associated with domain-wall motion and its interaction with structural defects [1–3]. The extrinsic contribution is considered responsible for nonlinearity in the macroscopic properties of piezoceramics [4–6]. The increase in permittivity, piezoelectric coefficient, and mechanical losses are commonly reported D.A. Ochoa ·R. Pérez ·J.E. García () Department of Applied Physics, Universitat Politècnica de Catalunya, 08034 Barcelona, Spain e-mail: [email protected] Fax: +34-934-016090 nonlinear effects. Therefore, control and modelling of the extrinsic response is an important issue in piezoceramics. Lead zirconate titanate (PZT) piezoceramics are widely used in sensor and actuator applications due to their high electromechanical coupling [1]. The addition of dopants severely modifies the PZT response, thereby improving their properties for specific applications [1,7]. Dopants are normally classified as donors (softeners) or acceptors (hardeners), depending on the valence of the ion substituted. The substitution of Zr4+or Ti4+by trivalent cations (i.e. Fe3+) produces oxygen vacancies forming the so-called complex defects, which present dipolar characteristics [8–10]. These defects operate as pinning centres, decreasing the domainwall motion [8,9]. Moreover, the substitution of Zr4+or Ti4+by pentavalent cations (i.e. Nb5+) produces lead vacancies and reduces oxygen vacancies, thus improving domainwall mobility [1,8]. As a consequence of compositional modification, many differences between hard and soft PZT properties are observed. For instance, compared with hard PZT properties, soft PZT exhibits higher dielectric and piezoelectric constants but higher losses and lower quality factor. Furthermore, when properties of these materials are compared, not only their values are different but also their behaviour [11]. The hysteresis loop is a well-known characteristic of ferroelectric materials that is unequal in soft and hard PZT, as is shown in Fig. 1. As can be observed, the shapes of the hysteresis curves are notably different. While soft PZT shows a ‘normal’ hysteresis loop, hard PZT shows a pinched one. The differences between macroscopic properties of both PZT-based materials are assumed to be due mainly to the presence of complex defects, which in polarized hard PZT present a dipolar moment oriented in the polarization direction [8–10]. Those complex defects are in charge of the well-known domain wall pinning effect [8–11]. 1082 D.A. Ochoa et al. Fig. 1 Typical hard and soft ferroelectric hysteresis loops, at room temperature. Curves are plotted from real hard and soft PZT data but different polarization scales have been used (actually, the maximum polarization of hard PZT is four times lower than soft PZT) In many applications, PZT ceramics are often subject to a bias electric field that influences domain-wall motion by changing the electromechanical properties of these materials [12–14]. As a consequence, nonlinear dielectric behaviour (i.e. permittivity as a function of applied ac electric field amplitude) is also modified. Understanding and modelling the influence of a bias electric field is an important question that must be solved in order to improve piezoelectric devices. In this work, the Preisach model is used in order to explain the extrinsic response behaviour of both soft and hard PZT piezoceramics when they are subjected to a bias electric field. The aim is to find the distribution function that correctly describes experimental results when the model, which considers only irreversible mechanisms of domainwall movement, is applied without any modification. 2 Preisach model The Preisach model is a classical method to describe the nonlinear behaviour in hysteretic systems and has been used satisfactorily in ferroelectrics [15–18]. According to this model, a hysteretic system contains a collection of simple bistable units (Fig. 2a) with both states contributing independently to the total response by the same amount ±P0. These units are characterized by an internal field Ei, and a coercive field Ec, so that the set of bistable units exhibits a statistical distribution f(E i,Ec)that depends on the parameters −∞<E i<∞and 0 <E c<∞. The response to an applied external field depends on both the amount of bistable units switched and the direction of applied field change. For increasing fields (Einc), when Einc >E i+Ec, the bistable units in the negative state of polarization (−P0)switch to the positive state (+P0), each one contributing 2P0to the total response (Fig. 2b). For decreasing fields, when Edec < Ei−Ec, each unit in the state +P0contributing −2P0to the total response switches to the state −P0(Fig. 2b). So, if the applied field is an alternating field E(t) =Eb+E0cosωt, where Ebis the applied bias electric field and E0is the amplitude of the applied alternating electric field, then according to whether the applied field increases (Fig. 2c) or decreases (Fig. 2d) the bistable units inside the grey triangular regions of the Preisach plane contributing to the total response switch to their opposite state [19]. The branches of the hysteresis loop can be obtained from the field-dependent integrals of the distribution function over the grey triangular regions given in Fig. 2c and d: D−(E) =DNL(E0,Eb) 2 −2P0[(Eb+E0−E)/2] 0 ×[Eb+E0−Ec] E+Ec f(E i,Ec)dEidEc,(1) D+(E) =−DNL(E0,Eb) 2 +2P0[(−Eb+E0+E)/2] 0 ×[E−Ec] Eb−E0+Ec f(E i,Ec)dEidEc;(2) which is DNL(E0,Eb)=2P0E0 0[Eb+E0−Ec] [Eb−E0+Ec]f(E i,Ec)dEidEc. (3) Equation (1) corresponds to the descending branch of the loop (Fig. 2c), where the bistable units change from a positive to a negative state of polarization. Equation (2)is the ascending branch (Fig. 2d), where the bistable units change from a negative to a positive state. The quantity DNL(E0,Eb), which represents the total extrinsic (nonlinear) contribution, depends on the number of switched units inside the triangular region of the Preisach plane. The term ±DNL/2 is added to the integrals in Eqs. (1) and (2) because the extreme values of each loop branch must coincide. In addition, it is also useful to centre the hysteresis loops. Bearing in mind the time dependence of the electric field E=E(ωt), we can use the phase angle ϕ=ωt to obtain the instantaneous electric displacement Di NL(ϕ) from Eqs. (1) and (2): Di NL(ϕ) =⎧ ⎨ ⎩ D−(E(ϕ)) when dE dt<0, D+(E(ϕ)) when dE dt>0. (4) Preisach modelling of nonlinear response in electrically biased lead zirconate titanate-based piezoceramics 1083 Fig. 2 Schematic representation of the bistable unit and the Preisach plane: (a) the bistable unit showing its characteristic parameters; (b) conditions fulfilled for the up or down switching of the bistable unit; (c) bistable elements that are switched when the field goes down; (d) bistable elements that are switched when the field goes up Taking the first term of the Fourier transform of Di NL(ϕ), and considering that the nonlinear instantaneous electric displacement can be expressed as Di NL(ϕ) =(ε−iε)× ε0E0eiϕ, it is possible to obtain the variation of the real and imaginary parts of the dielectric permittivity from the components in phase and quadrature phase with the applied electric field [20]: ε=1 ε0E01 π2π 0Di NL(ϕ) cos(ϕ)d(ϕ),(5) ε =1 ε0E01 π2π 0Di NL(ϕ) sin(ϕ) d(ϕ).(6) Nonlinear response depends on the specific distribution function associated with the material. In soft PZT, at subswitching regime, it can be described with a flat uniform distribution function f(E i,Ec)=g0, which results in the Rayleigh model equations [19]. This model is widely accepted in the literature for describing soft PZT behaviour [21–23]. Hard PZT, on the other hand, shows a different behaviour. The main observable characteristic is that the major hysteresis loop appears pinched (Fig. 1), which leads to a drastic decrease in the remanent polarization. Hard PZT behaviour has been described in the literature by a fourth-order polynomial distribution function, symmetric in terms of Ei, combined with a cut off at high coercive fields [15], although the same reference discusses that a more realistic distribution may be obtained by two superimposed Gaussians, symmetrical with respect to Ei. More general distribution functions, based on the shape of the hysteresis loop and its evolution under external actions (e.g. mechanical stress), have been proposed by other authors [24–27]. In any case, when the applied fields are far away from the coercive field (subswitching regime, i.e. low-field zone in Fig. 1), the distribution function may be expanded in a Maclaurin series [19]: f(E i,Ec)=g0+g2E2 i+gC 1Ec+gC 2E2 c.(7) The Preisach model allows further analysis of material response. So, the aim of this work is to study nonlinear (extrinsic) response of soft and hard PZT under the influence of a bias electric field. An attempt is made to achieve a distribution function able to describe the nonlinear behaviour of both soft and hard PZT, while at the same time describing the influence of a bias electric field on the nonlinear response. In a broader perspective, nonlinear response also depends on frequency [22] and temperature [28],soafrequency of 1 kHz and room temperature are taken as fixed parameters in this work. The results obtained could be extended to other piezoelectric systems. 3 Experimental procedure Pb(Zr1−xTix)O3(x=0.40) samples were prepared by solid-state reaction. Precursors were mixed in a ball mill for 2 h, dried at 80 °C, and calcined at 850 °C for 3.5 h. Nbdoped (soft PZT, hereafter labelled PNZT) and Fe-doped (hard PZT, hereafter labelled PFZT) powders were prepared by the addition of 1 wt.% of Nb2O5or Fe2O3to the green mixture before homogenization by a second ball milling. Powders were uniaxially pressed at 100 MPa into pellets of 20 mm in diameter and about 5 mm in length, and then isostatically pressed at 120 MPa. Sintering was carried out in a saturated PbO atmosphere. Samples were cut into discs of 15–16 mm in diameter and 0.8–0.9 mm in thickness, avoiding areas of inhomogeneous concentration of Pb. X-ray diffraction analysis revealed the expected rhombohedral crystallographic phase. Neither spurious phases nor contaminants were detected. After polishing, organic solvents were removed by heat treatment at 600 °C for 30 min, which also released the mechanical stress induced during polishing. Silver electrodes were painted on both sides of the discs. Nonlinear dielectric measurements were performed at room temperature by applying a 1-kHz driving ac electric field (sub-switching regime) for different bias electric fields. 1084 D.A. Ochoa et al. The electrical response was measured by using a capacitance comparator bridge, especially designed for this type of characterization [29]. The instantaneous nonlinear current was measured as a function of the instantaneous applied electric field. The electrical displacement D(E) was obtained by numerical integration of the current for each value of the voltage amplitude. The values of the increments in the dielectric permittivity and losses were computed from the components of the D(t) response in phase and in quadrature of phase with the electric field E(t), as reported elsewhere [29]. Steps were taken to ensure that the absolute value of the total applied field E(t) =Eb+E0cosωt was maintained in the low-field zone (lower than the coercive field), thereby preventing changes in the domain structure. Additionally, in order to verify that no permanent changes occurred in the domain structure (preserving the same domain configuration), weak-field linear properties were measured before and after each experiment by using an impedance analyzer. No significant differences in the linear properties were detected. 4 Experimental results Dielectric permittivity of PNZT versus the amplitude of the applied ac field with bias applied field Ebas parameter is shown in Fig. 3. It is observed that nonlinear response, understood as the variation of the permittivity with the applied ac field, remains constant for any value of Eb. The variation of the permittivity with the amplitude of the applied ac field is ε≈900 for every Eb. Furthermore, the variation of the permittivity is linear with the amplitude of the alternating applied field E0, which is a typical soft PZT behaviour. In this case, the direction of the bias applied field has no influence on the nonlinear response, irrespective of whether the Fig. 3 Dielectric constant εof PNZT as a function of the amplitude of the ac applied electric field E0, for different bias electric fields Eb.The sign of Ebindicates that it is antiparallel to the sample polarization. The values of the εversus E0slope are computed and plotted in the inset for Ebboth parallel and antiparallel to the sample polarization Ebfield is parallel or antiparallel to the sample polarization. Dielectric losses (not shown) have a similar behaviour to the dielectric constant and are analyzed below from the Preisach model. For each Eb, the slope of the εversus E0curve is computed from Fig. 2data and the values are plotted in the inset. As expected, similar slope values are obtained, showing a similar nonlinear response magnitude. The slope is described in the Rayleigh model, and for soft PZT the reported value in the literature is α=2400 (MV/m)−1[29]. The nonlinear response behaviour of PFZT (Fig. 4)is notably different. The first difference observed is that the variation of the permittivity is not linear with the amplitude of the applied field, as has been reported in previous work [30]. Another substantial difference is the strong dependence of the dielectric permittivity on Ebfor both the parallel (Fig. 4a) and antiparallel (Fig. 4b) directions to the sample polarization. The variation of the permittivity ε with the applied field changes from 15 to 35 depending on the bias applied. A more accurate analysis can be made from Fig. 4c and d. Irrespective of whether the Ebfield is parallel (Fig. 4c) or antiparallel (Fig. 4d) to the sample polarization, the nonlinear response increases with the increase in |Eb|, this rise being slightly higher when the Ebdirection is parallel to the sample polarization. Dielectric losses have the same behaviour as the dielectric constant and, as in PNZT, they are analyzed below from the Preisach model. 5 Discussion The simplest development of the Preisach model corresponds to a flat uniform distribution function f(E i,Ec)= g0. By substituting this distribution function into Eq. (3), and assuming that the applied electric field has the form E(t) =Eb+E0cosφ, it can be seen that the nonlinear contribution is proportional to the area of the Preisach plane swept by the field: DNL(E0)=2P0g0E2 0.(8) The nonlinear contribution obtained by the Preisach model for a flat distribution does not depend on the applied bias field Eb, that is, according to the results obtained for PNZT (Fig. 3). This can be explained by analyzing the area swept over the Preisach plane. For Eb=0, the nonlinear contribution is proportional to the area (E2 0) of the triangle in Fig. 5a. By applying a bias field Eb, this triangle moves up and down along the Eiaxis even though both the area and the number of bistable units swept in a cycle will remain constant, thus providing the same contribution to the nonlinear response. Continuing with the graphical analysis, one may observe that an increment of E0increases the swept area, and consequently the nonlinear contribution enhances (Eq. (8)). Preisach modelling of nonlinear response in electrically biased lead zirconate titanate-based piezoceramics 1085 Fig. 4 Dielectric constant, ε, and its increment, ε,asa function of E0for different values of bias electric field Eb, for Ebboth parallel (a)and(c) and antiparallel (b)and(d)to the sample polarization Fig. 5 Representation of a flat uniform distribution function in the Preisach plane, for: (a)Eb=0and(b)Eb=0. The area of the large triangle is proportional to the number of bistable elements swept in each cycle In order to analyze the dependence of the response on the amplitude E0, and also to compare this dependence with the results, it is necessary to analyze the dependence on E0 of both the dielectric constant and the losses. From Eqs. (5) and (6), for a flat distribution we obtain ε=P0g0E0,(9) ε =4 3πP0g0E0.(10) A linear dependence on E0is obtained for both εand ε, and the ratio ε/ε=4/3π=0.42, as foreseen by the Rayleigh model [30]. Additionally, there is no theoretical dependence on the applied bias field Eb. The branches of the sub-switching hysteresis loop can be obtained from Eqs. (1) and (2). The result for a flat distribution, Di NL(t) =P0g0E2 0cosϕ±1 2P0g0E2 01−cos2ϕ,(11) coincides with the expression given by the Rayleigh model by identifying the Rayleigh coefficient αas α=P0g0[30]. Therefore, the slope P0g0in Eq. (9) is the parameter α, which denotes the density of bistable units in the Preisach plane. A description of the more complex response of PFZT requires a deep analysis of the distribution function proposed in the literature (Eq. (7)). Here, the aim is to describe the observed effect of a bias field Ebon the PFZT nonlinear response (Fig. 4). According to the Preisach model, a bias field only implies a displacement of the sweeping zone up or down along the Preisach plane, as can be seen in Fig. 5b for Eb>0. This implies that the polynomial terms of the distribution function that do not fulfil Eimake no contribution to the change in the nonlinear response for the different Ebapplied. It is therefore possible to dispense with Ec terms, since they have no influence on the variation of the nonlinear response for different Eb. The distribution func- 1086 D.A. Ochoa et al. Fig. 6 Representation of a (a) symmetric and (b) non-symmetric quadratic distribution function in the Preisach plane. Two Gaussianlike functions (left) can describe the distribution function. Details of the zone corresponding to a sub-switching regime for (a) non-poled and (b) poled samples are shown on the right tion (Eq. (7)) is then reduced to f(E i,Ec)=g0+g2E2 i.(12) We can assume that Eq. (12) contains the first terms of a series development of a distribution function consisting of two Gaussian-like functions, symmetric with respect to Ei (Fig. 6a, left). The projection of the distribution function, givenbyEq.(12), in the Preisach plane is represented at the centre of Fig. 6a. The shaded zone, which corresponds to the Gaussian peaks, represents the increase in the density of bistable units with higher |Ei|due to the presence of complex defects with polar character [15]. Equation (12) may describe nonlinear response for a non-poled sample, where any direction of Ebproduces the same material response. However, experimental results (Fig. 4) show some differences, depending on whether Ebis applied parallel or antiparallel to the sample polarization. A more realistic description of a poled sample must take into account the asymmetry of the nonlinear response for the different orientations of Eb. To this end, it is necessary to break the symmetry that appears in the distribution function given by Eq. (12). A linear term in Eicould be added in this case, resulting in a new distribution function f(E i,Ec)=g0+g1Ei+g2E2 i.(13) This new function contains the first terms of the Maclaurin development series of a distribution function consisting of two Gaussian functions of different amplitudes, as represented on the left of Fig. 6b. In poled samples, complex defects tend to be oriented parallel to the direction of polarization [31], breaking the symmetry and giving preference to one of the directions, as can be seen in the shaded zones at the centre of Fig. 6b. It is necessary to point out that, although we do not consider the terms that depend on Ec, a more realistic distribution function for the description of the losses should contain terms in Ec. However, their inclusion in this work is regarded as irrelevant, because we focus our attention on the dependence with Ebof the nonlinear dielectric response. For poled PFZT, the total nonlinear contribution DNL is obtained by replacing the distribution function (Eq. (13)) in Eq. (3): DNL(E0,Eb)=2P0g0E2 01+g1 g0Eb+g2 g0E2 0 6+E2 b. (14) Moreover, the increments of the dielectric constant and losses are obtained from Eqs. (5) and (6), as follows: ε=P0g0E01+g1 g0Eb+g2 g07 48E2 0+E2 b,(15) ε =4 3πP0g0E01+g1 g0Eb+g2 g01 10E2 0+E2 b.(16) As can be observed, neither εnor ε depend linearly on the field amplitude E0. This result, which is far removed from the Rayleigh behaviour, is in accordance with the experimental results in PFZT and has been reported elsewhere [6]. Apart from the dependence on E0, a theoretical dependence on the bias field Ebis observed, which does not appear in the results obtained from a flat distribution. This dependence of the nonlinear response on Ebis consistent with the experimental results for PFZT (Fig. 4). The dependence of nonlinear dielectric response on Eb can be explained from an analysis of the distribution function (Eq. (13)) in the Preisach plane, as indicated in the zoom of Fig. 6(on the right) for sub-switching conditions. When a field Ebis applied, the triangular sweeping area is displaced along the Eiaxis. Although in PNZT this displacement does not change the number of switchable bistable elements, in PFZT we find more elements in this area when it is displaced towards the maximum concentration area (shaded area in Fig. 6). So, nonlinearity increases when the absolute value of Ebincreases. The differences observed in nonlinear response depending on the Ebdirection can be explained taking into account the asymmetry of the density of bistable units. Preisach modelling of nonlinear response in electrically biased lead zirconate titanate-based piezoceramics 1087 Fig. 7 Dependence of the coefficient A, defined in Eq. (18), on the bias electric field Eb. The polynomial fit of the coefficient Aversus Eb is also displayed As can be observed in Fig. 3, in PNZT the dielectric constant increases linearly with E0. However, Fig. 4shows that this dependence is not linear in PFZT, a fact that corresponds to the inclusion of the term E3 0in Eqs. (15) and (16). Notice that when the maximum applied electric field (Emax)is sufficiently low to verify g0g1Emax and g0g2E2 max, Eqs. (14), (15), and (16) are reduced to the form of the expressions (3), (5), and (6). Thus, the Rayleigh model can still be applied in this condition. In order to study the nonlinear response quantitatively, the coefficients that appear in the distribution function are obtained from the experimental data. For easy visualization of the fit, the Eq. (15) is rearranged as follows: ε=P0g0+P0g1Eb+P0g2E2 bE0+7 48P0g2E3 0.(17) The above expression can be rewritten as ε=AE0+ BE3 0, where A(Eb)=P0g0+P0g1Eb+P0g2E2 b(18) and B=7 48P0g2.(19) The value of B, which is independent of Eb, is obtained by fitting the experimental results when Eb=0, which yields B=(145 ±10)(MV/m)−3. Figure 7shows the dependence of Aon the bias electric field Eb.The values of Aare obtained by making a fit of the plot ε(E0)for each Eb.P0g0=(26.4±0.5)(MV/m)−1 and P0g2=(1110 ±30)(MV/m)−3are obtained from the fit of the function A(Eb). In addition, the value P0g1= (−33 ±4)(MV/m)−2is also obtained. From the Bvalue, it is possible to evaluate the P0g2value, which is P0g2= (1000 ±70)(MV/m)−3, and is very similar to the term obtained from the fit of A(Eb). The values of the P0g2coefficient, which are determined in two different ways, are similar and show the consistency of the distribution function taken into account for the description of the nonlinear behaviour in PFZT. The cubic dependence on E0that appears in the equation of the dielectric permittivity due to the term g2E2 ican also be obtained by the addition of a term gC 2E2 cinto the distribution function, thereby modifying the value of the coefficient B, as follows: B=7 48P0g2+gC 2.(20) The difference between the terms in g2and gC 2is that the term in g2also modifies the linear term in E0by adding a quadratic term in the equation for A(Eb). The fact that the two values obtained for P0g2are similar implies that the contribution of the term gC 2E2 cto Bis negligible with respect to the contribution of the term g2E2 i, which may justify the elimination of this term in the distribution function. Taking into account that the term P0g0is related to the parameter αof the Rayleigh model, and also that it is related to the irreversible motion of the domain walls, we suspect that the lowering of the variation of the permittivity in PFZT it is due to a decrease in the irreversible motion of the domain walls. This makes the reversible contribution more apparent, and it acquires greater importance in the total nonlinear response, although it is not considered in the present analysis. It is reasonable to assume that, in both soft and hard materials, there exist both reversible and irreversible contributions of domain-wall motion to the nonlinear response. While in the case of soft PZT the irreversible contribution is so large that it completely masks the reversible contribution, this phenomenon it is not so clear in the case of hard PZT. 6 Conclusions The dielectric behaviour of PZT ceramics when an electric field is applied strongly depends on the defects created by dopants. In soft PZT, there exists a linear dependence of the permittivity on the amplitude of the alternating field, while nonlinear response remains constant for any value of bias electric field. Otherwise, the dependence on the amplitude of the field is not linear in hard PZT and is notably dependent on the applied bias field. The Preisach model is satisfactorily used to describe the observed behaviour by taking an adequate distribution function of density of the bistable units. For the sub-switching regime, it is possible to consistently obtain a distribution function in the form of a series development. The model fits with the experimental measurements, thereby explaining 1088 D.A. Ochoa et al. the dependence on both the field amplitude, E0, and the bias field, Eb. A flat distribution, which explains the behaviour of soft PZT, and which is coherent with the Rayleigh model, can be seen as an approximation of a more complex distribution function. It is therefore possible to describe the behaviour of both soft and hard PZT by using the same form of the distribution function, changing exclusively the coefficient relations. Consequently, nonlinear response appears to be due to both reversible and irreversible contributions of domain-wall motion in both materials, although the irreversible contribution is very large in soft PZT. In fact, the term P0g0, which is related to the parameter αof the Rayleigh model, is two orders of magnitude greater in soft PZT than in hard PZT. It is further shown that the poling of hard PZT breaks the symmetry of the density distribution, so that some differences are observed depending on the direction of the bias field. Acknowledgements This work is supported by project MAT201021088-C03-02 of the Spanish Government. The authors wish to thank Prof. Jose A. Eiras, at the Universidade Federal de São Carlos, Brazil, for providing the samples. References 1. B. Jaffe, W.R. Cook, H. Jaffe, Piezoelectric Ceramics (Academic Press, New York, 1971) 2. Q.M. Zhang, W.Y. Pan, S.J. Jang, L.E. Cross, J. Appl. Phys. 64, 6445 (1988) 3. Q.M. Zhang, H. Wang, N. Kim, L.E. Cross, J. Appl. Phys. 75, 454 (1994) 4. B. Lewis, Proc. Phys. Soc. Lond. 73, 17 (1959) 5. D. Berlincourt, H.H. Krüger, J. Appl. Phys. 30, 1804 (1959) 6. J.E. Garcia, R. Perez, D.A. Ochoa, A. Albareda, M.H. Lente, J.A. Eiras, J. Appl. Phys. 103, 054108 (2008) 7. J.C. Burfoot, G.W. Taylor, Polar Dielectrics and Their Applications (University of California Press, Berkeley, 1979) 8. S. Takahashi, Ferroelectrics 41, 277 (1982) 9. W.L. Warren, G.E. Pike, K. Vanheusden, D. Dimos, B.A. Tuttle, J. Robertson, J. Appl. Phys. 79, 9250 (1996) 10. M.H. Lente, J.A. Eiras, J. Appl. Phys. 94, 2112 (2002) 11. D.A. Ochoa, J.E. García, R. Pérez, A. Albareda, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 55, 2732 (2008) 12. P.M. Chaplya, G.P. Carman, J. Appl. Phys. 92, 1504 (2002) 13. A. Pramanick, A.D. Prewitt, M.A. Cottrell, W. Lee, A.J. Studer, K. An, C.R. Hubbard, J.L. Jones, Appl. Phys. A 99, 557 (2010) 14. M. Marsilius, T. Granzow, J.L. Jones, J. Mater. Res. 26, 1126 (2011) 15. G. Robert, D. Damjanovic, N. Setter, Appl. Phys. Lett. 77, 4413 (2000) 16. S.A. Turik, L.A. Reznitchenko, A.N. Rybjanets, S.I. Dudkina, A.V. Turik, A.A. Yesis, J. Appl. Phys. 97, 64102 (2005) 17. G. Robert, D. Damjanovic, N. Setter, A.V. Turik, J. Appl. Phys. 89, 5067 (2001) 18. G. Robert, D. Damjanovic, N. Setter, J. Appl. Phys. 90, 2459 (2001) 19. G. Bertotti, I. Mayergoyz, The Science of Hysteresis,vol.3(Elsevier, Oxford, 2006) 20. H. Frohlich, Theory of Dielectrics (Clarendon, Oxford, 1958) 21. A.V. Turik, Sov. Phys., Solid State 5, 885 (1963) 22. D.V. Taylor, D. Damjanovic, J. Appl. Phys. 82, 1973 (1997) 23. D.A. Hall, J. Mater. Sci. 36, 4575 (2001) 24. A. Sutor, S.J. Rupitsch, R. Lerch, Appl. Phys. A 100, 425 (2010) 25. F. Wolf, A. Sutor, S.J. Rupitsch, R. Lerch, Sens. Actuators A 172, 245 (2011) 26. F. Wolf, A. Sutor, S.J. Rupitsch, R. Lerch, Sens. Actuators A 186, 223 (2012) 27. S.J. Rupitsch, F. Wolf, A. Sutor, R. Lerch, Acta Mech. 223, 1809 (2012) 28. J.E. Garcia, D.A. Ochoa, V. Gomis, J.A. Eiras, R. Perez, J. Appl. Phys. 112, 014113 (2012) 29. J.E. Garcia, R. Perez, A. Albareda, J. Phys. D, Appl. Phys. 34, 3279 (2001) 30. J.E. Garcia, R. Perez, A. Albareda, J. Phys. Condens. Matter 17, 7143 (2005) 31. W.L. Warren, D. Dimos, G.E. Pike, K. Vanheusden, R. Armes, Appl. Phys. Lett. 67, 1689 (1995)