# A Boltzmann statistical approach for the analysis of polarization states in mixed phase ferroelectric materials: Application to morphological phase boundary

A. Pramanick, L. Daniel

► To cite this version:

A. Pramanick, L. Daniel. A Boltzmann statistical approach for the analysis of polarization states in mixed phase ferroelectric materials: Application to morphological phase boundary. Journal of Applied Physics, 2025, 137 (20), 10.1063/5.0258087. hal-05288552

HAL Id: hal-05288552

https://hal.science/hal-05288552v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Boltzmann statistical approach for the analysis of polarization states in mixed phase ferroelectric materials: Application to morphological phase boundary

A. Pramanick ; L. Daniel

<figure><img src="image_2.png" /></figure>

Check for updates

J. Appl. Phys. 137, 204101 (2025)

https://doi.org/10.1063/5.0258087

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

## Articles You May Be Interested In

* Temperature-dependent anhysteretic behavior of co-doped PZT

* *J. Appl. Phys.* (September 2018)

* Revisiting the blocking force test on ferroelectric ceramics using high energy x-ray diffraction

* *J. Appl. Phys.* (May 2015)

* Three dimensional phase-field simulations on the frequency dependence of polarization vectors and hysteresis loops in ferroelectric crystals

* *J. Appl. Phys.* (February 2019)

<figure><img src="image_5.png" /></figure>

Nanotechnology &
Materials Science

<figure><img src="image_6.png" /></figure>

Optics &
Photonics

<figure><img src="image_7.png" /></figure>

Impedance
Analysis

<figure><img src="image_8.png" /></figure>

Scanning Probe
Microscopy

<figure><img src="image_9.png" /></figure>

Sensors

<figure><img src="image_10.png" /></figure>

Failure Analysis &
Semiconductors

<figure><img src="image_11.png" /></figure>

Unlock the Full Spectrum.
From DC to 8.5 GHz.

Your Application. Measured.

Find out more

<figure><img src="image_12.png" /></figure>

Zurich
Instruments

<figure><img src="image_13.png" /></figure>

AIP
Publishing
# A Boltzmann statistical approach for the analysis of polarization states in mixed phase ferroelectric materials: Application to morphological phase boundary

Cite as: J. Appl. Phys. 137, 204101 (2025); doi: 10.1063/5.0258087

Submitted: 14 January 2025 · Accepted: 8 May 2025 ·

Published Online: 23 May 2025

<figure><img src="image_14.png" /></figure>

View Online

<figure><img src="image_15.png" /></figure>

Export Citation

<figure><img src="image_16.png" /></figure>

CrossMark

A. Pramanick<sup>1,2,a)</sup> ID and L. Daniel<sup>1,2,a)</sup> ID

## AFFILIATIONS

## ABSTRACT

Ferroelectrics are widely used for a broad array of technological applications due to their attractive electrical and electromechanical properties. In order to obtain large functional properties, material compositions are often designed to favor the coexistence of multiple ferroelectric phases. For such compositions, the macroscopically observed properties are variously attributed to easier domain switching and/or phase transition. Nevertheless, modeling of concurrent domain switching and phase transition in mixed phase ferroelectrics remains a challenging task. Here, a methodology is presented to quantitatively evaluate the volume fractions of different domain variants in a mixed phase ferroelectric under complex electromechanical loading. The methodology combines the phenomenology of Landau free energy of ferroelectric phases with Boltzmann statistical analysis and is presented for Pb(Zr,Ti)O₃ near the morphotropic phase boundary (MPB). It is shown that specific grain orientation has a significant effect on how proximity to phase boundary affects microscopic phenomena at the single-crystal level. An estimate of phase and domain switching behavior in a polycrystalline aggregate is subsequently obtained, and the resultant polarization and strain responses at the macroscopic level are computed for a material with random texture. The results indicate the progressive evolution of domain and phase fractions in a material near MPB with mixed ferroelectric phases. We show that in polycrystalline materials with compositions slightly on the tetragonal side of the MPB, grains that exhibit large 90° domain switching have a larger contribution to the macroscopic strain response as compared to grains that undergo tetragonal-to-rhombohedral phase switching.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0258087

## I. INTRODUCTION

Ferroelectrics are widely used in many modern technologies including high-power capacitors, memories, precision actuators and sensors, MEMS, pyroelectric energy harvesters, and electrocaloric cooling.<sup>1-4</sup> In order to undertake a precise design of engineering components for these applications, it is important to develop appropriate modeling tools for predicting the performance of ferroelectric materials under different operating conditions. Although considerable works have been undertaken toward this end in the past, improvement over past approaches is desired for various reasons, as follows. First, as the scope of applications of ferroelectric materials

broadens, they are subjected to more extreme conditions, including higher electric fields, large multiaxial stress fields, or wider temperature variations.<sup>5-9</sup> Second, with increased complexity of engineering components that is facilitated by the recent developments in additive manufacturing technologies, different parts of even the same material component can experience widely different external stimuli, for example, stress or electric-field concentration, at nodes in an architected ferroelectric component.<sup>10,11</sup> The combination of microscopic mechanisms over such extreme conditions can deviate from the more well-known mechanisms based on which the current models are proposed. Third, driven by a desire to achieve higher engineering performance, ferroelectric compositions are often placed in close

<footnote n="1"><sup>1</sup>Laboratoire de Génie Électrique et Électronique de Paris, Université Paris-Saclay, CentraleSupélec, CNRS, Gif-sur-Yvette 91192, France</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Génie Électrique et Électronique de Paris, Sorbonne Université, CNRS, Paris 75252, France</footnote>

<footnote><sup>a)</sup>Authors to whom correspondence should be addressed: abhijit.pramanick@centralesupelec.fr and laurent.daniel@centralesupelec.fr</footnote>
proximity to phase boundaries for both Pb-based Pb(Zr,Ti)O₃ (PZT) as well as Pb-free ferroelectrics, such as solid solutions of BaTiO₃, (K,Na)NbO₃, and (Na,Bi)TiO₃.¹²,¹³ The microscopic mechanisms operational near the phase boundaries can be diverse and may also change with small changes in composition or temperature, which have not been explored fully. All of the above factors necessitate the development of new modeling tools.

The functional properties of ferroelectrics fundamentally arise from their spontaneous electrical polarization, whose state can be changed by the application of an electric field, stress, or temperature variations.¹⁴⁻¹⁶ The development of spontaneous polarization in a ferroelectric material below its Curie temperature ($T_C$) is, in turn, accompanied by its transformation from a high-temperature cubic phase to one or more non-centrosymmetric low-temperature polar crystal phase(s). For temperatures lower than $T_C$, the material can exist in either single or multiple ferroelectric phases as defined by their crystallography symmetries, such as tetragonal, orthorhombic, or rhombohedral. For each of these phases, there can exist multiple domain variants in which the polar vectors are misoriented by either 180° or non-180° depending on specific crystallography symmetry. For example, the polarization vectors in adjacent domains can be misoriented by 90° in the tetragonal phase and 71°/109° in the rhombohedral phase, in addition to the possible misorientation of 180° for both these phases. Under the application of an external stimulus, such as an electric field or stress, the polarization vectors can be reoriented in a reversible or irreversible manner into directions that are allowed by crystallography symmetry, which is defined as domain switching. In cases where multiple phases co-exist, or where it is energetically favorable, the application of an external stimulus may also trigger a broader realignment of polarization vectors through phase transformation. The contributions from domain switching and phase transformation to the functional response of a ferroelectric are generally categorized as the *extrinsic response*. In addition, a part of the macroscopic response of a ferroelectric also comes from lattice strain caused by the piezoelectric effect, which is categorized as the *intrinsic response*. In reality, many of these mechanisms can simultaneously occur, and all of them should ideally be considered while evaluating the macroscopic response of a ferroelectric material. However, for computational ease, in some cases, only the major contributions may be considered to obtain approximate solutions.

From a thermodynamic viewpoint, the overall behavior of a ferroelectric can be divided into a rate-independent reversible component and a rate-dependent irreversible component.¹⁶ Different modeling approaches have been developed to model these components independently. Within the framework of Landau–Ginzburg–Devonshire (LGD) theory, the rate-independent reversible behavior of a single-crystal ferroelectric can be predicted from thermodynamic models in which the Gibbs free energy of a material is expressed as a Taylor series expansion of appropriate order parameters, such as spontaneous polarization or spontaneous strain.¹⁷⁻²⁴ In principle, minimization of the free energy function can then be used to predict the stability of a particular phase or possible distortions within a single phase under different conditions of temperature, stress, or electric field. Note that minimization of the free energy function can result in multiple solutions for minima within the three-dimensional orientation space, which can be

interpreted as different polar domains within a single-crystal material. Similar models based on strain spinoidal were also presented, which can determine the possible number of different strain states under thermodynamic equilibrium.²⁵ In general, these models have not been used to predict the evolution of the volume fractions of different polar states in a polycrystalline material, especially under different conditions of the applied electric field or stress.

In contrast to lattice distortion, domain switching is considered as the major contributor to the rate-dependent irreversible component of the total response of a ferroelectric material. To deal with the rate-dependent irreversible component, mainly, two different approaches have been developed. In one approach, the switching between different domain variants is predicted based on the rate of dissipation of the Helmholtz free energy term due to the changes in remanent strain and polarization states.²⁶⁻³⁰ This approach is based on irreversible thermodynamic models that were originally developed to describe crystal plasticity. A switching criterion can then be developed, which requires that this energy dissipation under external driving forces must be non-negative. An alternate approach is to use phase-field models in which rate-dependent Landau equations are applied to model the evolution of domain and phase microstructures.³¹⁻³³ Both of these approaches can be computationally intensive for polycrystalline samples that take into account all possible grain orientations, especially if one considers the possibility of a wide range of polarization states, such as is the case for co-existing ferroelectric phases. In order to expedite the computation process, Daniel et al. proposed an alternative approach in which the population of different domain variants is calculated based on Boltzmann statistical distribution.³⁴ Inspired by similar formulations for ferromagnetic materials,³⁵,³⁶ Daniel et al. proposed that the relative volume fractions of the different domain variants are related to their Helmholtz free energies, that is, the useful work performed by the applied electric field and/or stress.³⁶ This method allows for the calculation of anhysteretic reversible changes in domain populations under an applied electric field or stress. The rate-dependent irreversible part to the domain switching fraction can then be calculated by secondary phenomenological approaches, such as in Refs. 37–40.

The method described in Ref. 34 is appealing because of its simplicity and computational ease; however, it requires modification for application to more complex scenarios, such as co-existing multiple phases, temperature changes, or induced phase transitions. This is because the Helmholtz free energy terms used in the formulation of Boltzmann statistical distribution in that work does not incorporate the effect of temperature changes or free energy changes due to phase transition. In order to take into effect these later aspects, we propose here an extension of the approach described in Refs. 34–36. Specifically, we propose incorporation of the Landau free energy for the application of Boltzmann statistical law for the computation of the volume fractions of the different domain or phase variants. The use of the Landau free energy, instead of only the Helmholtz free energy, terms, broadens the applicability of the model proposed in Ref. 34 to more complex scenarios for actual materials in service, which often encounter electric fields or stress in combination with temperature changes and can experience multiple microscopic phenomena including domain switching and/or phase transitions.
The objective here is to quantitatively evaluate the volume fractions of domain variants with different polarization orientations in a material with mixed ferroelectric phases. We propose a methodology that combines the phenomenology of Landau free energy of ferroelectric phases with Boltzmann statistical analysis. We present the application of this methodology for the case of the coexisting tetragonal and rhombohedral ferroelectric phases near a morphotropic phase boundary (MPB) of the solid-solution of Pb(Zr,Ti)O₃. The methodology presented is generic and has the potential to be applied for more complex systems, such as more than two coexisting phases or even highly disordered systems such as relaxors. First, we lay down the basic phenomenology for the application of Landau free energy within the framework of Boltzmann statistical distribution to predict the reversible changes in domain and phase fractions in a ferroelectric material near the morphotropic phase boundary (MPB) at the single-crystal level. Second, we apply the methodology to predict the collective behavior of a polycrystalline aggregate with random orientation distribution. The modeling results are compared with the reported experimental measurements of domain and phase switching in PZT ceramics near MPB composition.

## II. BASIC PHENOMENOLOGY

We first examine how the proportion of the various possible polar distortions changes while traversing across a phase transition point. We start with the thermodynamic expression of the Gibbs free energy, ΔG, in the framework of the Landau-Ginsburg-Devonshire (LGD) phenomenological model, under zero electric field and zero stress, and expanded to order 6,¹⁷⁻²⁴

$$ \Delta G = \alpha_1 [P_1^2 + P_2^2 + P_3^2] + \alpha_{11} [P_1^4 + P_2^4 + P_3^4] \\ + \alpha_{12} [P_1^2 P_2^2 + P_2^2 P_3^2 + P_3^2 P_1^2] + \alpha_{111} [P_1^6 + P_2^6 + P_3^6] \\ + \alpha_{112} [P_1^4 (P_2^2 + P_3^2) + P_2^4 (P_3^2 + P_1^2) + P_3^4 (P_1^2 + P_2^2)] \\ + \alpha_{123} P_1^2 P_2^2 P_3^2, \quad (1) $$

where P₁, P₂, and P₃ refer to the components of polarization P of a ferroelectric phase along the three axes of the coordinate system of its parent cubic paraelectric phase. In Eq. (1), coefficients α refer to the respective dielectric stiffnesses at constant dielectric displacement and constant stress. The stability of the various ferroelectric phases can be determined by equating $\frac{\partial\Delta G}{\partial P_i} = 0$, while applying specific symmetry conditions for the different polarization components; for example, P₁ = P₂ = 0 and P₃ ≠ 0 for the tetragonal phase, P₁ = 0 and P₂ = P₃ ≠ 0 for the orthorhombic phase, and so on. Also, in Eq. (1), α₁ is assumed to be temperature dependent, while all other coefficients are considered to be temperature independent.

LGD phenomenological theory has been successfully used to describe composition-dependent phase transition across a morphological phase boundary (MPB), such as in lead zirconate titanate Pb(Zr,Ti)O₃,¹⁹⁻²⁴ as well as temperature-dependent phase transition across a polymorphic phase boundary, such as in BaTiO₃.⁴¹,⁴² However, in many practical cases, one encounters the coexistence of different phases and, therefore, an important objective is to determine how the relative volume fractions of domains with

different polar states change under applied electric fields and/or stress—this is the problem we address below.

## III. COMPUTATION AT THE SINGLE-CRYSTAL LEVEL

We can apply Eq. (1) to evaluate the three-dimensional Landau energy surface for different compositions of Pb(Zr,Ti)O₃. The schematic in Fig. 1(a) illustrates the possible deviation of the polarization vector from the 001 axis, which is defined in terms of angles ψ and θ, where ψ represents the rotation angle from the 100 axis and θ represents the azimuthal angle from the 001 axis. The free energy for each orientation of the polarization vector is calculated following Eq. (1). The values of the coefficients α are obtained from Haun et al.²¹⁻²³ (see the table in Appendix C). Figure 1(b) shows the profile of Landau free energy G (in J/m³) as a function of different Zr/Ti ratios in Pb(Zr,Ti)O₃ solid-solution for the (010) plane, that is φ = 0°. For this plane, the minima in G are observed at θ = 0°, θ = ±90°, and θ = ±180°, which represent the orientations of the polarization vector within the different domain variants of the tetragonal phase. Note that the minima become sharper with increasing Ti ratio, which indicates that the tetragonal polar states become highly stable against possible deviations for a higher Ti content. Figure 1(c) shows the profile of Landau free energy G as a function of different Zr/Ti ratios in the Pb(Zr,Ti)O₃ solid-solution for the (110) plane, that is φ = 45°. For the (110) plane, the minima in G are observed at θ = 0° and θ = ±180°, which represent the 180° domain variants of the tetragonal-like polar states, as well as at θ = ±54.7° and θ = ±125.3°, which represent the orientation of the polarization vectors within the different domain variants of the rhombohedral-like polar states that are oriented along the <111> crystallographic directions. The energy for the minima with tetragonal and rhombohedral-like polar distortions has the same value for a Zr/Ti ratio of 52/48; however, the minima for the rhombohedral polar distortions increase relative to the tetragonal polar distortions with increasing Ti content. For intermediate values of φ between 0° and 45°, the free energy profile shows a gradual transition between these two extremes. Figures 1(d) and 1(e) show the complete energy surfaces for Pb(Zr,Ti)O₃ with reference to the Cartesian coordinate system of the parent cubic phase for compositions of Zr/Ti ratios of 52/48 and 46/54, respectively.

In order to determine the probability of the various orientations of the polarization vector, we introduce an internal variable at the single-crystal level, fₐ. Here, fₐ represents the volume fraction of a family of domains, which has the same orientation of the polarization vector that is defined by a. The value of fₐ under an equilibrium condition can be explicitly calculated following the Boltzmann probability function, Eq. (2), such as described earlier for both ferromagnetic and ferroelectric systems,³⁴⁻³⁶

$$ f_a = \frac{\exp(-A \cdot G_{i=a})}{\sum_{i=1}^{K} \exp(-A \cdot G_i)}. \quad (2) $$

In Eq. (2), Gᵢ is the free energy for a probable orientation of the polarization vector, which is based on one of the K minima in the free energy surface. A is an adjustment parameter, which can be related to the range of compositions over which the phase
<figure><img src="image_17.png" /><figcaption>FIG. 1. (a) Illustration of a modeling setup for the calculation of Landau free energy, G, corresponding to certain orientation of the polarization vector, which is defined by two angles φ and θ. (b) G for polarization vector orientations within the (010) crystallographic plane, as a function of angle to the 001 axis, θ, and for various Zr/Ti ratios in Pb(Zr,Ti)O₃. (c) Same as in (b), for the (110) crystallographic plane. In (b) and (c), the letters T and R in bold refer to the particular polar orientations corresponding to the domains in the tetragonal and rhombohedral phases, respectively, (d) and (e) Landau free energy surface for Zr/Ti ratios of 52/48 and 46/54, respectively. (f) Volume fraction of polar states corresponding to the tetragonal and rhombohedral distortions, as a function of Ti content.</figcaption></figure>

transition occurs as well as the slope of the dielectric susceptibility of a particular phase near the zero electric field, as we will see below. In Eq. (2), the denominator is obtained by the summation of the factor $\exp(-A.G_i)$ over all probable orientations of the polarization vector that are defined by the minima in the free energy surface.

For example, for $\phi = 45^\circ$, there are six minima (counting one of either $+180^\circ$ or $-180^\circ$ orientations), which have the same value of G for the composition Zr/Ti = 52/48. Hypothetically, if all the polarization vectors within the single-crystal are to reside within the crystallographic plane defined by $\phi = 45^\circ$, then the application of Eq. (2) predicts equal volume fractions for each of the six orientations of the polarization vector. Considering the polarization vectors with $\theta = \pm 54.7^\circ$ and $\theta = \pm 125.3^\circ$ to belong to the rhombohedral-like states and $\theta = 0^\circ$ and $\theta = \pm 180^\circ$ to belong to the tetragonal-like states, the relative volume fractions of the rhombohedral and tetragonal phases are 2/3 and 1/3, respectively. However, in order to calculate the relative proportions of all the probable orientations of the polarization vector, we must identify all the local minima from the three-dimensional energy surface, such as shown in Figs. 1(d) and 1(e), and subsequently apply Eq. (2) over all the probable polar states. For the Zr/Ti = 52/48 composition, there are six possible tetragonal-like

polar distortions and eight possible rhombohedral-like polar distortions, all of which have the same minimum energy value. Accordingly, the application of Eq. (2) predicts the volume fraction of the tetragonal domain variants to be equal to 6/14 and that for the rhombohedral domain variants to be equal to 8/14. With increasing Ti content, the minima for the rhombohedral states are higher as compared to that of the tetragonal states, and accordingly, as per Eq. (2), the relative volume fraction of the rhombohedral domains will decrease. The rate of change in volume fractions of the tetragonal and rhombohedral states as a function of composition depends on the value of A. The fractions of the tetragonal and rhombohedral domain variants as a function of composition, calculated for $A = 1 \times 10^{-5} \text{ m}^3/\text{J}$, near the MPB are shown in Fig. 1(f). The exact value of A can be adjusted to match the temperature range over which the coexistence of multiple polar phases is observed experimentally. Alternately, A can be adjusted to match the slope of the dielectric response of the material at low electric fields.³⁴ Ideally, the value of A can be adjusted to consistently satisfy both the above criteria. For the current chosen value of A, the relative volume fraction of the rhombohedral variants asymptotically approaches zero for a Zr/Ti ratio of ~46/54.
A. Application of electric-field

As shown by Damjanovic,¹⁶ the LGD phenomenological model
can be used to predict the susceptibility of polar distortions from the
ideal symmetry of a specific ferroelectric phase under small distur-
bances. For example, in the tetragonal phase, one can characterize the
susceptibility for rotation of the polarization vector from the ideal
polar 001 axis by evaluating the complete energy surface in the cubic
coordinate system and considering the flatness of the energy surface in
the three-dimensional orientation space. Here, we explore the use of
Boltzmann statistics to predict the change in fractions of different
polar states under the application of an electric field.

In order to identify the possible polar distortions under the
application of an external electric field, Eq. (1) can be further
expanded as follows:

$$
\begin{align}
\Delta G ={}& \alpha_1 [P_1^2 + P_2^2 + P_3^2] + \alpha_{11} [P_1^4 + P_2^4 + P_3^4] \nonumber \\
& + \alpha_{12} [P_1^2 P_2^2 + P_2^2 P_3^2 + P_3^2 P_1^2] + \alpha_{111} [P_1^6 + P_2^6 + P_3^6] \nonumber \\
& + \alpha_{112} [P_1^4 (P_2^2 + P_3^2) + P_2^4 (P_3^2 + P_1^2) + P_3^4 (P_1^2 + P_2^2)] \nonumber \\
& + \alpha_{123} P_1^2 P_2^2 P_3^2 - P_1 E_1 - P_2 E_2 - P_3 E_3, \tag{3}
\end{align}
$$

where $E_i$ refers to the electric-field components along the three axes of the coordinate system of the parent cubic paraelectric phase.

Note that in a polycrystalline material, the orientation of the electric field for each single crystallite will be different depending on its relative orientation with respect to the macroscopic electric-field direction. This is illustrated in Fig. 2(a). For Grain 1, the electric field is applied parallel to the 001 crystallographic axis, while for Grain 2, the electric field is applied parallel to the 111

<figure><img src="image_18.png" /><figcaption>FIG. 2. (a) Illustrative example of the orientation of two different grains with respect to the applied electric field. (b) Definition of the electric field direction with respect to the crystallographic directions, as defined by angles $\phi_E$ and $\theta_E$. (c) Illustration of the distortion of the Landau free energy surface for two different directions of an applied electric field, E.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 3. (a) and (b) Line profiles of G for polarization vectors residing within the (010) and (110) atomic planes, respectively, and for different electric-field magnitudes applied parallel to the 001 direction. In (b) and (c), the letters T and R in bold refer to particular polar orientations corresponding to the domains in the tetragonal and rhombohedral phases, respectively. (c) Type of polar domains as defined by the orientation of their polarization vectors with respect to the electric field direction. For example, domains marked R1 have their polarization vector oriented 54.7° with respect to the electric-field direction. (d) and (e) Change in the volume fraction of different domain variants as a function of the applied electric-field magnitude for two different values of the Zr/Ti ratio. (f) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of the applied electric field and for two different values of the Zr/Ti ratio, for the specific orientation of the grain with respect to the electric field shown in Fig. 3(c).</figcaption></figure>
<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

FIG. 4. (a) and (b) Line profiles of G for polarization vectors residing within the (010) and (110) atomic planes, respectively, and for different electric-field magnitudes applied parallel to the 111 direction. In (b) and (c), the letters T and R in bold refer to the particular polar orientations corresponding to the domains in the tetragonal and rhombohedral phases, respectively. (c) Type of polar domains as defined by the orientation of their polarization vectors with respect to the electric field direction. For example, domains marked T1 have their polarization vector oriented 54.7° with respect to the electric-field direction. (d) and (e) Change in the volume fraction of different domain variants as a function of the applied electric-field magnitude for two different values of the Zr/Ti ratio. (f) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of the applied electric field and for two different values of the Zr/Ti ratio.
tetragonal phases as a function of electric field and for two differ-
ent compositions.

The rate of different microscopic processes is expected to influ-
ence the electric-field induced response of a material. For example, as
shown in Fig. 3, the switching between tetragonal and rhombohedral
ferroelectric phases plays a dominant role for the 52/48 composition,
while the switching between different domain variants within the
tetragonal phase, instead, plays a more dominant role for the 50/50
composition. The consequence of this is illustrated in Fig. 5(a), which

shows the different polarization responses of grains with their
001 axes oriented parallel to the electric-field directions, for two
different compositions of 52/48 and 50/50. Notably, grains of this
orientation show a more gradual increase in polarization with
increasing electric field for the 52/48 composition than that for
the 50/50 composition. However, this situation changes for the
grains with their 111 axis parallel to the electric-field direction.
As shown in Fig. 4, for this specific orientation, switching
between the tetragonal and the rhombohedral phase variants is

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>FIG. 5. Polarization response of grains with two different orientations, (a) electric field parallel to the 001 direction, (b) electric field parallel to the 111 direction, plotted for two different values of the Zr/Ti ratio. Strain response of grains with two different orientations, (c) electric field parallel to the 001 direction, (d) electric field parallel to the 111 direction, plotted for two different values of the Zr/Ti ratio.</figcaption></figure>
the predominant mechanism for both 52/48 and 50/50 compositions; however, this process progresses more gradually with increasing electric field in the latter case. The consequence of this is illustrated in Fig. 5(b), which shows for grains with their 111 axis parallel to the electric-field direction; polarization increases more gradually with increasing electric field for the 50/50 composition as compared to the 52/48 composition.

Figure 5(c) shows the corresponding normal strain $\epsilon_{33}$ response of grains with their 001 axes oriented parallel to the electric-field direction, for two different compositions of 52/48 and 50/50. The method for calculating the total spontaneous strain matrix for a mixture of different domain variants is described in Appendix A. The values plotted in Fig. 5(c) indicate the strain for a certain grain relative to its state at zero electric field. Figure 5(d) similarly shows the corresponding normal strain $\epsilon_{33}$ response of grains with their 111 axis oriented parallel to the electric-field direction, for two different compositions of 52/48 and 50/50. Note that for the above calculations, we considered polarization and strain changes as a result of domain or phase changes alone. For clarity, the induced polarization and strain responses due to linear dielectric and piezoelectric effects are not included, which are, at least, one order of magnitude lower than the corresponding effects due to the domain/phase changes. Switching among the tetragonal domain variants produces a larger strain as compared to the switching among the rhombohedral domain variants, due to a lower value of spontaneous strain $S_0$ for the latter (see Appendix A). A comparison of the results shown in Figs. 5(c) and 5(d) indicates that the $\epsilon_{33}$ strain response measured along the 001 axis is larger for the 50/50 composition, while the $\epsilon_{33}$ strain response measured along the 111 axis is larger for the 52/48 composition.

## B. Application of mechanical stress

In many different applications, ferroelectric ceramics are also subjected to mechanical stress, either solely or in conjunction with electric field.43–45 For example, actuator components are likely to experience some pre-existing mechanical loads, even if unintentional. In addition, intergranular stresses are known to exist in polycrystalline ceramics.46 Therefore, it will be interesting to know how much such stresses affect the volume fractions of the different domain/phase variants.

In order to identify the possible polar distortions under the application of an external stress, Eq. (1) can be further expanded as follows:

$$ \begin{aligned} \Delta G ={}& \alpha_1 [P_1^2 + P_2^2 + P_3^2] + \alpha_{11} [P_1^4 + P_2^4 + P_3^4] \\ & + \alpha_{12} [P_1^2 P_2^2 + P_2^2 P_3^2 + P_3^2 P_1^2] + \alpha_{111} [P_1^6 + P_2^6 + P_3^6] \\ & + \alpha_{112} [P_1^4 (P_2^2 + P_3^2) + P_2^4 (P_3^2 + P_1^2) + P_3^4 (P_1^2 + P_2^2)] \\ & + \alpha_{123} P_1^2 P_2^2 P_3^2 - \frac{1}{2} S_{11} [X_1^2 + X_2^2 + X_3^2] \\ & - S_{12} [X_1 X_2 + X_2 X_3 + X_3 X_1] - \frac{1}{2} S_{44} [X_4^2 + X_5^2 + X_6^2] \\ & - Q_{11} [X_1 P_1^2 + X_2 P_2^2 + X_3 P_3^2] \\ & - Q_{12} [X_1 (P_2^2 + P_3^2) + X_2 (P_3^2 + P_1^2) + X_3 (P_1^2 + P_2^2)] \\ & - Q_{44} [X_4 P_2 P_3 + X_5 P_3 P_1 + X_6 P_1 P_2], \end{aligned} \quad (4) $$

where $X_1, \dots, X_6$ refer to the components of the stress tensor as $X_1 = \sigma_{11}$, $X_2 = \sigma_{22}$, $X_3 = \sigma_{33}$, $X_4 = \sigma_{23}$, $X_5 = \sigma_{13}$, and $X_6 = \sigma_{12}$, $S_{ij}$ refers to the elastic compliances at constant polarization, $Q_{ij}$ refers to the electrostrictive coupling between ferroelectric polarization and stress, and all other variables have the same meaning as described for Eq. (1). The volume fractions of the different domain/phase variants are then calculated based on the value of $\Delta G$ under applied stress and by applying Eq. (2) in the same manner as described in Sec. III A.

We evaluated the effect of normal stress, $X_3$, applied parallel to either 001 or 111 crystal axes, and for two different compositions 52/48 and 50/50, as illustrated in Figs. 6(a) and 6(d). Figure 6(b) shows changes in the volume fractions of the tetragonal and rhombohedral domain variants under tensile and compressive stresses. It shows that the application of either tensile or compressive stress favors an increase in the tetragonal domain variants at the expense of rhombohedral domain variants. Similarly, Fig. 6(e) shows that the application of normal tensile or compressive stress, $X_3$, along the 111 crystal axis favors rhombohedral domain variants at the expense of tetragonal domain variants. In both cases, there exists asymmetry whereby phase transition is stronger under tension than the same under compression. The resulting strains arising from ferroelastic domain-switching/phase transition under tensile/compressive stress applied along the 001 and 111 crystal axes are shown in Figs. 6(c) and 6(f), respectively. Note that linear elastic strain resulting from the application of stress is not included in this calculation for clarity, but they can be superimposed to the corresponding ferroelastic and phase transition strains. For stress applied parallel to the 001 crystal axis, the ferroelastic strains under both tensile and compressive stress are larger for the composition 50/50 than the same for composition 52/48. In comparison, for stress applied parallel to the 111 crystal axis, the strains under both tensile and compressive stresses are larger for the composition 52/48 than the same for composition 50/50. In addition, the saturation strain is also larger for tension as compared to that for compression. This reflects the fact that for the stress applied along 001 and for a positive value of spontaneous strain, a maximum of 2/3 of the tetragonal non-180° domain variants can be reoriented under tensile stress, while a maximum of 1/3 of the tetragonal non-180° domain variants can be reoriented under compressive stress, relative to the zero stress state. Similarly, for the stress applied along 111 and for a positive value of spontaneous strain, a maximum of 3/4 of the rhombohedral non-180° domain variants can be reoriented under tensile stress, while a maximum of 1/4 of the rhombohedral non-180° domain variants can be reoriented under compressive stress, relative to the zero stress state.

## C. Simultaneous application of mechanical stress and electric field

In reality, ferroelectric devices often experience a combination of mechanical loads and electric field.43–45 Even when only an electric field is applied to a ferroelectric ceramic component, differently oriented grains can experience intergranular stresses.46 It is, therefore, of interest to examine how domain switching and phase transition proceed under such a combined electromechanical loading scenario. We examine here a particular case of application of the incremental electric field under a constant compressive stress for
<figure><img src="image_35.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_40.png" /></figure>

FIG. 6. (a) Illustration of loading scenario for a grain with its 001 crystal axis parallel to that of the application of a tensile/compression stress. (b) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of applied stress along 001 and for two different values of the Zr/Ti ratio. (c) Resultant strain as a result of domain switching and phase transitions for the loading scenario illustrated in (a). (d) Illustration of the loading scenario for a grain with its 111 crystal axis parallel to that of the application of a tensile/compression stress. (e) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of applied stress and for two different values of the Zr/Ti ratio. (f) Resultant strain as a result of domain switching and phase transitions for the loading scenario illustrated in (d).
the MPB composition 52/48, and for two different crystal axes, as depicted in Figs. 7(a) and 7(e). For the evaluation of the Landau free energy under combined electromechanical loading, we combine Eqs. (2) and (3). Thereafter, the procedure followed to evaluate the equilibrium volume fractions of the different domain variants is the same as described above.

For the simultaneous application of mechanical stress and electric field along the 001 crystal axes, the change in the volume fraction of the tetragonal and rhombohedral phases is shown in Fig. 7(b). Under increasing magnitudes of compressive stress and for electric-field amplitudes greater than 1 kV/mm, the volume fraction of the tetragonal phase decreases, while that of the rhombohedral phase increases. Interestingly, the volume fraction of the tetragonal (rhombohedral) phase initially decreases (increases) with the application of smaller electric fields, before eventually increasing (decreasing) under higher electric fields. The resultant polarization and strain responses are shown in Figs. 7(c) and 7(d), respectively. The electric-field-induced polarization is lower for increasing magnitudes of compressive stress under intermediate electric-field values, although saturation polarization is approached under higher electric fields. A similar behavior is observed for the induced strain due to combined domain-switching and phase transition mechanisms, as shown in Fig. 7(d). Note that the apparent strain is slightly larger at higher electric fields under compressive stress, which is due to the fact that the overall strain is computed relative to the initial state and the material is initially compressed under zero electric field.

Figures 7(f)–7(h) show the changes in phase volume fractions, as well as polarization and strain responses under the simultaneous application of compressive stress and electric field parallel to the 111 crystal axis. In this case, as shown in Fig. 7(f), the volume fraction of the tetragonal phase increases, while that of the rhombohedral phase decreases under increasing magnitudes of compressive stress. Additionally, both electric-field-induced polarization and strain values are lower under intermediate electric-field values and for increasing magnitudes of compressive stress, as shown in Figs. 7(g) and 7(h); however, saturation is approached in both cases under higher electric fields. Note that the effect of a compressive stress on electric-field-induced polarization/strain is relatively weaker for the situation depicted in Fig. 7(b), as compared to the same for what is shown in Fig. 7(a). This is caused by how the differences between the energy minima of the different domain variants are relatively affected under different electromechanical loading scenarios.

### D. Generalization to all orientation states

At this point, we can generalize the methodology to compute the extent of phase transition for any relative misorientation between the crystal axes and the axes of the applied electric field and/or stress. A convenient way to represent the results is through the use of inverse pole figures, as shown in Fig. 8. The stereographic projections are made with respect to the (001) cubic axes as indicated in Fig. 8. Various points on the stereographic projection indicate the directions of the applied electric field or the axis of the applied uniaxial stress with respect to the crystallographic pseudo-cubic 001 crystal axes.

The top two panels in Fig. 8(a) indicate changes in the tetragonal and rhombohedral volume fractions for different directions of the applied electric fields and for different field magnitudes. The changes in phase volume fractions become more prominent for applied electric fields higher than 1 kV/mm. This is consistent with the general experimental observations, such as in Ref. 47. The plots indicate the different regions within the orientation space that favors an increase/decrease in either the tetragonal or rhombohedral phase volume fractions. The tetragonal phase fraction increases for electric field directions closer to the (001) crystal axes, while it decreases for electric field directions closer to the (111) crystal axes. Conversely, the rhombohedral phase fraction increases for electric field directions closer to the (111) crystal axes, while it decreases for electric field directions closer to the (001) crystal axes.

Figure 8(b) indicates the extent of changes in relative phase fractions for applied tensile/compressive stresses and for the combined application of an electric field and stress. In the case of applied compressive stress, the separation between the regions with preference for either the tetragonal or the rhombohedral phase is more diffused as compared to what is observed for an applied electric field or tensile stress alone. In the case that an electric field and stress are applied simultaneously, a convolution of the effects observed for each of these cases can be expected, as depicted in Fig. 8(c). For the simultaneous application of a positive electric field and compressive stress, the boundaries between the regions with different phase preferences are very much diffused and the rhombohedral phase has a higher preference, in general. In contrast, for simultaneously applied positive electric field and tensile stress, the different regions in the orientation space with preference for tetragonal and rhombohedral phases are more clearly defined.

### IV. APPLICATION TO POLYCRYSTALLINE AGGREGATES

Once the behavior of single crystallites for all different orientations (with respect to macroscopic electric field and/or stress) are determined, we can proceed to estimate the functional response of a polycrystalline aggregate by averaging over all orientation states. In a simple approach, the electrics field and stress can be assumed to be uniformly applied to each single crystallite in a polycrystalline ensemble. In this case, one can obtain the response of the polycrystalline aggregate by simply averaging over all the orientation states, such as those adopted in Refs. 48 and 49. In a further refinement, one can take into account the grain-to-grain interactions and re-optimize the level of the localized electric field and stress that are internally experienced by each single crystallite; such models have been used, for example, in Refs. 28, 34, and 50. Here, we use the former approach, that is, we consider the electric field and stress to be uniform for all single-crystal grains, as a first-approximation to characterize the progressive evolution of domain and phase volume fractions in a material. The effect of field and stress localization will be provided in a subsequent article.

#### A. Domain and phase switching

First, we simulate the microstructure of a polycrystalline ceramic with 20 000 grains with a random texture. This is done by assigning random values to the three Euler angles $\phi_1, \psi$, and $\phi_2$, as described in Appendix B. The computed results are
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>FIG. 7. (a) Illustration of the loading scenario for a grain with its 001 crystal axis parallel to that of the application of simultaneous electric field and compressive stress. (b) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of incremental electric field and constant compressive stress applied along the 001 crystal axis. (c) and (d) Resultant polarization and strain as a result of domain switching and phase transitions for the loading scenario illustrated in (a). (e) Illustration of the loading scenario for a grain with its 111 crystal axis parallel to that of the application of simultaneous electric field and compressive stress. (f) Change in the overall volume fractions of the tetragonal and rhombohedral domains as a function of incremental electric field and constant compressive stress applied along the 111 crystal axis. (g) and (h) Resultant polarization and strain as a result of domain switching and phase transitions for the loading scenario illustrated in (e).</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>FIG. 8. Inverse pole figures showing the degree of phase change for the applied electric field and/or uniaxial stress parallel to different crystallographic directions. All axes refer to the pseudo-cubic crystal lattice. Each point in the inverse pole figures refer to the electric-field direction with respect to the crystal axes. (a) Change in tetragonal and rhombohedral phase fractions for different electric field directions and magnitudes. (b) Same as in (a), for applied uniaxial tensile and compressive stresses. (c) Same as in (a), for the simultaneous application of the co-aligned electric field and uniaxial stress.</figcaption></figure>

presented by grouping different orientations in terms of their azimuthal angles $\psi$ into bins of 10° intervals. This is done to correlate the modeling results with the results of *in situ* x-ray diffraction measurements, which typically follow this convention, such as in Refs. 46 and 51–53.

Figure 9 shows the changes in volume fractions of the c- and a-domains of the tetragonal phase, as well as the rhombohedral phase volume fraction, for various azimuthal ($\psi$) angle bins, for composition with a Zr/Ti ratio of 52/48. The various azimuthal angular bins with respect to the electric field direction are
<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /><figcaption>FIG. 9. (a) Definition of different azimuthal angle bins with respect to the electric-field direction. Note that the material response is considered to be transversely isotropic with respect to the electric field direction. (b)–(j) Electric field dependent changes in volume fractions of the a- and c-domains of the tetragonal phase and the rhombohedral phase, for polycrystalline PZT with Zr/Ti = 52/48. (k) $\Delta\eta$, which signifies the extent of 90° domain switching in the tetragonal phase, is presented for various azimuthal angular bins. (l) $\Delta R$, which signifies the extent of tetragonal-to-rhombohedral phase transition, is presented for various azimuthal angular bins.</figcaption></figure>
indicated in Fig. 9(a). Within the azimuthal angular range of 0°–40°, the volume fraction of the c-domains increases while that of the a-domains decreases; the opposite is observed for the angular range of 70°–90°. For the azimuthal angular range of 30°–70°, the relative volume fraction of the a- and c-domains exhibit relatively smaller changes. Additionally, for the angular range of 30°–70°, an increase in the rhombohedral phase volume fraction with an applied electric field can be observed, especially for electric fields higher than 1 kV/mm. In order to quantitatively represent the microstructural changes, we further compute the following two quantities:

$$ \eta = \frac{V_{002}}{V_{002} + V_{200}}, \qquad (5) $$

$$ R = \frac{V_R}{V_T + V_R}, \qquad (6) $$

where $V_{002}$ and $V_{200}$ represent the c- and a-domain volume fractions of the tetragonal phase, respectively; $V_T$ and $V_R$ represent the volume fractions of the tetragonal and rhombohedral phases, respectively. $\Delta\eta$ and $\Delta R$ represent changes in the values of these variables, with respect to their original values at zero applied electric fields. In effect, $\Delta\eta$ represents the degree of 90° domain switching and $\Delta R$ represents tetragonal ↔ rhombohedral phase switching under the application of an electric field and/or stress. Figure 9(k) shows $\Delta\eta$ as a function of $\psi$ for various magnitudes of applied electric fields. The dotted line indicates the baseline of $\Delta\eta = 0$ for zero applied electric fields. The computed angular distribution for $\Delta\eta$ is consistent with that of the experimental measurements reported for PZT ceramics near MPB composition.46,51,52

Figure 9(l) similarly shows the angular distribution of $\Delta R$ for various applied electric field magnitudes. The results provide a phenomenological understanding of the experimental reports of electric-field-induced phase transformations for materials near MPB. For example, Fig. 10(a) shows changes in the {002} diffraction peak profile in PZT ceramics with composition near the MPB under an applied electric field along the direction that is rotated 45° from the 002 diffraction wavevector.47 This is made clear by the illustration on the right (added by us), which indicates that the 002 diffraction peak in Fig. 10(a) arises from specific sets of grain families for which the pseudocubic 001 crystal orientation is rotated 45° with respect to the electric-field direction, without any distinction for the crystal orientation transverse to the diffraction wavevector $k$. The result in Fig. 10(a) is correlated to the predicted response shown in Fig. 9(f), which shows an increase in the rhombohedral phase fraction with a corresponding decrease in the tetragonal phase fraction. Similar results were also reported for the solid-solution system of PbTiO$_3$-BiScO$_3$,54 which exhibits an MPB similar to PZT. Note that, in our modeling, we have considered the polarization states near (111) directions as “rhombohedral;” however, the flatness of the energy profile means that small deviations from the ideal rhombohedral states with (111) polarization directions can occur, which can be interpreted as “monoclinic,” such as reported in Refs. 47 and 54. For $\psi = 0°–30°$, our modeling indicates large negative values of $\Delta R$, that is, a rhombohedral-to-tetragonal phase

<figure><img src="image_62.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>(b)</figcaption></figure>

FIG. 10. (a) Change in the 002 diffraction peak profile from a PZT ceramic of MPB composition under the application of an electric field, $E$, which is rotated 45° from 002 diffraction wavevector; reproduced with permission from Hinterstein et al., Phys. Rev. Lett. 107, 077602 (2011). Copyright 2011 American Physical Society. (b) Change in the 002 diffraction peak profile from a PZT ceramic of MPB composition under the application of an electric field, $E$, which is parallel to the 002 diffraction wavevector; reproduced with permission from Liu et al., J. Am. Ceram. Soc. 88, 210–215 (2005). Copyright 2005 American Ceramic Society. The crystallographic illustrations are added by us to indicate the diffraction geometry. 001$_{pc}$ refers to the pseudocubic 001 crystal axes, while $k$ indicates the diffraction wavevector.

transition for grains of these orientations. This is also validated by the results from in situ XRD experiments, such as reported by Ref. 55 and shown in Fig. 10(b). Note that, for this case, the 002 diffraction peak in Fig. 10(b) arises from specific sets of grain families, for which the pseudocubic 001 crystal orientations are
approximately parallel with respect to the electric-field direction, without any distinction for the crystal orientation transverse to the diffraction wavevector **k**, such as illustrated in the bottom panel (added by us). As shown in Fig. 10(b), for a PZT ceramic of MPB composition, the rhombohedral domain variants have a larger contribution to the overall 002 diffraction peak intensity at zero electric field, as compared to that under electric field **E** that is applied approximately parallel to the 001 crystallographic poles. These results can be compared to the predicted response shown in Fig. 9(b), which shows an increase in the tetragonal phase fraction with a corresponding decrease in the rhombohedral phase fraction.

The domain and phase switching behavior in polycrystalline PZT with a Zr/Ti ratio of 50/50 is shown in Fig. 11. For Zr/Ti = 50/50, the volume fraction of the rhombohedral phase is much lower as compared to that observed for Zr/Ti = 52/48. However, the angular distribution of $\Delta\eta$ observed for Zr/Ti = 50/50 is similar to what is observed for Zr/Ti = 52/48. This is not surprising since the *change* in free energies of the polarization states for the *a*- and *c*-domains remains unchanged with composition, although the absolute values of their volume fractions are changed. In addition, we can note from Fig. 11(l) that the extent of rhombohedral-to-tetragonal phase transition is reduced as compared to Zr/Ti = 52/48, while the angular range for tetragonal-to-rhombohedral phase transition is broadened.

The domain and phase switching behavior can be also subject to change under the simultaneous application of an electric field and stress. An example of this is illustrated in Fig. 12, which shows electric-field-induced evolution of various domain volume fractions for different azimuthal angular ranges under the simultaneous application of an electric field and a constant compressive stress of -10 MPa. Interestingly, for $\psi = 0°-30°$ and $\psi = 70°-90°$, one can observe an initial increase in the rhombohedral phase fraction for electric-field magnitudes lower than 1 kV/mm, followed by a subsequent decrease of the same under higher electric fields. The angular distribution of $\Delta\eta$ and $\Delta R$ are shown in Figs. 12(k) and 12(l), respectively. It shows that $\Delta\eta$ or domain switching is initially suppressed for field magnitudes below 1 kV/mm, but then accelerates under higher electric fields. Also, for angles closer to the electric-field direction, there is an increase in the rhombohedral phase fraction for electric fields lower than 1 kV/mm, which is subsequently reversed for higher electric-field magnitudes.

## B. Contribution to macroscopic response

The modeling results presented above correlates with the earlier experimental reports and indicate that grains with different orientations within a polycrystalline aggregate can exhibit different domain/phase switching behaviors depending on their specific orientations with respect to the macroscopic field/stress direction. Depending on specific grain orientations, the dominant microscopic mechanism can be either 90° domain switching or rhombohedral-to-tetragonal or tetragonal-to-rhombohedral phase transition. Therefore, a pertinent question is which one of these mechanisms is the most significant contributor to the macroscopic functional response of a polycrystalline ceramic. In order to answer this, we

next calculate the relative contributions of grains with different orientations to the macroscopic material response, based on the results of $\Delta\eta$ and $\Delta R$, such as presented in Figs. 9, 11, and 12. Once the values of $\Delta\eta$ and $\Delta R$ are calculated for a specific grain orientation, the corresponding polarization and strain contributions from domain/phase switching are obtained through a tensorial transformation as explained in Appendixes A and B.

Figure 13 shows the contributions to the macroscopic polarization and strain from single-crystal grains of different orientations, which are defined with respect to the range of azimuthal $\psi$ angles, as defined in Appendix B. It is convenient to define the relative contributions from grains of different orientations based on their relative azimuthal orientation angle $\psi$, since a polycrystalline ceramic can be assumed to be transversely isotropic with respect to an applied electric field or uniaxial stress direction. The azimuthal sectors represented in Fig. 9 also correlate well with the convention typically used for in situ x-ray diffraction measurements, such as in Ref. 51. Figures 13(a) and 13(b) compare the polarization and strain contributions from grains of different azimuthal orientations for the two different compositions of Zr/Ti = 52/48 and Zr/Ti = 50/50. Figure 13(c) compares the contributions to polarization and strain from grains of different azimuthal orientations for the composition 52/48 under applied electric fields and a constant applied stress of -10 MPa.

The macroscopic polarization and strain curves shown in Fig. 13 can be compared with the corresponding, experimentally determined *anhysteretic* curves, such as described in Refs. 56 and 57. The *anhysteretic* curve refers to the predicted material response in the stable equilibrium state for certain values of the applied electric field and/or stress, without the effects of dissipative microscopic phenomena such as defect pinning of domain or phase boundary motion. Further details for the determination of anhysteretic polarization/strain curves can be found in Refs. 56 and 57. The computed values for macroscopic polarization and strain (Fig. 13) match reasonably well with those of anhysteretic curves measured for commercial PZT ceramics with composition near the MPB (52/48 ratio for Zr/Ti) (Ref. 57), as shown in Fig. 14. The experimentally measured maximum electric-field-induced anhysteretic strain is slightly lower than that obtained from the current model, which is likely due to the effect of intergranular constraints.³⁴,⁵² The effect of intergranular constraints on macroscopic response will be considered in future studies.

Furthermore, we can separate the relative contributions from grains with different orientations into three separate groups: (a) group I for $\psi = 0°-30°$ (~13% of grains for random texture) that exhibits positive(negative) values for $\Delta\eta$ ($\Delta R$), (b) group II for $\psi = 30°-70°$ (~52% of grains for random texture) that exhibits positive values for $\Delta R$ and relatively smaller $\Delta\eta$, and (c) group III for $\psi = 70°-90°$ (~35% of grains for random texture) that exhibits negative values for both $\Delta\eta$ and $\Delta R$. Groups I and III represent grain orientations in which 90° domain switching and rhombohedral-to-tetragonal phase transition are the dominant microscopic mechanisms. Group II represents grain orientations in which tetragonal-to-rhombohedral is the dominant microscopic mechanism. In all cases, the relative contributions to macroscopic polarization from grains of Groups I, II, and III roughly scale with their relative volume fractions. However, it is not the case for their strain contributions.
<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /></figure>

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /><figcaption>FIG. 11. Same as in Fig. 9, but for a Zr/Ti ratio equal to 50/50.</figcaption></figure>

The relative contributions to strain from grains of Groups I, II, and III are depicted as a function of electric field in Fig. 15. The contributions to macroscopic strain from grains of Group II for small electric fields are ~0.41 and ~0.32 for compositions of

52/48 and 50/50, respectively; however, both these values approach ~0.37 for an electric field of 2 kV/mm. This provides an estimate of the relative contributions to the overall strain from grains in which the tetragonal-to-rhombohedral phase transition constitute a major
<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

<figure><img src="image_77.png" /></figure>

<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>

<figure><img src="image_84.png" /></figure>

<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /><figcaption>FIG. 12. Same as in Fig. 9, but for the simultaneous application of an electric field and a constant compressive stress of -10 MPa.</figcaption></figure>

electric-field-induced microscopic phenomenon. The remaining
contribution comes from grains in which 90° domain switching
and rhombohedral-to-tetragonal phase transition are the dominant
microscopic mechanisms, that is Group I and Group II. It is notable

that the relative contribution from Group II grains increases under a
compressive stress as compared to the zero stress condition; however,
it also decreases with increasing electric fields. Overall, for PZT
ceramics of compositions slightly on the tetragonal side of the MPB,
<figure><img src="image_87.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_91.png" /></figure>

<figure><img src="image_92.png" /></figure>

FIG. 13. Contributions to macroscopic electric-field-induced polarization and strain due to phase and domain switching phenomena from grains with different azimuthal orientations with respect to the electric field or stress direction, as defined in Fig. 9. (a) For Zr/Ti = 52/48 under applied electric fields, (b) for Zr/Ti = 50/50 under applied electric fields, and (c) for Zr/Ti = 52/48 under applied electric fields and constant compressive stress of –10 MPa.
<figure><img src="image_93.png" /><figcaption>**FIG. 14.** Experimentally determined anhysteretic curves for electric-field-induced (a) polarization, and (b) strain for commercial PZT with composition near the MPB, from Ref. 57. Continuous curves denote the hysteretic polarization and strain responses, while the corresponding anhysteretic responses are indicated by blue dots.</figcaption></figure>

we can estimate that (a) grains with net tetragonal-to-rhombohedral phase transition (~52% of grains for random texture) contribute roughly 30–40% of the overall strain, and (b) grains with net rhombohedral-to-tetragonal phase transition and large 90° domain switching (~48% of grains for random texture) contribute roughly 60%–70% of the overall strain.

<figure><img src="image_94.png" /><figcaption>**FIG. 15.** Comparison of the relative contributions to macroscopic strain from grains with different azimuthal orientation ranges for (a) Zr/Ti ratios of 52/48 and 50/50, (b) a Zr/Ti ratio of 52/48 under combined applied electric fields and constant compressive stress.</figcaption></figure>

V. CONCLUDING REMARKS

In summary, we presented here a generic model based on a combination of LGD phenomenological theory and Boltzmann statistics, which can be used to predict domain switching and/or phase transition in complex ferroelectrics with multiple coexisting polar phases. The methodology allows for the examination of material behavior under different loading conditions, such as electric field and/or stress, as well as taking into account the relative grain orientations. Additionally, it is possible to incorporate more complex scenarios such as temperature variations, as long as the temperature dependence of the Landau coefficients is known.

The results presented here are noteworthy in underlining the significance of specific grain orientations while discussing the relative importance of either domain switching or phase transition phenomenon toward the functional response of a ferroelectric material. The results indicate that generalizations relating enhanced properties in the ferroelectrics of certain composition to solely domain switching or phase transition need to be properly examined. In order to interpret the macroscopic response of a polycrystalline material, it is necessary to consider grains of all possible orientations and weighing
their relative contributions to the macroscopic behavior. The methodology for the same is provided here in the framework of a polycrystalline aggregate with random texture, which allows for the estimation of relative contributions from the families of grains that undergo different degrees of phase transition and domain switching phenomena. An important result from our model is that, for PZT composition near the MPB, grain families that exhibit a larger 90° domain switching and a net rhombohedral-to-tetragonal phase switching have a larger contribution to macroscopic electric-field-induced strain, as compared to the grains that exhibit a net tetragonal-to-rhombohedral phase switching.

One limitation of the current model is that the effect of intergranular constraints on polarization and strain responses has not been included. This will be important for comparative studies of material behavior over a broader range of crystallographic distortions or to consider the effects of different grain sizes. A possible means for the inclusion of intergranular interaction effects has been described in Ref. 34; the application of which to the current model will be explored in future studies. In addition, the expansion of the current model to include hysteresis effects will also be undertaken in the future.

## ACKNOWLEDGMENTS

Support from CentraleSupelec and University Paris-Saclay is gratefully acknowledged. This work has been supported as part of France 2030 program ANR-11-IDEX-0003.

<figure><img src="image_95.png" /><figcaption>FIG. 16. Definition of the polarization directions for (a) tetragonal domains, (b) rhombohedral domains, with respect to pseudocubic crystallographic coordinates.</figcaption></figure>

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**A. Pramanick:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **L. Daniel:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

## APPENDIX A: CALCULATION OF DOMAIN-SWITCHING STRAINS

In a tetragonal crystal, the six possible polarization directions are along the (001) pseudo-cubic crystal axes, which are indicated in Fig. 16(a). The spontaneous strains corresponding to each of these six possible polarization directions are as follows:³⁸
$$ \epsilon^{T1} = \epsilon^{T2} = \frac{S_0}{3} \begin{bmatrix} -1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 2 \end{bmatrix}, $$

$$ \epsilon^{T3} = \epsilon^{T4} = \frac{S_0}{3} \begin{bmatrix} 2 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{bmatrix}, $$

$$ \epsilon^{T5} = \epsilon^{T6} = \frac{S_0}{3} \begin{bmatrix} -1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & -1 \end{bmatrix}, \qquad (A1) $$

where $S_0 = (\frac{d_{001}}{d_{100}} - 1)$, $d_{001}$, and $d_{100}$ correspond to 001 and 100 plane spacings, respectively. For tetragonal domains, we used $S_0 = 2.15\%$. T1 and T2 are referred to as c-domains, while T3, T4, T5 and T6 are referred to as a-domains.

In a rhombohedral crystal, the eight possible polarization directions are along the (111) pseudo-cubic crystal axes, which are indicated in Fig. 16(b). The spontaneous strain corresponding to each of these eight possible polarization directions are as follows:⁵⁹

$$ \epsilon^{R1} = \epsilon^{R2} = \frac{S_0}{3} \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix}, $$

$$ \epsilon^{R3} = \epsilon^{R4} = \frac{S_0}{3} \begin{bmatrix} 0 & -1 & 1 \\ -1 & 0 & -1 \\ 1 & -1 & 0 \end{bmatrix}, $$

$$ \epsilon^{R5} = \epsilon^{R6} = \frac{S_0}{3} \begin{bmatrix} 0 & -1 & 1 \\ -1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix}, $$

$$ \epsilon^{R7} = \epsilon^{R8} = \frac{S_0}{3} \begin{bmatrix} 0 & 1 & -1 \\ 1 & 0 & -1 \\ 1 & -1 & 0 \end{bmatrix}, \qquad (A2) $$

where $S_0 = (\frac{9}{8})(\frac{d_{111}}{d_{111}} - 1)$, $d_{111}$, and $d_{111}$ correspond to 111 and 111 plane spacings, respectively. For rhombohedral domains, we used $S_0 = 0.63\%$.

For a given condition, the total spontaneous strain matrix of a crystallite is obtained by the summation of the spontaneous strain matrices of domains corresponding to each of these fourteen different polarization states (six for tetragonal domains and eight for rhombohedral domains), weighted by their respective volume fractions.

## APPENDIX B: TENSORIAL TRANSFORMATION FOR POLYCRYSTALLINE ENSEMBLE

A polycrystalline microstructure is simulated by introducing relative orientations between the crystallographic axes (or the LGD frame) and the macroscopic frame that is used to define the direction of the applied electric field/stress. To accomplish the above, we define the three Eulerian angles that define the mutual orientation between the two set of axes: (1) one that is defined by the pseudocubic crystallographic axes (or the LGD frame), and (2) the macroscopic frame that is convenient to define the applied electric field/stress or induced polarization/strain response. The three Eulerian angles are defined as follows:

In the illustration shown in Fig. 17, the crystal axes are kept fixed while the macroscopic frame is rotated with respect to the fixed crystal axes. The first rotation is a counterclockwise rotation around axis $Z_3$ by angle $\phi_1$, which is defined by the transformation matrix

$$ a_I = \begin{bmatrix} \cos \phi_1 & \sin \phi_1 & 0 \\ -\sin \phi_1 & \cos \phi_1 & 0 \\ 0 & 0 & 1 \end{bmatrix}. $$

<figure><img src="image_96.png" /><figcaption>FIG. 17. Definition of Euler angles that define the mutual orientation between the crystal axes (fixed) and the macroscopic sample coordinates.</figcaption></figure>
The second rotation is a counterclockwise rotation around a new axis Z'₁ by angle ψ, which is defined by the transformation matrix

$$a_{II} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \psi & \sin \psi \\ 0 & -\sin \psi & \cos \psi \end{bmatrix}.$$

The third rotation is a counterclockwise rotation around a new axis Z''₃ by angle φ₂, which is defined by the transformation matrix

$$a_{III} = \begin{bmatrix} \cos \phi_2 & \sin \phi_2 & 0 \\ -\sin \phi_2 & \cos \phi_2 & 0 \\ 0 & 0 & 1 \end{bmatrix}.$$

The product of the above three rotations is given by

$$a = \begin{bmatrix} \cos \phi_2 & \sin \phi_2 & 0 \\ -\sin \phi_2 & \cos \phi_2 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \psi & \sin \psi \\ 0 & -\sin \psi & \cos \psi \end{bmatrix} \times \begin{bmatrix} \cos \phi_1 & \sin \phi_1 & 0 \\ -\sin \phi_1 & \cos \phi_1 & 0 \\ 0 & 0 & 1 \end{bmatrix}.$$

The relationships defining the electric field (E)/polarization (P)/stress (σ)/strain (ε) state in the coordinate system of the crystal axes (LGD frame) and the electric field (E'), polarization (P')/stress (σ')/strain (ε') state in the rotated coordinate system, as depicted in the above figure, are given by

$$E' = aE,$$

$$P' = aP,$$

$$E = a^T E',$$

$$P = a^T P',$$

$$\sigma' = a\sigma a^T,$$

$$\epsilon' = a\epsilon a^T,$$

$$\sigma = a^T \sigma' a,$$

$$\epsilon = a^T \epsilon' a.$$

<figure><img src="image_97.png" /><figcaption>FIG. 18. Inverse pole figure depicting the relative orientation between the crystal axes and the orientation of the applied electric field.</figcaption></figure>

The above relationships are reversed if the crystal axes are rotated with respect to a fixed macroscopic frame. For convenience, one can choose whether to fix the orientation of the crystal axes (such as for the calculation of inverse pole figures, viz., Fig. 8) or to fix the macroscopic frame (such as for the calculation of material response, viz., Figs. 9, 11, and 12). Both approaches provide identical results.

The inverse pole figure in Fig. 18 depicts the relative orientation between the crystal axes and the macroscopic coordinates for a simulated microstructure with 20 000 grains. Each dot represents the angular orientation of the applied electric field with respect to the 〈100〉 crystal axes.

The simulated random texture is checked for isotropy by computing the polarization responses for different directions of the applied electric fields for a polycrystalline material, as shown in Fig. 19.

<figure><img src="image_98.png" /><figcaption>FIG. 19. Polarization response of a polycrystalline material with random texture for different directions of collinear electric fields. The maximum deviation in the overall polarization response is within ±0.5%, which confirms the isotropic characteristic of the simulated random texture.</figcaption></figure>
APPENDIX C: MATERIAL PARAMETERS

The value of the Landau coefficients for certain compositions of PZT is given below,²¹ which can be used to compute stiffness values for intermediate values of the Zr/Ti ratio using regression analysis within definite composition ranges.

<table><thead><tr><th></th><th colspan="10">Mole fraction of PbTiO<sub>3</sub> in PZT<sup>a</sup></th></tr><tr><th></th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th><th>1.0</th></tr></thead><tbody><tr><th>P<sub>s</sub> (C/m<sup>2</sup>)</th><td>0.57</td><td>0.66</td><td>0.65</td><td>0.50</td><td>0.50</td><td>0.57</td><td>0.64</td><td>0.70</td><td>0.74</td><td>0.75</td></tr><tr><th>T<sub>c</sub> (°C)</th><td>256.5</td><td>300.6</td><td>334.4</td><td>364.4</td><td>392.6</td><td>418.4</td><td>440.2</td><td>459.1</td><td>477.1</td><td>492.1</td></tr><tr><th>C (105 °C)</th><td>2.050</td><td>2.083</td><td>2.153</td><td>2.424</td><td>4.247</td><td>2.664</td><td>1.881</td><td>1.642</td><td>1.547</td><td>1.500</td></tr><tr><th>α<sub>11</sub> (10<sup>7</sup> m<sup>5</sup>/C<sup>2</sup> F)</th><td>41.25</td><td>31.29</td><td>22.30</td><td>13.62</td><td>4.764</td><td>3.614</td><td>0.6458</td><td>-3.050</td><td>-5.845</td><td>-7.253</td></tr><tr><th>α<sub>12</sub> (10<sup>8</sup> m<sup>5</sup>/C<sup>2</sup> F)</th><td>-4.222</td><td>-0.0345</td><td>1.688</td><td>2.391</td><td>1.735</td><td>3.233</td><td>5.109</td><td>6.320</td><td>7.063</td><td>7.500</td></tr><tr><th>α<sub>111</sub> (10<sup>8</sup> m<sup>9</sup>/C<sup>4</sup> F)</th><td>5.068</td><td>4.288</td><td>3.560</td><td>2.713</td><td>1.336</td><td>1.859</td><td>2.348</td><td>2.475</td><td>2.518</td><td>2.606</td></tr><tr><th>α<sub>112</sub> (10<sup>8</sup> m<sup>9</sup>/C<sup>4</sup> F)</th><td>34.45</td><td>18.14</td><td>15.27</td><td>12.13</td><td>6.128</td><td>8.503</td><td>10.25</td><td>9.684</td><td>8.099</td><td>6.100</td></tr><tr><th>α<sub>123</sub> (10<sup>8</sup> m<sup>9</sup>/C<sup>4</sup> F)</th><td>-8.797</td><td>-7.545</td><td>-7.052</td><td>-5.690</td><td>-2.894</td><td>-4.063</td><td>-5.003</td><td>-4.901</td><td>-4.359</td><td>-3.660</td></tr></tbody></table>

ªAntiferroelectric composition is not included.

$$ \alpha_1 = (T - T_c)/(2 \times C \times \epsilon_0). $$

The elastic and electrostrictive coefficients used in the calculations are obtained from Refs. 21 and 60.

REFERENCES

¹G. E. Haertling, “Ferroelectric ceramics history and technology,” *J. Am. Ceram. Soc.* **82**, 797–818 (1999).

²K. Uchino, *Ferroelectric Devices* (Marcel Dekker, Inc., New York, 2000).

³C.-B. Eom and S. Trolier-McKinstry, “Thin-film piezoelectric MEMS,” *MRS Bull.* **37**, 1007–1017 (2012).

⁴C. R. Bowen, H. A. Kim, P. M. Weaver, and S. Dunn, “Piezoelectric and ferroelectric materials and structures for energy harvesting applications,” *Energy Environ. Sci.* **7**, 25–44 (2014).

⁵S. I. Shkuratov and C. S. Lynch, “A review of ferroelectric materials for high power devices,” *J. Materiomics* **8**, 739–752 (2022).

⁶Y. Liu, J. F. Scott, and B. Dkhil, “Direct and indirect measurements on electrocaloric effect: Recent developments and perspectives,” *App. Phys. Rev.* **3**, 031102 (2016).

⁷G. Huangfu, K. Zeng, B. Wang, J. Wang, Z. Fu, F. Xu, S. Zhang, H. Luo, D. Viehland, and Y. Guo, “Giant electric field-induced strain in lead-free piezoelectrics,” *Science* **378**, 1125–1130 (2022).

⁸H. Jin, X. Gao, K. Ren, J. Liu, L. Qiao, M. Liu, W. Chen, Y. He, S. Dong, Z. Xu, and F. Li, “Review on piezoelectric actuators based on high-performance piezoelectric materials,” *IEEE Trans. Ultrason., Ferroelectr., Freq. Control* **69**, 3057–3069 (2022).

⁹T. E. Hooper, J. I. Roscow, A. Mathieson, H. Khanbareh, A. J. Goetzee-Barral, and A. J. Bell, “High voltage coefficient piezoelectric materials and their applications,” *J. Eur. Ceram. Soc.* **41**, 6115–6129 (2021).

¹⁰J. Shi and A. H. Akbarzadeh, “3D hierarchical lattice ferroelectric metamaterials,” *Int. J. Eng. Sci.* **149**, 103247 (2020).

¹¹H. Cui, D. Yao, R. Hensleigh, H. Lu, A. Calderon, Z. Xu, S. Davaria, Z. Wang, P. Mercier, P. Tarazaga, and X. Zheng, “Design and printing of proprioceptive three-dimensional architected robotic metamaterials,” *Science* **376**, 1287–1293 (2022).

¹²S. Zhang, B. Malić, J.-F. Li, and J. Rödel, “Lead-free ferroelectric materials: Prospective applications,” *J. Mater. Res.* **36**, 985–995 (2021).

¹³D. Damjanovic and G. A. Rossetti, Jr., “Strain generation and energy-conversion mechanisms in lead-based and lead-free piezoceramics,” *MRS Bull.* **43**, 588–594 (2018).

¹⁴M. E. Lines and A. M. Glass, *Principles and Applications of Ferroelectrics and Related Materials, Classic Text in the Physical Sciences* (Oxford University Press, 2021).

¹⁵D. Damjanovic, “Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films and ceramics,” *Rep. Prog. Phys.* **61**, 1267–1324 (1998).

¹⁶D. Damjanovic, “Contributions to the piezoelectric effect in ferroelectric single crystals and ceramics,” *J. Am. Ceram. Soc.* **88**, 2663–2676 (2005).

¹⁷A. F. Devonshire, “XCVI. Theory of barium titanate,” *London, Edinburgh Dublin Philos. Mag. J. Sci.* **40**, 1040–1063 (1949).

¹⁸A. F. Devonshire, “CIX. Theory of barium titanate—Part II,” *London, Edinburgh Dublin Philos. Mag. J. Sci.* **40**, 1065–1079 (1951).

¹⁹M. J. Haun, E. Furman, S. J. Jang, and L. E. Cross, “Thermodynamic theory of the lead zirconate-titanate solid solution system, part I: Phenomenology,” *Ferroelectrics* **99**, 13–25 (1989).

²⁰M. J. Haun, E. Furman, H. A. McKinstry, and L. E. Cross, “Thermodynamic theory of the lead zirconate-titanate solid solution system, part II: Tricritical behavior,” *Ferroelectrics* **99**, 27–44 (1989).

²¹M. J. Haun, Z. Q. Zhuang, E. Furman, S. J. Jang, and L. E. Cross, “Thermodynamic theory of the lead zirconate-titanate solid solution system, part III: Curie constant and sixth order polarization interaction dielectric stiffness coefficients,” *Ferroelectrics* **99**, 45–54 (1989).

²²M. J. Haun, E. Furman, T. R. Halemane, and L. E. Cross, “Thermodynamic theory of the lead zirconate-titanate solid solution system. Part IV: Tilting of the oxygen octahedra,” *Ferroelectrics* **99**, 55–62 (1989).

²³M. J. Haun, E. Furman, S. J. Jang, and L. E. Cross, “Thermodynamic theory of the lead zirconate-titanate solid solution system, part V: Theoretical calculations,” *Ferroelectrics* **99**, 63–86 (1989).

²⁴M. J. Haun, E. Furman, S. J. Jang, H. A. McKinstry, and L. E. Cross, “Thermodynamics theory of PbTiO₃,” *J. Appl. Phys.* **62**, 3331–3338 (1987).

²⁵F. Xue, Y. Ji, and L.-Q. Chen, “Theory of strain phase separation and strain spinoidal: Applications to ferroelectric and ferroelastic systems,” *Acta Mater.* **133**, 147–159 (2017).
<footnote n="26">²⁶A. C. F. Cocks and R. M. McMeeking, "A phenomenological constitutive law for the behavior of ferroelectric ceramics," *Ferroelectrics* **228**, 219–228 (1999).</footnote>

<footnote n="27">²⁷J. E. Huber, N. A. Fleck, and R. M. McMeeking, "A crystal plasticity model for ferroelectrics," *Ferroelectrics* **228**, 39–52 (1999).</footnote>

<footnote n="28">²⁸J. E. Huber, N. A. Fleck, C. M. Landis, and R. M. McMeeking, "A constitutive model for ferroelectric polycrystals," *J. Mech. Phys. Solids* **47**, 1663–1697 (1999).</footnote>

<footnote n="29">²⁹C. M. Landis, "Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics," *J. Mech. Phys. Solids* **50**, 127–152 (2002).</footnote>

<footnote n="30">³⁰C. M. Landis, "Non-linear constitutive modelling of ferroelectrics," *Curr. Opin. Solid State Mater. Sci.* **8**, 59–69 (2004).</footnote>

<footnote n="31">³¹L.-Q. Chen, "Phase-field models for microstructure evolution," *Annu. Rev. Mater. Sci.* **32**, 113–140 (2002).</footnote>

<footnote n="32">³²M. R. Tonks and L. K. Aagesen, "The phase field method: Mesoscale simulation aiding material discovery," *Annu. Rev. Mater. Sci.* **49**, 79–102 (2019).</footnote>

<footnote n="33">³³L. Guin and D. M. Kochmann, "A phase-field model for ferroelectrics with general kinetics, part 1 model formulation," *J. Mech. Phys. Solids* **176**, 105301 (2023).</footnote>

<footnote n="34">³⁴L. Daniel, D. A. Hall, and P. J. Withers, "A multiscale model for reversible ferroelectric behavior of polycrystalline ceramics," *Mech. Mater.* **71**, 85–100 (2014).</footnote>

<footnote n="35">³⁵N. Buiron, L. Hirsinger, and R. Billardon, "A multiscale model for magnetoelastic couplings," *J. Phys. IV* **9**, 187–196 (1999).</footnote>

<footnote n="36">³⁶L. Daniel, O. Hubert, N. Buiron, and R. Billardon, "Reversible magnetoelastic behavior: A multiscale approach," *J. Mech. Phys. Solids* **56**, 1018–1042 (2008).</footnote>

<footnote n="37">³⁷S. Klinkel, "A phenomenological constitutive model for ferroelastic and ferroelectric hysteresis effects in ferroelectric ceramics," *Int. J. Sol. Struct.* **43**, 7197–7222 (2006).</footnote>

<footnote n="38">³⁸C. Miehe, D. Rosato, and B. Kiefer, "Variational principles in dissipative electro-magneto-mechanics: A framework for the macro-modeling of functional materials," *Int. J. Numer. Methods Eng.* **86**, 1225–1276 (2011).</footnote>

<footnote n="39">³⁹M. Meindlhumer, A. Pechstein, and A. Humer, "Variational inequalities for ferroelectric constitutive modelling," *J. Intell. Mater. Syst. Struct.* **32**, 317–330 (2021).</footnote>

<footnote n="40">⁴⁰C. Babori, M. Barati, and L. Daniel, "An energy-based model for ferroelectric ceramics," *Eur. J. Mech: A/Solids* **103**, 105151 (2024).</footnote>

<footnote n="41">⁴¹D. Vanderbilt and M. H. Cohen, "Monoclinic and triclinic phases in higher-order Devonshire theory," *Phys. Rev. B* **63**, 094108 (2001).</footnote>

<footnote n="42">⁴²Y. L. Wang, A. K. Tagantsev, D. Damjanovic, N. Setter, V. K. Yarmarkin, A. I. Sokolov, and I. A. Lukyanchuk, *J. Appl. Phys.* **101**, 104115 (2007).</footnote>

<footnote n="43">⁴³K. Uchino and S. Takahashi, "Multilayer ceramic actuators," *Curr. Opin. Solid State Mater. Sci.* **1**, 698–705 (1996).</footnote>

<footnote n="44">⁴⁴P. M. Chaplya, M. Mitrovic, G. P. Carman, and F. K. Straub, "Durability properties of piezoelectric stack actuators under combined electromechanical loading," *J. Appl. Phys.* **100**, 124111 (2006).</footnote>

<footnote n="45">⁴⁵J. Koruza, D. J. Franzbach, F. Schader, V. Rojas, and K. G. Webber, "Enhancing the operational range of piezoelectric actuators by uniaxial compressive preloading," *J. Phys. D: App. Phys.* **48**, 215302 (2015).</footnote>

<footnote n="46">⁴⁶D. A. Hall, A. steuwer, B. Cherdhirunkorn, T. Mori, and P. J. Withers, "A high-energy synchrotron x-ray study of crystallographic texture and lattice strain in soft lead zirconate titanate ceramics," *J. Appl. Phys.* **96**, 4245 (2004).</footnote>

<footnote n="47">⁴⁷M. Hinterstein, J. Rouquette, J. Haines, P. Papet, M. Knapp, J. Glaum, and H. Fuess, "Structural description of the morphotropic piezo- and ferroelectric properties of lead zirconate titanate," *Phys. Rev. Lett.* **107**, 077602 (2011).</footnote>

<footnote n="48">⁴⁸S. C. Hwang, C. S. Lynch, and R. M. McMeeking, "Ferroelectric/ferroeleastic interactions and a polarization switching model," *Acta Metall. Mater.* **43**, 2073 (1995).</footnote>

<footnote n="49">⁴⁹J. L. Jones, M. Hoffman, and K. J. Bowman, "Saturated domain switching textures and strains in ferroelastic ceramics," *J. Appl. Phys.* **98**, 024115 (2005).</footnote>

<footnote n="50">⁵⁰J. Rödel and W. S. Kreher, "Self-consistent modelling of non-linear effective properties of polycrystalline ferroelectric ceramics," *Comp. Mater. Sci.* **19**, 123 (2000).</footnote>

<footnote n="51">⁵¹J. L. Jones, E. B. Slamovich, and K. J. Bowman, "Domain texture distributions in tetragonal lead zirconate titanate by x-ray and neutron diffraction," *J. Appl. Phys.* **97**, 034113 (2005).</footnote>

<footnote n="52">⁵²A. Pramanick, D. Damjanovic, J. E. Daniels, J. C. Nino, and J. L. Jones, "Origins of electro-mechanical coupling in polycrystalline ferroelectrics during subcoercive electrical loading," *J. Am. Ceram. Soc.* **94**, 293 (2011).</footnote>

<footnote n="53">⁵³D. A. Hall, L. Daniel, M. Watson, A. Condie, T. P. Comyn, A. M. Kleppe, and P. J. Withers, "Domain switching and shear-mode piezoelectric response induced by cross-poling in polycrystalline ferroelectrics," *J. Appl. Phys.* **136**, 194101 (2024).</footnote>

<footnote n="54">⁵⁴H. Liu, J. Chen, H. Huang, L. Fan, Y. Ren, Z. Pan, J. Deng, L.-Q. Chen, and X. Xing, "Role of reversible phase transformation for strong piezoelectric performance at the morphotropic phase boundary," *Phys. Rev. Lett.* **120**, 055501 (2018).</footnote>

<footnote n="55">⁵⁵M. Liu, K. J. Hsia, and M. R. Sardela, Jr., "In situ diffraction study of electric-field-induced domain switching and phase transition in PZT-5H," *J. Am. Ceram. Soc.* **88**, 210–215 (2005).</footnote>

<footnote n="56">⁵⁶B. Kaeswurm, V. Segouin, L. Daniel, and K. G. Webber, "The anhysteretic polarization of ferroelectrics," *J. Phys. D: Appl. Phys.* **51**, 075305 (2018).</footnote>

<footnote n="57">⁵⁷C. babori, M. Barati, V. Segouin, R. Corcolle, and L. Daniel, "Anhysteretic strain in ferroelectric ceramics under electromechanical loading," *J. Phys. D: Appl. Phys.* **57**, 445001 (2024).</footnote>

<footnote n="58">⁵⁸F. X. Li and R. K. N. D. Rajapakse, "A constrained domain-switching model for polycrystalline ferroelectric ceramics. Part I: Model formulation and application to tetragonal materials," *Acta Mater.* **55**, 6472–6480 (2007).</footnote>

<footnote n="59">⁵⁹F. X. Li and R. Rajapakse, "A constrained domain-switching model for polycrystalline ferroelectric ceramics. Part II: Combined switching and application to rhombohedral materials," *Acta Mater.* **55**, 6481–6488 (2007).</footnote>

<footnote n="60">⁶⁰A. Amin, R. E. Newnham, and L. E. Cross, "Effect of elastic boundary conditions on morphotropic Pb(Zr,Ti)O₃," *Phys. Rev. B* **34**, 1595 (1986).</footnote>