Pyroelectric properties of Pb[Zr_0.5Ti_0.5]O₃ studied with a hybrid density functional method

Abstract

The primary and secondary pyroelectric coefficients of four different ordered models of lead zirconate titanate Pb[Zr_0.5Ti_0.5]O₃ solid solution have been investigated with a hybrid density functional method. Phonon anharmonicity and finite-temperature phonon properties necessary for the study of primary pyroelectricity are determined with the use of self-consistent phonon theory. Secondary pyroelectricity and lattice thermal expansion are investigated with quasi-harmonic approximation. The lattice thermal conductivity and other physical properties are also reported at the same level of theory. The largest absolute value of the pyroelectric coefficient is obtained for a rocksalt-type ordering. The results show that the pyroelectric properties of B-site solid-solution perovskites can be studied with the computational methodology previously used for simple titanate compounds BaTiO₃ and PbTiO₃, enabling further theoretical insights into computational screening of prospective new ferroelectric pyroelectrics.

---

1 Introduction

Heat fluctuation can be extracted as electrical energy with polar crystalline materials that exhibit the pyroelectric effect.¹ The strength of the pyroelectric effect is determined by the pyroelectric coefficient p, defined as the ratio of change in spontaneous polarization P_s and change in temperature T, p = ΔP_s/ΔT.² The coefficient consists of fixed-volume primary p⁽¹⁾ and secondary p⁽²⁾ parts, the secondary effect arising from the piezoelectric response to the thermal deformation of the unit cell. In ferroelectrics, a subclass of pyroelectrics, the spontanous polarization can be reversed. In ferroelectric materials, the pyroelectric coefficients are typically larger than in non-ferroelectrics, making them important for practical applications.³ In antiferroelectrics, the dipole polarizations are ordered antiparallel, so that they cancel out in the unit cell, resulting in zero net polarization.⁴ Thermal and waste energy harvesting is a prospective novel application for pyroelectrics, requiring further development of pyroelectric materials and enhancement of the pyroelectric coefficient for higher output.^3,5–7

The perovskite-type lead titanate PbTiO₃ was discovered to be a ferroelectric in 1951, with a single phase transition at 763 K from the tetragonal ferroelectric to a cubic paraelectric phase.^8,9 Around the same time, the lead zirconate PbZrO₃ with an orthorhombic crystal structure was found to be an antiferroelectric with a transition from the antiferroelectric to the cubic paraelectric phase at 503 K.^10–12 Further theoretical insights into the complex phases and phase transitions of PbZrO₃ require taking into account finite-temperature effects.¹³ The B-site solid solution of PbTiO₃ and PbZrO₃, Pb[Zr_1−xTi_x]O₃ (PZT), was soon found to exhibit ferroelectric behavior.^14,15 The phase diagram of the PZT solid solution is more complex than that of the parent phases.^16,17 The morphotropic phase boundary (MPB) at the composition Pb[Zr_0.52Ti_0.48]O₃, between the rhombohedral and tetragonal phases, is the optimal composition for performance in most technological applications due to the high piezoelectric response of this composition.¹⁸ In addition to technological importance, PZT is frequently used to further study ferroelectric and pyroelectric phenomena, such as ferroelectric domain structure effects on pyroelectricity.^19,20 Ferroelectricity can also be induced in PbZrO₃ in the polar R3c phase by high electric fields²¹ and in epitaxial thin films,²² and has computationally been predicted to be achievable at the extreme thickness limit.²³

Due to the technological importance of the PZT system, its piezoelectric properties have been investigated computationally for decades, starting with a GGA/LAPW+LO study for ordered structures of Pb[Zr_0.5Ti_0.5]O₃ in P4mm and I4mm space groups.²⁴ Further pioneering work has been done on the finite-temperature properties of PZT,²⁵ along with further studies on the piezoelectric properties with LDA/LAPW-LO,²⁶ GGA/FLAPW,²⁷ and DFPT.²⁸ Elastic properties of tetragonal PbZr_0.5Ti_0.5O₃ in planar, columnar, and rocksalt orderings has been investigated with LDA/mixed-basis pseudopotential.²⁹ The local atomic structure together with the dielectric and piezoelectric properties of tetragonal PZT have recently been investigated using 2 × 2 × 2 supercells in various configurations and with compositions x = 0.25, x = 0.50, and x = 0.75.³⁰ Furthermore, intrinsic piezoelectricity in the monoclinic PZT phase has been investigated with virtual crystal approximation.³¹ Soft-mode lattice dynamics and phonon dispersions of cubic PZT models have been investigated using primitive 2 × 2 × 2 supercells combined with virtual crystal approximation.³²

Despite the considerable previous theoretical studies on PZT, to the best of our knowledge, both the pyroelectric properties and lattice thermal conductivity remain entirely unexplored with first-principles methods. We have recently reported the calculation of pyroelectric properties of tetragonal ferroelectrics BaTiO₃, KNbO₃, and PbTiO₃ from first principles.^33,34 Here, we study the pyroelectricity, lattice thermal conductivity, and other physical properties of ordered models of Pb[Zr_0.5Ti_0.5]O₃ (P4mm) solid solution with DFT-PBEsol0 hybrid density functional method. Ordered structure models with planar, columnar, and rocksalt ordering of the B-site cations are studied. Extending the computational methodology for predicting pyroelectric properties is important for deepening the atomic-level understanding of pyroelectricity, enabling further insights into screening of novel solid solution compositions.

2 Computational methods

2.1 Crystal structure and physical properties

Quantum chemical calculations have been carried out using the CRYSTAL23³⁵ program package and hybrid PBEsol0^36–38 density functional method (DFT-PBEsol0, 25 exact exchange). Gaussian-type triple-ζ-valence + polarization (TZVP) basis sets, derived from Karlsruhe def2 basis,³⁹ were used for all elements. The basis sets for Ti, Zr, and O are described in ref. 40 and the basis set for Pb in ref. 34. Tight tolerance factors (TOLINTEG) of 8, 8, 8, 8, and 16 were used for the evaluation of the Coulomb and exchange integrals in the geometry optimization and property calculations. Default CRYSTAL23 structural optimization convergence criteria and DFT integration grids were used.

Spontaneous polarization was calculated with the Berry phase approach implemented in CRYSTAL (SPOLB keyword).^41,42 The piezoelectric and second-order elastic tensors at 0 K were obtained with the ELAPIEZO keyword implemented in CRYSTAL.^42,43 The Monkhorst–Pack-type k-meshes⁴⁴ used in the calculations for different structures and different steps are described in the ESI† (Table S1).

2.2 Phonon properties

Self-consistent phonon (SCPH) calculations⁴⁵ were carried out with the ALAMODE code^46,47 (version 1.5.0) employing a recently developed CRYSTAL interface.³³ To obtain finite-temperature phonons and phonon properties with ALAMODE, a number of randomly displaced or MD-based structures are first used for fitting force constants up to the quartic or higher order. For optimization of force constants, the elastic net linear model was used for all PZT structures except for the planar P4mm structure, for which the ordinary least squares model was used. These force constants are then used in the SCPH procedure for the calculation of the temperature-dependent phonon frequencies and eigenvectors. A tight SCF convergence criterion (TOLDEE keyword) of 10⁻¹⁰ au was applied in force calculations used for fitting the force constants. Fitting errors of quartic force constants, 2.38% for planar P4mm, 3.15% for planar Pmm2, 3.64% for rocksalt I4mm, and 2.41% for columnar P4mm, were slightly higher than in previous work (1.62% for PbTiO₃³⁴). Full details of the random displacements, force constant fitting and SCPH calculations are given in the ESI.† The results from the SCPH calculations employing force constants up to quartic order were used to plot the finite-temperature phonon dispersions and to obtain the temperature correction to second-order force constants for the calculation of the primary pyroelectric coefficient with the Berry phase approach (see below).

Harmonic frequencies obtained at the Γ point with the approach implemented in CRYSTAL^48–50 are reported in the ESI† (Tables S2–S6). Lattice thermal conductivity was evaluated within the relaxation-time approximation (RTA) for the Boltzmann transport equation as implemented in ALAMODE.^47,51 The quartic force constants were used. Isotope effects have been included as implemented in ALAMODE.⁵² Quasiharmonic approximation (QHA) implemented in CRYSTAL^53–55 was used to obtain the thermal expansion coefficients, necessary for calculating the secondary pyroelectric coefficient (see below). For the columnar P4mm structure, the QHA calculation was carried out in the P4 space group due to symmetry-related convergence issues.

Example ALM input (for fitting cubic and quartic force constants) and ANPHON input (for SCPH calculations) are given in the ESI.† The reciprocal space path used in plotting the phonon dispersion relations is also included in the ANPHON input examples given in the ESI.†

2.3 Pyroelectric properties

Calculation of the pyroelectric properties follows the scheme developed and described in ref. 33 for tetragonal BaTiO₃ and KNbO₃ and in ref. 34 for PbTiO₃. Because pyroelectricity is a tensor property, it is only necessary to evaluate the phonon modes affecting the polar c axis, as in the point group 4mm the pyroelectric coefficients in other directions are zero.² Compared to PbTiO₃, the number of phonon modes is different for ordered supercells of PZT, as the unit cell is doubled to accomodate for the Zr substitution and ordering. For structures that retain the point group 4mm (planar P4mm, rocksalt I4mm, columnar P4mm), there are now six optical modes at the Γ point with the irreducible representation A₁, as opposed to three in PbTiO₃.⁵⁶ These A₁ modes correspond to atomic displacements in the direction of the polar c axis (Cartesian z direction). The primary pyroelectric coefficient was thus calculated for these six A₁ modes of tetragonal PZT modifications. For the columnar structure in the space group Pmm2 there are nine A₁ modes.

Temperature-corrected interatomic force constants from the SCPH yield finite-temperature phonon eigenvectors. They are used for finite-temperature displacements, representing finite-temperature and mode-specific behavior of ordered PZT structures at elevated temperatures. These displacements were created using the normal coordinate Q basis displacements in ALAMODE, modified so that the displacements are not drawn randomly, but instead are carried out at the Γ point only, and their magnitude is calculated from the phonon frequency ω and input temperature. Direction of displacement was kept uniform for each A₁ mode. The primary pyroelectric coefficient p⁽¹⁾ for temperature T is obtained with the central difference method as p⁽¹⁾ = (P_T2 − P_T1)/(T₂ − T₁), where T₂ = T + 20 K and T₁ = T − 20 K. The choice of the temperature step of ΔT = ±20 K is based on observations in previous work on BaTiO₃ and PbTiO₃.^33,34 The temperature-dependence of the finite-temperature displacements obtained from the phonon eigenvectors at a certain temperature was found to be linear and a temperature step of 20 K was found to be sufficiently large to yield numerically stable results with the Berry phase method.

As in the previous work,^33,34 the thermal expansion of the unit cell is not taken into account when calculating the primary pyroelectric coefficient. The secondary pyroelectric effect p⁽²⁾, which is a piezoelectric effect driven by thermal expansion, is taken into account separately. This is done by calculating the piezoelectric coefficients d (in CN), second-order elastic constants C (in N m⁻²) and thermal expansion coefficients α (in K⁻¹). In the ∞mm point group the secondary pyroelectric coefficient is given as p⁽²⁾ = 2d₃₁(c₁₁α₁ + C₁₂α₁ + C₁₃α₃) + d₃₃(2C₁₃α₁ + C₃₃α₃).⁵⁷ In the mm2 point group the secondary pyroelectric coefficient is given as p⁽²⁾ = d₃₁(c₁₁α₁ + C₁₂α₂ + C₁₃α₃) + d₃₂(c₁₂α₁ + C₂₂α₂ + C₂₃α₃) + d₃₃(2C₁₃α₁ + C₃₃α₃).⁵⁷

As in the previous work, the tertiary pyroelectric effect⁵⁸ and macroscopic effects caused by orientation, the structure of the ferroelectric domain, and domain walls^19,20 affecting the total pyroelectric coefficient are not included here.

3 Results

3.1 Structural properties and spontaneous polarization

Typical B-site orderings for A[BB′]O₃ perovskites are planar, columnar, and rocksalt configurations.⁵⁹ Five different ordered models for Pb[Zr_0.5Ti_0.5]O₃ solid solution were investigated: two planar (P4mm and Pmm2), two columnar (P4mm and Amm2), and rocksalt (I4mm). For the tetragonal P4mm space group of the parent PZT structure, there are two variations of the planar and columnar orderings, depending on whether the ordering is along the ab plane or along the c axis. Stacking the B and B′ cation layers on top of each other in the c axis direction results in the P4mm space group, while stacking them along the a (or b) axis results in Pmm2 space group. In the columnar ordering, having the columns run along the c axis results in space group P4mm, while having them run along the a (or b) axis results in space group Amm2. The ordered models are illustrated in Fig. 1, and given in CIF format in the ESI.† It should be noted that these ordered models are just the simplest possible approximations, and do not fully scope the set of orderings that would require larger supercells to accommodate more complex orderings. Experimental investigations of the local order reveal that monoclinic short range order can persist even in the tetragonal or rhombohedral phases with long range order.^60–62

Fig. 1 Illustration of the optimized structures of the five different ordered Pb[Zr_0.5Ti_0.5]O₃ models used in this work. Ti and its octahedra in blue, Zr and its octahedra in green, O in red, and Pb in grey. All structures are viewed from the same direction. Visualization made with the VESTA software.⁶³

The relative electronic energies per formula unit (ΔE, kJ mol⁻¹ Z⁻¹) together with the absolute values of spontaneous polarization (P_s, C m⁻²) are given in Table 1 for all ordered models. The higher symmetry variants for the ordered columnar and planar models are energetically more favorable. The energetics of the different ordered models were further investigated on the basis of the QHA results, as discussed below. The columnar-ordered P4mm structure can be identified as the lowest-energy structure, while the highest-energy columnar structure Amm2 was omitted from the calculation of the pyroelectric properties. The values of P_s for the ordered models range from 0.41 to 0.48 C m⁻². For comparison, a P_s value of 0.51 C m⁻² was obtained for PbTiO₃ at the same level of theory.³⁴ Composition, sample type, as well as domain effects can explain differences in the spontaneous polarization between the experimental measurements and calculated values. To evaluate P_s of a truly disordered PZT, a thermodynamic average would be needed over a larger number of ordered models with more complex compositions.

Table 1 Relative electronic energies ΔE, and the absolute values of spontaneous polarization P_s for the ordered PZT models. The relative energies are given per one formula unit (Z). Experimental values of spontaneous polarization are reported in ref. 64 calculated from saturation polarization and in ref. 19 from room temperature ferroelectric hysteresis response

PZT structure	ΔE, kJ mol^−1 Z^−1	P s, C m^−2
Planar P4mm	6.9	0.45
Planar Pmm2	5.1	0.41
Rocksalt I4mm	5.0	0.48
Columnar P4mm	0.0	0.48
columnar Amm2	7.5	0.47
Exp. PbZr0.4Ti0.6O3^64		0.90
Exp. PbZr0.52Ti0.48O3^19		0.30

---

The Ti–O and Zr–O bond lenghts and the shortest Pb–O, Pb–Ti, and Pb–Zr distances of the studied ordered PZT models are listed in Table 2. In the planar P4mm structure, the Zr–O bond length in the polar c direction is clearly elongated, while in the Pmm2 structure the Ti–O length is elongated. In the columnar P4mm, the Ti–O lengths along the c axis are elongated, while in the columnar Amm2 structure the longer Ti–O and Zr–O bond lengths are closest to each other for the investigated PZT models. Notably in the planar P4mm, rocksalt I4mm, and columnar Amm2 structures, the longer Ti–O bond length is shorter than in PbTiO₃. In the lowest-energy columnar P4mm structure, this Ti–O bond length is longer than in PbTiO₃. For the three energetically most unfavorable structures (planar P4mm, planar Pmm2, and columnar Amm2), all or some Ti–O bond lengths in the ab direction are larger compared to PbTiO₃.

Table 2 Ti–O and Zr–O bond lengths and the shortest Pb–O, Pb–Ti, and Pb–Zr distances in the studied ordered Pb[Zr_0.5Ti_0.5]O₃ models (all distances in Å units)

System	Ti–O (c direction)	Ti–O (ab direction)	Zr–O (c direction)	Zr–O (ab direction)	Pb–O (shortest)	Pb–Ti (shortest)	Pb–Zr (shortest)
PbTiO3^34	1.76, 2.38	1.96	—	—	2.49	3.34	—
PZT planar P4mm	1.75, 2.19	2.01	1.97, 2.46	2.05	2.52	3.31	3.47
PZT planar Pmm2	1.74, 2.51	1.92, 2.04	1.97, 2.29	2.04, 2.11	2.53	3.37	3.42
PZT rocksalt I4mm	1.81, 2.17	1.96	1.94, 2.36	2.08	2.53	3.37	3.39
PZT columnar P4mm	1.74, 2.52	1.96	1.98, 2.27	2.08	2.50	3.40	3.38
PZT columnar Amm2	1.76, 2.27	1.95, 2.03	1.97, 2.34	2.04, 2.09	2.44	3.33	3.43

---

3.2 Elastic and piezoelectrics properties

The calculated bulk moduli of the PZT structures are 122 GPa for columnar P4mm, 131 GPa for rocksalt I4mm, 127 GPa for planar P4mm, 120 GPa for planar Pmm2, and 120 GPa for columnar Amm2. Values from 128 GPa up to 186 GPa have been reported in earlier computational work for various ordered models.²⁹ The full elastic constants are given in ESI† (Table S7).

The calculated piezoelectric strain coefficients d (pC N⁻¹ or pm V⁻¹) and the stress constant ε₃₃ (C m⁻²) of the PZT ordered models are presented in Table 3 along with available experimental data.^65,66 Values for Pmm2 and Amm2 are averaged as d̄₃₁ = (d₃₁ + d₃₂)/2 and d̄₁₅ = (d₂₄ + d₁₅)/2. For the stress constant, calculated values previously reported for ordered Pb[Zr_0.5Ti_0.5]O₃ models are also included. The difference between our calculated values of d and ε and the experimentally reported ones is considerable. The values of ε₃₃ are in line with previous theoretical studies, but are much lower than the value based on experimental data.^24,65 The largest value of ε₃₃ is obtained for the rocksalt I4mm structure.

Table 3 Calculated piezoelectric strain constants d (pC N⁻¹) and stress constants ε₃₃ (C m⁻²) of ordered structures of tetragonal Pb[Zr_0.5Ti_0.5]O₃ compared with available experimental and previous computational data. Abbreviation pl. stands for planar and rs. for rocksalt ordering

	d 31 (pC N^−1)	d 33 (pC N^−1)	d 15 (pC N^−1)	ε 33 (C m^−2)
Planar P4mm	−5.0	28.6	70.8	3.40
Planar Pmm2	−6.8	37.6	129.9	3.13
Rocksalt I4mm	−7.8	34.5	155.1	4.40
Columnar P4mm	−6.4	32.0	139.4	3.09
Columnar Amm2	−6.4	33.9	132.9	3.50
Exp. PZT 50/50^65	43	110	166	27.0
Exp. PZT^66	−120 to −170	60–130
Calc. pl., P4mm^24				4.81
Calc. rs., I4mm^24				3.60
Calc. pl., P4mm^27				3.42
Calc. rs., I4mm^27				5.10
Calc. pl., P4mm^28				ca. 4

---

One possible explanation for the discrepancy is the proximity of the computational composition to the MPB, which is known to significantly enhance the piezoelectric response. The change in the range of structural order has been proposed as the reason for the enhanced properties close to the MPB,⁶⁰ something that our ordered models cannot replicate. The fact that elastic and piezoelectric constants are calculated for perfect crystals at 0 K, omitting defects, temperature effects, and possible anharmonicity, can also contribute to the difference.⁶⁷ Furthermore, the experimental values also differ greatly, while many of the measured and reported values are for commercial ceramics, such as PZT-5H, that have compositions different from pure Pb[Zr_1−xTi_x]O₃.⁶⁸

3.3 Phonon properties

The finite-temperature phonon dispersion of the lowest-energy columnar P4mm ordering of PZT is illustrated in Fig. 2. While experimental measurements on the lattice dynamics of PZT near the MPB have been carried out with various spectroscopic methods,⁶⁹ full phonon dispersion relations have not been reported. Based on the measurements of the phonon frequencies, the disordered PZT solid solution shows one-mode behavior corresponding to mixed B–O modes, while our calculated phonon dispersions show two-mode behavior, where the Ti–O and Zr–O modes are distinguishable. The phonon dispersions and the related reciprocal space paths for the other PZT structures are included in the ESI† (Fig. S1–S4). Differences between the harmonic and anharmonic dispersions are relatively small for the studied models, as the parent tetragonal PZT phase is known to be thermodynamically stable down to 0 K in the PbTiO₃–PbZrO₃ phase diagram.¹⁶ The planar ordered models show a few small imaginary phonon modes outside the Γ point, around the S and M points. These imaginary modes might arise from the fact that ordered models do not represent true minima for PZT, which in reality is a disordered B-site solid solution. Furthermore, in the rocksalt I4mm structure, small numerical inaccuracies persist along the Γ–M path. The pyroelectric properties of the ordered models are calculated using the optical A₁ phonon modes at the Γ point, and thus the slight imaginary frequencies for the acoustic phonon modes are not likely to affect the calculated pyroelectric coefficients.

Fig. 2 Anharmonic 300 K (red, solid) and harmonic (black, dotted) phonon dispersion of columnar P4mm PZT. Phonon dispersions of other ordered models are included in the (ESI†).

The quasiharmonic approximation (QHA) implemented in CRYSTAL was used to investigate thermal expansion of ordered PZT models and to obtain thermal expansion coefficients for the calculation of secondary pyroelectric coefficient. Lattice thermal expansion was evaluated from 0 to 500 K. Experimentally, an average thermal expansion coefficient α = (α_a + α_b + α_c)/3 of ca. 2 × 10⁻⁶ K⁻¹ has been reported at 300 K.⁷⁰ The calculated values of α_c are presented in Fig. 3. The calculated average values of α are included in the ESI† (Fig. S5). Negative thermal expansion along the c axis is observed for multiple PZT structures, most notably planar P4mm and rocksalt I4mm orderings at the low temperature region. The thermal expansion of the columnar P4mm stays uniformly positive throughout the investigated tempreature range. As the value of α is obtained with QHA, anharmonicity is not fully included, possibly affecting the results.

Fig. 3 Thermal expansion coefficients in the c direction, α_c (K⁻¹), calculated with quasiharmonic approximation in a temperature range of 10 to 500 K. Data for PbTiO₃ from ref. 34.

QHA can also be used to evaluate the relative energetics of the structures at finite temperatures. The calculated relative Helmholtz free energies per formula unit ΔF (kJ mol⁻¹ Z⁻¹) of the studied PZT models are presented in Fig. 4. Noticeably, the energy difference between the phases become even smaller compared to the electronic energies in Table 1, with the difference between columnar P4mm and rocksalt I4mm reduced to 0.8 kJ mol⁻¹ Z⁻¹ at 700 K. Gibbs free energies ΔG (kJ mol⁻¹ Z⁻¹) were also investigated at 300 K with the same phonon sampling as in the QHA results. The energy order remains the same, with the columnar P4mm structure at 0 kJ mol⁻¹ Z⁻¹, followed by rocksalt I4mm (1.1 kJ mol⁻¹ Z⁻¹), planar P4mm (3.8 kJ mol⁻¹ Z⁻¹), planar Pmm2 (3.8 kJ mol⁻¹ Z⁻¹), and columnar Amm2 (6.6 kJ mol⁻¹ Z⁻¹).

Fig. 4 Relative Helmholtz free energies ΔF (kJ mol⁻¹ Z⁻¹) calculated with quasiharmonic approximation in a temperature range of 100 to 700 K.

3.4 Lattice thermal conductivity

The lattice thermal conductivity of various PZT compositions, including Pb[Zr_0.5Ti_0.5]O₃, has been measured by Tachibana et al.⁷¹ The experimental and our calculated total lattice thermal conductivity κ (W m⁻¹ K⁻¹) values are presented in Fig. 5. Quartic force constants were employed in the lattice thermal conductivity calculations. The results obtained with DFT-PBEsol0/TZVP are in agreement with experimental data, with PZT having a smaller lattice thermal conductivity than the parent phase PbTiO₃, as also observed experimentally. It should also be noted that the experimentally measured values are strongly affected by the sample type, as observed in Fig. 5 for PbTiO₃ data using a single crystal⁷² and polycrystalline sample,⁷¹ the polycrystalline sample having a lattice thermal conductivity ca. 2 W m⁻¹ K⁻¹ smaller than the single crystal sample around 300 K. The presented calculated values are an average of contributions in the three Cartesian directions.

Fig. 5 Average lattice thermal conductivity of ordered PZT models (red, blue, green, and yellow, dashed) and PbTiO₃ (black, solid line)³⁴ calculated with RTA, compared with experimental data of Tachibana et al.⁷² (grey, dotted) for PbTiO₃ single crystal (low temperature data up to 400 K) and Tachibana et al.⁷¹ (grey, stars) for polycrystalline PbTiO₃ (high temperature data from 300 K) and Pb[Zr_0.5Ti_0.5]O₃ (grey, crosses).

3.5 Pyroelectric properties

Fig. 6 shows the primary p⁽¹⁾ and secondary p⁽²⁾ pyroelectric coefficient of the studied ordered PZT models from 100 K to 500 K. The largest absolute value of the total pyroelectric coefficient is found for the rocksalt ordering, −549 μC m⁻² K⁻¹ at 200 K. At 300 K, the pyroelectric coeffients are −506 μC m⁻² K⁻¹ for rocksalt I4mm, −346 μC m⁻² K⁻¹ for columnar P4mm, −184 μC m⁻² K⁻¹ for planar Pmm2, and −110 for planar P4mm. The coefficients at each temperature point are also given in the ESI† (Tables S8–S11).

Fig. 6 Calculated pyroelectric coefficients for planar P4mm (red), planar Pmm2 (blue), columnar P4mm (yellow) and rocksalt I4mm (green) orderings. Primary coefficient p⁽¹⁾ in solid lines and total p (p⁽¹⁾ + p⁽²⁾) in dashed lines. Experimental data in diamonds for PZT-5H⁷⁶ and PZNTU⁷⁷ ceramics, and with crosses for thin film,¹⁹ sol–gel thick film,⁷⁵ and thick film.⁷³

Experimental values of the total pyroelectric coefficient p of PZT vary greatly depending on the sample type and composition. Reported values range from −110 μC m⁻² K⁻¹ (ref. 73) to −1000 μC m⁻² K⁻¹ (ref. 74) for some multilayer thick films. A value of ca. −200 μC m⁻² K⁻¹ has been reported for a thin film with composition Pb[Zr_0.52Ti_0.48]O₃¹⁹ and a value of −352 μC m⁻² K⁻¹ for a thick film deposited with sol–gel method.⁷⁵ Ceramic samples PZT-5H (Pb_1.0[Zr_0.49Ti_0.46 (Nb_0.25Sb_0.75)_0.05]_1.0O₃)⁷⁶ and PZFNTU (Pb[Zr_1−(2x+y)Fe_xNb_xTi_y]_1−zU_zO₃),⁷⁷ for which the coefficient has been measured as −416 μC m⁻² K⁻¹ and −380 μC m⁻² K⁻¹, respectively, are far from the Pb[Zr_0.5Ti_0.5]O₃ composition used in the calculations. On the other hand, our calculations are for a bulk system, and thus fail to reproduce effects, such as strain, present in thin and thick films.

Due to the negative thermal expansion in the low temperature regions of the rocksalt I4mm and planar P4mm structures, their secondary coefficients are also negative close to 100 K. For columnar P4mm and Pmm2 the secondary contribution is positive throughout the temperature range. This starkly contrasts with the negative and much larger contribution observed recently in thin film PZT samples,¹⁹ with a value of ca. −75 μC m⁻² K⁻¹ for Pb[Zr_0.52Ti_0.48]O₃. However, the secondary contribution in Zr-rich PbZr_0.95Ti_0.05O₃ has been evaluated as positive with a value of 37.7 μC m⁻² K⁻¹.⁷⁸ The major difference compared to experiment in the piezoelectric strain constants near the MPB is the likely reason for the smaller values of the secondary effect. Noticeably it is also smaller in magnitude than what we have calculated for PbTiO₃, for which we obtained a value of p⁽²⁾ of −43 μC m⁻² K⁻¹ at 300 K.³⁴ A separate plot of p⁽²⁾ is available in the ESI† (Fig. S6).

The effect of the discrepancy of the calculated and experimental values of d₃₁ and d₃₃ can further be evaluated by using the experimental values for the calculation of p⁽²⁾. The p⁽²⁾ of rocksalt I4mm is 5.2 μC m⁻² K⁻¹ at 300 K according to DFT. Using the smaller and larger experimental values of d₃₁ and d₃₃ reported in ref. 66, respectively, values of 40 μC m⁻² K⁻¹ and 46 μC m⁻² K⁻¹ are obtained. This difference is rather significant for the secondary effect. Combined with the primary effect, however, the obtained values of p are not considerably changed, nor does the difference affect the relative ordering of the ordered structures according to the value of p.

The value of p⁽¹⁾ is obtained by summing contributions over modes contributing to pyroelectricity. This includes the A₁ modes with atomic displacements in the direction of the polar c axis (Cartesian z direction). The mode-specific contributions for the structure with the largest absolute value of the coefficient, rocksalt I4mm, are given in Fig. 7. As previously found for BaTiO₃, KNbO₃,³³ and PbTiO₃,³⁴ the major contributions arise from the lower frequency phonon modes, corresponding to the modes with indices 6 and 12 in those simple perovskites. In the solid-solution model here, modes 12 and 15 arise from the splitting of mode 6 in the simple perovskite, while 21 and 28 arise from the splitting of mode 12, rendering the contributions arising from the splitting of mode 15 into modes 29 and 30 negligible. Just like for BaTiO₃, KNbO₃, and PbTiO₃, a positive, canceling contribution to the total pyroelectric effect arises from higher frequency mode 28. Mode-specific contributions for other PZT ordered models are given in ESI† (Fig. S7–S9).

Fig. 7 Mode-specific contributions (colored dashed and dotted lines) to the primary pyroelectric coefficient p⁽¹⁾ (black solid line) of rocksalt-ordered I4mm PZT structure. Due to doubling of the cell size, the number of A₁ modes is doubled compared to tetragonal PbTiO₃.³⁴ Modes 12 and 15 (green and olive circles) originate from A₁ mode 6 of PbTiO₃, while modes 21 and 28 (red and pink triangles) originate from PbTiO₃ A₁ mode 12, and modes 29 and 30 (blue and cyan squares) from PbTiO₃ A₁ mode 15.

4 Conclusions

We have demonstrated that our computational methodology^33,34 for determining the pyroelectric coefficients of ferroelectric perovskites can be applied for solid solutions and beyond the tetragonal crystal system. Despite the lowering of the space group symmetry and the increase in the number of A₁ phonon modes, the methodology based on summing contributions to pyroelectricity over phonon modes yields consistent results. The pyroelectric coefficients of the four lowest-energy ordered models of PZT have been reported. The energetically most favorable PZT structure investigated, columnar P4mm, has the most elongated Ti–O bond length of 2.52 Å. The largest absolute values of the pyroelectric coefficient are found for ordered models that also have the largest value of spontaneous polarization.

The primary pyroelectric coefficient of rocksalt I4mm PZT is larger than that of tetragonal P4mm PbTiO₃. The rocksalt I4mm structure also exhibits the largest value of the piezoelectric stress constant ε₃₃ of the studied ordered models. The columnar P4mm ordering has a value of the primary pyroelectric coefficient closer to that of PbTiO₃, while the primary pyroelectric coefficients of the planar P4mm and Pmm2 are smaller than in PbTiO₃. The calculated secondary contribution to pyroelectricity is not as pronounced for PZT structures as it is for PbTiO₃. Noticeably, unlike for PbTiO₃, for PZT structures negative thermal expansion is only observed at the lower temperature range, and not for all ordered models. The piezoelectric properties calculated for PZT differ from the experimental results, which may explain the differences in secondary pyroelectric contributions.

While PZT is a solid solution without B-site cation ordering, in true Pb[BB′]O₃ compounds ordering depends on the synthesis conditions in a controllable way.⁵⁹ As the ordered PZT models undergo considerable distortion in structural optimizations, the investigations could be expanded by structural optimization at finite temperatures.⁷⁹ Simultaneous investigation of structural properties, phonon properties, piezoelectricity, and pyroelectric coefficients at the same computational level of theory can enable further insight into the discovery and design of new ferroelectric and pyroelectric materials via computational screening and prediction.

Author contributions

Kim Eklund: investigation, visualization, writing – original draft preparation. Antti. J. Karttunen: conceptualization, writing – review and editing, supervision, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

The authors thank the Research Council of Finland for funding (grant no. 363733), and the Finnish IT Center for Science (CSC) for computational resources. K. E. thanks the Aalto University School of Chemical Engineering for funding and the Walter Ahlström Foundation, the Magnus Ehrnrooth Foundation, and the Finnish Foundation for Technology Promotion (TES) for grant support.

Notes and references

1. S. B. Lang, Sourcebook of pyroelectricity, Gordon and Breach Science Publishers, London, New York, 1974.
2. J. F. Nye, Physical properties of crystals: their representation by tensors and matrices, Clarendon, Oxford, 1957.
3. C. R. Bowen, J. Taylor, E. LeBoulbar, D. Zabek, A. Chauhan and R. Vaish, Energy Environ. Sci., 2014, 7, 3836–3856.
4. C. A. Randall, Z. Fan, I. Reaney, L. Chen and S. Trolier-McKinstry, J. Am. Ceram. Soc., 2021, 104, 3775–3810.
5. Y. Zhang, M. Xie, J. Roscow, Y. Bao, K. Zhou, D. Zhang and C. R. Bowen, J. Mater. Chem. A, 2017, 5, 6569–6580.
6. Y. Zhang, P. T. T. Phuong, E. Roake, H. Khanbareh, Y. Wang, S. Dunn and C. Bowen, Joule, 2020, 4, 301–309.
7. Y. Zhou and G. W. Ho, Next Energy, 2023, 1, 100026.
8. G. Shirane and S. Hoshino, J. Phys. Soc. Jpn., 1951, 6, 265–270.
9. R. Nelmes and W. Kuhs, Solid State Commun., 1985, 54, 721–723.
10. E. Sawaguchi, G. Shirane and Y. Takagi, J. Phys. Soc. Jpn., 1951, 6, 333–339.
11. E. Sawaguchi, H. Maniwa and S. Hoshino, Phys. Rev., 1951, 83, 1078.
12. F. Jona, G. Shirane, F. Mazzi and R. Pepinsky, Phys. Rev., 1957, 105, 849–856.
13. H. Zhang, H.-C. Thong, L. Bastogne, C. Gui, X. He and P. Ghosez, Phys. Rev. B, 2024, 110, 054109.
14. E. Sawaguchi, J. Phys. Soc. Jpn., 1953, 8, 615–629.
15. B. Jaffe, R. S. Roth and S. Marzullo, J. Appl. Phys., 1954, 25, 809–810.
16. J. Frantti, S. Ivanov, S. Eriksson, H. Rundlöf, V. Lantto, J. Lappalainen and M. Kakihana, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 064108.
17. D. I. Woodward, J. Knudsen and I. M. Reaney, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 104110.
18. A. J. Bell, J. Mater. Sci., 2006, 41, 13–25.
19. G. Velarde, S. Pandya, L. Zhang, D. Garcia, E. Lupi, R. Gao, J. D. Wilbur, C. Dames and L. W. Martin, ACS Appl. Mater. Interfaces, 2019, 11, 35146–35154.
20. C.-C. Lin, Y. Hu, J. Kim, D. Lou, A. Bhat, P. Kavle, T. Y. Kim, C. Dames, S. Liu and L. W. Martin, Phys. Rev. X, 2025, 15, 011063.
21. O. E. Fesenko, R. V. Kolesova and Y. G. Sindeyev, Ferroelectrics, 1978, 20, 177–178.
22. P. Dufour, T. Maroutian, M. Vallet, K. Patel, A. Chanthbouala, C. Jacquemont, L. Yedra, V. Humbert, F. Godel, B. Xu, S. Prosandeev, L. Bellaiche, M. Otoničar, S. Fusil, B. Dkhil and V. Garcia, Appl. Phys. Rev., 2023, 10, 021405.
23. N. Maity, M. Haddad, N. Bassiri-Gharb, A. Kumar, L. Jones, S. Lisenkov and I. Ponomareva, npj Comput. Mater., 2025, 11, 48.
24. G. Sághi-Szabó, R. E. Cohen and H. Krakauer, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 12771–12776.
25. L. Bellaiche, A. García and D. Vanderbilt, Phys. Rev. Lett., 2000, 84, 5427–5430.
26. Z. Wu and H. Krakauer, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 014112.
27. Y. Okuno, K. Kawato, M. Suzuki, A. Harada and T. Oguchi, Jpn. J. Appl. Phys., 2007, 46, 5199.
28. S. Kim, W.-J. Lee, Y.-H. Cho, M. Shim and S. Kim, Jpn. J. Appl. Phys., 2013, 52, 091101.
29. P. Marton and C. Elsässer, Phys. Status Solidi B, 2011, 248, 2222–2228.
30. J. Niu, C. Li and Z. Xi, Phys. Chem. Chem. Phys., 2024, 26, 29457–29465.
31. J. Niu, C. Li and Z. Xi, RSC Adv., 2024, 14, 38245–38252.
32. J. S. Baker and D. R. Bowler, Phys. Rev. B, 2019, 100, 224305.
33. K. Eklund and A. J. Karttunen, J. Phys. Chem. C, 2023, 127, 21806–21815.
34. K. Eklund and A. J. Karttunen, J. Phys. Chem. C, 2024, 128, 16199–16207.
35. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Flament and B. Kirtman, J. Chem. Theory Comput., 2023, 19, 6891–6932.
36. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
37. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
38. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
39. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
40. M. Kuklin, K. Eklund, J. Linnera, A. Ropponen, N. Tolvanen and A. Karttunen, Molecules, 2022, 27, 874.
41. S. DallOlio, R. Dovesi and R. Resta, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 10105–10114.
42. Y. Noel, C. M. Zicovich-Wilson, B. Civalleri, P. DArco and R. Dovesi, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 65, 014111.
43. A. Erba, K. E. El-Kelany, M. Ferrero, I. Baraille and M. Rérat, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 035102.
44. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
45. A. Togo, J. Phys. Soc. Jpn., 2023, 92, 012001.
46. T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys.: Condens. Matter, 2014, 26, 225402.
47. T. Tadano and S. Tsuneyuki, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 054301.
48. F. Pascale, C. M. Zicovich-Wilson, F. López Gejo, B. Civalleri, R. Orlando and R. Dovesi, J. Comput. Chem., 2004, 25, 888–897.
49. C. M. Zicovich-Wilson, F. Pascale, C. Roetti, V. R. Saunders, R. Orlando and R. Dovesi, J. Comput. Chem., 2004, 25, 1873–1881.
50. L. Maschio, B. Kirtman, M. Rérat, R. Orlando and R. Dovesi, J. Chem. Phys., 2013, 139, 164102.
51. R. Peierls, Ann. Phys., 1929, 395, 1055–1101.
52. S.-i Tamura, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 27, 858–866.
53. A. Erba, J. Chem. Phys., 2014, 141, 124115.
54. A. Erba, M. Shahrokhi, R. Moradian and R. Dovesi, J. Chem. Phys., 2015, 142, 044114.
55. A. Erba, J. Maul, M. De La Pierre and R. Dovesi, J. Chem. Phys., 2015, 142, 204502.
56. J. D. Freire and R. S. Katiyar, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 2074–2085.
57. A. S. Bhalla and R. E. Newnham, Phys. Status Solidi A, 1980, 58, K19–K24.
58. R. W. Whatmore, Rep. Prog. Phys., 1986, 49, 1335–1386.
59. S. Vasala and M. Karppinen, Prog. Solid State Chem., 2015, 43, 1–36.
60. A. M. Glazer, P. A. Thomas, K. Z. Baba-Kishi, G. K. H. Pang and C. W. Tai, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 184123.
61. N. Zhang, H. Yokota, A. M. Glazer, Z. Ren, D. A. Keen, D. S. Keeble, P. A. Thomas and Z.-G. Ye, Nat. Commun., 2014, 5, 5231.
62. K. Hayashi, C. Lu, A. K. R. Ang, K. Ohwada, Y. Xie, W. Hu, T. Matsushita, N. Happo, S. Hosokawa, A. A. Bokov and Z.-G. Ye, Phys. Status Solidi B, 2020, 257, 2000191.
63. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
64. H. Morioka, S. Yokoyama, T. Oikawa, H. Funakubo and K. Saito, Appl. Phys. Lett., 2004, 85, 3516–3518.
65. D. Berlincourt, C. Cmolik and H. Jaffe, Proc. IRE, 1960, 48, 220–229.
66. M. A. Fraga, H. Furlan, R. S. Pessoa and M. Massi, Microsyst. Technol., 2014, 20, 9–21.
67. K. Eklund, J. Alajoki and A. J. Karttunen, Cryst. Growth Des., 2023, 23, 3427–3436.
68. M. W. Hooker, Properties of PZT-Based Piezoelectric Ceramics between −150 and 250 °C, National Aeronautics and Space Administration Technical Report NASA/: CR-1998-208708;, 1998.
69. E. Buixaderas, V. Bovtun, M. Kempa, D. Nuzhnyy, M. Savinov, P. Vanek, I. Gregora and B. Malic, Phys. Rev. B, 2016, 94, 054315.
70. T. Peters, C. Cheng, G. A. Rossetti and S. Trolier-McKinstry, J. Am. Ceram. Soc., 2022, 105, 4058–4070.
71. M. Tachibana, C. Bourgès and T. Mori, Appl. Phys. Express, 2023, 16, 101002.
72. M. Tachibana, T. Kolodiazhnyi and E. Takayama-Muromachi, Appl. Phys. Lett., 2008, 93, 092902.
73. J. More-Chevalier, P. V. Yudin, C. Cibert, P. Bednyakov, P. Fitl, J. Valenta, M. Novotný, M. Savinov, M. Poupon, T. Zikmund, G. Poullain and J. Lančok, J. Appl. Phys., 2019, 126, 214501.
74. M. Botea, L. Hrib, I. Pasuk, A. Iuga, L. Trupina, R. Negrea, N. Becherescu and L. Pintilie, Curr. Appl. Phys., 2019, 19, 804–810.
75. M. Es-Souni, M. Kuhnke, A. Piorra and C.-H. Solterbeck, J. Eur. Ceram. Soc., 2005, 25, 2499–2503.
76. E. Tang, Y. Wang, R. Wang, Y. Han, M. Chang, C. Chen, K. Guo and L. He, Mater. Chem. Phys., 2023, 307, 128109.
77. R. W. Whatmore and F. W. Ainger, Proc. Soc. Photo-Opt. Instrum. Eng., 1983, 395, 261–266.
78. J. L. Wentz and L. Z. Kennedy, J. Appl. Phys., 1964, 35, 1767–1770.
79. R. Masuki, T. Nomoto, R. Arita and T. Tadano, Phys. Rev. B, 2022, 106, 224104.

---

Footnote

† Electronic supplementary information (ESI) available: Additional computational details, Γ-point harmonic frequencies, phonon dispersion relations, elastic constants, and average thermal expansion coefficients for ordered PZT models, tabulated primary, secondary, and total pyroelectric coefficients, mode-specific contributions to pyroelectricity, example inputs of the ALAMODE ALM and ANPHON SCPH calculations, and optimized PZT structures both in CIF and CRYSTAL format. See DOI: https://doi.org/10.1039/d5cp01655j

---