## Abstract

Most key elements of ferroelectric properties are defined through the hysteresis loops. For a ferroelectric ceramic, its loop is contributed collectively by its constituent grains, each having its own hysteresis loop when the ceramic polycrystal is under a cyclic electric field. In this paper, we propose a polycrystal hysteresis model so that the hysteresis loop of a ceramic can be calculated from the loops of its constituent grains. In this model a micromechanics-based thermodynamic approach is developed to determine the hysteresis behaviour of the constituent grains, and a self-consistent scheme is introduced to translate these behaviours to the polycrystal level. This theory differs from the classical phenomenological ones in that it is a micromechanics-based thermodynamic approach and it can provide the evolution of new domain concentration among the constituent grains. It also differs from some recent micromechanics studies in its secant form of self-consistent formulation and in its application of irreversible thermodynamics to derive the kinetic equation of domain growth. To put this two-level micromechanics theory in perspective, it is applied to a ceramic PLZT 8/65/35, to calculate its hysteresis loop between the electric displacement and the electric field (*D* versus *E*), and the butterfly-shaped longitudinal strain versus the electric field relation (*ϵ* versus *E*). The calculated results are found to be in good quantitative agreement with the test data. The corresponding evolution of new domain concentration *c*_{1} and the individual hysteresis loops of several selected grains—along with those of the overall polycrystal—are also illustrated.

## 1. Introduction

The hysteresis loops of a ferroelectric ceramic define many key elements of its ferroelectric characteristics. In the electric displacement versus electric field (*D* versus *E*) plot, it provides the dielectric constant, saturation and remanent polarization, and coercive field. In the longitudinal strain versus electric field (*ϵ* versus *E*) plot, it gives the piezoelectric constant and the maximum strain upon reversal of the electric field. These quantities are significant for its functional applications such as sensors, actuators and transducers. These two seemingly single-lined loops, however, are collectively contributed by the loops of all constituent grains, for each has its own loop when the ceramic is under a cyclic electric loading. This is due to the different orientations of the constituent grains in relation to the externally applied field, and it is through the orientational average over the properties of these differently oriented grains that these two loops arise.

It is in this spirit that the present work was undertaken. Since the hysteresis behaviour of the constituent grain is caused by domain switch and a crystallite at its generic state consists of the parent and switched domains, we will first develop a micromechanics-based thermodynamic approach to determine the evolution of new domain concentration, *c*_{1}, as the electric field increases. Then by means of a dual-phase composite theory the effective polarization and effective strain of each crystallite at a given *c*_{1} will be determined. The evolution of *c*_{1} is driven by the thermodynamic driving force, which we will also derive. This force depends on the local electric field and local stress of the crystallite, which differ from grain to grain, and also from the externally applied one. As a consequence the evolution of the new domain among the constituent grains is heterogeneous, and so are their effective polarization and electromechanical moduli. After these ferroelectric properties are determined for all constituent grains, their transition to the ceramic polycrystal level can be calculated through an electromechanical self-consistent scheme.

Classical ferroelectric theories—such as those of Devonshire (1949, 1951, 1954) and Slater (1950)—are phenomenological and are often based on the functional representations of Gibbs or Helmholtz free energy of the system. These theories are of course useful in their own right, for they have provided explanations for many experimentally observed behaviours, including the hysteresis loop, but they could not provide information regarding the evolution of microstructure. For a general account of these phenomenological theories and pertinent ferroelectric phenomena, one may refer to Jaffe *et al*. (1971), Lines & Glass (1977), Jona & Shirane (1993), Ikeda (1996) and Strukov & Levanyuk (1998), among others. The proposed theory differs from these classical ones in some significant ways—it is a micromechanics-based thermodynamic approach, and it could provide the evolution of microstructures in terms of the volume concentration of switched domain in each constituent grain.

We have seen several important studies on the micromechanical modelling of domain switch in recent years. Hwang *et al*. (1995) have used a polycrystal model and an energy criterion to determine the condition of switch in each grain, assuming that all grains are in a uniform state of stress and electric field. The heterogeneous nature of the grain fields was subsequently accounted for by Arlt (1996) in a dielectric hysteresis analysis. Subsequently, mean-field approaches using spherical inclusions to represent the switched domains were developed by Chen *et al*. (1997), Chen & Lynch (1998) and Hwang *et al*. (1998). A self-consistent polycrystal model in the spirit of crystallographic slips in crystal plasticity was developed by Huber *et al*. (1999), and a rate-dependent model—patterned after the rate-dependent polycrystal viscoplasticity—was introduced by Huber & Fleck (2001). A micromechanics-based thermodynamic approach making use of an effective domain to represent the collective effects of all active variants was developed by Li & Weng (1999) and Li (2003*a*). In addition, finite element simulations have also been carried out by Hwang & Waser (2000) again using the energy criterion, and by Fotinich & Carman (2000) using the critical electric displacement criterion. On the phenomenological front, continuum theories have also been proposed by Chen & Peercy (1979), Chen (1980), Chen & Marsden (1981), Bassiouny *et al*. (1988*a*,*b*), Bassiouny & Maugin (1989*a*,*b*), Lynch (1998), Cocks & McMeeking (1999), Kamlah & Tsakmakis (1999), Landis & McMeeking (1999) and Landis (2002), among others.

In light of the significant amount of work cited above, it is necessary to point out the new features of this model before we proceed to present it. This is a two-level, micromechanics-based self-consistent polycrystal model: the first level exists on the grain level that is treated as a dual-phase crystallite involving the parent domain and the switched domain, and the second one involves the transition form the constituent grains to the transversely isotropic polycrystal. Inside each grain the dual-phase crystallite is taken to form lamellar domains. The thermodynamic driving force for domain switch will be derived from a micromechanics-based irreversible thermodynamic principle, and the effective ferroelectric behaviour at a given level of new domain concentration will be calculated by means of a homogenization theory. Finally, the transition of the ferroelectric response from the local grains to the overall polycrystal will be carried out through a self-consistent scheme. Among the work cited above, only Huber *et al*. (1999) and Huber & Fleck (2001) have treated the grain–grain interactions through a self-consistent scheme. Their formulations were incremental, involving the tangent moduli of the crystallites and the overall polycrystal. Our formulation instead will be based on a *total* form involving the secant moduli. While the incremental formulation in general has a wider applicability (e.g. it can cover the non-radial loading), this total scheme is substantially simpler for the calculation of hysteresis loops. At the grain level, Huber *et al*. made a direct application of crystal plasticity with various slip systems—or ferroelectric domains now—to capture the ferroelectric response of each crystallite. They have included the contributions of all potential variants (six in tetragonal crystallites), but their formulation does not account for the domain morphology nor does it consider the difference in the electromechanical moduli of the parent and switched domains along common axes. In our formulation all these factors will be fully accounted for. In addition, their transformation criterion was established in a way analogous to the yield criterion for the operation of a slip system. In our approach, the thermodynamic driving force for domain switch will be derived from the reduction of Gibbs free energy of the crystallite. This driving force is exact when the domain morphology is lamellar. Of course, it should be recognized that lamellar domains tend to develop only under simple proportional or cyclic loading, so they are not recommended for calculations involving a non-radial loading. But despite such a restriction the two-domain lamellar morphology does carry the virtue that it has both the original parent domain and the evolving switched domain, and by this model the evolution of new domain concentration *c*_{1} for all grain orientations and for the overall polycrystal, and the derivation of the hysteresis loops of the polycrystal from those of its constituent grains, can all be addressed.

## 2. Constitutive equations of the parent and switched domains inside the crystallites

In order to pave the way for the development of a micromechanics-based thermodynamic approach to domain switch in the constituent grain, we first briefly recapitulate the constitutive equations of a piezoelectric crystal. In anticipation that the switched domain grows from the parent phase during the switching process, we will take the parent domain as reference. Then its linear electromechanical constitutive relations can be written by taking stress *σ*_{ij} and electric displacement *D*_{m} as one pair and strain *ϵ*_{kl} and electric field *E*_{n} as another, as(2.1)in tensorial form, where ^{(E)} and ^{(ϵ)} are the elastic modulus and dielectric permittivity tensors measured at constant electric field and constant strain, respectively, and is the piezoelectric tensor. The superscript *T* stands for transpose. Such a pairing has the advantage that, in the absence of body force and free-charge density, both *σ* and *D* satisfy the divergence theorem(2.2)whereas *ϵ* and *E* are derivable from the potential-like displacement *u*_{i} and electric potential *ϕ*,(2.3)Such a pairing is particularly useful to address the boundary value problems in piezoelectricity.

There are, however, other ways to write the constitutive equations, and the choice largely depends on the issue at hand. To study domain switch under a mechanical stress and electric field, it is better to put them as a pair, and strain and electric displacement as another. This type of pairing is convenient for the formulation of Gibbs free energy which we will undertake later. In this way, we have(2.4)where *s*^{(E)} and (*σ*) are the elastic compliance and dielectric permittivity tensors measured at constant electric field and constant stress, respectively, and is the piezoelectric compliance tensor (sometimes called piezoelectric moduli tensor, as it represents the direct piezoelectric effect under application of stress).

The same material constants still hold for the reoriented domain when measured along its local oriented axes. But due to domain switch, a change of unit cell orientation and dipole moment also occurs, resulting in the eigenstrain and eigen-polarization inside the reoriented domain. The constitutive relations of the switched domain then carry these additional terms, as(2.5)when expressed in its local oriented axes.

These coupled relations can be cast in a unified fashion for computational purpose, as(2.6)with *X* serving as the load, *Y* the response and the compliance. Specifically in Nye's (1979) contracted notations,(2.7)where , is the eigenfield that automatically vanishes in the parent domain. We further write the moduli tensor by so that =^{−1}. Since at a generic state the parent and switched domain may coexist in a crystallite, we shall henceforth refer the parent domain as phase 0 and the switched domain as phase 1, and denote the volume concentration of phase *r* as *c*_{r}. We will also write their respective moduli in the common parent coordinate as _{0} and _{1}, respectively, with the difference(2.8)For later calculations the fully poled crystallite will be taken to be transversely isotropic, with the 6*mm* symmetry, so that the components of its compliance tensor are(2.9)These components still hold for the reoriented domain in its local axes, but require the usual fourth-order tensor transformation to convert them to the parent coordinates for _{1} and for Δ above.

The eigenstrain and eigen-polarization of the switched domain result from the reorientation process. As an illustration, figure 1*a*,*b* shows the change of crystal structure during spontaneous polarization of BaTiO_{3} at room temperature. For the tetragonal unit cell in figure 1*b*, there are four potential 90° variants as depicted in figure 1*c–f*, and a 180° one, for reorientation. By taking the poled state with an upward dipole as reference, the eigenfields for the four 90° switches are given by(2.10)and(2.11)where *c* and *a* are the dimensions of the tetragonal cell, and *P*_{s} is the saturation polarization. The 180° switch produces a negative polarization of −2*P*_{s} in the 3-direction but it produce no net strain with such a direct switch. In order to capture the sharp tails of the butterfly-shaped longitudinal strain versus electric field relation, two consecutive turns—from 0 to 90°, and then from 90 to 180°—are required.

Of course not all ferroelectric domains are in the tetragonal state at room temperature. For the PLZT to which we will make calculations later, it is in a rhombohedral state. Figure 2 shows the four commonly observed crystal structures: (i) cubic, (ii) tetragonal, (iii) orthorhombic and (iv) rhombohedral. The tetragonal structure can be viewed as a stretch of the cubic one along its edge, whereas the orthorhombic and rhombohedral structures can be viewed as a stretch along the face diagonal and the body diagonal, respectively. The eigenstrain and eigen-polarization associated with a rhombohedral structure can be similarly determined as with the tetragonal one, but in this case there are three 70.5° potential domains and another three 109.5° ones, in addition to the remaining 180° switch.

## 3. The polycrystal model

The proposed polycrystal model is motivated by the observation that ferroelectric ceramics possess the polycrystalline structure and domain switch takes place inside the constituent grains. As a consequence we believe it is more physically consistent to address the issue of ferroelectric hysteresis first from domain switch inside the constituent grains and then translate the obtained grain properties to the polycrystal level through an averaging process.

A typical morphology of a ferroelectric ceramic—taken from DeVries & Burke (1957)—is shown in figure 3*a*, where domains with lamellar structure inside the grains are clearly visible. While sub-domains with different orientations can also be seen in some grains, the morphology is predominantly lamellar and, for this reason, we model the polycrystal with the sketch of figure 3*b*. In this diagram, the white regions in each grain represents the parent domain and the grey regions the switched domain. In each crystallite, the active—or switched—domain is chosen from all potential variants based on the criterion of maximum thermodynamic driving force (to be given in §4) in relation to the parent one. The constitutive equations of the parent phase are then given by equation (2.4), and those of the switched domain by equation (2.5) but transformed to the coordinate system defined by the parent domain. Each grain is now a dual-phase crystallite consisting of the parent domain and the switched one with a volume concentration *c*_{1} that differs from grain to grain. The eigenstrain and eigen-polarization exist only inside the switched domains. While it is appreciated that this two-variant representation could be an over idealization for general applications, it does possess the virtue that both the parent and the switched domain are present, and it could yield the major features of ferroelectric phenomena.

In order to provide a global view of the polycrystal model first, we will outline the mechanics of transition from the grain to the polycrystal level in this section and leave the calculation of new domain inside the constituent grains to §4.

Now under a given external electromechanical load , let the average stress and electric field of a grain (i.e. crystallite) be denoted by , where an overbar in signifies that it is a volume-averaged quantity over the entire polycrystal while in , it is averaged over the considered crystallite. This local differs from the global , and it also differs from grain to grain due to their orientations. For this reason, the is to be understood as , where (*θ*, *ϕ*) are the Euler angles of a grain with respect to the axial direction of the electric field to be applied to the polycrystal. (Note: since we will make specific calculations only for the axial response, two Euler angles are sufficient.) The average of over all grain orientations, however, must give rise to for consistency. We further denote the overall electromechanical response of the polycrystal by and that of a generic crystallite by or simply . The objective of the polycrystal model is to calculate at a given level of through of all constituent grains, which basically is controlled by , which in turn depends on the applied . In a polycrystalline aggregate these quantities are inter-related to each other through the self-consistent relation.

To establish such a relation for a ferroelectric ceramic, we may recall Hill's self-consistent relation in polycrystal plasticity (Hill 1965). In the ferroelectric context it can be generalized to(3.1)in the incremental form, where ^{*} now stands for the electromechanical constraint tensor of the matrix and is related to its *tangent* moduli. Under monotonic, proportional loading it can be modified to a total form, as (see Berveiller & Zaoui 1979 in the context of polycrystal plasticity)(3.2)where now depends on the *secant* moduli of the matrix, _{s} (defined as ), through(3.3)The secant **S**–tensor here is the Eshelby-type**S**–tensor (1957) but cast in the electromechanical context. For a spherical inclusion embedded in a transversely isotropic matrix, its components can be calculated from Huang & Yu (1994) and Dunn & Wienecke (1997). Both derivations, however, made use of (*ϵ*^{*}, *E*^{*}) as the pair of eigenfields, so they need to be converted to the (*ϵ*^{*}, *D*^{*}) pair for the present application. A direct conversion procedure can be found in Li & Dunn (2001).

Equation (3.2) provides the average electromechanical load for a generic crystallite whose orientation is defined by the Euler angles (*θ*, *ϕ*). It follows that, from the orientational average over all orientations of the constituent grains, the overall electromechanical load and response of the polycrystal must satisfy(3.4)where the brackets 〈.〉 stand for the said orientational average. Since in equation (3.2) depends only on the effective properties of the polycrystal and is independent of the crystallite properties, application of the orientational average to equation (3.2) automatically leads to equation (3.4), making this local-to-global transition truly ‘self-consistent’.

It is this —not the applied —that controls the evolution of new domain concentration *c*_{1} in each crystallite. After *c*_{1} is determined, the electromechanical response can be calculated. Once is known for all grain orientations, the overall of the polycrystal follow from equation (3.4). However, the corresponding , , and must satisfy the self-consistent relation (3.2), so these quantities are intertwined and some iteration procedure is often required to find the solution.

We now turn to the domain switch process—that is, the kinetic equation—and the overall response of the crystallites.

## 4. Evolution of new domain concentration and calculation of electromechanical response of the crystallites

At a given level of , the next issue is the determination of for each crystallite. These terms now will be simply written as and for brevity, with the understanding that the calculation procedure must be performed for every grain orientation. In each grain, we take the fully polarized state as reference and write the fields of both parent and switched domains in the common parent coordinate, so that there is the electromechanical heterogeneity Δ=_{1}−_{0} as stated in equation (2.8), and eigenfield in the switched domain. In order to calculate , it is noted that, as the magnitude of increases, domain switch will continue to take place and the volume concentration of the switched domain, *c*_{1}(*θ*, *ϕ*)—or simply *c*_{1}—will continue to increase. Such an increase is the source of nonlinearity in ferroelectric response. If *c*_{1} is known, the dual-phase homogenization theory will yield in terms of and *c*_{1}, as(4.1)or in terms of the mechanical and electrical components,(4.2)where *Y*^{*}=(*ϵ*^{*}, *P*^{*}), is the Eshelby-type equivalent eigenfield introduced into the regions of the switched domain to join as the total eigenfield, so that its _{1} could be replaced by _{0} to yield the same electromechanical field in this phase. To plot the left branch of the hysteresis loop, the eigenfield *Y*_{sp}=(ε_{sp}, *P*_{sp}), resulting from the spontaneous polarization—such as from the cubic to tetragonal phases in figure 1*a*,*b*—should be added to equations (4.1) and (4.2).

In view of equation (4.1) or (4.2), it is evident that the evolution of *c*_{1} inside the crystallite must be established now. This is the issue of kinetic equation for domain growth. To obtain the kinetic equation, we may appeal to the principle of irreversible thermodynamics, from which the thermodynamic driving force for domain switch can be derived. In this context we note that, under an external on the crystallite, the driving force arises from the reduction in its Gibbs free energy from the switched state to the un-switched, referenced state. The former is given by (Stratton 1941; Eshelby 1957)(4.3)where , are the mechanical traction and electric potential on the boundary surface *S* of the grain with an outward normal *n*_{j} at position *x*_{i}, to give rise to a uniform stress and electric field . Gibbs free energy without domain switch is simply(4.4)with a *unit* volume of appropriate scale for the crystallite, where and are strain and electric displacement of the parent domain without switch, given by in tensorial notations. The difference in Gibbs free energy then follows as(4.5)Following the principle of irreversible thermodynamics, we have(4.6)where equality holds only for the reversible process (i.e. d*c*_{1}=0). Then by taking *c*_{1} as the microstructural parameter, the conjugate thermodynamic driving force for domain switch can be calculated from (Rice 1975)(4.7)To derive this driving force explicitly, it is important to find the analytical expression of Δ*G*. This has been done in the context of martensitic transformation in shape-memory alloys (Lu & Weng 1997, 1998). Now with the help of lamellar morphology as already illustrated in figure 3*a* for a polycrystal, and by Merz (1954) and Hooton & Merz (1955) for many BaTiO_{3} single crystals, it can be extended to the present electromechanical context. Indeed, lamellar structure is morphologically favoured over other ones for it represents the minimum energy associated with the martensitic-type, diffusionless transformation (Khachaturyan 1983; Wayman 1992). Since lamellar structure is a special form of aligned ellipsoidal inclusions, the change Δ*G* can be derived explicitly by means of the Eshelby (1957) and Mori–Tanaka (1973) method. This approach is similar to the Maxwell–Garnett approximation in transport phenomena (Nan 1993) and it has also been successfully applied to study the electrostriction of composite materials (Nan & Weng 2000). In the pure elastic context it has the virtue of being related to the Hashin & Shtrikman (1963) and Willis' (1977) bounds with aligned ellipsoidal inclusions and, with the lamellar configuration at hand, it becomes exact (see Walpole (1969) and Weng (1984, 1990, 1992) for further discussions).

With this approach, it can be shown that(4.8)where, as stated for equations (4.1) and (4.2), *Y*^{*}=(*ϵ*^{*}, *P*^{*}) is Eshelby's-type equivalent electromechanical eigenfield of the switched domain, and , which carries a subscript 1, is its average field (all quantities are tensorial). In addition,(4.9)where is the electromechanical field concentration tensor of a single ellipsoidal inclusion with _{1} embedded in an infinite matrix with _{0}, and *Q*=(*c*_{1}*B*_{1}+*c*_{0}*I*)^{−1}, *I* being the fourth-order symmetric identity tensor. The Eshelby-type -tensor here is associated with the lamellar inclusion that is generally oriented with a certain angle from the poling direction of the parent domain inside the considered crystallite; explicit forms of its non-vanishing components can be found in Li (2003*b*). In the thermoelastic context, this energy returns to the mechanical potential energy for martensitic transformation in shape-memory alloys.

With equations (4.8) and (4.9), the thermodynamic driving force can be derived from equation (4.7), and put into the form(4.10)after some algebra, where(4.11)and . We reiterate that, with lamellar morphology, this driving force is exact.

In contrast to the driving force, there is a resistance force to domain switch that arises primarily from the energy dissipation associated with the domain wall movement. This force tends to increase with increasing domain concentration *c*_{1}. Such an increase is required because if it stays constant, the switching process would be instantly complete and the nonlinear portion of the *D* versus *E* loop of a ferroelectric single crystal would be solely vertical, but experiment on BaTiO_{3} crystal clearly indicates that it is not the case (Merz 1952). In the study of martensitic transformation in shape-memory alloys, Sun & Hwang (1993*a*,*b*) suggested a linear energy function to derive the thermodynamic resistance force. We have found that, in order to describe the general hysteresis loop for ferroelectrics, the double-exponential function,(4.12)has a wider applicability over the entire range of *c*_{1}. Here *h*, *b*_{0} and *b*_{1} are material constants. When the driving force is sufficient to overcome the resistance force, that is, when *f*_{drv}=*f*_{res}, the switching condition is reached, i.e.(4.13)This represents the kinetic equation of domain growth. Solving for *c*_{1} at a given , the evolution of new domain concentration can be determined as increases. Since *A*, *B* and *C* all contain *c*_{1}, it is easier to use *c*_{1} as input in real calculations and treat as output. This calculated *c*_{1} then can be substituted into equation (4.1) or (4.2) to find for the considered crystallite. After it has been calculated for grains of all orientations, the effective electromechanical response of the polycrystal can be evaluated from the orientational average of equation (3.4). For the calculation of the left branch of the hysteresis loop, the applied —or under an ac field—can be increased in the reverse direction. The right branch follows simply as a symmetric counterpart of the left one. Together they form the hysteresis loop.

## 5. Hysteresis loops of a 8/65/35 PLZT under a cyclic electric field

The established framework allows us to determine the hysteresis behaviour of a ceramic polycrystal from those of its constituent grains in a self-consistent manner. It also allows us to see the corresponding evolution of new domain concentration among grains of different orientations. As an application we now use it to examine the *D* versus *E* hysteresis loop of a ceramic 8/65/35 PLZT and the associated butterfly-shaped *ϵ* versus *E* relation. This material has been tested by Hwang *et al*. (1995) and Lynch (1996). Hwang *et al*. (1998) further pointed out that this ceramic at room temperature is composed mostly of rhombohedral crystallites. For this reason, we have adopted the rhombohedral instead of the tetragonal structure to model the domain switch. The material constants of the crystallites adopted in computations arewhere direction 3 is taken to be the poling direction of the crystallite, and *E*_{11} and *E*_{33} are Young's moduli. In addition, we further took *b*_{0}=1.72, *b*_{1}=4.12, *h*=0.29 MPa.

Before we proceed to present the results, it is fitting to give a brief description of the iteration process in our computational procedure. We started out by setting initially for each grain to compute the according to its specific orientation. After was determined for every grain orientation, we then proceeded to find the corresponding from equation (3.4). With these , and , we went back to equation (3.2) to compute . If this newly computed was equal—or very close—to the originally assumed , the solution was taken to have been found; otherwise a new reflecting this newly calculated one was assumed for another round of iteration, until the solution was found.

In the computation, the polycrystal was modelled with 155 discrete grain orientations. It is known that during polarization of an initially isotropic ferroelectric ceramic, the polar orientation within each grain undergoes switch upon application of an electric field. To simulate the initial isotropy, these 155 grain orientations were taken to be uniformly distributed in their Euler angles, with *θ* indicated in figure 4*a*. Here, the vertical bars signify that grains with the *θ* values on the horizontal axis exist (the height has no meaning). Such an initially isotropic state is represented by point *A* in figure 5. Then the switching process under continues through *B* until the applied electric field reaches a certain level at which the polycrystal is fully poled as indicated by point *C*. This is the end state of the initial polarization and it also serves as the starting point for calculation of the hysteresis loop. The 155-grain orientations at this starting point are represented in figure 4*b*, which are highly confined in a cone of about 55° with respect to the original applied electric field direction (direction 3). The hysteresis computation then started with this transversely isotropic polycrystal. We then let the applied electric field decrease (or apply a reversed electric field) to zero, passing the coercive field at point *D*, and then increase further in this opposite direction to reach another saturation state at point *E*, where the polycrystal is again fully poled in the opposite direction. After that, if the electric field is reversed, it will follow the right branch of figure 5, back to *C*. As the right branch is a mirror image of the left one, it suffices to compute only the left branch. During this half loading cycle, the longitudinal electric displacement *D*_{3} versus *E*_{3} and the longitudinal strain *ϵ*_{33} versus *E*_{3} (overbar removed for simplicity) of the polycrystal were calculated with the self-consistent scheme outline in §3, and the evolution of *c*_{1}(*θ*, *ϕ*) and of the crystallites were determined with the micromechanics-based thermodynamic approach in §4. The choice of the active variant-pair in each crystallite was based on the maximum thermodynamic driving force *f*_{drv} given by equation (4.10).

The calculated hysteresis loop between *D*_{3} and *E*_{3}, and the butterfly-shaped *ϵ*_{33} versus *E*_{3} relations are shown in figures 6 and 7, respectively. The solid lines represent the theoretical results while the solid-dotted lines are measured data taken from Lynch (1996). The computed loop is seen to capture the fundamental characteristics of ferroelectric hysteresis, and the calculated remanent polarization, coercive field, and longitudinal strain, are all in good accord with the experimental data. The amplification of electric field in computation is a little larger than the one in experiment to guarantee a fully poled state at the end of each loading cycle. These theoretical results serve to substantiate the validity of the developed polycrystal model.

A closer look at figure 6 indicates that initially the electric displacement *D*_{3} decreases linearly with the electric field. After the field reduces to the vicinity of zero, the electric displacement *D*_{3} begins to decrease nonlinearly, marking the onset of domain switch. The evolution is slow in the beginning but becomes very fast near the coercive field, at 0.37 MV m^{−1}, after that it again proceeds slowly toward completion. A slight departure between the calculated results and test data is observed at the turning points. It reflects that domain reorientation in the theoretical model initially takes place a bit faster than the experimental one and proceeds to the saturation slower at the end, leading to a less abrupt theoretical curve (it may be noted that this departure could be reduced with the choice of a different resistance function—such as the double tangent function—but we decided not to pursue such a fine tuning because our objective is to bring about the essential features of the theory).

The *ϵ*_{33} versus *E*_{3} loop in figure 7 is seen to display the typical butterfly-shaped longitudinal strain versus electric field relation with sharp tails. As it is known that a direct 180°-switch does not generate any longitudinal strain during the evolution, and that the sharp tails of the butterfly-shaped loop reflect a non-stopping 180°-switch that, in the case of tetragonal crystals, can be viewed as a result of two consecutive 90° switches: the first part results in a reduction of the strain and the subsequent one leads to an increase. For the rhombohedral crystallites in our calculation, each grain undergoing the 180° switch was also taken to consist of two consecutive switching steps: the 70.5° switch and the 109.5° switch. The sharp tails displayed in figure 7 are a consequence of such a non-stopping two-step process. For the evolution of new domains in a polycrystal, it is also known that the majority of the constituent grains will undergo a 180° switch at the end. This is verified in figure 8, showing that, after completion of domain switch, all the grains have reoriented to a cone making about 55° (or spanning from the original 125° to 180°) along the (reversed) loading axis. This corresponds to the end state at point *E* of figure 5. Figure 8 in essence is identical to figure 4*b* when viewed from the opposite direction. It is also worth noting that, consistent with the measured data, the polycrystal model can generate a negative strain near the coercive field. This is due to the fact that the strain of the initially isotropic ceramic at point *A*—when the grain orientations were uniformly distributed—was taken to be zero, but at the coercive field the poling directions of all constituent grains are somewhat flattened out. That is, when the electric field drops to the coercive field, the majority of the grains are half way through their non-stopping 180° switch. Thus, the average strain of the polycrystal in the global longitudinal direction is smaller than in its isotropic, random state, resulting in a negative strain.

Finally, in order to disclose the heterogeneous nature of domain evolution among the differently oriented constituent grains, we provide the hysteresis loops and the evolution of *c*_{1} for four selected grain orientations in figures 9 and 10, respectively. These four grain orientations are denoted as *G*_{33}, *G*_{112}, *G*_{104} and *G*_{145} among the 155 orientations, with their respective Euler angles (*θ*, *ϕ*) indicated in the parentheses. The averaged loop of the polycrystal over the total 155 orientations is also depicted as a solid curve for comparison. The first two grain orientations have smaller *θ* values than the latter two, indicating that they are more opposed to the reversed applied field and thus experiencing a faster turn than the latter. But the average of these two groups—one more opposed and the other less so—must give rise to the averaged response of the polycrystal for self-consistency. The evolution process of each grain is seen to be highly dependent upon its Euler angle *θ*; a grain with a smaller *θ* evolves faster than one with a larger value.

## 6. Concluding remarks

Based on the observation that domain switch takes place inside the constituent grains and that it is their averaged behaviour that gives rise to the overall response of the polycrystal, a two-level micromechanical model has been developed to calculate the hysteresis loops of a ferroelectric ceramic. This involves a thermodynamic-driving switching process inside the grain, and a self-consistent transition from the constituent grains to the polycrystal level. The overall electromechanical response of the ceramic polycrystal has been calculated by the orientational average over its constituent grains. These grains are taken to be fully poled initially by an external electric field, with the grain orientations closely aligned within a cone of this field direction. The hysteresis loops of the polycrystal are then calculated by reversing the applied field until another fully poled state has been reached. The calculated hysteresis loop between the electric displacement and applied electric field, and between the longitudinal strain and electric field, are found to display the essential features of ferroelectric response, and are in good agreement with the measured data. In particular, the magnitudes of the remanent polarization and coercive field that define the size of the loop, and the butterfly-shaped strain with its sharp tails, are well captured. This theory differs from the classical phenomenological one in that it is a micromechanics-based thermodynamic approach, and that it could deliver the evolution of new domain concentration among the constituent grains. As an illustration of this special feature the hysteresis loops of several selected grains and evolution of their new domain concentration *c*_{1} have been displayed along with those of the overall polycrystal. The final grain orientations upon complete reversal of the applied field are also illustrated to demonstrate the reversal of polar orientations from the initial, fully poled state.

## Acknowledgments

This work was supported by the National Science Foundation, under grants CMS-0114801 and 0510409.

## Footnotes

- Received June 15, 2005.
- Accepted November 17, 2005.

- © 2006 The Royal Society