Electrothermal Phenomena in Ferroelectrics

Abstract

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.

1. 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-compressor-embedded 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(Zn_1/3Nb_2/3)O₃-PbTiO₃ (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³ [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 off-resonance, 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 $v_{max}^{RMS}=0.6 \mathrm{m}/\mathrm{s}$ in a rectangular plate PZT’s, which corresponds to the energy density of 40 W/cm³. 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” Q_m 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.

2. 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 Section 2.1.3 and Section 4.1, (3) harmonic dynamic in Section 2.2 and Section 4.2, and (4) anharmonic dynamic in Section 4.3.

2.1. Static/Equilibrium Ferroelectric Phenomenology

2.1.1. 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:

dU = dQ + Xdx + EdD
     ≦ TdS + Xdx + EdD

(1)

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 G₁ = U − TS − Xx as examples in this paper; that is,

$$dA=-SdT+Xdx+Ed$$

(2a)

$$dG=-SdT-xdX-DdE$$

(2b)

$$d{G}_1=-SdT-xdX+EdD$$

(2c)

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 G₁ 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 ${\left(\frac{\partial G_1}{\partial P}\right)}_{T,X}$.

We can obtain the famous “Maxwell relationship” examples from Equation (2) as follows, taking the second partial derivative:

$${\left(\frac{\partial D}{\partial X}\right)}_{T,E}={\left(\frac{\partial }{\partial X}{\left(\frac{-\partial G}{\partial E}\right)}_{T,X}\right)}_{T,E}={\left(\frac{\partial }{\partial E}{\left(\frac{-\partial G}{\partial X}\right)}_{x,T}\right)}_{T,D}={\left(\frac{\partial x}{\partial E}\right)}_{T,X}$$

(3a)

$${\left(\frac{\partial X}{\partial D}\right)}_{T,x}={\left(\frac{\partial }{\partial D}{\left(\frac{\partial A}{\partial x}\right)}_{X,D}\right)}_{T,x}={\left(\frac{\partial }{\partial x}{\left(\frac{\partial A}{\partial D}\right)}_{x,T}\right)}_{T,D}={\left(\frac{\partial E}{\partial x}\right)}_{T,D}$$

(3b)

$${\left(\frac{\partial S}{\partial E}\right)}_{T,X}={\left(\frac{\partial }{\partial E}{\left(\frac{-\partial G}{\partial T}\right)}_{X,E}\right)}_{T,X}={\left(\frac{\partial }{\partial T}{\left(\frac{-\partial G}{\partial E}\right)}_{X,T}\right)}_{T,X}={\left(\frac{\partial D}{\partial T}\right)}_{X,E}$$

(3c)

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.

2.1.2. 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 − T_R (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]:

$$\begin{gathered} G\left(T,X,E\right)=G_0+\left(\frac{\partial G}{\partial T}\right)\theta +\left(\frac{\partial G}{\partial X}\right)X+\left(\frac{\partial G}{\partial E}\right)E \\ \frac{+1}{2}\left(\frac{{\partial }^{2}G}{\partial T^2}\right){\theta }^2+\frac{1}{2}\left(\frac{{\partial }^{2}G}{\partial X^2}\right)X^2+\frac{1}{2}\left(\frac{{\partial }^{2}G}{\partial E^2}\right)E^2 \\ +\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)\theta X+\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)\theta E+\left(\frac{{\partial }^{2}G}{\partial X\partial E}\right)X{E}^{ } \end{gathered}$$

(4)

Taking into account $dG=-SdT-xdX-DdE$, we obtain first the relations, ${\left(\frac{\partial G}{\partial T}\right)}_{\theta ,X,E=0}=-S_0$, ${\left(\frac{\partial G}{\partial X}\right)}_{\theta ,X,E=0}=-x_0$ and ${\left(\frac{\partial G}{\partial E}\right)}_{\theta ,X,E=0}=-D_0$. Since $S_0$ is the entropy density at $\theta ,X,E=0$, we take this as the original value and set it to $S_0=0$. The values $x_0$ and $D_0\left(\approx P_0\right)$ 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

$$S=-\left(\frac{\partial G}{\partial T}\right)=-\left(\frac{{\partial }^{2}G}{\partial T^2}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)X-\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)E$$

(5a)

$$x=-\left(\frac{\partial G}{\partial X}\right)=-\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial X^2}\right)X-\left(\frac{{\partial }^{2}G}{\partial X\partial E}\right)E$$

(5b)

$$D=-\left(\frac{\partial G}{\partial E}\right)=-\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial X\partial E}\right)X-\left(\frac{{\partial }^{2}G}{\partial E^2}\right)E$$

(5c)

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 G₁ for the phase transition discussion, as in Section 2.1.1.

2.1.3. 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:

G₁(P) = G₁(0) + a₁P + a₂P² + a₃P³ + a₄P⁴ + a₅P⁵ + a₆P⁶ + ⋯

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, G₁(P) = G₁(−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.

2.1.3.1. Second-Order Phase Transition

We start from the simplest nonlinear phenomenon “second-order phase transition” model [14]:

G₁(P,T) = (1/2)αP²+ (1/4)β P⁴

(6)

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² 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 T₀, called Curie–Weiss temperature. From this concept, we assume α as a linear relation:

$$\alpha = (T - T_0)/{\epsilon }_0 C [{\epsilon }_0 = 8.854 \times {10}^{-12} \mathrm{F}/\mathrm{m}]$$

(7)

where C is taken as a positive constant called the Curie–Weiss constant, and T₀ is equal to or lower than the actual transition temperature T_C (Curie temperature), in general. This temperature dependence can also be interpreted as the adoption of the “TP²” electrothermal coupling term [$\left(\frac{{\partial }^3{G}_1}{\partial T{\partial }^{2}P}\right)\theta P^2$] in the Taylor expansion of the G₁ energy. The temperature parameter θ here is the deviation from the Curie–Weiss temperature T₀. From $d{G}_1=-SdT-xdX+EdP$ [$P\approx D$ in our discussion], the equilibrium polarization under an electric field E should satisfy the condition:

(∂G₁/∂P) = E = α P + β P³ = P∙(α + β P²)

(8)

With no electric field applied, Equation (8) provides two cases: (i) P = 0 and (ii) P² = −α/β. The trivial solution in (i) case corresponds to a paraelectric (high temperature) state, while the polarization solutions (P_S = $\pm \sqrt{-\alpha /\beta }$ 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 T_C (in this case also T₀). 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 T_C. 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 +P_S (=+$\sqrt{-\alpha /\beta }$ or −P_S (=$-\sqrt{-\alpha /\beta }$ should be equal. Thus, with decreasing the temperature, passing through T_C, +P_S domains and −P_S 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 $E_C=\sqrt{\frac{-4}{27}{\alpha }^3/\beta }$ 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.

2.1.3.2. 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:

G₁(P,X,T) = (1/2)α P² + (1/4)β P⁴ + (1/6)γ P⁶
    − (1/2)s X² − Q P²X, (α = (T − T₀)/ε₀ C)

(9)

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, $d{G}_1=-SdT-xdX+EdD$, leads to

E = (∂G₁/∂P) = α P + β P³ + γ P⁵ − 2Q PX minus sign

(10a)

x = − (∂G₁/∂X) = sX + QP²

(10b)

When the external stress X is zero, the following equations are derived:

E = α P + β P³ + γ P⁵

(11a)

x = QP²

(11b)

1/ε₀ ε = (∂E/∂P) = α + 3 β P² + 5 γ P⁴

(12)

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² = (−β +$\sqrt{{\beta }^2-4\alpha \gamma }$)/2γ.

(i) Paraelectric phase: P_S = 0 or P = ε₀ε E (under small E)

Permittivity: ε = C/(T − T₀) (Curie-Weiss law)

(13)

Electrostriction: x = Q (ε₀ε)² E²

(14)

(ii) Ferroelectric phase: P_S² = (−β + $\sqrt{{\beta }^2-4\alpha \gamma }$)/2γ or P = P_S + ε₀ε E (under small E)

x = Q(P_S + ε₀ε E)² = QP_S² + 2 ε₀ε QP_S E (+Q (ε₀ε)²E²)

(15)

where we define the spontaneous strain x_S and the piezoelectric constant d as:

Spontaneous strain: x_S = QP_S²

(16)

Piezoelectric constant: d = 2 ε₀ ε QP_S

(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 > T₀

P_S = 0

(18a)

1/ε₀ ε = α, then ε = C/(T − T₀) (Curie-Weiss Law)

(18b)

(ii) Ferroelectric phase—T < T₀

$$P_S= \pm \sqrt{\frac{-\alpha }{\beta }}=\pm \sqrt{\left(T_0-T\right)/{\epsilon }_{0}C\beta }$$

(19a)

1/ε₀ ε = α + 3 β P² = −2α, then ε = C/2(T₀ − T)

(19b)

x_S = QP_S² = Q (T₀ − T)/ε₀Cβ

(20)

From Equations (18b) and (19b), the piezoelectric constant is obtained as

$$\begin{gathered} d= 2 {\epsilon }_0\epsilon Q P_S \\ = Q\sqrt{{\epsilon }_{0}C/\beta }\sqrt{T_0-T} \end{gathered}$$

(21)

2.1.4. Isothermal Process—Piezoelectric Coupling

2.1.4.1. Thermodynamical Meaning of Piezoelectric Constant

We continue the “linear energy handling” in Section 2.1.2. When the temperature is constant (i.e., isothermal), $\theta =0$ in Equation (5b,c), we can obtain the “intensive” parameter-based piezoelectric constitutive equations:

$$x=s^{E}X+dE$$

(22a)

$$D=dX+{\epsilon }_0{\epsilon }^{X}E$$

(22b)

where the following notations are used, and denoted as $s^E$ elastic compliance under constant E, ${\epsilon }_0{\epsilon }^X$ dielectric permittivity under stress free:

$$\{\begin{array}{c}{s}^E=-\left(\frac{{\partial }^{2}G}{\partial X^2}\right)\\ {\epsilon }_0{\epsilon }^X=-\left(\frac{{\partial }^{2}G}{\partial E^2}\right)\\ d=-\left(\frac{{\partial }^{2}G}{\partial X\partial E}\right)\end{array}$$

(23)

The Maxwell relation, Equation (3a), ${\left(\frac{\partial D}{\partial X}\right)}_{T,E}={\left(\frac{\partial x}{\partial E}\right)}_{T,X}$ verifies that the piezoelectric coefficient d in Equation (22a,b) are thermodynamically the same.

When we start from the Helmholtz free energy A ($dA=-SdT+Xdx+EdD$), by taking a similar Taylor expansion approach, we obtain another set of piezoelectric constitutive equations in terms of “extensive” parameters, x and D:

X = c^D x − h D,

(24a)

E = −h x + κ₀κ^x D,

(24b)

where c^D is elastic stiffness under constant D, and κ₀κ^x is inverse permittivity (${\kappa }_0=1/{\epsilon }_0$) under strain free condition, and these coefficients are expressed by:

$$\{\begin{array}{c}{c}^D=\left(\frac{{\partial }^{2}A}{\partial x^2}\right)\\ {\kappa }_0{\kappa }^x=\left(\frac{{\partial }^{2}A}{\partial D^2}\right)\\ h=-\left(\frac{{\partial }^{2}A}{\partial x\partial D}\right)\end{array}$$

(25)

The Maxwell relation, ${\left(\frac{\partial X}{\partial D}\right)}_{T,x}={\left(\frac{\partial E}{\partial x}\right)}_{T,D}$ in Equation (3b), verifies that the inverse piezoelectric coefficients h in Equation (24a,b) are thermodynamically the same.

2.1.4.2. Electromechanical Coupling Factor

The term, electromechanical coupling factor k, is defined as the square value k² being the ratio of the converted energy over the input energy: when electric to mechanical

k² = (Stored mechanical energy/Input electrical energy),

(26a)

and when mechanical to electric,

k² = (Stored electrical energy/Input mechanical energy)

(26b)

Let us derive Equation (26a) first practically, when an external electric field E₃ 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)ε₀ε₃^XE₃² from Equation (4), and the output strain generated by E₃ should be d₃₃E₃ from Equation (22a). Since the converted/stored mechanical energy is obtained as (1/2 s₃₃^E) x₃², we obtain:

k₃₃² = [(1/2)−(d₃₃ E₃)²/s₃₃^E]/[(1/2) ε₀ε₃^XE₃²]
= d₃₃²/ε₀ε₃^X∙s₃₃^E

(27a)

Let us now consider Equation (26b), when an external stress X₃ is applied to a piezoelectric material in a pseudo-static process. Refer to Figure 4b. Under short-circuit condition (E₃ = 0), the input mechanical energy must be equal to (1/2)s₃₃^EX₃² from Equation (4), and the electric displacement D₃ (or polarization P₃) generated by X₃ should be equal to d₃₃X₃ from Equation (22b). This D₃ 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 ε₀ ε₃^X) D₃², we obtain:

k₃₃² = [(1/2 ε₀ε₃^X) (d₃₃X₃)²]/[(1/2)s₃^EX₃²]
= d₃₃²/ε₀ε₃^X∙s₃₃^E      

(27b)

It is essential to understand that the electromechanical coupling factor k (or k², which has a physical meaning of energy transduction/conversion rate), can be exactly the same for both converse (27a) and direct (27b) piezoelectric effects. The conditions under constant X (free stress) or constant E (short-circuit) are considered to be non-constrained.

2.1.4.3. 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:

$$\left(\frac{1}{2}\right){\epsilon }^X{\epsilon }_0{E}_0^2=\left(\frac{1}{2}\right){\epsilon }^x{\epsilon }_0{E}_0^2+\left(\frac{1}{2s^E}\right)x^2=\left(\frac{1}{2}\right){\epsilon }^x{\epsilon }_0{E}_0^2+\left(\frac{1}{2s^E}\right){\left(d{E}_0\right)}^2$$

such that:

$$\begin{gathered} {\epsilon }^X{\epsilon }_0={\epsilon }^x{\epsilon }_0+\left(\frac{d^2}{s^E}\right) \mathrm{or} \\ {\epsilon }^x={\epsilon }^X\left(1-k^2\right)\hspace{1em}\hspace{1em} [k^2=\frac{d^2}{{\epsilon }^X{\epsilon }_0{s}^E}] \end{gathered}$$

(28a)

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:

$$\left(\frac{1}{2}\right)s^E{X}_0^2=\left(\frac{1}{2}\right)s^D{X}_0^2+\left(\frac{1}{2}\right){\epsilon }^X{\epsilon }_0{E}^2=\left(\frac{1}{2}\right)s^D{X}_0^2+\left(\frac{1}{2}\right){\epsilon }^X{\epsilon }_0{(\mathrm{d}/{\epsilon }_0{\epsilon }^{\mathrm{X}})}^2{X}_0^2$$

which leads to:

$$s^E=s^D+\left(\frac{d^2}{{\epsilon }^X{\epsilon }_0}\right) \mathrm{or} s^D=s^E\left(1-k^2\right) [k^2=\frac{d^2}{{\epsilon }^X{\epsilon }_0{s}^E}]$$

(28b)

In principle, if we measure the permittivity in a piezoelectric specimen under stress-free and completely-clamped conditions, we can obtain ${\epsilon }^X$ and ${\epsilon }^x$, respectively. However, in practice, ${\epsilon }^x$ 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 $s^E$ and $s^D$, respectively. However, in practice, $s^D$ 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 ${\epsilon }^x$ or open-circuit D-constant $s^D$ 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:

$${\epsilon }^x/{\epsilon }^X=\left(1-k^2\right)$$

(29a)

$$s^D/s^E=\left(1-k^2\right)$$

(29b)

where

$$k^2=\frac{d^2}{s^E{\epsilon }^X{\epsilon }_0}$$

(30)

We can also write equations of a similar form for the corresponding reciprocal quantities:

$${\kappa }^X{\kappa }^x=\left(1-k^2\right)$$

(31a)

$$c^E/c^D=\left(1-k^2\right)$$

(31b)

where, in this context,

$$k^2=\frac{h^2}{c^D({\kappa }_0{\kappa }^x)}$$

(32)

This new parameter k in Equation (32) is also the electromechanical coupling factor in the “extensive” parameter description, and identical to the k in Equation (30). Note the k expression derivation from the piezoelectric constitutive equations, Equations (22a,b) and (24a,b):

$$k^2=\frac{{\left(Couplingfactor\right)}^2}{\left(Productofthediagonalparameters\right)}=\frac{d^2}{s^E{\epsilon }^X{\epsilon }_0}=\frac{h^2}{c^D{(}^x/{\epsilon }_{0 })}$$

2.1.5. 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, $E=0$ in Equation (5a,b), we can obtain the following equations:

$$S=-\left(\frac{{\partial }^{2}G}{\partial T^2}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)X$$

(33a)

$$x=-\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial X^2}\right)X$$

(33b)

or

$$S=\frac{c_p}{T}\theta +{\alpha }_{L}X$$

(34a)

$$x={\alpha }_L\theta +s^{E}X$$

(34b)

where the following notations are used, and denoted as $c_p$ specific heat capacity under X = 0 and E = 0, $s^E$ elastic compliance under constant E:

$$\{\begin{array}{c}{c}_p=-T{\left(\frac{{\partial }^{2}G}{\partial T^2}\right)}_{X,E}\\ s^E=-{\left(\frac{{\partial }^{2}G}{\partial X^2}\right)}_{E,T}\\ {\alpha }_L=-{\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)}_E\end{array}$$

(35)

Let us discuss here the diagonal expansion coefficient $-\left(\frac{{\partial }^{2}G}{\partial T^2}\right)$ 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 $c_p$ is defined by

$$c_p=\frac{\partial q}{\partial T}=T{\left(\frac{\partial S}{\partial T}\right)}_{X,E}=-T{\left(\frac{{\partial }^{2}G}{\partial T^2}\right)}_{X,E}$$

(36)

It is noteworthy that the piezothermal coefficient ${\alpha }_L$ is usually called “linear thermal expansion coefficient” because $-{\left(\frac{{\partial }^{2}G}{\partial T\partial X}\right)}_E=\left(\frac{\partial x}{\partial T}\right)$. This piezothermal coefficient ${\alpha }_L$ 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 ${\alpha }_L$ 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 (34a,b) by

$${k^{PT}}^2=\frac{{\left(Couplingfactor\right)}^2}{\left(Productofthediagonalparameters\right)}=\frac{{{\alpha }_L}^2}{\left(\frac{c_p}{T}\right)s^E}$$

(37)

We can derive the completely “clamped” (strain x free) specific heat capacity and the “adiabatic” elastic compliance in terms of ${k^{PT}}^2$ theoretically as follows, though these may not be useful physical parameters: Under strain-free condition (x = 0), we obtain $X=\frac{-{\alpha }_L}{s^E}\theta$ from Equation (34b). Since Equation (34a) gives $S=\frac{c_p}{T}\theta +{\alpha }_L\left(\frac{-{\alpha }_L}{s^E}\theta \right)=\frac{c_p}{T}\left(1-{k^{PT}}^2\right)\theta$, we finally obtain $c_p^x=c_p^X\left(1-{k^{PT}}^2\right)$. On the other hand, under the adiabatic condition (S = 0), we obtain $\theta =\frac{-T}{c_p}{\alpha }_{L}X$ from Equation (34a). Since Equation (34b) gives $x={\alpha }_L\left(\frac{-T}{c_p}{\alpha }_{L}X\right)+s^{E}X=s^E\left(1-{k^{PT}}^2\right)X$, we obtain $s_{adia}^{E,S}{s}_{iso}^{E,T}\left(1-{k^{PT}}^2\right)$.

2.1.6. 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, $X=0$ in Equation (5a,c), we can obtain the following equations:

$$S=-\left(\frac{{\partial }^{2}G}{\partial T^2}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)E$$

(38a)

$$D=-\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)\theta -\left(\frac{{\partial }^{2}G}{\partial E^2}\right)E$$

(38b)

or

$$S=\frac{c_p}{T}\theta -pE$$

(39a)

$$D=-p\theta +{\epsilon }_0{\epsilon }^{X}E$$

(39b)

where the following notations are used, and denoted as $c_p^E$ specific heat capacity (per unit mass) under X = 0 and E = 0, ${\epsilon }_0{\epsilon }^X$ permittivity under constant stress X:

$$\{\begin{array}{c}{c}_p=-T{\left(\frac{{\partial }^{2}G}{\partial T^2}\right)}_{X,E}\\ {\epsilon }_0{\epsilon }^X=-{\left(\frac{{\partial }^{2}G}{\partial E^2}\right)}_{T,X}\\ p={\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)}_X\end{array}$$

(40)

The primary electrothermal coupling coefficient $p$ is usually called “pyroelectric coefficient”, defined by

$$p={\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)}_X=-{\left(\frac{\partial P}{\partial T}\right)}_X,$$

(41)

where we used the relation ${\left(\frac{\partial G}{\partial E}\right)}_X=D\left(\approx P\right)$.

2.1.6.1. Constraint Specific Heat Capacity

In Equation (36), we introduced $c_p$ 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 (39a) with respect to T by keeping X = D = 0, and $c_p=c_p^E$,

$${\left(\frac{\partial S}{\partial T}\right)}_{X,D}=\frac{c_p^E}{T}{\left(\frac{\partial \theta }{\partial T}\right)}_{X,D}-p{\left(\frac{\partial E}{\partial T}\right)}_{X,D}$$

(42)

From Equation (38b) at D = 0, we obtain

$$\begin{gathered} E=\frac{p}{{\epsilon }_0{\epsilon }^X}\theta \\ {\left(\frac{\partial E}{\partial T}\right)}_{X,D}=\frac{p}{{\epsilon }_0{\epsilon }^X} \end{gathered}$$

If we denote $c_p^D=T{\left(\frac{\partial S}{\partial T}\right)}_{X,D}$ and ${\left(\frac{\partial \theta }{\partial T}\right)}_{X,D}=1$, we can obtain

$$c_p^D=c_p^E-\frac{T{p}^2}{{\epsilon }_0{\epsilon }^X}=c_p^E\left[1-\frac{p^2}{\left(\frac{c_p^E}{T}\right){\epsilon }_0{\epsilon }^X}\right]=c_p^E\left(1-{k^{ET}}^2\right)$$

(43)

Here, we denoted the primary “electrothermal coupling factor” $k^{ET}$ from Equation (39a,b) as

$${k^{ET}}^2=\frac{{\left(Couplingfactor\right)}^2}{\left(Productofthediagonalparameters\right)}=\frac{p^2}{\left(\frac{c_p^E}{T}\right){\epsilon }_0{\epsilon }^X}$$

(44)

It is important to note that Equation (43) is analogous to Equation (29b) to correlate the D-constant and E-constant parameters in terms of “coupling factors”, ${k^{ET}}^2$ and $k^2$.

2.1.6.2. 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 (39b), isothermal permittivity ($\theta =0$) is given by

$${\epsilon }_0{\epsilon }^{X,T}={\left(\frac{\partial D}{\partial E}\right)}_{X,T}$$

(45)

Under adiabatic condition, by putting S = 0 in Equation (39a), we obtain

$$\theta =\frac{pT}{c_p^E}E$$

(46)

Inserting Equation (46) into Equation (39b),

$$D=-p\frac{pT}{c_p^E}E+{\epsilon }_0{\epsilon }^{X,T}E={\epsilon }_0{\epsilon }^{X,T}\left[1-\frac{T{p}^2}{c_p^E{\epsilon }_0{\epsilon }^X}\right]E$$

(47)

Thus, “adiabatic (S =constant) permittivity” is related with “isothermal (T = constant) permittivity” again by using primary electrothermal coupling factor $k^{ET}$ as:

$${\epsilon }_0{\epsilon }^{X,S}={\epsilon }_0{\epsilon }^{X,T}\left(1-{k^{ET}}^2\right)$$

(48)

2.1.6.3. Electrocaloric Effect

Equation (39a) gives the necessary formula for the electrocaloric effect. When we take “adiabatic” condition; that is, constant entropy dS = 0, and Equation (46)

$$\theta =\frac{pT}{c_p^E}E$$

Here, $c_p^E$ is specific heat capacity (per unit mass) under E = 0, and p is pyroelectric coefficient given by $p=-{\left(\frac{\partial P}{\partial T}\right)}_X$. The material’s development strategy should be higher p and lower $c_p^E$ 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, $\frac{pT}{c_p^E}$. Recall the following equations obtained in the ferroelectric phase (T < T_C = T₀) in Section 2.1.3.2:

$$\begin{gathered} G_1(P,T=(1/2)\alpha P^2 + (1/4)\beta P^4 \\ \alpha = (T - T_0)/{\epsilon }_{0}C\hspace{1em}\hspace{1em}[{\epsilon }_0{\epsilon }_0=8.854\times {10}^{-12} \mathrm{F}/\mathrm{m}] \\ {P}_S= \pm \sqrt{\frac{-\alpha }{\beta }}= \pm \sqrt{\left(T_0-T\right)/{\epsilon }_{0}C\beta } \\ {P_S}^2=\frac{-\alpha }{\beta }=\left(T_0-T\right)/{\epsilon }_{0}C\beta \\ {G}_1\left(P,T\right)=\frac{1}{2}\alpha {P_S}^2+\frac{1}{4}\beta {P_S}^4=\frac{-1}{4}\frac{{\alpha }^2}{\beta }=\frac{-1}{4\beta {\left({\epsilon }_{0}C\right)}^2}{\left(T_0-T\right)}^2 \end{gathered}$$

(49)

(a)

Specific Heat Capacity

$$c_p^E=-T{\left(\frac{{\partial }^{2}G}{\partial T^2}\right)}_{X,E}=\frac{T}{2\beta {\left({\epsilon }_{0}C\right)}^2}$$

(50)

(b)

Pyroelectric Coefficient

Noting here $P_S=\sqrt{\left(T_0-T\right)/{\epsilon }_{0}C\beta }$, the pyroelectric coefficient is obtained as

$$p=-{\left(\frac{\partial P}{\partial T}\right)}_X=\frac{1}{2\sqrt{{\epsilon }_{0}C\beta }}\frac{1}{\sqrt{\left(T_0-T\right)}}$$

(51)

(c)

Figure of Merit of Electrocaloric Effect

$$\frac{pT}{c_p^E}=\sqrt{\beta }{\left({\epsilon }_{0}C\right)}^{3/2}\frac{1}{\sqrt{\left(T_0-T\right)}}$$

(52)

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

$$E_C=\sqrt{\frac{-4}{27}{\alpha }^3/\beta }=\sqrt{\frac{4}{27}{\left(T_0-T\right)}^3/\beta {\left({\epsilon }_{0}C\right)}^3},$$

(53)

the E_C 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.

2.1.6.4. “Electrothermal Coupling Factor” $k^{ET}$

Let us obtain a practical number of “electrothermal coupling factor” $k^{ET}$ Equation (44) in PZTs at room temperature: ${k^{ET}}^2=\frac{p^2}{\left(\frac{c^E}{T}\right){\epsilon }_0{\epsilon }^X}$. Initially, all parameters are listed below:

• Pyroelectric coefficient p
  Pb(Zr_0.52Ti_0.48)O₃—$0.3 0.7\times {10}^4$ C/m²K [15]
  PZT-5A—$0.3 0.5\times {10}^4$ C/m²K [16]
  Pb(Zr_0.74Ti_0.36)O₃—$2.7\times {10}^4$ C/m²K [17]
• Specific heat capacity $c_p^E$
  PZT-APC841—340 J/KgK [9]
  PLZT(6/80/20)—338 J/KgK [18]

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, $c^E$ and $c^D$. 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 ${\epsilon }^X$

PZT-APC841—1375 [19]

Finally, the electrothermal coupling factor is calculated as follows:

$${k^{ET}}^2=\frac{p^2}{\left(\frac{c^E}{T}\right){\epsilon }_0{\epsilon }^X}=\frac{{(0.5\times {10}^{-4} \mathrm{C}/{\mathrm{m}}^2\mathrm{K})}^2}{\left(\frac{340 \mathrm{J}/\mathrm{KgK}}{300 K}\right)\left(8.854\times {10}^{-12} \frac{F}{m}\right)\times 1375}=0.18.$$

(54a)

$$k^{ET}=43\%$$

(54b)

We can conclude that the electrothermal coupling factor $k^{ET}$ in PZTs has a comparable order of magnitude to the electromechanical coupling factor k of 30–70%. It is interesting to calculate the $k^{ET}$ value for the second-order phase transition case: using Equations (49)–(51), $k^{ET}=100\%$, irrelevant to any temperature, very idealistically.

2.2. 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.

2.2.1. 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 $\tau \propto 1/\left(T_C-T\right)$ near a second-order phase transision point, T_C. The key assumption is that the change of polarization P with time (i.e., $\frac{\partial P}{\partial t}$) is proportional to internal energy decrease rate with the unit order parameter change (i.e., $\frac{-\partial U}{\partial P}$). By introducing “relaxation time” $\tau$, the inverse of which is taken as a proportional constant:

$$\frac{\partial P}{\partial t}=\frac{-1}{\tau }\left(\frac{\partial U}{\partial P}\right)$$

(55)

We use again the simple “second-order phase transition” model here for a practical calculation:

G₁(P,T) = (1/2)α P² + (1/4)β P⁴  [α = (T − T₀)/ε₀ C]

Recall the fact that, in the ferroelectric phase (T < T₀), the spontaneous polarization is expressed as

$${\mathrm{P}}_{\mathrm{S}} = \pm \sqrt{\frac{-\alpha }{\beta }}= \pm \sqrt{\left(T_0-T\right)/{\epsilon }_{0}C\beta }$$

Equation (55) is expressed in practice by

$$\frac{\partial P}{\partial t}=\frac{-1}{{\tau }^{\prime }}\frac{\partial G_1}{\partial P}=\frac{-1}{{\tau }^{\prime }}\left(\alpha P+\beta P^3\right)$$

(56)

We consider the initial polarization $P_0$ ($P_S+{\epsilon }_0{\epsilon }^X{E}_0$ by applying small electric field $E_0$ on a ferroelectric specimen with $P_S=\sqrt{\frac{-\alpha }{\beta }}$. The meaning of Equation (56) is to provide how quickly the P should return to the equilibrium state $P_S$ after sudden termination of the field. Equation (56) can be solved as

$$\begin{gathered} {{\displaystyle \int }}_0^t\left(\frac{-1}{{\tau }^{\prime }}\right)dt={{\displaystyle \int }}_{P_0}^P\frac{dP}{\left(\alpha P+\beta P^3\right)}=\frac{1}{\beta {P_S}^2}\underset{P_0}{\overset{P}{{\displaystyle \int }}}[-\frac{1}{P}+\frac{1}{2}\frac{1}{\left(P-P_S\right)}+\frac{1}{2}\frac{1}{\left(P+P_S\right)}]dP \\ \frac{-t}{{\tau }^{\prime }}=\frac{1}{-\alpha }\left[-ln\left(\frac{P}{P_0}\right)+\frac{1}{2}ln\left(\frac{P-P_S}{P_0-P_S}\right)+\frac{1}{2}ln\left(\frac{P+P_S}{P_0+P_S}\right)\right] \end{gathered}$$

Thus, we obtain

$$\sqrt{\left[1-{\left(\frac{P_S}{P}\right)}^2\right]/\left[1-{\left(\frac{P_S}{P_0}\right)}^2\right]}=e^{-t/\left(\frac{{\tau }^{\prime }}{-\alpha }\right)}$$

(57)

The actual “relaxation time constant” τ is provided by

$$\tau =\frac{{\tau }^{\prime }}{-\alpha }=\frac{{\tau }^{\prime }{\epsilon }_{0}C}{\left(T_0-T\right)}\propto 1/\left(T_0-T\right)\hspace{1em}[\alpha = (T-T_0)/{\epsilon }_0 C]$$

(58)

From Equations (57) and (58), we can conclude the following items:

• Time constant is proportional to “permittivity” ${\epsilon }_0{\epsilon }^X$.
• t → $\infty$ gives P → P_S.
• With approaching T → T₀, $\tau$ → $\infty$. The recovery time of P is very slow around T_C = T₀.

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.

2.2.2. 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:

$$\frac{\partial \theta }{\partial t}=\frac{-1}{\tau }\left(\frac{\partial U}{\partial \theta }\right)$$

(59)

If we take the Gibbs energy transformed from Equation (4) under stress-free condition (X = 0),

$$\begin{gathered} G\left(T,X,E\right)=\frac{1}{2}\left(\frac{{\partial }^{2}G}{\partial T^2}\right){\theta }^2+\frac{1}{2}\left(\frac{{\partial }^{2}G}{\partial E^2}\right)E^2+\left(\frac{{\partial }^{2}G}{\partial T\partial E}\right)\theta E \\ \hspace{1em}\hspace{1em}-\frac{1}{2}\left(\frac{c_p^E}{T}\right){\theta }^2-\frac{1}{2}\left({\epsilon }_0{\epsilon }^X\right)E^2+p\theta E \end{gathered}$$

(60)

Equation (59) is transformed into

$$\frac{\partial \theta }{\partial t}=\frac{1}{{\tau }^{\prime }}\left[\left(\frac{c_p^E}{T}\right)\theta -pE\right]$$

(61)

The meaning of Equation (61) is to provide how quickly the θ should approach the equilibrium status ${\theta }_s$ 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 (46) as

$${\theta }_s=\frac{pT}{c_p^E}E.$$

Equation (61) can be solved as

$$-{{\displaystyle \int }}_0^t\frac{1}{{\tau }^{\prime }}\left(\frac{c^E}{T}\right)dt={{\displaystyle \int }}_0^{\theta }\frac{d\theta }{\left({\theta }_s-\theta \right)}$$

(62a)

$$\frac{-1}{{\tau }^{\prime }}\left(\frac{c_p^E}{T}\right)t=ln\left({\theta }_s-\theta \right)-ln\left({\theta }_s\right)$$

(62b)

$$\left(1-\theta /{\theta }_s\right)=e^{-t/{\tau }^{\prime }\left(\frac{T}{c_p^E}\right)}$$

(63)

The actual relaxation time constant τ is expressed by

$$\tau ={\tau }^{\prime }\left(\frac{T}{c_p^E}\right)$$

(64)

From Equations (63) and (64), we can conclude the following items:

• In a uniform (no space gradient) specimen, temperature change follows an exponential trend with time: $\left(1-\theta /{\theta }_s\right)=e^{-t/\tau }$.
• t → $\infty$ gives $\theta$ → ${\theta }_s$. ${\theta }_s$ corresponds to the temperature change by the electrocaloric effect.
• Time constant is proportional to temperature, and inversely proportional to the specific heat capacity $c_p^E$. The larger the specific heat capacity $c_p^E$, and the lower the temperature, the lower the time constant $\tau$. 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 (64) is very relevant to the time constant in thermal diffusivity measurement discussed in Section 3.2, where $\tau \propto 1/{\alpha }_T=\frac{c_p\rho }{\lambda }$. Note the opposite effect of the specific heat capacity to the above relaxation time constant, and the heat flow relaxation time.

3. 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 $\left(1-\theta /{\theta }_s\right)=e^{-t/\tau }$ is not sustained correctly.

3.1. 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) ($\left(T-T_R\right)$ difference from the T_R) at a point a distance x from one end of rod. We introduce the following three parameters:

• q—Heat flux = quantity of heat passing through cross section of rod per unit area per unit time. The unit of heat quantity is calorie (cal), and 1 cal corresponds to 4.186 Joule in energy.
• λ—Thermal conductivity = thermal or heat conductance per unit length of material. The inverse 1/λ is the resistance which measures the temperature drop per unit length when heat flux is unity.
• C—Heat capacitance of material = specific heat capacity c_p ($c_p^E$, $c_p^D$) $\times$ density ρ, number of heat units to raise block of unit area and unit length (i.e., per unit volume) 1 °C in temperature. Note that the specific heat capacity [J/(kg·K)] is determined per unit mass (kg), which must be multiplied by density (kg/m³) to obtain the heat capacitance per unit volume [J/(m³·K)].

First, if $\theta \left(x+dx,t\right)$ is the temperature at a point at a distance of $\left(x+dx\right)$ from the end of the rod, we have (considering the unit cross section area tube)

$$\left(\frac{1}{\lambda }\right)qdx=\theta \left(x,t\right)-\theta \left(x+dx,t\right)=\theta \left(x,t\right)-\left[\theta \left(x,t\right)+\left(\frac{\partial \theta }{\partial x}\right)dx\right]$$

(65a)

Thus, we obtain

$$-\left(\frac{\partial \theta }{\partial x}\right)=\left(\frac{1}{\lambda }\right)q$$

(65b)

Next, the heat flux $q\left(x,t\right)$ at the point x raises the temperature of the lamina of thickness dx at the rate $\left(\frac{\partial \theta }{\partial t}\right)$, and heat flux of magnitude $q\left(x+dx,t\right)$ emerges from the lamina of width dx. We thus have the equation

$$q\left(x,t\right)=Cdx\left(\frac{\partial \theta }{\partial t}\right)+q\left(x+dx,t\right)$$

(66a)

Again expanding $q\left(x+dx,t\right)$ by Taylor’s expansion, we obtain

$$-\left(\frac{\partial q}{\partial x}\right)=C\left(\frac{\partial \theta }{\partial t}\right)$$

(66b)

From Equations (65b) and (66b) to eliminate the parameter q, we obtain 1D heat transfer equation as

$$\left(\frac{\partial \theta }{\partial t}\right)=\frac{\lambda }{C}\left(\frac{{\partial }^2\theta }{\partial x^2}\right)=\frac{\lambda }{c_p\rho }\left(\frac{{\partial }^2\theta }{\partial x^2}\right)$$

(67)

The proportional constant $\frac{\lambda }{c_p\rho }$ is called “thermal diffusivity” ${\alpha }_T$, which is practically measured by thermal experiments. Units and typical numbers of the necessary parameters in PZT are summarized:

• Thermal conductivity $\lambda$—$0.6 1.3$ [W/(m·K)]
• Specific heat capacity $c_p$—$340 420$ [J/(kg·K)]
• Mass density $\rho$—7600 [kg/m³]
• Thermal diffusivity ${\alpha }_T$—$2.7 5.1\times {10}^{-7}$ [m²/s]

3.2. Solution of 1D Heat Transfer Equation

Let us solve the heat equation Equation (67) using the Laplace operator:

$$\left(\frac{\partial \theta }{\partial t}\right)={\alpha }_T\left(\frac{{\partial }^2\theta }{\partial x^2}\right)$$

(68)

Taking $\Theta$ as the Laplace transform of $\theta \left(x,t\right)$ ($\Theta =L\left[\theta \left(x,t\right)\right]$), we have

$$\left(\frac{{\partial }^2\Theta }{\partial x^2}\right)=\frac{s}{{\alpha }_T}\Theta$$

(69)

We may assume the solution form as

$$\Theta =A{e}^{\sqrt{\frac{s}{{\alpha }_T}}x}+B{e}^{-\sqrt{\frac{s}{{\alpha }_T}}x}$$

(70)

Knowing the inverse Laplace transform,²²⁾

$$L^{-1}\left[e^{-a\sqrt{s}}\right]=\frac{a}{2\sqrt{\pi t^3}}{e}^{-a^2/4t}$$

(71)

and $a=\mp \frac{x}{\sqrt{{\alpha }_T}}$, Equation (70) can be transformed into

$$\theta \left(x,t\right)=-A\frac{\frac{x}{\sqrt{{\alpha }_T}}}{2\sqrt{\pi t^3}}{e}^{-x^2/4t{\alpha }_T}+B\frac{\frac{x}{\sqrt{{\alpha }_T}}}{2\sqrt{\pi t^3}}{e}^{-x^2/4t{\alpha }_T}$$

(72a)

We consider the boundary condition: Large heat source ($\theta$ = $\theta$₀) 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:

$$\theta \left(x,t\right)={\theta }_0{e}^{-x^2/4t{\alpha }_T}$$

(72b)

Equation (72b) is a sort of “Error function” in terms of coordinate x, and indicates that

• At x = 0, $\theta \left(x=0,t\right)={\theta }_0$ —Temperature is constant at ${\theta }_0$, irrelevant to time after attaching to heat source.
• At x = L, $\theta \left(L,t\right)={\theta }_0{e}^{-L^2/4t{\alpha }_T}$ —By measuring temperature change with t at x = L, we can obtain the thermal diffusivity ${\alpha }_T$.

Figure 7a shows temperature distribution profiles as a function of position x for various time lapses, and Figure 7b shows the temperature change with time at the position x = L. This simulation used the following values: ${\theta }_0=15 °\mathrm{C}$, L = $2\times {10}^{-3}\left[\mathrm{m}\right],{\alpha }_T$ = $8\times {10}^{-7}$ [m²/s]. As can been seen, heat flux penetration speed for 2 mm thick sample does not take longer than 0.5 s, and the temperature monitoring finishes in less than 3 s. Even for a thicker specimen with 10 mm, the measuring time is still shorter than 30–40 s.

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]:

$$\theta \left(L,t\right)={\theta }_0\left[1-e^{-t/\left(\frac{L^2}{{\alpha }_T{\left(\pi /2\right)}^2}\right)}\right]$$

(73a)

We verify this handling equivalency with the exact formula,

$$\theta \left(L,t\right)={\theta }_0{e}^{-L^2/4t{\alpha }_T}$$

(73b)

By putting $A=e^{-\left(\frac{L^2}{4{\alpha }_T}\right)/t}$, we obtain $\frac{1}{lnA}=\frac{-t}{\left(\frac{L^2}{4{\alpha }_T}\right)}$, then $e^{-t/\left(\frac{L^2}{4{\alpha }_T}\right)}=e^{\frac{1}{lnA}}$. Taking Taylor expansion,

$$e^{\frac{1}{lnA}}\approx 1-A+\cdots$$

we verify the relation

$$e^{-L^2/4t{\alpha }_T}\approx \left[1-e^{-t/\left(\frac{L^2}{4{\alpha }_T}\right)}\right]$$

(74)

However, note the numerical factor difference; that is, ${\left(\pi /2\right)}^2$ in Equation (73a), and 4 in Equation (74). This means that the experimentally determined ${\alpha }_T$ with Equation (73a) was overestimated from the ${\alpha }_T$ value with Equation (73b) or (74) by the factor of $4/{\left(\pi /2\right)}^2$. The absolute ${\alpha }_T$ values may need to be modified by this factor from the previously reported values in Ref. [9], if required.

3.3. 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 (74). 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 $c_p^E$ (short-circuit condition) and $c_p^D$(open-circuit condition). The accuracy of our setup has already been demonstrated in several materials with low to medium thermal diffusivity (0.1–3 $\times {10}^{-6}$ m²/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 ${\alpha }_T$ 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 d₃₃ 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 c_p = 340 J/kgK and ρ = 7600 kg/ m³. The specific heat capacity difference between $c_p^E$ and $c_p^D$ 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. Thermal diffusivity, specific heat capacity, and thermal conductivity of a Hard PZT, for unpoled and poled under different electrical boundary conditions.

Hard PZT	Thermal Diffusivity ${\mathit{\alpha }}_{\mathit{T}}$ (10−7 m2/s)	+/−	Specific Heat Capacity cp (J/kg K)	Thermal Conductivity $\mathit{\lambda }$ (W/m K)	+/−
Unpoled	4.32	0.34	$c_p^D$ = 279	0.91	0.10
Open Circuit	5.02	0.23	$\frac{1}{3}\left(c_p^D+2c_p^E\right)$ = 320	1.22	0.06
Short Circuit	8.25	0.78	$c_p^E$ = 340	2.13	0.23

summarizes thermal diffusivity, specific heat capacity, and thermal conductivity of a Hard PZT (APC 841), for unpoled and poled under different electrical boundary conditions (open- and 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.

3.4. 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

$${\alpha }_T=\frac{\lambda }{c_p\rho }$$

(75)

where $\lambda$, c_p, ρ are thermal conductivity, specific heat capacity, and mass density, respectively. The difference in ${\alpha }_T$ must be the combination effect of the $\lambda$ and c_p differences (the density ρ may not be changed by the electric constraint condition).

3.4.1. Specific Heat Capacity—Scalar Parameter

Recall the two types of specific heat capacity, $c_p^E$ and $c_p^D$, which are correlated by

$$c_p^D=c_p^E\left(1-{k^{ET}}^2\right)$$

Here, the primary “electrothermal coupling factor” $k^{ET}$ is given by ${k^{ET}}^2=\frac{p^2}{\left(\frac{c_p^E}{T}\right){\epsilon }_0{\epsilon }^X}$.

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 $E_{dep}=-\left(\frac{\Delta P}{{\epsilon }_0{\epsilon }^X}\right)$ 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, $c_p^E$ 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, $c_p^D$ 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, $c_p^D$ may approach $c_p^E$ (slightly larger). The complex situation is found in a specimen of poled under open-circuit condition. Because of the surface electrode, E-constant two-dimensionally 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 $\frac{1}{3}\left(c_p^D+2c_p^E\right)$ because of scalar characteristics. We decided first $c_p^E=340J/kg·K$ from the DSC equipment. Then, using ${k^{ET}}^2=$, $c_p^D=c_p^E\left(1-{k^{ET}}^2\right)=279J/kg·K$. Finally, $\frac{1}{3}\left(c_p^D+2c_p^E\right)=320J/kg·K$ 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 k_t or k₃₃ piezoelectric resonance modes, where the thickness or length stress distribution generates the depolarization field $E_{dep}=-\left(\frac{\Delta P}{{\epsilon }_0{\epsilon }^X}\right)$ anti-parallel to the induced polarization (∆P = d₃₃X₃) in order to satisfy the D—constant (constraint) condition. The elastic compliance $s_{33}^D$ and sound propagation speed ($v_{33}^D=1/\sqrt{\rho s_{33}^D}$ are modified by the electromechanical coupling factor ${k_t}^2$ or ${k_{33}}^2$ as $s_{33}^D=s_{33}^E\left(1-{k_t}^2\right)$. 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 ${k^{ET}}^2$ as $c^D=c^E\left(1-{k^{ET}}^2\right)$ depending on the electrical constraint. We may also imagine the heat-flow (i.e., thermal conductivity) modulation according to the electrical constraint condition.

3.4.2. Thermal Conductivity—Tensor Parameter

Using ${\alpha }_T=\frac{\lambda }{c_p\rho }$ with $\rho =7600\mathrm{kg}/{\mathrm{m}}^3$, we calculated thermal conductivity λ as listed in the last column of Table 1. Note first that the thermal conductivity ${\lambda }_{ij}$ is a second-rank tensor (similar to the electrical conductivity ${\sigma }_{ij}$), different from the scalar “specific heat capacity”.

Let us discuss the difference among the measured thermal conductivity values for poled short-circuit, poled open-circuit, and unpoled specimens. Note that, once poled, the piezo-ceramic crystallographic symmetry becomes $\infty mm$ (equivalent symmetry to $6mm$ in terms of tensor notations). First, the open-circuit thermal conductivity λ^D₃₃ is 0.57 times smaller than the short-circuit one λ^E₃₃. 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 q₃ and x-component of the temperature gradient grad (θ) (i.e., vector) along the spontaneous polarization P_S direction. Second, the unpoled specimen exhibits a crystallographic symmetry of isotropic cubic $m3m$ (different from the poled specimens), and the thermal conductivity λ^u₁₁ 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^λ₃₃ below, in the relationship between the open-circuit λ^D₃₃ and short-circuit thermal conductivity λ^E₃₃ such as:

$${\lambda }_{33}^D={\lambda }_{33}^E\left(1-{k_{33}^{\lambda }}^2\right)$$

(76)

Let us consider a derivation of this “secondary coupling coefficient k^λ₃₃” with using a “space-gradient model” illustrated in Figure 9. Analogous to Equation (39a,b),

$$\begin{gathered} S=\frac{c_p}{T}\theta -pE \\ D=-p\theta +{\epsilon }_0{\epsilon }^{X}E \end{gathered}$$

We define “heat flux” per unit time per unit area $q=-T\nabla \left(S\right)$ [W/m²] and $\nabla \left(\theta \right)=\frac{\partial \theta }{\partial x}$ [K/m], and gradient of vector $\nabla \left(E\right)$ [V/m²]. Nabra $\nabla$ is the gradient, $\nabla =i\frac{d}{dx}+j\frac{d}{dy}+k\frac{d}{dz}$. In our simple model, we focus only on the x-direction component. Knowing the facts in 1D model: (1) $q_x=-\lambda \nabla \left(\theta \right)$, and (2) heat flux seems to be controlled by the electric field E (because of the difference between ${\lambda }_{33}^D$ and ${\lambda }_{33}^E$), we assume the following constitutive equations, resembling Equation (39a,b) above, but taking the gradient with respect to the heat flux direction x:

$$\begin{gathered} -T\nabla \left(S\right)=\lambda \left(-\nabla \left(\theta \right)\right)-T{E}_x\nabla \left(p^{\prime }\right)-T{p}^{\prime }\left(-\nabla \left(E_x\right)\right) \\ \left(-\nabla \left(D_x\right)\right)=\theta \nabla \left(p^{\prime }\right)-p^{\prime }\left(-\nabla \left(\theta \right)\right)-E_x\nabla \left({\epsilon }_0{\epsilon }^X\right)+{\epsilon }_0{\epsilon }^X\left(-\nabla \left(E_x\right)\right) \end{gathered}$$

If we assume that the room temperature is reasonably low from the Curie temperature, enough to satisfy $\nabla \left(p^{\prime }\right)$, $\nabla \left({\epsilon }_0{\epsilon }^X\right)\ll p^{\prime },{\epsilon }_0{\epsilon }^X$ (that is, the pyroelectric coefficient and permittivity change with temperature (or space coordinate x) is reasonably small), we can approximate the following equations:

$$-\nabla \left(S\right)=\left(\lambda /T\right)\left(-\nabla \left(\theta \right)\right)-p^{\prime }\left(-\nabla \left(E_x\right)\right)$$

(77a)

$$-\nabla \left(D_x\right)=-p^{\prime }\left(-\nabla \left(\theta \right)\right)+{\epsilon }_0{\epsilon }^X\left(-\nabla \left(E_x\right)\right)$$

(77b)

We adopt another assumption, that is, D-constant along the thermal flow. From $dD=\left(\frac{\partial D}{\partial x}\right)dx+\left(\frac{\partial D}{\partial t}\right)dt=0$, we obtain $-\nabla \left(D_x\right)=\left(\frac{\partial D}{\partial t}\right)/\left(\frac{\partial x}{\partial t}\right)=\frac{1}{v}\left(\frac{\partial D}{\partial t}\right)$, where $v=\left(\frac{\partial x}{\partial t}\right)$, drift velocity of electric-displacement-charge. We can transform the gradient equations into time dependence:

$$-\nabla \left(S\right)=\left(\lambda /T\right)\left(-\nabla \left(\theta \right)\right)-p^{\prime }\left(\frac{1}{v}\dot{E_x}\right)$$

(78a)

$$\left(\frac{1}{v}\dot{D_x}\right)=-p^{\prime }\left(-\nabla \left(\theta \right)\right)+{\epsilon }_0{\epsilon }^X\left(\frac{1}{v}\dot{E_x}\right)$$

(78b)

Since $\left(\frac{\partial D}{\partial t}\right)$ means the electric-displacement-current, Equation (78b) stands for the dynamic (transient) current density and electric field relation, modified by additional current originating from the temperature gradient via pyroelectric effect. The constitutive equations are valid theoretically only under the equilibrium conditions. However, even a transient and irreversible status, if the phenomenon is in a pseudo-steady state like electric current or heat flow, we may assume similar constitutive equations: Equation (78a,b) are with respect to $-\nabla \left(S\right)$, $\left(\frac{1}{v}\dot{D_x}\right)$, $-\nabla \left(\theta \right)$, and $\left(\frac{1}{v}\dot{E_x}\right)$, instead of S, D, θ, and E in Equation (39a,b). Recall that $\theta$ and $E_x$ are intensive (irrelevant to volume) and $S$ and $D_x$ are extensive (size dependent) parameters. Let us confirm the units of all equation parameters: since the unit of entropy S is [J/K], we may assume the unit of $\nabla \left(S\right)$ as flow per second per unit area (normalize by size) [W/K·m²]. Because T and $\nabla \left(\theta \right)$ have [K] and [K/m], respectively, $\lambda$ should possess the unit of [W/K·m] (normalized by time). While $\left(\frac{1}{v}\dot{E_x}\right)$ has [V/m²], $p^{\prime }$ should have [Cm/K·s]. On the contrary, Equation (78b) requires a tricky explanation: we assume $\left(\frac{1}{v}\dot{D_x}\right)$ is total charge flow (multiplied by cross-section area) per second with the unit [C/s]. In this case, taking $\nabla \left(\theta \right)$ as [K/m], $p^{\prime }$ should have [Cm/K·s]. The remaining problem is the last term: ${\epsilon }_0{\epsilon }^X$ has [C/Vm] and $\left(\frac{1}{v}\dot{E_x}\right)$ has [V/m²], thus the product should have the unit of [C/m³]. We need to translate this term by multiplying the volume, such as ${\epsilon }_0{\epsilon }^X\dot{E_x}\to C\dot{V}$ (here C: capacitance, V: voltage). The new coupling parameter p′ seems to be analogous to pyroelectric coefficient p [C/K·m²] but is the total charge flow (multiplied by volume) per second because we are discussing basically space-gradient temperature profile. From Equation (78a,b), we may define a new “electrothermal coupling coefficient” k^λ₃₃ as

$${k_{33}^{\lambda }}^2=\frac{{\left(Couplingfactor\right)}^2}{\left(Productofthediagonalparameters\right)}=\frac{T{p\prime }^2}{\lambda {\epsilon }_0{\epsilon }^X}$$

(79)

though Equation (77a,b) are not in an equilibrium state, but in an irreversible state. When we consider D₃ = constant (i.e., open-circuit), we get $\left(-\nabla \left(E_x\right)\right)=\frac{p^{\prime }}{{\epsilon }_0{\epsilon }^X}\left(-\nabla \left(\theta \right)\right)$ from Equation (77b), then Equation (77a) is transformed into

$$q_x={\lambda }_{33}^E\left(1-{k_{33}^{\lambda }}^2\right)\left(-\nabla \left(\theta \right)\right)$$

(80)

The reader can understand that Equation (80) provides the support for the Equation (76) formula. The experimental result, λ^D₃₃ = 0.57 λ^E₃₃, gives ${k_{33}^{\lambda }}^2=\frac{Tp{\prime }^2}{{\lambda }_{33}^E{\epsilon }_0{\epsilon }^X}=0.43$. Knowing ${\lambda }_{33}^E=2.13W/K·m$ and ${\epsilon }^X=1375$, T = 300 K, we can obtain $p^{\prime }=6.1\times {10}^{-6}\left[\mathrm{Cm}/\mathrm{s}·\mathrm{K}\right]$. The remaining problem is how to integrate the volume and unit time calibration into the ${\epsilon }_0{\epsilon }^X$ 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.

3.5. 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 k₃₁ type plate specimens.

Table 2. Thermal properties of NKN-Cu in comparison with Hard PZT.

Thermal Properties	cp (J/kg K)	λ (W/m/K)
Hard-PZT	420	1.25
NKN-Cu	580	3.10

summarized the preliminary data on the thermal conductivity in the NKN-based material with the PZT’s value (both unpoled samples). Pb-free piezoceramics such as (Na,K)NbO₃− and (Bi,Na)TiO₃− 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.

4. Nonlinear Elastic Performances and “Phonons”

We discussed the phenomenological discussions in ferroelectrics on (1) linear equilibrium in Section 2.1, (2) nonlinear equilibrium in Section 2.1.3, and (3) Harmonic dynamic in Section 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”.

4.1. 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:

$$U=\frac{-M{q}^2}{r}+\frac{Nb}{r^n}$$

(81)

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 ($r_0^{n-1}=nNb/M{q}^2$), we obtain the form as a function of $\Delta r$ (=r − r₀):

$$\Delta U=U\left(r\right)-U\left(r_0\right)=f{\left(\Delta r\right)}^2-g{\left(\Delta r\right)}^3$$

(82)

where $f=\left(n-1\right)M{q}^2/2r_0^3$ and $g=\left(n+4\right)\left(n-1\right)M{q}^2/6r_0^4$. Using the Boltzmann distribution for Gibbs energy $\Delta V_{\pm }=\Delta U\pm qE\Delta r$, the average equilibrium separation at high temperatures under an applied electric field E is as follows:

$$\begin{gathered} \Delta r_{\pm }> {{\displaystyle \int }}_{-\infty }^{\infty }\Delta rexp\left(\frac{-\Delta V_{\pm }}{kT}\right)d\Delta r/{{\displaystyle \int }}_{-\infty }^{\infty }exp\left(\frac{-\Delta V_{\pm }}{kT}\right)d\Delta r \\ \frac{3gk}{4f^2}T\pm \frac{q}{2f}E+\frac{3g{q}^2}{4f^3}{E}^2, \end{gathered}$$

(83)

where subscripts $\pm$ denote the ion pairs in the positive and negative directions of the electric field respectively. The strain is therefore given by:

$$\Delta r_{+>+\Delta r_{->\frac{}{2r_0}=\frac{3gk}{4f^2{r}_0}T+\frac{3g{q}^2}{4f^3{r}_0}{E}^2}}$$

(84)

where the first term represents thermal expansion and the second term, electrostriction. The reader can understand that both coefficients are originated from “$g$”, which is an anharmonic term, and when the crystal is purely harmonic, neither thermal expansion nor electrostriction occur. Knowing the polarization P, given by,

$$P=q\langle \Delta r\rangle \frac{}{2r_0^3}=\frac{q^2}{4f{r}_0^3}E$$

(85)

the thermal expansion coefficient ${\alpha }_L$ and electrostrictive Q coefficient can be calculated

$${\alpha }_L=\frac{3gk}{4f^2{r}_0}=\frac{\left(n+4\right)k}{2n\left(n-1\right)Nb}{r}_0^n$$

(86)

$$Q=\left(\frac{3g{q}^2}{4f^3{r}_0}\right)/{\left(\frac{q^2}{4f{r}_0^3}\right)}^2=\frac{4\left(n+4\right)M}{nNb}{r}_0^{n+3}$$

(87)

Finally, using the “First Law” of thermodynamic of $dU=-pdV$ (p: hydrostatic pressure, V: volume), the isothermal compressibility ${\chi }_T$ becomes:

$$\frac{1}{{\chi }_T}=-V{\left(\frac{\partial p}{\partial V}\right)}_T=V{\left(\frac{{\partial }^{2}U}{\partial V^2}\right)}_T$$

(88)

Changing the variable from V to r and developing the result leads to the equation

$${\chi }_T=\frac{18r_0^4}{\left(n-1\right)M{q}^2}=\frac{18}{n\left(n-1\right)Nb}{r}_0^{n+3}$$

(89)

Comparing Equations (87) and (89), we obtain the relationship, ${\chi }_T\propto Q$. The compressibility and electrostrictive coefficient are induced strain originated from the same crystal lattice anharmonic potential against the force and electric field.

4.2. 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].

4.2.1. 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 $u_n$, we can obtain the following set of Newton equations:

$$\begin{gathered} M{\ddot{u}}_1=-K\left(u_1-u_2\right) \\ M{\ddot{u}}_n=-K\left(u_n-u_{n-1}\right)-K\left(u_n-u_{n+1}\right) \left(\mathrm{n} = 2,3,\cdots ,\mathrm{N}-1\right) \\ M{\ddot{u}}_N=-K\left(u_N-u_{N-1}\right) \end{gathered}$$

(90a)

Assuming a harmonic vibration, $u_n\propto e^{j\omega t}$, where ω is the frequency,

$$\begin{gathered} -M{\omega }^2{u}_1=-K\left(u_1-u_2\right) \\ -M{\omega }^2{u}_n=-K\left(u_n-u_{n-1}\right)-K\left(u_n-u_{n+1}\right) \left(\mathrm{n} = 2,3,\cdots ,\mathrm{N}-1\right) \\ -M{\omega }^2{u}_N=-K\left(u_N-u_{N-1}\right) \end{gathered}$$

(90b)

Thus, we can write the following matrix formula:

$${\omega }^{2}u=\left[F\right]u$$

(91)

where u is an N-Dimensional vector with the components of $u_n$, and $\left[F\right]$ is a symmetric matrix with the components:

$$F_{nm}=\left(\frac{K}{M}\right)\left\{\left(2-{\delta }_{n1}-{\delta }_{nN}\right){\delta }_{nm}-{\delta }_{n,m-1}-{\delta }_{n,m+1}\right\}$$

(92)

Here, we used Kronecker ${\delta }_{nm}$, which means that ${\delta }_{nm}=1$ for $n\ne m$; ${\delta }_{nm}=0$ for $n=m$.

We will now solve the “eigenvalue equation” of Equation (91). Knowing the matrix component relation

$$\left(\frac{K}{M}\right)\left[\left(\frac{M{\omega }^2}{K}-2\right)u_n+u_{n-1}+u_{n+1}\right]=0$$

(93)

and the boundary condition for the “free-end” ($u_0-u_1=u_N-u_{N+1}=0$), assuming the a general solution of “difference equation” Equation (93) in the form of $u_n\propto P^n$, we obtain the following characteristic equation:

$$P^2+\left(\frac{M{\omega }^2}{K}-2\right)P+1=0$$

(94)

When $\left|\frac{M{\omega }^2}{K}-2\right|<2$, we can put

$$\left(\frac{M{\omega }^2}{K}-2\right)\equiv -2coska\hspace{1em}\hspace{1em}(0 < \mathit{ka} < \mathsf{\pi }),$$

(95)

where k, a are the wave number and the lattice unit size. The solution of Equation (94) is expressed as

$$P=exp\left(\pm jka\right)$$

(96)

Then, we can express a general solution of Equation (93) as

$$u_n=Aexp\left(jkan\right)+Bexp\left(-jkan\right)$$

(97)

The “free-end” boundary condition gives the relations:

$$\{\begin{array}{c}A+B=Aexp\left(jka\right)+Bexp\left(-jka\right)\\ Aexp\left(jka\left(N+1\right)\right)+Bexp\left(-jka\left(N+1\right)\right)=Aexp\left(jkaN\right)+Bexp\left(-jkaN\right)\end{array}$$

(98)

To realize the meaningful solution (or A and B values), at least the following “determinant = 0” should be satisfied:

$$\left|\begin{array}{cc}1-exp\left(jka\right)& 1-exp\left(-jka\right)\\ exp\left(jkaN\right)\left[1-exp\left(jka\right)\right]& exp\left(-jkaN\right)\left[1-exp\left(-jka\right)\right]\end{array}\right|=0$$

(99)

This requirement determines the solutions regarding the k value

$$ka=\left(\frac{\alpha }{N}\right)\pi (\mathsf{\alpha } = 1, 2, \cdots , \mathrm{N}-1)$$

(100)

and for $u_n$,

$$\begin{gathered} u_n\propto \left[1-exp\left(-jka\right)\right]exp\left(jkan\right)-\left[1-exp\left(jka\right)\right]exp\left(-jkan\right) \\ 4jsin\left(\frac{ka}{2}\right)coska\left(n-\frac{1}{2}\right) \end{gathered}$$

(101)

Taking the normalizing factor determined by ${\displaystyle \sum }_{n=1}^N{u_n}^2=1$, we obtain the eigenvectors and eigen frequencies as follows:

$$u_n\left(\alpha \right)=\sqrt{\frac{2}{N}}cos\frac{\alpha \pi }{N}\left(n-\frac{1}{2}\right)\hspace{1em}\hspace{1em}(\mathsf{\alpha } = 1, 2, \cdots , \mathrm{N}-1)$$

(102)

$${\omega }_{\alpha }=2\sqrt{\frac{K}{M}}sin\frac{\alpha \pi }{2N}=2\sqrt{\frac{K}{M}}sin\frac{k_{\alpha }a}{2}$$

(103)

When $\left|\frac{M{\omega }^2}{K}-2\right|=2$, unique solution ${\omega }_0=0$, $u_n\left(0\right)=\sqrt{\frac{1}{N}}$ exists. However, when $\left|\frac{M{\omega }^2}{K}-2\right|>2$, no solution exists to satisfy the “free-end” boundary condition.

$k_{\alpha }a\equiv \left(\frac{\alpha }{N}\right)\pi$ is called “wave number”, and the relationship between the frequency and wave number (Equation (103)) is called “dispersion relation”. As shown in Figure 11a, $\omega \propto k$ 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 “$\omega \left(k\right)vs.k$” 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²³ per mol) of atoms of N, the energy gap among phonons is quantum-mechanically small. According to the Gibbs statistical dynamics for the equilibrium status, thermo-dynamical 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 $\hslash$ as

$$U=\frac{1}{2N}{\displaystyle \sum }_{\alpha =1}^N\hslash {\omega }_{\alpha } \cot \mathrm{h}\left(\frac{\hslash {\omega }_{\alpha }}{2k_{B}T}\right)$$

(104)

Here, $\hslash$ is the Planck constant, $k_B$ the Boltzmann constant, and T is temperature. The reader remembers in your high school physics that Einstein introduced “$\hslash \omega$” quantum energy unit for “photon” introduction. Equation (104) comes from the same philosophical concept.

We usually take N →$\infty$ in an ideal crystal, and ${\omega }_{\alpha }$ exists densely (N modes) just $0\le {\omega }_{\alpha }\le {\omega }_L=2\sqrt{K/M}$. In the case of large crystal (N → $\infty$), 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 (92) is modified as

$$F_{nm}=\left(\frac{K}{M}\right)\left(2{\delta }_{nm}-{\delta }_{n,m-1}-{\delta }_{n,m+1}-{\delta }_{n1}{\delta }_{mN}-{\delta }_{nN}{\delta }_{m1}\right)$$

(105)

We solve the above difference equation under the boundary condition $u_0-u_N=u_{N+1}-u_1=0$.

The results are slightly different from Equations (102) and (103):

$$u_n\left(\alpha \right)=\{\begin{array}{c}\sqrt{\frac{2-{\delta }_{\alpha ,0}-{\delta }_{2\left[\frac{\alpha }{2}\right],N}}{N}}cos\frac{2\pi n}{N}\left[\frac{\alpha }{2}\right]\left(\alpha =evennumber\right)\\ \sqrt{\frac{2}{N}}sin\frac{2\pi n}{N}\left[\frac{\alpha }{2}\right]\left(\alpha =oddnumber\right)\end{array}$$

(106)

$${\omega }_{\alpha }=2\sqrt{\frac{K}{M}}sin\frac{\pi }{N}\left[\frac{\alpha }{2}\right]\hspace{1em}\hspace{1em}(\mathsf{\alpha } = 0, 2, 3, \cdots , \mathrm{N})$$

(107)

Here, $\left[\frac{\alpha }{2}\right]$ is the Gauss symbol, which means the maximum integer number $\le \alpha /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 $k_{\alpha }a\equiv \left(\frac{2\pi }{N}\right)\left[\frac{\alpha }{2}\right]$ exhibits double-degenerated eigenvalues, which correspond to the following two eigen functions:

$$u_n\left(t\right)=Acos\left(k_{\alpha }an\pm {\omega }_{\alpha }t\right)$$

(108)

These waves correspond to left-ward and right-ward progressive waves. Taking an atomic distance “a”, the wavelength $\Lambda$ and wave velocity v are given by

$$\Lambda =2\pi a/k_{\alpha }$$

(109)

$$v=\Lambda \left(\frac{{\omega }_{\alpha }}{2\pi }\right)=2\sqrt{\frac{K}{M}}\frac{sin\left(k_{\alpha }/2\right)}{k_{\alpha }}a$$

(110)

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 ($k_{\alpha }a\ll 1$), the wave velocity approaches $v_0=2\sqrt{\frac{K}{M}}a$ and $\omega \propto k$ is valid.

4.2.2. 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 $u_n$, the dynamic equations are expressed as

$$\{\begin{array}{c}-m{\omega }^2{u}_{2n}=-K\left(u_{2n}-u_{2n+1}\right)-K\left(u_{2n}-u_{2n-1}\right)\\ -M{\omega }^2{u}_{2n+1}=-K\left(u_{2n+1}-u_{2n}\right)-K\left(u_{2n+1}-u_{2n+2}\right)\end{array}$$

(111)

Mass m (or M) is located at an even (or odd) number lattice position, and assuming the solutions

$$\{\begin{array}{c}\sqrt{m}{u}_{2n}=Q{e}_{1}exp\left(jkan\right)\\ \sqrt{M}{u}_{2n+1}=Q{e}_{2}exp\left(jkan\right)\end{array}$$

(112)

Equation (111) can be calculated as

$$\{\begin{array}{c}{\omega }^2{e}_1=\frac{2K}{m}{e}_1-\frac{K}{\sqrt{Mm}}\left[1+exp\left(-jka\right)\right]e_2\\ {\omega }^2{e}_2=\frac{-K}{\sqrt{Mm}}\left[1+exp\left(jka\right)\right]e_1+\frac{2K}{M}{e}_2\end{array}$$

(113)

In the matrix notation,

$${\omega }^{2}e=D\left(k\right)e$$

(114)

where

$$e\equiv \left(\begin{array}{c}{e}_1\\ e_2\end{array}\right)\left({e_1}^2+{e_2}^2=1\right)$$

(115)

$$D\left(k\right)=\left(\begin{array}{cc}\frac{2K}{m}& \frac{-K}{\sqrt{Mm}}\left[1+exp\left(-jka\right)\right]\\ \frac{-K}{\sqrt{Mm}}\left[1+exp\left(jka\right)\right]& \frac{2K}{M}\end{array}\right)$$

(116)

The lattice chain has a unit cell of size 2a with M and m atoms, and $e_1$ and $e_2$ are the “polarization (displacement) vectors” of m (even number position) and M (odd number position), respectively. From Equation (114), the eigenvalue ω should satisfy

$${\omega }^4-2K\left(\frac{1}{M}+\frac{1}{m}\right){\omega }^2+\frac{4K^2}{Mm}si{n}^2\left(\frac{ka}{2}\right)=0$$

(117a)

The solution of the above equation is given by

$${{\omega }_{\pm }}^2\left(k\right)=K\{\frac{1}{m}+\frac{1}{M}\pm \sqrt{{\left(\frac{1}{m}+\frac{1}{M}\right)}^2-\frac{4}{Mm}si{n}^2\left(\frac{ka}{2}\right)}\}$$

(117b)

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 → $\infty$, k will vary from 0 to 2π. At $k=0$,

$${{\omega }_{\pm }}^2\left(0\right)=\{\begin{array}{c}0\\ 2K\left(\frac{1}{m}+\frac{1}{M}\right)\end{array},$$

(118a)

while, at $ka=\pi$, from ${{\omega }_{\pm }}^2\left(0\right)=K\left(\frac{1}{m}+\frac{1}{M}\right)\pm \left(\frac{1}{m}-\frac{1}{M}\right)$ [supposing M > m]

$${{\omega }_{\pm }}^2\left(0\right)=\{\begin{array}{c}2K/M\\ 2K/m\end{array}$$

(118b)

Figure 11b illustrates the dispersion curves, plotted in the ka range of $-\pi \pi$ (because of the period symmetry). You can initially find frequency gap between $\sqrt{2K/M}$ and $\sqrt{2K/m}$, 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 $ka=\pi$. ${\omega }_{+}^2\left(k\right)$ corresponds to the upper curve, which is called “optical branch”. This branch approaches to ${\omega }_L$ with $k\to 0$:

$${\omega }_L=\sqrt{2K\left(\frac{1}{m}+\frac{1}{M}\right)}$$

(119)

Around $k=0$, the “polarization vector” e satisfies

$$\sqrt{m}{e}_1+\sqrt{M}{e}_2=0$$

(120)

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 Cl^- ions 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 ${\omega }_{-}^2\left(k\right)$ of in Figure 11b is called “acoustic branch”. This branch approaches $\omega \propto k$ with $k\to 0$; a long-wavelength vibration transfer with a constant wave speed. The polarization vector e satisfies the following condition at k = 0:

$$\sqrt{m}{e}_1=\sqrt{M}{e}_2$$

(121)

Two species of atoms displace in the same direction with the same amplitude.

4.2.3. Phonon Spectral Density

As expressed in Equation (100), 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, $\omega \propto sin\frac{ka}{2}$, 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 $\omega$ and $\left(\omega +d\omega \right)$ is denoted as $Ng\left(\omega \right)d\omega$, where N is the total number of atoms and $g\left(\omega \right)$ is called “spectral density”. Figure 11a,b inserted a modified phonon spectral density (i.e., frequency square spectral density), $\frac{g\left(\omega \right)}{2\omega }=\frac{1}{\pi \omega \sqrt{{{\omega }_L}^2-{\omega }^2}}$, for the reader’s reference.

4.3. Anharmonic Phonon Modes and Electrothermal Coupling

4.3.1. 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 (82) or the quartic term in Equation (6), the resonance frequency should include at least $2{\omega }_0$ or $3{\omega }_0$, higher order harmonic modes. According to the uncertainty of the resonance frequency from the definite $\omega ={\omega }_0$ to $\Delta \omega =\omega -{\omega }_0\vee$, the phonon distribution changes from $\infty$ 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

$$u_{\alpha }\left(t\right)\propto sin{k}_{\alpha }a·cos{\omega }_{\alpha }t$$

(122)

Taking into account a modified “phonon spectral density” (i.e., frequency square spectral density), $\frac{g\left(\omega \right)}{2\omega }$ as the summation weight, and the low frequency acoustic mode dispersion relation $k_{\alpha }\propto {\omega }_{\alpha }$, we calculate the sum of $u_n\left(t\right)$’s, which exist in $\Delta \omega$

$${{\displaystyle \sum }}_{\alpha =1}^n{u}_{\alpha }\left(t\right)\propto {{\displaystyle \sum }}_{\alpha =1}^n\frac{sin{\omega }_{\alpha }acos{\omega }_{\alpha }t}{{\omega }_{\alpha }}$$

(123)

Individual wave forms (time dependence) $u_{\alpha }\left(t\right)$ are illustrated for ${\omega }_{\alpha }$= 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 (${\omega }_{\alpha }$ = 0.2) are observed. When we sum up wave forms up to ${\omega }_{\alpha }=3$, the top wave packet can be obtained. Note that, if we add ${\omega }_{\alpha }$ up to $\infty$ 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).

4.3.2. 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 $\acute{n}$ deviates from the equilibrium value ${\acute{n}}_0$, thermal current arises as expressed by

$$q_x=\frac{1}{V}{\displaystyle \sum }\hslash \omega \left(\acute{n}-{\acute{n}}_0\right)v_x,$$

(124)

where $v_x$ is the energy transport velocity of phonons, equivalent to v in Equation (78a,b). There are two mechanisms to cause time variation of $\acute{n}$ 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

$$\frac{d\acute{n}}{dt}={\left(\frac{\partial \acute{n}}{\partial t}\right)}_{diff.}+{\left(\frac{\partial \acute{n}}{\partial t}\right)}_{decay}$$

(125)

Time variation in the latter term is described with a relaxation time (τ) approximation.²⁸⁾ Assuming the spirit of the Drude model, the phonon collision occurs in a distance $l=v_x\tau$ ($v_x$: 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:

$${\left(\frac{\partial \acute{n}}{\partial t}\right)}_{decay}=-\left(\acute{n}-{\acute{n}}_0\right)/\tau$$

(126)

When the “steady state” condition is assured, that is, $\frac{d\acute{n}}{dt}=0$, the local thermal equilibrium (i.e., dynamic equilibrium) is assumed, leading to the equation

$${\left(\frac{\partial \acute{n}}{\partial t}\right)}_{diff}=-{\left(\frac{\partial \acute{n}}{\partial t}\right)}_{decay}=-v_x\frac{\partial {\acute{n}}_0}{\partial T}\frac{\partial T}{\partial x}$$

(127)

Using the relaxation time approximation for the Boltzmann equation and assuming steady-state conditions, the phonon thermal conductivity λ can be determined:

$$\lambda =c_p{v_x}^2\tau$$

(128)

Let us now discuss the thermal conductivity changes under different electrical constraints. Phonon transmission speed $v_x$ 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₃₃ and λ^D₃₃), according to the crystallographic symmetry change from $m3m$ to $\infty mm$. Third, because the elastic compliance under the short-circuit condition (s^E₃₃) is softer than the elastic compliance under open-circuit conditions (s^D₃₃) via the relation s^E₃₃(1 − k₃₃²) = s^D₃₃, and the sound velocity $v_{33}^E$ is smaller than $v_{33}^D$, 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., D-constant), 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₃₃). This speculation can also suggest the analogous relations among two coupling factors: electromechanical k₃₃ and secondary electrothermal k^λ₃₃.

5. 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.

5.1. Pseudo-DC Piezoelectric Actuators

5.1.1. 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 V_e/A, where V_e 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 V_e generates the heat, which is dissipated through the area A. Thus, if we need to suppress the temperature rise, a small V_e/A design is preferred. Instead of one ML, four (1/4) small-ML’s connected in parallel are preferred.

5.1.2. 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, q_g, and that at which the heat is dissipated, q_d, can be expressed as

q_g − q_d = V ρ c_p (dT/dt),

(129)

where it is assumed that a “uniform temperature distribution” exists throughout the sample and V is the total volume, ρ is the mass density, and c_p 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, q_g, is expressed as:

q_g = w f V_e,

(130)

where w is the loss per driving cycle per unit volume, f is the driving frequency, and V_e 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), w_e, which is expressed in terms of the intensive dielectric loss tan δ’ as:

w = w_e = π ε^Xε₀ E₀² tan δ’.

(131)

If we neglect the transfer of heat through conduction, the rate of heat dissipation (q_d) from the sample is the sum of the rates of heat flow by radiation (q_r) and by convection (q_c):

q_d = q_r + q_c = eAσ (T⁴ − T_o⁴) + h_c A(T − T_o),

(132)

where e is the emissivity of the sample, A is the sample surface area, σ is the Stefan–Boltzmann constant, T_o is the initial sample temperature, h_c is the average “convective heat transfer coefficient”. Thus, Equation (98) can be written in the form:

w f V_e − A k(T) (T − T_o) = V ρ C_p (dT/dt),

(133)

where the quantity

k(T) = σe (T² + T_o²)(T + T_o) + h_c

(134)

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:

T − T_o = [w f V_e/k(T) A] (1 − e^−t/τ),

(135)

where the time constant τ is expressed as:

τ = ρC_p V/k(T) A.

(136)

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:

ΔT = w f V_e/k(T) A,

(137a)

while, as t → 0, the initial rate of temperature rise is given by

dT/dt = (w_T f V_e/ρC_p V) = ΔT/τ,

(137b)

where w_T 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, w_T, as calculated from Equation (137b) is given for three multilayer specimens in

Table 3. Loss and overall heat transfer coefficient for PZT multilayer samples (under E = 3 kV/mm, f = 300 Hz) [34].

Actuator	4.5 × 3.5 × 2.0 mm3	7.0 × 7.0 × 2.0 mm3	17 × 3.5 × 1.0 mm3
wT (kJ/m3) [$\frac{\rho C_{p}V}{f{V}_e}{\left(\frac{dT}{dt}\right)}_{t\to 0}$]	19.2	19.9	19.7
P-E hysteresis loss (kJ/m3)	18.5	17.8	17.4
k(T) (W/m2 K)	38.4	39.2	34.1

, 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 w_e 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 δ’.

5.2. Resonance Drive Piezoelectric Transducers

5.2.1. Heat Generation from a Resonating Piezoelectric Specimen

Tashiro et al. observed the heat generation in a rectangular piezoelectric plate during a resonating drive [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 k₃₁ 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.

5.2.2. 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 Q_A at the resonance under constant voltage drive, while by Q_B at the antiresonance under constant current drive. In the k₃₁ mode, we derived the mechanical quality factors Q_A and Q_B as follows [37]:

$${\mathit{Q}}_{\mathit{A},31}=\frac{1}{\tan {{\mathit{\varphi }}_{11}}^{\prime }}$$

(138a)

$$\frac{1}{{\mathit{Q}}_{\mathit{B},31}}=\frac{1}{{\mathit{Q}}_{\mathit{A},31}}-\frac{2}{1+{(\frac{1}{{\mathit{k}}_{31}}-{\mathit{k}}_{31})}^2{{\Omega }_{\mathit{B},31}}^2}(2\tan {{\mathbf{\theta }}_{31}}^{\prime }-\tan {{\mathbf{\delta }}_{33}}^{\prime }-\tan {{\mathit{\varphi }}_{11}}^{\prime })$$

(138b)

Here, tanδ_33′, tanϕ_11′, and tanθ_31′ are intensive loss factors for ε₃₃^X, s₁₁^E, d₃₁, respectively, and Ω_B,31 is the normalized antiresonance frequency:

$${\Omega }_{\mathit{A},31}=\frac{{\mathbf{\omega }}_{\mathit{a}}\mathit{l}}{2{\mathit{v}}_{11}^{\mathit{E}}}=\frac{\mathbf{\pi }}{2}\hspace{1em}\hspace{1em}\left[v_{11}^E=1/\sqrt{\rho {s_{11}}^E}\right],$$

(139a)

$${\Omega }_{\mathit{B},31}=\frac{{\mathbf{\omega }}_{\mathit{b}}\mathit{l}}{2{\mathit{v}}_{11}^{\mathit{E}}},$$

(139b)

where the antiresonance (normalized) frequency Ω_B,31 is k₃₁ dependent, and can be obtained from the equation $1-{{\mathit{k}}_{31}}^2+{{\mathit{k}}_{31}}^2\frac{\tan {\Omega }_{\mathit{B}}}{{\Omega }_{\mathit{B}}}=0$. Because the intensive dielectric loss tanθ_31′ is larger than $(\tan {{\mathbf{\delta }}_{33}}^{\prime }+\tan {{\mathit{\varphi }}_{11}}^{\prime })/2$ in PZT piezoceramics, Q_B at antiresonance is higher than Q_A 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 k₃₁ = 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 k₃₁ 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

5.2.3. Thermal Analysis on the Resonance Mode

5.2.3.1. Heat Transfer Modeling

We develop a 1D heat transfer model for the k₃₁ 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:

$$\frac{\partial T\left(x,t\right)}{\partial t}=\frac{\lambda }{c_p\rho }\frac{{\partial }^{2}T\left(x,t\right)}{\partial x^2}+\frac{q_g\left(x\right)}{c_p\rho }-\frac{h_d}{c_p\rho }\left[T\left(x,t\right)-T_{air}\right]$$

(140)

where λ is thermal conductivity [unit: W/m K], c_p (c^E in the k₃₁ case) is specific heat (unit: J/kg K), ρ density, and $\frac{\lambda }{c_p\rho }={\alpha }_T$ is called thermal diffusivity [unit: m²/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 $c_p\rho$), which may exhibit a sinusoidal distribution. The third term indicates the heat dissipation proportional to the temperature difference ΔT from the ambient temperature $T_{air}$.

Regarding the heat generation, we further assume that $q_g\left(x\right)$ is expressed in proportion to the square of the sinusoidally-distributing strain, as visualized in Figure 19:

$$q_g\left(x\right)=g_{h}co{s}^2\left(\frac{\pi x}{L}\right)$$

(141)

The heat dissipation h_d is a proportional constant to [T − T_air], similar to “heat transfer coefficient” introduced in Equation (131).

5.2.3.2. Temperature Distribution Profile Change with Time

We prepared a piezoelectric k₃₁ rectangular plate specimen with $80\times 14\times 2{ \mathrm{mm}}^3$ 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 k₃₁ specimen is shown in Figure 20b, which can be used for calculating the total thermal dissipation energy g_h. 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

$$Q_m=2\pi \frac{EnergyStored/Cycle}{EnergyLost/Cycle}$$

(142)

we can obtain Q_m from

$$Q_m=2\pi f\left(\frac{\rho V_{rms}^2}{h_g}\right)$$

(143)

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\times 14\times 2 {\mathrm{mm}}^3$) with the resonance frequency at 20.04 kHz at room temperature, and obtained Q_m = 507 by a 3-dB-down 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². Figure 21 shows the Q_m obtained at three frequencies slightly above the resonance frequency (20.04 kHz) [7]. The Q_m values around 550 agree with the extraporated values from the above 507. Thus, we can conclude that this thermal method can determine the Q_m reasonably at any frequency around the resonance and anti-resonance region. An increase of Q_m with increasing the frequency suggests that a maximum-Q_m frequency (i.e., the highest efficiency) exists between the resonance and antiresonance frequencies [38].

6. 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.