Electrothermal Phenomena in Ferroelectrics

: Physical properties of lead-zirconate-titanate (PZT) ceramics change according to the initial electric poling process and electrical boundary conditions. This paper reports the electrothermal, piezothermal, and piezoelectric coupling phenomena in ferroelectrics from thermodynamics viewpoints, in particular, thermal property differences between unpoled and poled PZT’s in the poling direction for open circuit and short circuit conditions. We propose a new terminology, “secondary electrothermal” coupling factor k λ , which is analogous to the electromechanical coupling factor k , relating the elastic compliances under short- and open-circuit conditions, in order to explain the fact that the short-circuit condition exhibited the larger thermal diffusivity than the open-circuit condition. On the other hand, the unpoled specimen exhibits the lowest thermal diffusivity. This tutorial paper was authored for providing comprehensive knowledge on equilibrium and time-dependent thermodynamics in ferroelectrics.


Background
Thermoelectric devices generate voltage when each side of the device is set at different temperatures. Conversely, when voltage is applied to it, heat is transferred from one side to the other, creating a temperature gradient. This effect can be used to generate voltage in "thermocouple" to measure temperature of a specimen. Because heating or cooling is determined by the applied electric field direction, thermoelectric devices can also be utilized as temperature controllers. Popular thermoelectric phenomena include the Seebeck effect, Peltier effect, and Thomson effect in metallic or semiconductive materials based on the "electron" flow. On the contrary, this paper handles the "electrothermal" effect in resistive ferroelectric materials based on the "phonon" flow. The "pyroelectric effect" is a phenomenon in which a ferroelectric material shows voltage under a change in temperature. Its inverse effect, the so-called "electrocaloric effect", shows a reversible temperature (increase or decrease) change under an applied electric field polarity. It should not be confused in principle with the semiconductor "thermoelectric effect", in which a temperature difference occurs when a current is driven through an electric junction with two dissimilar conductors.
Electrocaloric materials were initially the focus of significant application interest in the 1940s (during World War II). The US and Japan accelerated the research for developing the air-cooling system particularly in submarines without generating mechanical noise (like air-compressorembedded refrigerators). Due to the confidential military research, not many publications have been disclosed on this research. The author personally learned from his former bosses in the university in the 1960s, when the research on this topic was fading out because the electrocaloric effect was insufficient for practical applications. However, in the beginning of 21st century, we restarted its development in parallel to the enthusiasm on new compact refrigeration technique in the "sustainable" society. Our report on giant response in Pb(Zn1/3Nb2/3)O3-PbTiO3 (PZN-PT) bulk single crystals in 2003 [1], being 0.3 °C under 1 kV/mm electric field applied, ignited the "renaissance" of electrocaloric effect. Successive reports on the high response, 12 °C under an applied electric field of 48 kV/mm in thin film (or ribbon) lead zirconate-titanate (PZT) by Mischenko et al. [2] (that is equivalent to 0.25 °C under 1 kV/mm, lower performance than PZN-PT single crystals), and the demonstration of 12 °C of cooling near room temperature with a ferroelectric polymer by Neese et al. in 2008 [3] also accelerated the research boom on electrocaloric devices.
On the other hand, our research group at the Penn State University has been concentrating on the development of high power piezoelectrics for 40 years, which can handle the power density as high as 40 W/cm 3 [4]. One of the Figures-of-Merit for the high-power density is the "maximum vibration velocity" in the same material species such as PZT's, which is defined by the vibration velocity generated at the sample edge under its resonance operation, when the maximum temperature rise (i.e., at the nodal point) reaches 20 °C above the room temperature [5,6]. When we operate the piezo-sample under higher vibration levels, even though highly resistive under the offresonance, the piezoelectric PZT starts generating heat via its dielectric, elastic and piezoelectric losses. Most of the additional input electric energy over the maximum vibration velocity level is converted to heat, like a "ceramic heater"! The current commercially available top data is around 0.6 m/s in a rectangular plate PZT's, which corresponds to the energy density of 40 W/cm 3 .
We found that the maximum vibration velocity strongly depends on the thermal conductivity. Figure  1 demonstrates the simulation results of the saturated temperature distribution profile difference between two compositions, PZT-5H and PZT-19, under the same vibration velocity operating condition. These two compositions exhibit significant difference in thermal conductivity: 0.14 W/m K versus 2 W/m K [7,8]. The profile curve of PZT-5H (low thermal conductivity) fits a sinusoidal line beautifully, while PZT-19 shows considerable edge temperature rise. You can understand easily that this profile difference is originated primarily from the thermal conductivity or diffusivity difference. Taking the total thermal energy dissipated from the specimen by integrating the temperature rise with respect to the position coordinate, we can expect a similar "mechanical quality factor" Qm value (i.e., similar intensive elastic loss tan ϕ') for these two specimens. However, you can notice that the peak temperature at the nodal point is significantly lower for the larger thermal diffusivity material, which means we can excite the vibration more under higher voltage for the PZT-19, since the maximum vibration velocity is defined by the highest temperature rise 20 °C above room temperature at the node (plate center) from the human safety viewpoint in the practical industrial standard (e.g., human fingers will be burned above 50 °C). This is our motivation on the thermal property measurements, in particular, conductivity/diffusivity on ferroelectric/piezoelectric ceramics from the practical application viewpoint. We propose a "secondary electrothermal" coupling factor k λ , which is applied for irreversible thermal flow status and explain the difference in the thermal conductivity between the electric constraint difference, that is, under short-and open-circuit conditions [9]. This tutorial paper was written for the encouragement to the researchers, who are currently studying electrocaloric or high-power piezoelectric materials, by providing them comprehensive and underlying thermodynamical backgrounds on the "electrothermal effect" in comparison with the "mechanothermal" and "electromechanical coupling". Suitable textbooks or review papers have not been found in this area, to the best of the author's knowledge.

Review of Phenomenology in Ferroelectrics
We discuss the phenomenology in ferroelectrics on the following four stages with increasing the complexity in analyses: (1) linear equilibrium in Section 2.1, (2) nonlinear equilibrium in Sections 2.1.3 and 4.1, (3) harmonic dynamic in Sections 2.2 and 4.2, and (4) anharmonic dynamic in Section 4.3.

Thermodynamical Functions
A thermodynamic phenomenological theory is discussed basically in the form of expansion series of the "internal energy" U as a function of the physical properties; the free energy is composed of "thermal energy" dQ, "mechanical energy" Xdx (stress X, strain x) and "electrical energy" EdD (dielectric displacement D (almost equal to polarization P in high permittivity materials), electric field E) in ferroelectrics. When a material is ferromagnetic, magnetic energy HdB is integrated. Further, thermal energy dQ is given by dQ ≦ TdS, where temperature is T and entropy is S. Equality "=" is for the reversible process, and "<" is true for the irreversible process. We can describe U as follows: According to the combination possibilities among (T, S), (X, x) and (E, D), there are eight types of energy. We will use Helmholtz free energy A = U − TS, Gibbs free energy G = U − TS − Xx − ED and elastic Gibbs energy G1 = U − TS − Xx as examples in this paper; that is, Internal energy U is described in terms of "extensive" physical parameters (S, x, D). Since the phase transition and experiments are normally conducted under external T = constant, X = 0 and E = 0, Gibbs free energy described in terms of "intensive" physical parameters (T, X, E) is the most popularly utilized. The elastic Gibbs energy G1 is most convenient for discussing the ferroelectric phase transition because of the two merits: (1) by describing higher order Taylor expansion terms of the "order parameter" D (or P), we can derive the spontaneous polarization, and (2) external X can be controlled explicitly, and E is related with easily with , .
We can obtain the famous "Maxwell relationship" examples from Equation (2) as follows, taking the second partial derivative: The meaning of Equation (3a,b) are the "piezoelectric coefficients", discussed in Section 2.1.4.1, while Equation (3c) is used in the "pyroelectric coefficient" p derivation in Section 2.1.6.

Linear Energy Handling
We consider a practical formula of the Gibbs free energy G(T, X, E) for the case of a small value change in temperature θ = T − TR (room temperature), external X and E (1D case). If the change of parameters is small, we may adopt the three-parameter Taylor expansion approximation up to second derivatives in order to discuss just the linear relationships, based on the description by Mitsui et al. [10]: Taking into account =--− , we obtain first the relations, = − . Since is the entropy density at , , = 0, we take this as the original value and set it to = 0. The values and ≈ are considered to be spontaneous strain and spontaneous polarization in the ferroelectric phase of this material (refer to the next Section 2.1.3), and set them as zero in the discussion merely in the ferroelectric phase. Now, Equation (4) can be transformed as Continued discussions on the linear handling are made after Section 2.1.4. Though nonlinear energy handling is not our direct purpose in this section, we describe Landau/Devonshire theory briefly in order to provide the reader the key idea on how the spontaneous polarization and strain are generated in a ferroelectric material [11][12][13][14]. We use the elastic Gibbs energy G1 for the phase transition discussion, as in Section 2.1.1.

Nonlinear Energy Handling
We first construct the electric energy formula, expanding the energy in terms of so-called order parameter, polarization P, including nonlinear terms as follows: Considering that the free energy form of the crystal should be applied for any temperature range below or above the phase transition temperature (which is called "Curie temperature"), and that if the mother phase (high temperature phase) has "cubic symmetry" (like perovskite PZT and PZN-PT), the Gibbs energy should not change with polarization reversal (P → −P), that is, G1(P) = G1(−P).
Thus, the energy expansion series should not contain terms in odd powers of a component, P, or the expansion series include only even powers of P.

Second-Order Phase Transition
We start from the simplest nonlinear phenomenon "second-order phase transition" model [14]: The coefficients α and β depend, in general, on the temperature, but only α is assumed to be temperature dependent and β is positive constant in the following discussion for the sake of minimizing the required parameters. The phenomenological formulation should be applied to the whole temperature range over which the material is in the paraelectric (high temperature) and ferroelectric (low temperature) states. Because the spontaneous polarization should be zero in the paraelectric state, the free energy should be zero in the paraelectric phase at any temperatures above its phase transition temperature (i.e., Curie temperature). To stabilize the ferroelectric state with a spontaneous polarization, the free energy for a certain polarization P should be lower than "zero".
Otherwise, the paraelectric state is maintained. Thus, at least, the coefficient α of the P 2 term must be negative for the spontaneously polarized state to be stable, while, in the paraelectric state, it must be positive passing through zero at a certain temperature T0, called Curie-Weiss temperature. From this concept, we assume α as a linear relation: where C is taken as a positive constant called the Curie-Weiss constant, and T0 is equal to or lower than the actual transition temperature TC (Curie temperature), in general. This temperature dependence can also be interpreted as the adoption of the "TP 2 " electrothermal coupling term [ ] in the Taylor expansion of the G1 energy. The temperature parameter θ here is the deviation from the Curie-Weiss temperature T0. From =-− + [ ≈ in our discussion], the equilibrium polarization under an electric field E should satisfy the condition: With no electric field applied, Equation (8) provides two cases: (i) P = 0 and (ii) P 2 = −α/β. The trivial solution in (i) case corresponds to a paraelectric (high temperature) state, while the polarization solutions (PS = ± − ⁄ in (ii) case correspond to a ferroelectric (low temperature, α < 0) state. Figure 2 shows the temperature dependence of the free energy curve. One potential minimum at a high temperature will split into two minimum-branches at a low temperature phase. This "Y"shape phase splitting is called "bifurcation" in terms of temperature. The critical point corresponds to the Curie temperature TC (in this case also T0). In a lattice vibration model, we can rephrase the above situation as follows: at a high temperature, the nearly "harmonic" atomic (or lattice) vibration (usually an "optical phonon" mode, which generates the polarization deviation) in a sharp convex "parabolic" potential at a high temperature will reduce the vibration frequency with decreasing the temperature approaching to the Curie temperature (i.e., critical slowing down) because of the flat curvature of the potential around TC. This "phonon softening" phenomenon is interpreted by the nonlinearity of the phonon mode [Refer to Section 4.2]. Because the free energy curve is symmetric with respect to the polarization P when the external field E = 0, the probability of the state +PS (=+ − ⁄ or −PS (= − − ⁄ should be equal. Thus, with decreasing the temperature, passing through TC, +PS domains and −PS domains may arise locally in a specimen with equal volumetric ratio, leading to the "multidomain states" in general. In order to produce the net spontaneous polarization, strain, and piezoelectricity, we need to apply large electric field (higher than "coercive field", which is given by = ⁄ in the second order phase transition), which is called the "electric poling" process to make a crystal with mono-domain-like status. The double-well potential change with the applied field is illustrated in the inserted figure. Note first that the minimum potential well energy is negative (that is, stabler than the paraelectric phase energy zero), and second that external electric field enforces one spontaneous state more stable. Though the poling process in polycrystalline ceramic samples, composed of many micro single crystals (i.e., "grain") with random crystal orientations, is not easy, a similar phenomenology is assumed to be adopted qualitatively in this paper.

First-Order Phase Transition
The theories for electrostriction in ferroelectrics were formulated in the 1950s by Devonshire [11,12] and Kay [13]. The elastic Gibbs energy is expanded in a 1D form: where P, X, T are the polarization, stress, and temperature, respectively, and s and Q are called the elastic compliance and the electrostrictive coefficient, respectively. The elastic Gibbs energy definition, =-− + , leads to When the external stress X is zero, the following equations are derived: When β < 0, γ > 0, the "first-order" phase transition comes out. If the external electric field is equal to zero (E = 0), two different states are derived; P = 0 and P 2 = (−β + − 4 )/2γ. Electrostriction: (ii) Ferroelectric phase: PS 2 = (−β + − 4 )/2 or P = PS + ε0ε E (under small E) where we define the spontaneous strain xS and the piezoelectric constant d as: Spontaneous strain: xS = QPS 2 (16) Piezoelectric constant: d = 2 ε0 ε QPS (17) Note here that the electrostriction is a nonlinear phenomenon originated from the anharmonic potential energy, which will be again discussed in Section 3.1. We see from Equation (17) that piezoelectricity is equivalent to the "electrostrictive phenomenon" biased by the spontaneous polarization", which was already mentioned in the derivation process of Equation (4). The temperature dependences of the spontaneous polarization and dielectric constant are plotted in Figure 3 by comparing the second-order (a) and first-order (b) phase transitions. Temperature dependence of various parameters are expressed in simple formulas in the case of second-order transition as follows: (i) Paraelectric phase -T > T0 From Equations (18b) and (19b), the piezoelectric constant is obtained as We continue the "linear energy handling" in Section 2.1.2. When the temperature is constant (i.e., isothermal), = 0 in Equation (5b,c), we can obtain the "intensive" parameter-based piezoelectric constitutive equations: where the following notations are used, and denoted as elastic compliance under constant E, dielectric permittivity under stress free: The Maxwell relation, Equation (3a), , = , verifies that the piezoelectric coefficient d in Equation (22a,b) are thermodynamically the same. When we start from the Helmholtz free energy A ( =-+ + ), by taking a similar Taylor expansion approach, we obtain another set of piezoelectric constitutive equations in terms of "extensive" parameters, x and D: where c D is elastic stiffness under constant D, and κ 0 κ x is inverse permittivity ( = 1 ⁄ ) under strain free condition, and these coefficients are expressed by: The Maxwell relation, , = , in Equation (3b), verifies that the inverse piezoelectric coefficients h in Equation (24a,b) are thermodynamically the same.

Electromechanical Coupling Factor
The term, electromechanical coupling factor k, is defined as the square value k 2 being the ratio of the converted energy over the input energy: when electric to mechanical k 2 = (Stored mechanical energy/Input electrical energy), and when mechanical to electric, k 2 = (Stored electrical energy/Input mechanical energy) (27a) Let us derive Equation (26a) first practically, when an external electric field E3 is applied to a piezoelectric material in a pseudo-static process. See Figure 4a first, when we apply the electric field on the top and bottom electrodes under stress free conditions (X = 0). Input electric energy must be equal to (1/2)ε0ε3 X E3 2 from Equation (4), and the output strain generated by E3 should be d33E3 from Equation (22a). Since the converted/stored mechanical energy is obtained as (1/2 s33 E ) x3 2 , we obtain: Let us now consider Equation (27a), when an external stress X3 is applied to a piezoelectric material in a pseudo-static process. Refer to Figure 4b. Under short-circuit condition (E3 = 0), the input mechanical energy must be equal to (1/2)s33 E X3 2 from Equation (4), and the electric displacement D3 (or polarization P3) generated by X3 should be equal to d33X3 from Equation (22b). This D3 can be obtained by integrating the short-circuit current in terms of time through the electric lead. Since the converted/stored electric energy is obtained as (1/2 ε0 ε3 X ) D3 2 , we obtain: It is essential to understand that the electromechanical coupling factor k (or k 2 , which has a physical meaning of energy transduction/conversion rate), can be exactly the same for both converse (26b) and direct (27b) piezoelectric effects. The conditions under constant X (free stress) or constant E (short-circuit) are considered to be non-constrained.

Constraint Physical Parameters
It is important to consider the conditions under which a material will be operated when characterizing the dielectric constant and elastic compliance of that material. When a constant electric field is applied to a piezoelectric sample as illustrated in Figure 5 Top, the total input electric energy (left) should be equal to a combination of the energies associated with two distinct mechanical conditions that may be applied to the material: (1) stored electric energy under the mechanically clamped state, where a constant strain (zero strain) is maintained and the specimen cannot deform, and (2) converted mechanical energy under the mechanically free state, in which the material is not constrained and is free to deform. This situation can be expressed by: such that: When a constant stress is applied to the piezoelectric as illustrated in Figure 5 Bottom, the total input mechanical energy will be a combination of the energies associated with two distinct electrical conditions that may be applied to the material: (1) stored mechanical energy under the open-circuit state, where a constant electric displacement is maintained, and (2) converted electric energy (i.e., "depolarization" field) under the short-circuit condition, in which the material is subject to a constant electric field. This can be expressed as: which leads to: In principle, if we measure the permittivity in a piezoelectric specimen under stress-free and completely-clamped conditions, we can obtain and , respectively. However, in practice, cannot be measured because of the experimental difficulty to maintain the ideal strain-free (clamped) condition for a long period. Similarly, if we measure the strain in a piezoelectric specimen as a function of applied stress pseudo-statically, under short-circuit and open-circuit conditions, we can obtain and , respectively. However, in practice, cannot be measured because of the induced bound charge (or polarization) "screening" by migrating charge in the electrode, specimen or surrounding atmosphere in a couple of minutes. Thus, the clamped permittivity or open-circuit D-constant can only be measured with high-frequency dynamical methods, such as an impedance analyzer around the resonance/antiresonance frequencies.
In conclusion, we obtain the following equations: We can also write equations of a similar form for the corresponding reciprocal quantities: where, in this context, This new parameter k in Equation (30b) is also the electromechanical coupling factor in the "extensive" parameter description, and identical to the k in Equation (30a). Note the k expression derivation from the piezoelectric constitutive equations, Equation (22a,b) and Equation (24a,b):

Adiabatic Process 1-Piezothermal Effect
When we discuss the piezoelectric coupling phenomena, we assume the isothermic condition; that is, the specimen temperature is maintained constant, even the energy conversion is conducted between the electrical and mechanical energy. We consider next the adiabatic process; that is, the specimen is isolated from the external heat source, and temperature may be changed during the energy conversion process by the external input. When the electric field is constant, = 0 in Equation (5a,b), we can obtain the following equations: where the following notations are used, and denoted as specific heat capacity under X = 0 and E = 0, elastic compliance under constant E: Let us discuss here the diagonal expansion coefficient − in terms of specific heat capacity.
Recall the relation dq = TdS in the "reversible" thermal process, where dq is the thermal energy flow per unit "mass", given by the total energy flow dQ = ρVdq (ρ: mass density, V: volume). The specific heat capacity is defined by It is noteworthy that the piezothermal coefficient is usually called "linear thermal expansion This piezothermal coefficient contributes to the converse effect; that is, temperature change under the stress application, which can be called "piezo-caloric" effect. The reader might have experienced in your elementary or middle school age that you felt "cool" temperature when you touched a thick rubber band on your lip immediately after expanding it. The piezothermal coefficient is originated from a nonlinear elastic vibration or the anharmonic phonon interaction, which will be discussed in Section 4.3.1. The "piezothermal coupling factor" k PT can be defined from Equation (33a,b) by We can derive the completely "clamped" (strain x free) specific heat capacity and the "adiabatic" elastic compliance in terms of theoretically as follows, though these may not be useful physical parameters: Under strain-free condition (x = 0), we obtain = from Equation (33b). Since the other hand, under the adiabatic condition (S = 0), we obtain = from Equation (33a). Since

Adiabatic Process 2-Electrothermal Effect
Direct effect is known as pyroelectricity, which is popularly used for temperature sensors and night-vision devices, while the converse effect is "electrocaloric" phenomenon, which is a recent hot topic for the application to a new refrigeration technology, as mentioned in Section 1 Background.
When the stress is constant, = 0 in Equation (5a,c), we can obtain the following equations: where the following notations are used, and denoted as specific heat capacity (per unit mass) under X = 0 and E = 0, permittivity under constant stress X: The primary electrothermal coupling coefficient is usually called "pyroelectric coefficient", defined by where we used the relation = (≈ ).

Constraint Specific Heat Capacity
In Equation (35), we introduced specific heat capacity (per unit mass) under X = 0 and E = 0; that is, under a short-circuit condition of a ferroelectric specimen's surface electrodes. We may consider a different specific heat capacity under an open-circuit condition (i.e., D = constant or zero) [14].
Taking the first derivative of Equation (38a) with respect to T by keeping X = D = 0, and = , Here, we denoted the primary "electrothermal coupling factor" from Equation (38a,b) as It is important to note that Equation (42) is analogous to Equation (29b) to correlate the Dconstant and E-constant parameters in terms of "coupling factors", and .

Constraint (Adiabatic) Permittivity
Permittivity was defined isothermally so far. However, we may consider adiabatic permittivity theoretically when no heat flow is hypothesized, such as the case where a ferroelectric specimen is suspended in a vacuum chamber [10].
From Equation (38b), isothermal permittivity ( = 0) is given by Thus, "adiabatic (S =constant) permittivity" is related with "isothermal (T = constant) permittivity" again by using primary electrothermal coupling factor as: 2.1.6.3. Electrocaloric Effect Equation (38a) gives the necessary formula for the electrocaloric effect. When we take "adiabatic" condition; that is, constant entropy dS = 0, and Equation (45) = Here, is specific heat capacity (per unit mass) under E = 0, and p is pyroelectric coefficient given by = − . The material's development strategy should be higher p and lower at room temperature (i.e., operation temperature 300 K).
Using a simple second-order phase transition, let us evaluate the temperature dependence of the electrocaloric figure of merit, . Recall the following equations obtained in the ferroelectric phase (T < TC = T0) in Section 2.1.3.2: (a) Specific Heat Capacity Since the pyroelectric coefficient increases significantly with approaching the phase transition temperature, in practical applications, the Curie temperature of the material is designed slightly higher than room temperature (i.e., operation temperature).
In order to realize the refrigeration function, a large "negative electric field" is required theoretically, so that the material's "depoling" is anticipated. Since the coercive field is evaluated in the second-order phase transition [4] as the EC reduces significantly with approaching the Curie temperature. Thus, in practice, by applying a high positive electric field gradually (pseudo-equilibrium or isothermally), then, a sudden decrease down to zero is conducted to lower the temperature to escape from the depoling problem in a low coercive field material. Specific heat capacity is measured by Differential Scanning Calorimetry (DSC), where endothermic and exothermic processes are monitored by keeping the input energy flow constant. Since the measuring time period is around 5-10 min, we cannot observe the difference between the short-circuited and open-circuited ferroelectric specimens, that is, and . This happens probably because the pyroelectrically induced charge is almost screened by the migrating free charge during a rather long measuring time (several minutes): • Permittivity Finally, the electrothermal coupling factor is calculated as follows: We can conclude that the electrothermal coupling factor in PZTs has a comparable order of magnitude to the electromechanical coupling factor k of 30-70%. It is interesting to calculate the value for the second-order phase transition case: using Equations (48)-(50), = 100% , irrelevant to any temperature, very idealistically.

Time-Dependent Ferroelectric Phenomenology
We handle transident reponse of a physical parameter after slightly deviated from the equilibrium state in terms of the "relaxation" time.

Polarization Relaxation
Landau and Khalatnikov developed a theory for the temperature dependence of the "relaxation time" of the order parameter, P, that is based on the Ginzburg-Landau theory of phase transitions [20,21]. This mean-field approach shows a divergence of ∝ 1 ( − ) ⁄ near a second-order phase transision point, TC. The key assumption is that the change of polarization P with time (i.e., ) is proportional to internal energy decrease rate with the unit order parameter change (i.e., ). By introducing "relaxation time" , the inverse of which is taken as a proportional constant: We use again the simple "second-order phase transition" model here for a practical calculation: Recall the fact that, in the ferroelectric phase (T < T0), the spontaneous polarization is expressed as We consider the initial polarization ( + by applying small electric field on a ferroelectric specimen with = . The meaning of Equation (55) is to provide how quickly the P should return to the equilibrium state after sudden termination of the field. Equation (55) can be solved as The actual "relaxation time constant" τ is provided by From Equations (56) and (57), we can conclude the following items: • Time constant is proportional to "permittivity" The relaxation time τ is a measure of time delay of P generation after the external E is applied, which corresponds to the "reaction resistance" or "dielectric loss". A similar concept is introduced in Section 4.3.2 for explaining the thermal resistivity/conductivity.

Temperature Relaxation
Time dependence of temperature (uniform temperature in a small volume) is considered in this section. We start from an equation similar to the previous section: If we take the Gibbs energy transformed from Equation (4) under stress-free condition (X = 0), The meaning of Equation (60) is to provide how quickly the θ should approach the equilibrium status after sudden application of the field E. As you realize, the electrocaloric effect under an adiabatic condition indicates the temperature change in the equilibrium status [Equation (45)] as = .
Equation (60) can be solved as The actual relaxation time constant τ is expressed by From Equations (62) and (63), we can conclude the following items: • In a uniform (no space gradient) specimen, temperature change follows an exponential trend with time: • t  ∞ gives  . corresponds to the temperature change by the electrocaloric effect. • Time constant is proportional to temperature, and inversely proportional to the specific heat capacity . The larger the specific heat capacity , and the lower the temperature, the lower the time constant . The recovery time of θ is quicker. This τ is roughly a suitable rise time period of an applying pseudo-step electric field.
The relaxation time constant τ in Equation (63) is very relevant to the time constant in thermal diffusivity measurement discussed in Section 3.2, where ∝ 1⁄ = . Note the opposite effect of the specific heat capacity to the above relaxation time constant, and the heat flow relaxation time.

Thermal Diffusivity and Conductivity
So far, we discussed the equilibrium and time-dependent (dynamic) phenomenology in the uniform ferroelectric specimen (that is, uniform temperature distribution in a specimen). We will discuss spatially non-uniform (i.e., "space-gradient") phenomenology in this chapter. The concept of "heat flow" is introduced, where the exponential trend with time such as (1 − ⁄ ) = ⁄ is not sustained correctly.

1D Heat Conduction Model
We will develop a 1D heat transfer model, as shown in Figure 6a, where a uniform disk or rod is thermally isolated so that no heat escapes from its sides [22]. Temperature gradient is introduced, and we define the temperature function θ(x,t) (( − ) difference from the TR) at a point a distance x from one end of rod. We introduce the following three parameters: Again expanding ( + , ) by Taylor's expansion, we obtain From Equations (64b) and (65b) to eliminate the parameter q, we obtain 1D heat transfer equation as The proportional constant is called "thermal diffusivity" , which is practically measured by thermal experiments. Units and typical numbers of the necessary parameters in PZT are summarized: • Thermal conductivity -0. 6

Solution of 1D Heat Transfer Equation
Let us solve the heat equation Equation (66) We consider the boundary condition: Large heat source ( = 0) is attached suddenly (t = 0) on the surface A (x = 0) in Figure 6a; that is, heat flux is given as a step-function with respect to time. Excluding immediately after t = 0 and adjacent to x = 0 point, we consider the region t > δ, x > δ', and assume the following equation in most of the rod material's length:  It is noteworthy that the temperature distribution profiles as a function of position x is not a simple exponential curve, but resembles a half positive part of so-called "error function", when the abrupt step-function heat flux is inserted from the position x = 0 at t = 0. Because the curve profile of Figure 7b seems to be a simple exponential curve, we previously used the approximate formula [9]: We verify this handling equivalency with the exact formula, By putting = , we obtain = , then = . Taking Taylor expansion, However, note the numerical factor difference; that is, ( 2 ⁄ ) in Equation (72a), and 4 in Equation (73). This means that the experimentally determined with Equation (72a) was overestimated from the value with Equation (72b) or (73) by the factor of 4 ( 2 ⁄ ) ⁄ . The absolute values may need to be modified by this factor from the previously reported values in Ref. [9], if required.

Thermal Diffusivity Measurements
We can find previous papers [23,24], which report on the thermal properties of PZT ceramics at low temperature (20 K-300 K), indicating a transition temperature between 50 K and 80 K, and at high temperature (300 K-800 K), characterizing the effect of phase transition on the thermal properties. However, to the author's knowledge, the relationship of thermal properties with electrical boundary conditions and poling status in the PZTs and other ferroelectric piezoelectric ceramics has not been systematically studied yet.
We utilized our new experimental setup for measuring the thermal diffusivity of PZT's. Though Ref. [25] reported the details of the experimental setup, we summarize it briefly for the reader's sake (see Figure 6b). Hot isothermal condition is imposed on a flat disk sample in our experiments by suddenly being applied on the top surface (x = 0) at t = 0. The transient temperature increase of the bottom side of the rod specimen is analytically similar to an exponential formulation with a classic time constant, as described in the previous section in Equation (73). The thermal diffusivity is given proportional to the "inverse time constant". This method hosts a variety of advantages over other methods such as (1) high accuracy, (2) low cost, (3) elimination of interface effects, and (4) small specimen size. Best of all, (5) the measuring time period is short (10-20 s) enough for measuring the difference between (short-circuit condition) and (open-circuit condition). The accuracy of our setup has already been demonstrated in several materials with low to medium thermal diffusivity (0.1-3 × 10 m 2 /s). The thermal diffusivity measurements in this experiment have an accuracy of 5% or better for some standard materials such as Fused Quartz and Pyrex 7740, in comparison to the literature values [25].
Using this method, the thermal diffusivity was measured on "poled" and "unpoled" commercially-available hard PZT ceramic disks, APC 841 (APC International, Mill Hall, PA, USA), of a diameter of 51 mm with thickness 2-5 mm. The 'unpoled' sample was "depoled" by heating the specimen, annealing, and then checking their no piezoelectric response on a conventional d33 m. The thermal diffusivity was measured in the direction of polarization for the poled samples. In parallel, the absolute value of heat capacity of a small specimen was measured using the DSC Q2000 (TA Instruments, New Castle, DE, USA), by comparing it to a sapphire reference sample. The "specific heat capacity (per mass)" and density are respectively cp = 340 J/kgK and ρ = 7600 kg/ m 3 . The specific heat capacity difference between and could not be detected under electrical boundary condition difference (short-and open-circuit) because the differential scanning calorimetry (DSC) took too long (5 min) to keep the "depolarization field". Table 1 summarizes thermal diffusivity, specific heat capacity, and thermal conductivity of a Hard PZT (APC 841), for unpoled and poled under different electrical boundary conditions (openand short-circuit) [9]. Note the significant difference in "thermal diffusivity" obtained experimentally in these three electrical constraint conditions. The origin of this difference is discussed in the following section.

Thermal Diffusivity under Different Electrical Constraints
We discovered the fact that experimentally obtained "thermal diffusivity" exhibits the largest in a poled specimen under short-circuit condition, second in a poled specimen under open-circuit condition, and the smallest in an unpoled specimen; all disk samples have electrode on both the top and bottom surfaces. Because the thermal diffusivity is expressed by where , cp, ρ are thermal conductivity, specific heat capacity, and mass density, respectively. The difference in must be the combination effect of the and cp differences (the density ρ may not be changed by the electric constraint condition).

Specific Heat Capacity-Scalar Parameter
Recall the two types of specific heat capacity, and , which are correlated by Here, the primary "electrothermal coupling factor" is given by = .
We propose an assumption on the "depolarization field" generation during temperature rise by pyroelectric effect, as shown in Figure 8: (a) poled specimen under short-circuit condition, (b) poled specimen under open-circuit condition, and (c) unpoled specimen. Because a polycrystalline PZT specimen is composed of multiple micro crystals, the situation is not actually so simple as described here. However, we assume that poled specimens are handled as a completely aligned single crystal with a monodomain; while the unpoled specimens have no net polarization, but each grain (i.e., microcrystal) generates the polarization change with temperature, which generates so-called "depolarization field" described as = − ∆ in each domain and grain [13]. Note first that the heat capacity is a "scalar" parameter, the value of which is not affected by the specimen orientation. When a specimen is poled and short-circuited with electrodes on the whole specimen (except for the narrow specimen sides) in Figure 8a, the depolarization field inside the material should be zero because the free charge in the electrode covers the bound charge on the surface and compensate the internal field; thus, is expected to be observed. On the contrary, in an unpoled specimen, though the polarization directions are randomly aligned in the specimen, since each grain experiences the polarization reduction with temperature, the internal depolarization field should be generated in each grain. We assume no free charge on the grain boundary region inside the specimen to compensate the depolarization field, though total macroscopic field may be zero (Figure 8c). Thus, is expected to be observed. We may need to modify this assumption slightly, if required, because the bound charge of the adjacent grains may partially compensate on the grain boundary; that is, may approach (slightly larger). The complex situation is found in a specimen of poled under open-circuit condition. Because of the surface electrode, E-constant twodimensionally on the disk surface, but along the heat flow direction, D-constant seems to be maintained by the depolarization field generation (Figure 8b). We assume the effective specific heat capacity be evaluated by the average of 3D such as + 2 because of scalar characteristics. We decided first = 340 ⁄ • from the DSC equipment. Then, using =, = 1 − = 279 ⁄ • . Finally, + 2 = 320 ⁄ • for an unpoled specimen was calculated. See the center column data of Table 1.
The "depolarization field" generation mechanism is very analogous to the piezoelectric effect. Recall kt or k33 piezoelectric resonance modes, where the thickness or length stress distribution generates the depolarization field = − ∆ anti-parallel to the induced polarization (∆P = d33X3) in order to satisfy the D-constant (constraint) condition. The elastic compliance and sound propagation speed ( = 1 ⁄ are modified by the electromechanical coupling factor or as = 1 − . While in the electrothermal material, the polarization modulation is made from the temperature change via the pyroelectric effect. The specific heat capacity is modified by the electrothermal coupling factor as = 1 − depending on the electrical constraint. We may also imagine the heat-flow (i.e., thermal conductivity) modulation according to the electrical constraint condition.

Thermal Conductivity-Tensor Parameter
Using = with = 7600 kg m ⁄ , we calculated thermal conductivity λ as listed in the last column of Table 1. Note first that the thermal conductivity is a second-rank tensor (similar to the electrical conductivity ), different from the scalar "specific heat capacity".
Let us discuss the difference among the measured thermal conductivity values for poled shortcircuit, poled open-circuit, and unpoled specimens. Note that, once poled, the piezo-ceramic crystallographic symmetry becomes ∞ (equivalent symmetry to 6 in terms of tensor notations). First, the open-circuit thermal conductivity λ D 33 is 0.57 times smaller than the short-circuit one λ E 33. We denote the superscripts E and D of the conductivity λ, supposing that the sample is roughly maintained "electric field E constant" and "electric displacement D (or P) constant" during a rather short experimental time period less than 10 s. Subscript "33" stands for the tensor component relating to heat flux vector q3 and x-component of the temperature gradient grad (θ) (i.e., vector) along the spontaneous polarization PS direction. Second, the unpoled specimen exhibits a crystallographic symmetry of isotropic cubic 3 (different from the poled specimens), and the thermal conductivity λ u 11 shows the smallest value, 25% less than that of the poled open circuit case (see the last column of Table 1). Different from the scalar analysis, the tensor analysis is directly coupled with the material's crystallographic symmetry. The unpoled specimen cannot be discussed in parallel with the poled specimens because of the crystallographic symmetry difference.
We would like to define a new "electrothermal coupling coefficient" k λ 33 below, in the relationship between the open-circuit λ D 33 and short-circuit thermal conductivity λ E 33 such as: Let us consider a derivation of this "secondary coupling coefficient k λ 33" with using a "spacegradient model" illustrated in  (2) heat flux seems to be controlled by the electric field E (because of the difference between and ), we assume the following constitutive equations, resembling Equation (38a,b) above, but taking the gradient with respect to the heat flux direction x: If we assume that the room temperature is reasonably low from the Curie temperature, enough to satisfy ( ), ( ) ≪ ′, (that is, the pyroelectric coefficient and permittivity change with temperature (or space coordinate x) is reasonably small), we can approximate the following equations: We adopt another assumption, that is, D-constant along the thermal flow. From = + = 0 , we obtain − ( ) = = , where = , drift velocity of electricdisplacement-charge. We can transform the gradient equations into time dependence: Since means the electric-displacement-current, Equation (77b) stands for the dynamic The reader can understand that Equation (79)  problem is how to integrate the volume and unit time calibration into the term to satisfy all units. We will discuss a more microscopic "phonon model" in Section 4 on the heat flow in conjunction with the above phenomenological explanation.

Thermal Conductivity in Pb-Free Piezoelectrics
In 2006, European Union started RoHS (Restrictions on the use of certain Hazardous Substances), which explicitly limits the usage of lead (Pb) in electronic equipment. Basically, we may need to regulate the usage of lead zirconate titanate (PZT), currently the most-widely used piezoelectric ceramics, in the future. Because of this background, the research on Pb-free piezoelectrics is now accelerated including their high-power performance. When we compare the maximum vibration velocity in Pb-free piezoelectric with PZT's, 1 m/s (rms) vs. 0.6 m/s in k31 type plate specimens. Table 2 summarized the preliminary data on the thermal conductivity in the NKNbased material with the PZT's value (both unpoled samples). Pb-free piezoceramics such as (Na,K)NbO3-and (Bi,Na)TiO3-based materials show much higher maximum vibration velocity than the PZT's [26,27]. Much larger thermal conductivity in the NKN-based material than the PZT's may also contribute to this good high-power performance in NKN-ceramics because of the definition of the maximum vibration velocity, as we have discussed in Figure 1 in Section 1.

Nonlinear Elastic performances and "Phonons"
We discussed the phenomenological discussions in ferroelectrics on (1) linear equilibrium in Section 2.1, (2) nonlinear equilibrium in Sections 2.1.3, and (3) Harmonic dynamic in Sections 2.2. We discuss in this section the remaining nonlinear equilibrium and (4) Anharmonic dynamic. Important nonlinear equilibrium effects include the "electrostriction", "thermal expansion", and the anharmonic dynamic lattice vibration creates the "phonon thermal conductivity".

Electrostriction and Thermal Expansion
Heat transport in dielectric solids is by way of elastic dynamic vibrations of the lattice (i.e., phonons) [28]. If the effective elastic energy potential for the crystal atoms possesses nonlinearity or anhormonicity, thermal expansion, electrostriction, and isothermal compressibility come out under equilibrium status. According to our previous paper Ref. [29], we introduce the effect derivation process briefly.
We shall use the theory of cohesive forces in ionic crystals here, proposed by Born [30]. In the interests of simplicity, we derive various formulae for two-ion rock-salt structure. We start from the potential function involving an inverse power type of repulsive energy of the form: where M is the Madelung constant for Coulombic energy, N the coordination number and b is a potential constant for the quantum mechanical energy. Expanding the potential function around the equilibrium position ( = ⁄ ), we obtain the form as a function of ∆ (=r − r0): . Using the Boltzmann distribution for Gibbs energy ∆ ± = ∆ ± ∆ , the average equilibrium separation at high temperatures under an applied electric field E is as follows: where subscripts ± denote the ion pairs in the positive and negative directions of the electric field respectively. The strain is therefore given by: where the first term represents thermal expansion and the second term, electrostriction. The reader can understand that both coefficients are originated from " ", which is an anharmonic term, and when the crystal is purely harmonic, neither thermal expansion nor electrostriction occur. Knowing the polarization P, given by, the thermal expansion coefficient and electrostrictive Q coefficient can be calculated Finally, using the "First Law" of thermodynamic of = − (p: hydrostatic pressure, V: volume), the isothermal compressibility becomes: Changing the variable from V to r and developing the result leads to the equation Comparing Equations (86) and (88), we obtain the relationship, ∝ . The compressibility and electrostrictive coefficient are induced strain originated from the same crystal lattice anharmonic potential against the force and electric field.

Lattice Vibration and Phonon
We may speculate the physical reasons of the thermal conductivity difference under short-and open-circuit conditions from the aspects of phonon mode scattering, orientation of domains, and elastic compliance in this section. Remember the fact that ferroelectric-originated piezoelectric materials exhibit the phase transition associated with phonon mode "softening" (in particular, optical phonon mode) [14]. Using the "first-principle" calculations, Zhang et al. obtained elastic compliances, heat capacity, and thermal expansion coefficient, based on phonon dispersion curves [31]. Though the detailed discussion is beyond our target in this paper, we introduce here briefly basic concepts on the lattice vibration and "phonon", based on books by Ashcroft and Mermin [28] and Ishii [32].

Lattice Vibration-One-Atom Chain Model
We consider the simplest one-atom (mass M) 1D chain (number of atoms N) connected by a harmonic spring (spring constant K), as shown in Figure 10a (Note two end atoms in a free condition). When we denote the displacement of the n-th atom from the equilibrium position as , we can obtain the following set of Newton equations: Assuming a harmonic vibration, ∝ , where ω is the frequency, Thus, we can write the following matrix formula: where u is an N-Dimensional vector with the components of , and [ ] is a symmetric matrix with the components: Here, we used Kronecker , which means that = 1 for ≠ ; = 0 for = .
We will now solve the "eigenvalue equation" of Equation (90). Knowing the matrix component relation where k, a are the wave number and the lattice unit size. The solution of Equation (93) is expressed as Then, we can express a general solution of Equation (92) as The "free-end" boundary condition gives the relations: To realize the meaningful solution (or A and B values), at least the following "determinant = 0" should be satisfied: This requirement determines the solutions regarding the k value and for , Taking the normalizing factor determined by ∑ = 1 , we obtain the eigenvectors and eigen frequencies as follows: When − 2 = 2, unique solution = 0, (0) = exists. However, when − 2 > 2, no solution exists to satisfy the "free-end" boundary condition. ≡ is called "wave number", and the relationship between the frequency and wave number (Equation (102)) is called "dispersion relation". As shown in Figure 11a, ∝ is valid only for a low k range (long wavelength λ); that is, a wave packet shape distorts significantly via frequency dispersion. This is the reason why " ( ) . " curve is called "dispersion relation".
The reader now understands that there are N lattice vibration eigenmodes for N eigenvalues in the N-chain crystal lattice. Each mode is identified as a "phonon". Because the actual crystal size includes Avogadro number (10 23 per mol) of atoms of N, the energy gap among phonons is quantummechanically small. According to the Gibbs statistical dynamics for the equilibrium status, thermodynamical properties in a crystal lattice is determined by these eigenvalues of the phonons. The internal energy U is expressed by using quantum-mechanically small unit ℏ as Here, ℏ is the Planck constant, the Boltzmann constant, and T is temperature. The reader remembers in your high school physics that Einstein introduced "ℏ " quantum energy unit for "photon" introduction. Equation (103) comes from the same philosophical concept.
We usually take N ∞ in an ideal crystal, and exists densely (N modes) just 0 ≤ ≤ = 2 ⁄ . In the case of large crystal (N  ∞), we occasionally use so-called "periodic boundary condition", instead of "end-free" condition; that is, the end atom (N-th) is connected with the first atom with a spring K to the lattice chain as a "ring". In this case, Equation (91) is modified as We solve the above difference equation under the boundary condition − = − = 0.
The results are slightly different from Equations (101) and (102): Here, is the Gauss symbol, which means the maximum integer number ≤ 2 ⁄ . The key difference due to the boundary condition change is the doubled unit cell; that is, even and odd number position atoms are not identical. The wave number ≡ exhibits doubledegenerated eigenvalues, which correspond to the following two eigen functions: These waves correspond to left-ward and right-ward progressive waves. Taking an atomic distance "a", the wavelength and wave velocity v are given by Equation (110) means that the wave velocity depends on the wave number; thus, a wave packet with different frequency mixture distorts. Just for a long wavelength range ( ≪ 1), the wave velocity approaches = 2 and ∝ is valid.

Lattice Vibration-Two-Atom Chain Model
We now consider the two-atom (mass M and m, rock-salt-like) 1D chain (total number of atoms 2N) connected by a harmonic spring (spring constant K), as shown in Figure 10b [28,32]. When we denote the displacement of the n-th atom from the equilibrium position as , the dynamic equations are expressed as Mass m (or M) is located at an even (or odd) number lattice position, and assuming the solutions Equation (111) can be calculated as In the matrix notation, where The lattice chain has a unit cell of size 2a with M and m atoms, and and are the "polarization (displacement) vectors" of m (even number position) and M (odd number position), respectively. From Equation (114), the eigenvalue ω should satisfy The solution of the above equation is given by We assume there are N unit cells with two atoms (M and m) in unit cells with a unit cell size of 2a, and the "period boundary condition" is adopted in this model. For N  ∞, k will vary from 0 to 2π. At = 0, while, at = , from ± (0) = + ± − [supposing M > m] Figure 11b illustrates the dispersion curves, plotted in the ka range of − (because of the period symmetry). You can initially find frequency gap between 2 ⁄ and 2 ⁄ , where no solution (i.e., phonon mode) exists. If we take the case of M = m, this gap disappears and the dispersion curve becomes equivalent to that of one-atom model in Figure 11a, by considering the unfold shape at = . ( ) corresponds to the upper curve, which is called "optical branch". This branch approaches to with → 0: Around = 0, the "polarization vector" e satisfies which indicates that the two species of atoms move in opposite directions from each other by keeping the center of gravity position. When we consider ionic crystal such as rock-salt (or PZT), Na + and Clions move in opposite directions, leading to the electric polarization deviation. This phonon mode can easily couple with the external electric field such as "optical light" electric wave; thus, this is called "optical branch." The reader can imagine that this interaction may be a part of the origin of "electrothermal coupling" in ferroelectric materials. On the contrary, the lower branch ( ) of in Figure 11b is called "acoustic branch". This branch approaches ∝ with → 0; a long-wavelength vibration transfer with a constant wave speed. The polarization vector e satisfies the following condition at k = 0: Two species of atoms displace in the same direction with the same amplitude.

Phonon Spectral Density
As expressed in Equation (99), the phonon modes are uniformly distributed in terms of wave number k. However, in terms of frequency ω, the distribution is not uniform, but non-uniform. Due to the sinusoidal relationship, ∝ , the flat parts of the dispersion curve in Figure 11a,b condense phonon modes significantly. The number of phonon eigenvalues which exist between the frequency and ( + ) is denoted as ( ) , where N is the total number of atoms and ( ) is called "spectral density". Figure 11a,b inserted a modified phonon spectral density (i.e., frequency square spectral density), ( ) = , for the reader's reference.

Anharmonic Phonon Modes and Wave Packet
Thermal energy can be stored in the lattice vibration normal phonon modes as a vibration amplitude. However, in a perfectly harmonic crystal, the phonon states are stationary, wide-spread in the crystal uniformly as a wave (not an isolated wave packet or particle-like phonon). Since phonon vibration is not disturbed, we can say the thermal resistance is zero, or thermal conductivity is infinite [28]. However, when the atomic energy potential includes nonlinear terms, the cubic anharmonic term just introduced in Equation (81) or the quartic term in Equation (6), the resonance frequency should include at least 2 or 3 , higher order harmonic modes. According to the uncertainty of the resonance frequency from the definite = to ∆ = − ∨, the phonon distribution changes from ∞ wide spread state to a sort of packet status. Figure 12 demonstrates a phonon wave packet generation from the anharmonic lattice vibration. From Equation (101), we assume the normal phonon eigenfunctions as Taking into account a modified "phonon spectral density" (i.e., frequency square spectral density), ( ) as the summation weight, and the low frequency acoustic mode dispersion relation ∝ , we calculate the sum of ( )'s, which exist in ∆ Individual wave forms (time dependence) ( ) are illustrated for = 1, 2, 3 and 0.2 in Figure  12, for the reader's reference. Here, we took the atomic distance a = 1 as a normalized unit. One, two, and three sinusoidal wave forms and a rather flat long-wavelength shape ( = 0.2) are observed. When we sum up wave forms up to = 3, the top wave packet can be obtained. Note that, if we add up to ∞ ideally, the wave shape should theoretically be a square pulse inserted in Figure 12. Once multiple phonons are generated synchronously via the anharmonic potential interactions, the phonon energy localization is expected (not uniformly distributed as a separated harmonic vibration), leading to the concept of wave-to-particle transition (i.e., the reason of the terminology introduction of "phonon" from the elastic lattice vibration wave). If we rephase the above argument in a reverse way, by having an excess of phonons with similar directed group velocities (i.e., by touching a heat source on the solid crystal), the anharmonic part of the ionic interaction plays the wave packet into creation, destruction or scattering of various frequency phonons, which corresponds to "thermal resistivity" (i.e., inverse of thermal conductivity).

Heat Flow and Electrothermal Coupling
First, the sound velocity of the longitudinal wave is much faster than the transverse waves in dielectric solid. Thus, thermal conductivity seems to be determined primarily by the speed of longitudinal phonons. Second, phonons in the acoustical branch dominate the phonon heat conduction as they have greater energy dispersion and a greater distribution of phonon velocities. Additional "optical modes" can also be caused by the presence of charge or mass at lattice points. Though it is implied that the sound velocity of the optical modes is low and their contribution to the lattice thermal conduction is small in general [33], the optical mode contribution seems to also be important in ferroelectric materials such as PZT's because of their low-frequency soft-phonon modes, which are highly interactive with the electric field. Third, when the phonon number ´ deviates from the equilibrium value ´ , thermal current arises as expressed by where is the energy transport velocity of phonons, equivalent to v in Equation (77a,b). There are two mechanisms to cause time variation of ´ in a particular region, the thin segmented disk in Figure  9: (1) The number of phonons that diffuse into and diffuse out between the considering region and neighboring regions, and (2) phonon decay inside the same region into other phonons. A special form of the Boltzmann equation provides Time variation in the latter term is described with a relaxation time (τ) approximation. 28) Assuming the spirit of the Drude model, the phonon collision occurs in a distance = ( : energy transport velocity) in terms of the relaxation time τ. Though the microscopic origin is not clear, this relaxation time concept has also been integrated in Section 2.2: When the "steady state" condition is assured, that is, ´ = 0, the local thermal equilibrium (i.e., dynamic equilibrium) is assumed, leading to the equation Using the relaxation time approximation for the Boltzmann equation and assuming steady-state conditions, the phonon thermal conductivity λ can be determined: Let us now discuss the thermal conductivity changes under different electrical constraints. Phonon transmission speed is proportional to the thermal conductivity λ. First, because of the random orientation of domains in the unpoled sample, it is expected that there will be the most phonon scattering at the domain boundaries in this specimen (i.e., acoustic mismatch at the grain boundaries). Therefore, it will have the lowest thermal conductivity. Second, by increasing the degree of domain orientation by the electric poling process, less scattering is expected in the poled specimens and thermal conductivity will be larger for the poled material (λ E 33 and λ D 33), according to the crystallographic symmetry change from 3 to ∞ . Third, because the elastic compliance under the short-circuit condition (s E 33) is softer than the elastic compliance under open-circuit conditions (s D 33) via the relation s E 33(1 − k33 2 ) = s D 33, and the sound velocity is smaller than , the lattice vibration and domain wall dynamic motion are expected to be larger in the polarization direction under the short-circuit condition. In the thermal conductivity measurement, the "depolarization field" in parallel to the spontaneous polarization direction under open-circuit condition (i.e., Dconstant), originated from the pyroelectric charge during the temperature rise, stabilizes the crystal vibration. The internal electric field in piezo-materials in general stabilizes the lattice vibration (in particular optical branch) and domain wall motion. The larger lattice vibration and domain wall dynamics in the short-circuit condition may also introduce a larger phonon transportation, i.e., larger thermal conductivity (λ E 33). This speculation can also suggest the analogous relations among two coupling factors: electromechanical k33 and secondary electrothermal k λ 33.

Thermal Analysis on Piezoelectric Transducers
Based on the aforementioned knowledge, we introduce two practical thermal analysis case studies in piezoelectric transducers in this section.
Heat generation in piezoelectric materials are originated from three losses: dielectric, elastic, and piezoelectric losses. The hysteresis curves: dielectric D vs. E, elastic x vs. X, and piezoelectric x vs. E under a stress-free condition and D vs. X under a short-circuit condition correspond dielectric, elastic, and piezoelectric coupling losses, respectively. We discuss the heat generation mechanisms of piezoelectric actuators under (1) off-resonance operation for actuator applications (under a large electric field, 1 kV/mm or higher), where the dielectric loss is the primary origin, and (2) resonance (or anti-resonance) operation for ultrasonic transducer applications (under a high vibration condition at low electric field, 100 V/mm or lower), where the intensive elastic loss is the primary heat generation factor.

Heat Generation from Multilayer Actuators
Zheng et al. reported the heat generation at an off-resonance frequency from various configurations of multilayer (ML) type piezoelectric ceramic (soft PZT) actuators [34]. Figure 13 shows a structure of the multilayer piezoelectric actuators. The temperature change with time in the actuators was monitored when driven at 3 kV/mm (high electric field) and 300 Hz (much lower frequency than the resonance frequency) (Figure 14a). The specimen temperature reached up to 140 °C, depending on the size, showing an exponential increase with the operation time lapse. Figure 14b plots the saturated temperature as a function of Ve/A, where Ve is the effective volume (electrode overlapped part, abL in the figure) and A is the all surface area. Because the temperature was uniformly generated in a bulk sample (no significant stress distribution, except for the small inactive portion of the external electrode sides), this linear relation is reasonable because the volume Ve generates the heat, which is dissipated through the area A. Thus, if we need to suppress the temperature rise, a small Ve/A design is preferred. Instead of one ML, four (1/4) small-ML's connected in parallel are preferred.

Thermal Analysis on ML Actuators
According to the law of energy conservation, the amount of heat stored in the piezoelectric, which is just the difference between the rate at which heat is generated, qg, and that at which the heat is dissipated, qd, can be expressed as where it is assumed that a "uniform temperature distribution" exists throughout the sample and V is the total volume, ρ is the mass density, and cp is the specific heat capacity (per mass) of the specimen. The heat generation in the piezoelectric is attributed to losses. Thus, the rate of heat generation, qg, is expressed as: where w is the loss per driving cycle per unit volume, f is the driving frequency, and Ve is the effective volume of active ceramic (no-electrode parts are omitted). According to the measurement conditions (no significant stress in the sample at off-resonance), this w may correspond primarily to the dielectric hysteresis loss (i.e., P-E hysteresis), we, which is expressed in terms of the intensive dielectric loss tan δ' as: If we neglect the transfer of heat through conduction, the rate of heat dissipation (qd) from the sample is the sum of the rates of heat flow by radiation (qr) and by convection (qc): where e is the emissivity of the sample, A is the sample surface area, σ is the Stefan-Boltzmann constant, To is the initial sample temperature, hc is the average "convective heat transfer coefficient". Thus, Equation (97) can be written in the form: where the quantity is the overall "heat transfer coefficient". If we assume that k(T) is relatively insensitive to temperature change (if the temperature rise is not large), solving Equation (130) for the rise in temperature of the piezoelectric sample yields: where the time constant τ is expressed as: We can understand that the specimen temperature rise follows only when the heat transfer coefficient k(T) is nearly constant.
As t  ∞, the maximum temperature rise in the sample becomes: while, as t  0, the initial rate of temperature rise is given by where wT can be regarded under these conditions as a measure of the total loss of the piezoelectric. The dependence of k(T) on applied electric field and frequency are shown in Figure 15a,b, respectively. Note that k(T) is almost constant, as long as the driving voltage or frequency is not very high (E < 1 kV/mm, f < 2 kHz). The total loss, wT, as calculated from Equation (137b) is given for three multilayer specimens in Table 3, three values of which are almost the same in less than 2% deviation. In parallel, we measured the P-E hysteresis losses under stress-free conditions. The we values obtained (Equation (130)) are also listed in Table 3 for comparison. It is intriguing that the extrinsic P-E hysteresis loss contributes more than 90% of the calculated total loss associated with the heat generated in the operating piezoelectric specimen [34,35]. We can conclude that the heat generation of the piezoelectric specimen under high-electric field off-resonance operation is primarily originated from the intensive dielectric loss factor, tan '.
(a) (b) Figure 15. Overall heat transfer coefficient, k(T), plotted as a function of applied electric field (a), and of frequency (b) for a PZT ML actuator with 7 × 7 × 2 mm 3 driven at 400 Hz.  [36]. Even though the electric field is not large, considerable heat is generated due to the large induced strain/stress at the resonance. The maximum heat generation was observed at the nodal regions for the resonance vibration, which correspond to the locations where the maximum strains/stresses are generated.
The ICAT at the Penn State University also worked on the heat generation comprehensively in rectangular piezoelectric k31 plates ( Figure 16) when driven at the resonance [37]. The temperature distribution profile in a PZT-based plate sample was observed with a pyroelectric infrared camera (FLIR Systems ThermaCAM S40 outfitted with a 200 mm lens), as shown in Figure 17, where the temperature variations are shown in the sample driven at the first (28.9 kHz) (a) and second resonance (89.7 kHz) (b) modes, respectively. The highest temperature (bright spot) is evident at the nodal line areas of the length resonance for the specimen in Figure 17a,b. This observation supports that the heat generated in a resonating sample is primarily originated from the intensive elastic loss, tanφ'. Remember that the maximum vibration velocity is defined as the velocity under which a 20 °C temperature increase occurs at the maximum temperature point (i.e., nodal line!) in the sample.

Heat Generation at the Antiresonance Mode
The resonance and antiresonance are both mechanical resonances with the electrical admittance maximum and almost zero (or maximum), respectively, when electrically excited. We can amplify the generating displacement (or vibration velocity) significantly by the factor of mechanical quality factor QA at the resonance under constant voltage drive, while by QB at the antiresonance under constant current drive. In the k31 mode, we derived the mechanical quality factors QA and QB as follows [37]: Here, tanδ33′, tanϕ11′, and tanθ31′ are intensive loss factors for ε33 X , s11 E , d31, respectively, and ΩB,31 is the normalized antiresonance frequency: . Because the intensive dielectric loss tanθ31′is larger than 33 11 (tan ' tan ') / 2 δ φ + in PZT piezoceramics, QB at antiresonance is higher than QA at resonance; that is, the antiresonance operation seems to be more efficient than the resonance drive. Figure 18a,b show temperature variations in a PZT-based plate specimen driven at the antiresoance (a) and resonance frequency (b) under the same vibration velocity (i.e, almost the same output mechanical energy), which clearly exhibit lower temperature rise in the antiresonance, than that in the resonance drive. Remember the similar resonance and antiresonance vibration modes for a reasonable electromechanical coupling factor k31 = 30% [37]. In comparison with the fundamental resonance mode, though the antiresonance mode exhibits the strain-zero lines (i.e., anti-node lines) a little inside from the plate edges for this small k31 case, the overall vibration configurations for resonance and antiresonance modes are rather close to each other. Namely, as long as the vibration velocity at the plate edge is the same, the total mechanical energy is assumed to be the same for both resonance and antiresonance operations. Numerical profiles of the temperature distribution for the A-and B-type resonance modes are shown in Figure 18c for various vibration velocity, which seems to be pseudo-sinusoidal curves in terms of the length position coordinate. Under the same vibration velocity, 550 mm/s RMS, the resonance nodal line area shows the temperature up to 100 °C, while the antiresonance mode shows the maximum around 60 °C; dramatic heat generation reduction under the antiresonance drive (a) (b) (c) Figure 18. Temperature variations in a PZT-based plate specimen observed with a pyroelectric infrared camera, when driven at the antiresoance (a) and resonance frequency (b); (c) numerical temperature profile for the A-and B-type resonance modes.

Heat Transfer Modeling
We develop a 1D heat transfer model for the k31 mode piezoelectric rectangular plate around the resonance/antiresonance frequency range. In comparison with the off-resonance model, where the heat is generated primarily from the dielectric loss and the uniform temperature distribution profile due to no particular stress in the specimen, the resonance case generates the heat originated from the intensive elastic loss on a sinusoidal stress distribution of a specimen. Because of the initial heat source being sinusoically distributed in the specimen, we need to integrate the thermal diffusivity of the piezo-ceramic in order to analyze the temperature distribution of the sample.
We set the following assumptions for developing the heat diffusion equation: (1) 1D heat conduction in the specimen.
(2) Heat generation is proportional to strain squared (i.e., elastic energy), distributed on the specimen. (3) Heat dissipation via convection (to air) and radiation. Conduction is neglected.
Using a temperature parameter T(x,t), which is defined as temperature of a sliced volume Δx from the position coordiate x to (x + Δx) at time t, the following equation is assumed: where λ is thermal conductivity [unit: W/m K], cp (c E in the k31 case) is specific heat (unit: J/kg K), ρ density, and = is called thermal diffusivity [unit: m 2 /s]. 4,7) The first term of the right-hand-side of Equation (140) describes the temperature distribution in respect of the position x, which changes the shape with time. The second term corresponds to the temperature increment caused by heat generation per unit volume (devided by ), which may exhibit a sinusoidal distribution. The third term indicates the heat dissipation proportional to the temperature difference ΔT from the ambient temperature . Regarding the heat generation, we further assume that ( ) is expressed in proportion to the square of the sinusoidally-distributing strain, as visualized in Figure 19: The heat dissipation hd is a proportional constant to [T − Tair], similar to "heat transfer coefficient" introduced in Equation (131).

Temperature Distribution Profile Change with Time
We prepared a piezoelectric k31 rectangular plate specimen with 80 × 14 × 2 mm in size with a hard PZT composition (APC 841, American Piezo Company, Mill Hall, PA, USA) for both admittance and thermal imaging observation purposes. We used a constant vibration velocity method for the measurements under 300 mm/s RMS. Figure 20a shows the temperature distribution profile change with time after driving. You can notice that the plate edge temperature increases significantly with time lapse, primarily due to the thermal diffusion in the PZT from the nodal highest point to the edge lowest temperature point. The saturated temperature distribution profile for the k31 specimen is shown in Figure 20b, which can be used for calculating the total thermal dissipation energy gh. A small temperature dent at x = 0 originated from the heat dissipation by the center sample-holding rod conduction, which is neglected in the simulation. If we adopt the most general definition of the mechanical quality factor as we can obtain Qm from Note that this is the unique approach to determine the mechanical quality factor with using merely the thermal data, without using any electrical energy information.
Shekhani et al. measured the admittance spectrum on the above PZT sample (80 × 14 × 2 mm ) with the resonance frequency at 20.04 kHz at room temperature, and obtained Qm = 507 by a 3-dBdown method on an admittance spectrum. Then, the sample was excited under the vibration velocity of 400 mm/s for 30 s, which corresponds to the heat dissipation of 11.6 W/m 2 . Figure 21 shows the Qm obtained at three frequencies slightly above the resonance frequency (20.04 kHz) [7]. The Qm values around 550 agree with the extraporated values from the above 507. Thus, we can conclude that this thermal method can determine the Qm reasonably at any frequency around the resonance and antiresonance region. An increase of Qm with increasing the frequency suggests that a maximum-Qm frequency (i.e., the highest efficiency) exists between the resonance and antiresonance frequencies [38].

Conclusions
The electrothermal, piezothermal, and piezoelectric coupling phenomena in ferroelectrics have been introduced with thermodynamics viewpoints-in particular, thermal property differences among unpoled and poled PZT's in the poling direction for open circuit and short circuit conditions. The derivation processes for piezoelectric, piezothermal and electrothermal "coupling factors" have been introduced comprehensively. We proposed a new "electrothermal" coupling factor k λ for the thermal conductivity analysis, relating the parameter change under short-and open-circuit conditions, in order to explain the fact that the short-circuit condition exhibited larger thermal diffusivity than the open-circuit condition. On the other hand, the unpoled specimen exhibits the lowest thermal conductivity. This tutorial paper was authored for providing comprehensive knowledge on equilibrium and time-dependent thermodynamics in ferroelectrics. Electrothermal coupling may be categorized as primary and secondary effects in general. The primary effect includes "pyroelectric" and "electro-caloric" phenomena, and the secondary effect includes the space-gradient (i.e., first-derivative) phenomena, such as "thermal conductivity" and local "electric displacement current". The author is expecting further theoretical expansion in this electrothermal coupling area, and industrial application products.