High-response piezoelectricity modeled quantitatively near a phase boundary

High-response piezoelectric materials belong to the perovskite solid solution family Pb((Zn or Mg)_1∕3Nb_2∕3)_1–xTi_xO₃ with a complex phase diagram as a function of composition x and T,^1–3 as illustrated in Fig. 1(a). At compositions lying close to a structural instability–termed the morphotropic phase boundary (MPB)–the electromechanical responses are maximal, far above the responses of conventional piezoelectric ceramics. These materials are attractive for technological applications, such as medical diagnostics, sonar, the Piezoelectronic transistor (PET) family of devices,⁴ among others.

The most responsive materials are located in composition close to the MPB, on the rhombohedral side¹ or just over it into the orthorhombic/monoclinic mixed phase labeled M⁵ in Fig. 1(a) (see Ref. 6 for a thorough discussion). Achieving an understanding of this composition range below the Curie temperature is a key objective to enable advances in materials through fundamental study.

Experiment has discovered, in the near-Curie temperature region, a line of critical points for the paraelectric to ferroelectric transformations in the composition-temperature-electric field phase diagram of PMN-PT.⁷ In the supercritical region, a continuous passage in the electric field-temperature plane from cubic to tetragonal crystal symmetry occurs, accompanied by a maximum in piezo coefficient as the Widom line is crossed. These results show the facile rotation of the polarization order parameter associated with near-criticality.

A corresponding extension of the facile rotation phenomenon occurs in the subcritical region. Here, the magnitude of the polarization vector is relatively well-defined, but the facile rotation of the polarization for compositions close to the MPB is the key phenomenon. This rotator mechanism, activated with a composition rather than the electric field (see Fig. 1), was proposed early on in Refs. 8 and 9.

The two main phases below the ferroelectric Curie point T_c are, on the low-PT side the rhombohedral phase R, with (1, 1, 1) electric polarization, and on the high-PT side, the tetragonal phase T with (0, 0, 1) polarization, separated by an intrusive monoclinic/orthorhombic^10,11 or mixed phase¹² (see Fig. 1(a)). Near the MPB polarization is labile, simplified here as lying between the set of (1, 1, 1) and (0, 0, 1) directions, allowing an unusually large response in polarization direction to an external perturbation such as electric field or stress.¹³ The result on the R side is a giant value of permittivity $ϵ33T$⁠, and, since the crystal morphology is coupled to electric polarization, compliance $s33E$ and piezo coefficient d₃₃ (superscripts T and E represent constant stress and electric field, respectively),¹³ although the largest $ϵ33T$ values are in the M phase.

A first-principles approach to piezoelectrics has been developed¹⁴ in systems with crystalline order, but in the case of the relaxor piezoelectrics, it is challenged by the local perturbation presented by frozen-in B-site cation disorder. A phenomenological method^1,9,15 based on the rotating crystal reference frames provides valuable insights on the orientation dependence of crystal responses to external perturbations but is not predictive (being based on linear combinations of experimentally measured responses). A molecular dynamics approach,¹⁶ finding an easy rotation of the polarization direction, has reproduced the phase diagram, and provides an understanding of the nature of the interleaved phase between R and T. A phase field modeling study¹⁷ agrees that the ultrahigh piezoelectricity of [001]-oriented ferroelectric single crystals at the MPB originates mainly from polarization rotation.

Although significant steps in fundamental understanding have emerged from the work of Refs. 16–17, these are inherently simulation-based. A predictive analytic theoretical formulation, where key phenomena can be connected by closed-form expressions would be an important synergetic contribution. We proceed through the Landau-Devonshire (LD) approach^9,15,18,19 concentrating on the phases at the phase boundary on either side of the MPB and show how the analytic structure of the model connects to the phenomenology.

We seek to construct the minimal model which will capture the essential physics of the high-response piezoelectric materials PMN-PT and PZN-PT at the two edges of the MPB (though the model is quantitative, points are missed like the permittivity maximum lying in the M phase⁵). As shown in a recent article,²⁰ a simplified description of the R/T phase transition, is possible using only up to fourth order terms in the L-D free energy, which in the present notation takes the form

where p_i are the components of polarization, the coefficient α of the quadratic term is positive, and the coefficient β of the quartic term is also positive. The coefficient k of the anisotropy term is a proxy for the titanium concentration x, with k < 0 defining the R state, k > 0 defining the T state, and k = 0 defining the MPB; the model of Eq. (1) has the property that on the k < 0 side, the R state p = p₀(±1, ±1, ±1) is a stable minimum, whilst on the k > 0 side the T state $p=p0′(0,0,±1)$⁠, etc., is a stable minimum (a more complicated description is needed to describe the M region^10,16). When k is small and negative, on the R side near the MPB, the anisotropy is small and the polarization nearly isotropic; in this region the polarization may be easily rotated⁸ and there is a hyper-response to external probes, e.g., in the dielectric polarizability χ.^19,20

However, the level of description provided in Eq. (1) is not adequate to fully account for the properties of the intensively studied³ and commercially important engineered-domain or R4 samples. The R4 system has global tetragonal symmetry, but is in fact a micro-to-nanoscale mixture of the four kinds of z-poled rhombohedral domains (±1, ±1, 1)⁵ (see Fig. 1(b), black arrows). Table I illustrates an example¹³ of the dramatic non-cubic deviations in piezo and dielectric properties presented by these hyper-responsive samples indicating that a successful model must break the cubic symmetry. Treating the R4 phase requires an averaging procedure described in supplementary material.

We see in Table I, the characteristic hyper-responsiveness of PMN-30%PT in the giant values of dielectric permittivity (⁠$ϵ33T=7800$⁠), piezo coefficient (d₃₃ = 1981 pm/V), and uniaxial compliance (⁠$s33E=67.7 10−12 m2/N$⁠). But a puzzling complexity is added by the anisotropy, e.g., in dielectric constant (⁠$ϵ11T=3600$⁠), whilst there is a dramatic absence of hyper-responsiveness in the shear quantity d₁₅⁶–why is this so small?

In terms of modeling, we note that LD is a global model. Hence if the global symmetry of the sample is tetragonal, despite a local composition of (±1, ±1, 1) domains, then the symmetry needs to be broken in a tetragonal sense. Thus we introduce the main development in this paper, a L-D model for the R4 samples, which omitting lattice terms is

The innovative step in this Eq. (2) is the addition of the tetragonal symmetry-breaking term σ, where σ, obtained by the fitting experimental data, is negative. The effect of σ is to rotate the rhombohedral polarization vector away from the corners of the cubic reference towards the xy plane, as shown in Fig. 1(b) (red arrows) (note that the σ-term is multiplied by k, so that it vanishes at the MPB), allowing anisotropy in the polarization rotation.

The nearly coplanar polarization vectors (Fig. 1(b), red arrows) are easy to rotate by an external E-field E_z, which is nearly orthogonal to the polarization vector, hence leading to hyper-active ϵ₃₃, but the vectors are hard to rotate by a x-field E_x, which is nearly parallel to polarization, hence the much lower value of ϵ₁₁ (Table I); a similar explanation accounts for the small d₁₅. In fact we shall see that the model (2) can quantitatively describe the properties of the hyperactive piezoelectric materials.

Another aspect of the data for these materials is that a set of sum rules approximately hold: for example, in Table I, note that, since d₃₁ = d₃₂, then d₃₃ + d₃₂ + d₃₁ ≃ 0, and other sum rules are also found (see below). The physical content is that in systems dominated by a near-critical behavior of the polarization rotational degree of freedom, there is little coupling to volume, hence the existence of the sum rules supports the rotational dominance originally proposed by Fu and Cohen⁸ and by Damjanovic,⁹ who also stresses the relevance of the extender mechanism. These sum rules are also exactly satisfied by the LD theory of (2)–extended to include the lattice coupling–as will be demonstrated below.

Starting from the simplest possible model, the fourth-order model considered in Ref. 20 which predicts the stability of the R and T phases by varying the sign of the anisotropy terms, we add three key extensions. First, as discussed above, we include the σ term. Second, in describing the phenomena such as piezoelectricity and compliance, the terms coupling the polarization to the lattice must be included. Third, a useful device, based on the fact that in the rotational-dominance paradigm the high-response phenomena are associated with proximity to the MPB (i.e., to k = 0 in Eq. (2)), is that in the theoretical analysis, we typically expand in powers of k keeping only the lowest order terms. Care is also taken (supplementary material) with treatment of the cancellation of the x and y polarization components.⁶ The resulting simple, intuitive expressions also will be shown to attain quantitative accuracy. In addition the treatment enables the discovery of previously unknown sum rules and explains their physical origin, and reveals a fresh understanding of film clamping by the substrate.

The complete LD model is

Here p_i is component i = x, y, z of polarization, and S_p is the strain component p in compressed (Voight) notation (indices follow the IEEE index convention²¹ 1 ≤ p, q ≤ 6; 1 ≤ i, j ≤ 3). The first two terms describe an isotropic polarization $p2=α/2β′$ (⁠$β′$ positive, α positive if, as we assume, T lies below the ferroelectric transition temperature T_c, see Fig. 1(a)). The third, the anisotropy term, determines the direction of polarization which controls the crystal structure (Fig. 1); here $σ′$ quantifies the key pseudo-tetragonal anisotropy of the R4 phase, and will later be shown to strongly impact the material response. The fourth and fifth terms in Eq. (3), defined in terms of the unpolarized (cubic) material, describe, respectively, the elastic energy, with elastic modulus matrix C⁰, and the strain-polarization electrostriction with coupling matrix Q⁰.¹⁸ The sixth and seventh terms describe the effects of external electric field E and stress T.

Minimizing the free energy, we obtain the polarization, (assuming −1 < σ < 0–see supplementary material Eq. (26)) in the R (k < 0) and T (k > 0) states

Here the unprimed notations α, β, and σ define renormalized quantities modified by coupling to strain (supplementary material). Since σ is found to be negative, it is seen that in R4, the polarization is canted away from the (1, 1, 1) direction towards the xy-plane (see Fig. 1(b) red arrows).

The free energy difference between the R and T states to leading order in k is

it is seen that for small k, the R state is stable for k negative, and the T state is stable for k positive, our approximation to Fig. 1(a). We have also shown that the states are a local minimum of F for appropriate k—no two minima exist simultaneously in this model.

Of great interest are the response functions (RF's) such as piezoelectric coefficient in the R4 state. The RFs (see supplementary material) are the elements in the inverse of the 9 × 9 matrix of second derivatives of Eq. (3), taken at the stationary point and at zero external field T, E. Of these elements only 11 are nonzero and independent, and only 10 are accounted for in our theory—s₆₆ is not coupled to the rotator mechanism we seek to understand. These RFs are all singular in k, going as 1/k, expressing the fact that at the MPB where k = 0, there is free reorientation of polarization. Our universal approach retains only the k⁻¹ term in these expressions. The results can be parameterized in terms of 4 coefficients

where the $q̂1$ are components of the electrostrictive coefficient in terms of stress: $q̂1=[ (C0)−1Q0 ]11, q̂2=[ (C0)−1Q0 ]12, q̂3=(C440)−1Q440$⁠. Key RF components are

Figure 2 shows a graphical comparison of the L-D formulae with the experimental data^22,23 for 5 materials, (numerical values are given in the supplementary material). The L-D approach is seen to be remarkably successful in reproducing the experimental data. The data peculiarities we drew attention to above are now explained in terms of angular effects, especially strong in PMN-33%PT, where the polarization is canted to within only 13° of the xy-plane (Fig. 1(b)). This orientation is very efficient for producing a large d₃₃, and especially ϵ₃₃,⁵ because a E-field almost transverse to the polarization vector produces an optimally large rotational effect. But in the case of ϵ₁₁ and d₁₅, the relevant E-field lies almost parallel to the polarization vector, with little effect in rotating it, so these quantities are relatively small, despite their formally divergent 1/k behavior.

Various known identities²¹ and previously unrecognized sum rules are obeyed by the L-D theory. The coupling factor $k33l$ is seen from Eq. (7) to be unity in the leading-k L-D model

The identity $d152/(ϵ11ϵ0)=s44E−s44D$ is satisfied. The sum rules

have not been previously reported and arise from the requirement that in a strict rotator model close to the MPB polarization rotation does not couple to volume. Hence any strain along say, the z-axis must be compensated by negative strains along the x- and y-axes so as to conserve volume. In Table II, we present a quantitative comparison of the accuracy of the identity Eq. (8) and the sum rules Eq. (9) with experiment. It is seen that both are accurate to within a few percent (the negative value in column 4 for PZN-8%PZT is an artefact of the experimental uncertainty). While non-rotationally dominant piezoelectric materials such as PMN-42%PT²² still exhibit signs of the sum rules, the sum rules are obeyed to less precision than in the mainly rotator piezoelectrics.

As an application of our results, consider the piezoelectric coefficient of a domain-engineered thin film of relaxor piezoelectric on an infinitely rigid substrate. The effective film piezo response d_33,f is given by the well-known Lefki-Dormans²⁴ result

Now if we insert the sum rules (9)$s11E+s12E+s13E=0, d33+2d31=0$⁠, Lefki-Dormans reduces to d_33,f ≡ 0! The clamping effect of adhesion to a rigid substrate completely eliminates the universal piezoelectric response, leaving experimentally a small residual from polarization extension and extrinsic effects;²⁵ pure rotators are perfectly clamped.

The underlying philosophy of our model is the rotator mechanism,^8,9 in which the polarization angle is subject to facile rotation. The result of this facile rotation is seen in several phenomena: multiple phases (Fig. 1) and super-high permittivities, piezo coefficients, and compliances (e.g., $ϵ33T, d33,s33E$ in Table I). In addition to these phenomena, data on the widely used R4 or engineered-domain samples reveals another mechanism at work: symmetry breaking caused by global tetragonality. For example we cite the strong anisotropy in the permittivity (⁠$ϵ33T$ vs. $ϵ11T$ in Table I), and the very small value of the shear coefficient d₁₅ in Table I. To capture these effects, we modify and extend the existing state-of-the-art LD^9,15,18–20 by adding an anisotropy term σ (Eqs. (2) and (3)) to treat the global tetragonality and the coupling to the lattice; it is important that in our formalism, σ is multiplied by k, so that its effect vanishes at the MPB and for simplicity, we show that a very low order polynomial in polarization is sufficient for accuracy.

The result is a set of simple formulae (e.g., Eq. (7)) which quite accurately describes the data with a few parameters. Our model predicts that the effect of global tetragonality is to shift the polarization vector away from the nominal rhombohedral (±1, ±1, 1) directions towards the xy-plane (Fig. 1(b), red arrows). We also discovered in the dataset, a hitherto ignored set of sum rules, e.g., $∑i=13d3i=0$⁠, an exact result of our theory (Table II), reflecting that 3-dimensionally isotropic perturbations do not impact the rotational degree of freedom.

The R-phase properties at the left side of the MPB have an analogue in the properties of the T-phase on the right side of the MPB. However since the polarization vector in the T-phase is along (0, 0, 1), this phase is insensitive to properties involving coupling to the z-polarization, such as $ϵ33T$⁠, and d₃₃. Instead the large responses in this phase are those involving the x- and y-components.

Finally, our theory has implications in the context of the Lefki-Dormans theory²⁴ for the effective d₃₃ piezo coefficient of a PE film clamped to an inelastic substrate (⁠$d33,f≡0$⁠), capturing the sensitivity of high-response piezoelectrics to clamping.

Our L-D model is constructed in a minimalist fashion and yet embodies the intrinsic physics. The physical origin of the form can be thought of as a coarse graining of a more complex many-body Hamiltonian. This includes reflecting the polarization rotation mechanism identified from atomistic simulations by Cohen et al. and the critical fluctuations extending down to lower temperatures as described by Blinc. The model describes well the behavior of the materials at small but finite “k” on either side of the transition but is too simple to encompass their coexistence—a reasonable compromise in order to achieve simple analytic forms with physical content. To make contact with the MPB, a simple analytically solvable two state model with a first order transition can be constructed following the prescription of Ref. 26 

where W > 2 for a first order transition, x = 1 → R, x = 0 → T, with x a renormalized concentration. We reserve this extended model for future work, and a multicomponent formalism could also be constructed, including, for instance, the monoclinic/orthorhombic phases near the MPB following the seminal work of Ref. 11.

See supplementary material for the calculation of the full 9 × 9 response function, including details of the averaging over parallel polarization components and numerical results.

This research was partially supported by the DARPA MESO (Mesodynamic Architectures) Program under Contract No. N66001–11-C-4109. The views expressed are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government. One of us (F.C.) was funded jointly by the Scottish Condensed Matter Doctoral Training Centre and the National Physical Laboratory under an EPSRC CASE scholarship and by the Scottish Universities Physics Alliance under an Early Career Research Visit scholarship.