Determination of the PIC700 Ceramic’s Complex Piezo-Dielectric and Elastic Matrices from Manageable Aspect Ratio Resonators

Achieving good piezoelectric properties, such as the widely reported d33 charge coefficient, is a good starting point in establishing the potential applicability of piezoceramics. However, piezoceramics are only completely characterized by consistent piezoelectric-elastic-dielectric material coefficient matrices in complex form, i.e., including all losses. These matrices, which define the various alternative forms of the constitutive equations of piezoelectricity, are required for reliable virtual prototyping in the design of new devices. To meet this need, ten precise and accurate piezoelectric dielectric and elastic coefficients of the material, including all losses, must be determined for each alternative. Due to the difficulties arising from the coupling of modes when using the resonance method, this complete set of parameters is scarcely reported. Bi0.5Na0.5TiO3-based solid solutions are already commercially available in Europe and Japan. Here, we report a case study of the determination of these sets of material coefficients (diα, giα, eiα and hiα; sE,Dαβ and cE,Dαβ; εTik and εSik; and βTik and βSik), including all losses, of the commercial PIC700 eco-piezoceramic. Plate, disk, and cylinder ceramic resonators of a manageable aspect ratio were used to obtain all the material coefficients. The validation procedure of the matrices is also given by FEA modeling of the considered resonators.


Introduction
Piezoceramics of Bi 0.5 Na 0.5 TiO 3 -based solid solutions are already commercially available in Europe under the denomination, PIC700 [1]. Moreover, there are commercial ceramics based on the same material in Japan [2]. Achieving good piezoelectric properties, such as the widely reported d 33 charge coefficient for sensors, is a good starting point in establishing their potential applicability. However, this knowledge is largely insufficient. Consistent piezoelectric-elastic-dielectric material coefficient matrices in complex form, i.e., including all losses, which define the various alternative forms of the constitutive equations of piezoelectricity, are required for reliable virtual prototyping in the design of new devices, which is a relatively inexpensive and time-saving procedure [3,4].
Finite Element Analysis (FEA) is widely used for device simulation. Nevertheless, its potential for modeling lead-free commercial piezoceramics is still insufficiently exploited. This is due to the lack of complete sets of data, including all losses, for many of these. The FEA modeling is often restricted to consider piezoceramics as lossless and even isotropic materials, which is just a rough approximation of their real performance. Piezoelectric ceramics have lossy properties, which are conveniently defined by complex values S ij = s D ijkl T kl + g kij D k (5) T ij = c D ijkl S kl − h kij D k (7) where e kij , d kij , g kij , and h kij are the piezoelectric coefficient matrices; s E,D ijkl and c E,D ijkl are the elastic compliances and elastic stiffness tensors, respectively, at a constant electric field (E), closed circuit, or constant electric displacement (D), open circuit; ε T ik and ε S ik , β T ik , and β S ik are the free (at constant, zero, stress (T)) and clamped (at constant, zero, strain (S)) relative dielectric permittivity matrices and dielectric impermeabilities or impermittivities (β = 1/ε), respectively. These equations are currently being simplified using Einstein s summation convention of repeated subscripts (i, k = 1, 2, 3; ij = α; and kl = β, where α, β = 1, 2, . . . , 6) and the subscript 3 to denote the polarization direction in the material. The mentioned matrices' inconsistency derives from the use of resonator modes at a relatively wide range of frequencies and whose microstructural characteristics or poling levels may not be identical. For this reason, it is a good practice to provide one or more validity and consistency criteria of the set of parameters determined [8].
To provide a better consistency between matrices, it was shown previously that all material coefficients can be determined from only 3 resonators and four modes of resonance [3,10]. Measurements of the radial and thickness extensional modes of thicknesspoled thin disk can be used to obtain those parameters derived from the Length Thickness Extensional (LTE) mode of a bar and the Thickness Extensional (LTE) mode of a thin plate, which is thickness-poled. The full set of matrices is completed using well-known relationships between coefficients [3,10]. However, the thin plate for shear coefficient determination (e.g., d 15 , s E 55 , and ε S 11 ) from the electrically-induced shear mode cannot be excluded from this reduced resonator set. Two types of difficulties arise when using the standard thin shear plate, which is length-poled and excited in thickness. They cause the associated parameters to be scarcely reported and affect the consistency of the matrices. On the one hand, longitudinal poling is a problem when dealing with materials of a high coercive field [11] or low dielectric breakdown [12]. On the other hand, the coupling of shear modes and other lateral and contour modes, which takes place in the shear resonance of standard length-poled shear plates, is unavoidable [3,4,13,14]. Unfortunately, this leads to the stringent requirement of aspect ratios between the length for poling (L), the width (w), and the thickness for electrical excitation (t) of the plate of L, w > 20 t [15] and to an underestimation of the piezoelectric coefficients [13].
To create an improved approach to the determination of the properties using the resonance method, we present, in this work, the determination of piezoceramic coefficients by the measurement of four resonance modes in three poled ceramic resonators, including a, non-standard, thickness-poled and shear plate [3]. An analysis of the complex impedance spectra was carried out using a widely tested iterative automatic method that overcomes the difficulties and limitations of the current standard procedures of measurement [7,8]. Noticeably, this method has been implemented in the analysis of all resonances under consideration [3,[16][17][18]. Besides, the method is self-sustainable, as it does not need an initial estimate of the parameters. Some issues to consider when attempting to obtain weakly coupled resonances and reliable coefficients are discussed. Commercial PIC700 ceramics were studied, and the piezoelectric-dielectric-elastic matrices of material coefficients in complex form were determined for the alternative types of constitutive sets of piezoelectricity (Equations (1)- (8)). Besides, the validity of the matrices was tested by using them in Finite Element Analysis modeling for the reproduction of the experimental measurements.

Materials
Bi 0.5 Na 0.5 TiO 3 -based commercial lead-free (PIC700; PI Ceramic GmbH, Lederhose, Germany) poled ceramic resonators were studied in this work [19]. The ceramics were prepared from high-purity oxide precursors by the conventional ceramic method. Samples were cut and lapped to the desired dimensions. Afterwards, they were electroded and poled in an oil bath, according to the company standards [19], whose details are protected knowledge. All the resonators under study were poled to saturation in the same conditions to ensure the consistency of the results obtained from them. All measurements were taken for at least two samples-although statistical analysis is not an issue discussed in this work-to verify that the results are not affected by microstructure accidents. Changes from batch to batch could take place. The density of the considered samples was 5.49 g·cm −3 .

Material Coefficient Determination
Piezoelectric, elastic, and dielectric complex coefficients in the linear range (d iα , g iα , e iα , and h iα ; s E,D αβ and c E,D αβ ; ε T ik and ε S ik ; and β T ik and β S ik ), accounting for all losses, were determined at room temperature from measurements of complex impedance curves. The data were automatically acquired from an HP-4192A LF precision impedance analyzer (Hewlett-Packard, Palo Alto, CA, USA) and controlled by a PC via a GPIB-PCIIA (National Instruments). An automatic iterative method explained elsewhere [3,16,17] was here used for the analysis of the curves for each resonance considered. Contrarily to standard methods, this method does not require additional measurements, besides those of the resonance. Four complex impedance spectra, Z*, of the electromechanical resonance modes of three poled ceramic resonators were analyzed. These modes are: (i) Shear mode of non-standard, thickness-poled, and longitudinally excited shear plates, (ii) Radial and (iii) Thickness extensional modes of thickness-poled thin disks; and (iv) Length extensional mode of a longitudinally poled cylinder. The software used is available to the scientific community [20].
Commonly, Z* is represented as the modulus (Z) and phase angle (θ) vs. frequency plots ( Figure 1a). Here, we also use an equivalent representation: the peak of the real part of Z*, the resistance (R), and the peak of the conductance (G), the real part of the complex admittance, which is the inverse of the impedance (Y* = 1/Z*). This is a convenient plot for using the automatic iterative method (Figure 1b). After plotting these peaks, the frequencies of the maximum G and R values, f s and f p , respectively, are automatically determined by the software. For each of the considered modes, some material coefficients were directly determined by the automatic iterative solution of the impedance/admittance expression as a function of these coefficients and the resonator density (ρ) and dimensions (thickness between electrodes for electrical excitation (t); electroded surface area (S)) [3,16,17]. In each iteration of the analysis, the software solves a system of four equations for the values of Z* or Y* measured at four frequencies, including f s and f p , by a numerical method, until a convergence criterion of the last determined coefficients in comparison with the previous ones is met.
Besides, using well-known expressions among coefficients of the constitutive equations [3,10], several other coefficients related to the studied mode were also determined after the directly calculated ones. These are the electromechanical coupling factors (k x , where x = 33, 31, 15, p (radial mode of the disc) and t (thickness mode of the disc)) and the corresponding frequency numbers (N x = f x s (kHz) · D(mm), where f x s is the series resonance of the x mode, and D is the leading dimension of the resonance), as well as the planar Poisson's ratio.
After the directly obtained coefficients are available, a reconstruction of the spectra (Figure 1b) is carried out. Calculation of the complex impedance as a function of the frequency is conducted using the analytical expression that was just solved and the calculated coefficients. The residuals for these Rand G-reconstructed peaks (a residual being the difference between the experimental value for a given frequency and the value of the reconstructed peaks at the same frequency) give us a regression factor for the iterative analysis ( 2 ). This parameter accounts for the validity of the analytical expression and these coefficients for 1D modeling of this mode of resonance. The higher the coupling of the resonance with other undesired resonances, the lower the 2 . The closer the experimental resonance to a monomodal resonance, the closer 2 is to 1.

Finite Element Modelling
The FEA modeling software (COMSOL Inc., Burlington, MA 01803, USA, COMSOL Multiphysics ® 4.3) was used. The piezoelectric device module (PZD) was used to simulate the piezoelectric response of the four resonance modes of the three resonators considered here. The mesh used typically has five nodes per wavelength, which allowed for simulating the modes of resonance that are excited, together with the fundamental resonance under study. The number of frequencies analyzed was chosen to obtain the needed resolution of each calculated spectra. The shear plate was simulated as a full 3D item. This typically results in a calculation time of up to 5 h for each sweep of 1000 frequencies. The cylinder and disk resonators were simulated using their rotational symmetry, which results in faster calculations. The 3D harmonic analysis used here provides the complex impedance values in a given interval of frequencies. The intervals of the frequency of interest in the present work were those of the experimental measurements. Thin PIC700 ceramic plates of a ratio between the length for electrical excitation to thickness for poling L/t = 3.75 were initially prepared with t = 2.00 mm. Figure 1a shows the experimental spectra in (Z, θ) plots for this resonator. In addition to the fundamental shear mode associated to the thickness of the plate, which takes place at f s = 797.7 kHz and is marked as 3, four other modes are measured. The resonances marked as 1 and 2 correspond to modes taking place at lower frequencies, which are associated with the larger lateral dimensions of the plate. Their overtones appear at periodically distributed frequencies. Those marked with 4 and 5 are coupled with the fundamental shear resonance, which is marked 3. This invalidates the sample for the accurate determination of complex material coefficients.
The major surfaces of the non-standard plate are not electroded for the electrical excitation of the resonance and the impedance measurement. It is therefore feasible to reduce the thickness, which increases the L/t ratio and the frequency of the fundamental resonance by polishing the major surfaces in steps of approximately 0.05 mm. Measurements were taken after each step for the control of the coupling. As the lateral dimensions remain unchanged, the frequencies marked 1 and 2 for the lateral modes are unaffected. The coupling diminishes as the shear frequency moves away from these frequencies to higher frequencies. Noticeably, these modes are weakened as their frequency is separated from the frequency for the fundamental shear. The calculations carried out in each step are not strictly valid for the accurate determination of the properties [21], but they allow for the quantitative evaluation of the coupling, using the regression factor for the iterative analysis (ℜ 2 ). Figure 1a also shows the resonance spectra, measured again after the thickness of the plate was reduced in steps to a plate with t = 1.66 mm (L/t = 4.52), which has a virtually uncoupled fundamental shear mode taking place at fs = 892.0 kHz. 3.0x10 -4 The major surfaces of the non-standard plate are not electroded for the electrical excitation of the resonance and the impedance measurement. It is therefore feasible to reduce the thickness, which increases the L/t ratio and the frequency of the fundamental resonance by polishing the major surfaces in steps of approximately 0.05 mm. Measurements were taken after each step for the control of the coupling. As the lateral dimensions remain unchanged, the frequencies marked 1 and 2 for the lateral modes are unaffected. The coupling diminishes as the shear frequency moves away from these frequencies to higher frequencies. Noticeably, these modes are weakened as their frequency is separated from the frequency for the fundamental shear. The calculations carried out in each step are not strictly valid for the accurate determination of the properties [21], but they allow for the quantitative evaluation of the coupling, using the regression factor for the iterative analysis ( 2 ). Figure 1a also shows the resonance spectra, measured again after the thickness of the plate was reduced in steps to a plate with t = 1.66 mm (L/t = 4.52), which has a virtually uncoupled fundamental shear mode taking place at f s = 892.0 kHz.  Figure 1b shows also the data for a plate with L/t = 3.89, after a further thickness reduction to t = 1.61 mm, which again shows a coupled shear mode due to interference, with modes 4 and 5 marked in the spectra of the t = 1.66 mm sample in Figure 1a. For this sample, 2 is reduced to 0.9560. This indicates that the maximum decoupling of modes was surpassed, and the resonators evolve away from a monomodal resonance once again.
The described procedure does not follow the standard procedure for preventing coupling by the separation of frequencies of the resonance modes of the resonator, associated with each one of its dimensions, by measurement in high-aspect-ratio resonators. It is assumed that, in this way, the overtones of low-frequency modes become weaker at higher frequencies, and the same occurs in the case of coupling with the high frequency resonance of interest.
Nevertheless, previous work on thickness-poled shear plates has shown that undesirable overtones of lateral modes are activated by coupling when the shear mode frequency is in coincidence with their frequency. This occurs periodically. This periodical sequence of coupling-decoupling processes of the non-standard shear resonator was observed for a soft lead titanate-zirconate (PZT; four periods studied from L/t = 7.5 to 15) [3] and for a barium calcium titanate-zirconate (BCZT; five periods studied from L/t = 6 to 12) [19]. Besides, it was modeled by FEA for bismuth sodium barium titanate (BNBT; three periods modeled from L/t = 7.5 to 14.5) [22].
According to the results shown in Figure 1, the plate used here for the coefficient determination was the t = 1.66 mm one. This was efficiently decoupled from lateral modes, despite its low aspect ratio, compared with the standard recommended procedure. Besides, this was achieved using a method that is not feasible for the standard resonator, because it requires multiple steps for the removal of electrodes and re-electroding at room temperature.
The directly calculated coefficients in this mode are ε S 11 , e 15 , and s E 55 , and, from these,

Resonances of a Thickness-Poled and Excited Thin Disk
(a) Extensional radial resonance of the disk This is the only resonance mode for which it is necessary to conduct additional measurement outside the fundamental resonance spectrum to apply the iterative method. The series frequency of the first overtone of this mode (f 2s ) is used to calculate the first estimate of the planar Poissón s ratio (σ P ) from the ratio, f 2s /f s [17].
Radial resonance-derived parameters are the most commonly reported, since it is the mode of the lower frequency of the disk, and it is not affected by coupling. Figure 2 shows the monomodal impedance spectrum of the radial mode as an (R, G) plot of a thicknesspoled and excited thin disk with a diameter of 12 mm and t = 1.15 mm. The experimental and reconstructed spectra are shown. The coefficients from the iterative method lead to a a high regression factor 2 = 0.9999.
The directly calculated coefficients in this mode are ε T 33 , d 31 , and s E 11 , which are the coefficients that can otherwise be obtained by the LTE of a thickness-poled bar. The coefficient s E 12 is also determined in the radial mode from the Poisson s ratio and given the expression σ P = (−s 12 E ⁄s 11 E ). Besides, g 31 , s D 11 , and s D 12 are also calculated from the previous values, along with s D,E 66 = 2(s D,E 11 − s D,E 12 ) and c D,E 66 = 1/s D,E 66 . Due to the symmetry d 32 = d 31, g 32 = g 31, s 21 = s 12 , and s 22 = s 11 are determined. Additionally, the electromechanical coupling factors and frequency numbers, k p , N p , and k 31 , were obtained. These are shown in Tables 1-3.

(b) Extensional thickness resonance of the disk
In comparison with the relatively easy to handle monomodal resonances of bars and plates, great care must be taken with the disc resonator and its dimensional ratio to avoid complex high-frequency resonance modes. These are the overtones of the planar mode or the contour modes, which are coupled with the high-frequency fundamental thickness mode of the disk. For this purpose, the recommended standard procedure is to have a diameter to thickness ratio (D/t) > 10. However, it is also widely reported that coupling takes place also under this condition [9,[23][24][25]. The directly calculated coefficients in this mode are ε T 33, d31, and s E 11, which are the coefficients that can otherwise be obtained by the LTE of a thickness-poled bar. The coefficient s E 12 is also determined in the radial mode from the Poisson´s ratio and given the expression σ P = (−s12 E ⁄s11 E ). Besides, g31, s D 11, and s D 12 are also calculated from the previous values, along with s D,E 66 = 2(s D,E 11 − s D, E 12) and c D,E 66 = 1/s D,E 66. Due to the symmetry d32 = d31, g32 = g31, s21 = s12, and s22 = s11 are determined. Additionally, the electromechanical coupling factors and frequency numbers, kp, Np, and k31, were obtained. These are shown in Tables 1-3.

(b) Extensional thickness resonance of the disk
In comparison with the relatively easy to handle monomodal resonances of bars and plates, great care must be taken with the disc resonator and its dimensional ratio to avoid complex high-frequency resonance modes. These are the overtones of the planar mode or the contour modes, which are coupled with the high-frequency fundamental thickness mode of the disk. For this purpose, the recommended standard procedure is to have a diameter to thickness ratio (D/t) > 10. However, it is also widely reported that coupling takes place also under this condition [9,[23][24][25]. Figure 3 shows the impedance curve in an (R, G) plot for the thickness extensional mode used for the evaluation of the parameters in this mode. A distorted spectrum is shown, which is the result of coupling. This deficiency would pose some problems with the accuracy of the material coefficients obtained, particularly the losses, as reflected by the regression factor of the calculation: ℜ 2 = 0.8720. Undoubtedly, research must still be conducted with the aim of enhancing this experimental result, which compromises the overall determination of matrices. However, this type of spectra is what is commonly accepted as a manageable aspect ratio and impedance curve for analysis when working with the thickness modes of disks [9,[23][24][25].  Figure 2. Fundamental radial extensional resonance mode of a thickness-poled and excited thin disk. Impedance spectra in an (R, G) plot of a disk with a diameter of 12 mm and t = 1.15 mm. The symbols are the experimental data, and the lines are the reconstructed peaks after the calculation of the coefficients. The value of the regression factor for this calculation ( 2 ) is shown. The black symbols and lines (G data) can be read in the left-y-axis, and the blue ones (R data) are in the right-y-axis. Figure 3 shows the impedance curve in an (R, G) plot for the thickness extensional mode used for the evaluation of the parameters in this mode. A distorted spectrum is shown, which is the result of coupling. This deficiency would pose some problems with the accuracy of the material coefficients obtained, particularly the losses, as reflected by the regression factor of the calculation: 2 = 0.8720. Undoubtedly, research must still be conducted with the aim of enhancing this experimental result, which compromises the overall determination of matrices. However, this type of spectra is what is commonly accepted as a manageable aspect ratio and impedance curve for analysis when working with the thickness modes of disks [9,[23][24][25].  The directly calculated coefficients in this mode are ε S 33, h33, and c D 33. They can also be obtained from the LE mode of a thickness-poled thin plate. From these coefficients, e33 and c E 33 are also calculated. At this stage of the process, the ε S ij matrix is complete, which allows for the calculation of the β S ij matrix. Additionally, the electromechanical coupling factor and frequency number, kt and Nt, are obtained. These are shown in Tables 1-3 Figure 3. Fundamental thickness extensional resonance mode of a thin disk. The impedance spectra in an (R, G) plot of a disk with a 12 mm diameter and 1.00 mm thickness. The symbols are the experimental data, and the lines are the reconstructed peaks after the calculation of the coefficients. The value of the regression factor for this calculation ( 2 ) is shown. The black symbols and lines (G data) can be read in the left-y-axis, and the blue ones (R data) are in the right-y-axis.
The directly calculated coefficients in this mode are ε S 33 , h 33 , and c D 33 . They can also be obtained from the LE mode of a thickness-poled thin plate. From these coefficients, e 33 and c E 33 are also calculated. At this stage of the process, the ε S ij matrix is complete, which allows for the calculation of the β S ij matrix. Additionally, the electromechanical coupling factor and frequency number, k t and N t , are obtained. These are shown in Tables 1-3. Figure 4 shows the impedance curve in an (R, G) plot for the length mode used for the evaluation of the parameters. As it corresponds to the mode of the resonator with the lowest frequency, this is uncoupled, and the regression factor from the calculation is high ( 2 = 0.9998). The directly calculated coefficients in this mode are ε T 33 , g 33 , and s D 33 . Due to the frequency dependence of the permittivity, the ε T 33 determined in this mode is higher than the one determined in the thickness mode taking place at a higher frequency ( Figure 3). Consequently, at this point, there is a need for an expert selection of the value that is more amenable for the validity of the matrices for FEA modelling [10].

Combined Determination of the Remaining Material Coefficients
The combination of the coefficients obtained from the radial and thickness extensional resonances of a thin disk and the length extensional resonance of the cylinder, together with the well-known relationships among these coefficients [3,7,8,10], also gives the following: (a) from the inversion of the matrices, we have the following expression: and, hence, we can obtain s E 13. From s E 13, the symmetry considerations provide s E 23 = s E 13, s E 31 = s E 13, and s E 32 = s E 13. This completes the s E αβ matrix and allows the calculation of its inverse, c E αβ, which at this point was only lacking the values of c E 11, c E 12, and c E 13, to be completed. This also ensures consistency between the stiffness and compliance matrices.
(b) knowing c E 11, c E 12, and c E 13, we can make use of Equation (62) in [8] to obtain e31 from the following: This completes the eiα matrix.
(c) by the relationships between the coefficients, we can calculate h31 using the following expression:  At this stage of the process, with the coefficients determined for this resonance mode, the g iα and ε T ij matrices are complete, which allows the β T ij matrix to be calculated. From these coefficients, d 33, which completes the d iα matrix as well, and s E 33 are also calculated. Due to the frequency dependence of the material coefficients, this d 33 value at resonance is lower than the quasi-static one measured at 100 Hz. Additionally, the electromechanical coupling factor and frequency number, k t and N t , are obtained. These are shown in Tables 1-3.

Combined Determination of the Remaining Material Coefficients
The combination of the coefficients obtained from the radial and thickness extensional resonances of a thin disk and the length extensional resonance of the cylinder, together with the well-known relationships among these coefficients [3,7,8,10], also gives the following: (a) from the inversion of the matrices, we have the following expression:  [8] to obtain e 31 from the following: This completes the e iα matrix.
(c) by the relationships between the coefficients, we can calculate h 31 using the following expression: This completes the h iα matrix.
(d) making use of Equation (29) in [8], we can obtain s D 13 from the following: (12) and, similarly, by the symmetry considerations, we can also determine s D 23 = s D 13 , s D 31 = s D 13 , and s D 32 = s D 13 . This completes the s D αβ matrix and allows the calculation of its inverse, c D αβ , which at this point was lacking the values of c D 11 , c D 12 , and c D 13 , to be completed. This completes all the matrix determinations, as shown in Tables 1-3, and determines the consistency between the c D αβ and s D αβ matrices. The PIC700 eco-piezoceramic has been considered a good replacement for PIC 255 (PI Ceramic GmbH, Lederhose, Germany) material [1]. PIC255 is a soft-type modified lead zirconate-titanate, which is characterized by a moderate permittivity, high thickness coupling factor, high charge coefficient, and low mechanical quality factor (ratio P /P" for elastic coefficients) [19]. Applications of this ceramic type are in low-power ultrasonic transducers, non-resonant broadband systems, and force and acoustic sensors. PIC700 is additionally characterized by a lower planar coupling factor, higher electromechanical anisotropy [26], and lower density (Table 3).
When standard procedures are used, this is not a sufficiently precise method for determining a physically meaningful and mathematically coherent set of properties. To obtain a coherent set of matrices suitable for FEA analysis, these results must be optimized, using a suitable algorithm to minimize the overall error.
Here, we have used a validity check at each considered resonance mode through the reconstruction of the measured spectrum and considering only results with high enough 2 values. Additionally, we have used a, non-standard, thickness-poled and longitudinally excited shear plate producing accurate complex shear coefficients, as the shear mode can be efficiently decoupled from other modes. The procedure reported here for the calculation of the coefficients by combining parameters that were obtained from previously analyzed spectra is not a universal solution, as there are options for calculating the indirectly determined coefficients [8,10].

Validation of the Piezoelectric, Elastic, and Dielectric Material Coefficients
The validity of the coefficients shown in Tables 1-3 will be evaluated in the following.

Meaningful Losses
Holland's criteria [26] for the material to be passive, in the sense of not being an energy generator, are given for the matrices of complex coefficients by the following expressions: These criteria were checked, and they are fulfilled by the set of data in Tables 1-3, which provide the first consistency criteria of the accomplished material coefficients. Figure 5a shows the FEA modelled (/Z/, θ) spectra, including the lateral and the fundamental shear modes, for the coupled and decoupled resonances of the nonstandard shear plates of a thickness t = 2.00 mm and t = 1.66 mm, respectively, whose experimental spectra are shown in Figure 1a. Both were calculated using the values of the complex coefficients shown in Tables 1-3. Figure 5b shows the (R, G) spectra for the plate used for the material coefficient determination (t = 1.66 mm). The height and frequency, where the main resonances are positioned in the FEA modeled spectra, are in good agreement with their respective experimental data, marked as 3 for the fundamental shear mode. It is worth noting that the agreement is also good for the secondary modes, marked as 1, 2, 4, and 5. The calculated parameters show the validity of the reproduction of all details of the spectra. Figures 6-8 show the FEA modeled spectra in (R, G) plots in the fundamental radial and thickness extensional modes of the thin disk and length extensional mode of the cylinder. The monomodal resonances are well modeled (Figures 6 and 8). Figure 7 shows a good reproduction of the positions for the minimum and maximum impedance modulus (Figure 7a) and the maxima of the main R and G peaks (Figure 7b) of the coupled resonance. Besides, all secondary features of the experimental spectrum in the distorted thickness extensional mode are also reproduced in the FEA-modeled curves.

Finite Element Analysis
The overall comparison between experimental and 3D FEA modeled spectra shown in Figures 5-8 is a positive criterion for the validity of the complex material coefficients shown in Tables 1-3 in a frequency interval from 100 kHz to at least 3 MHz. main resonances are positioned in the FEA modeled spectra, are in good agreement with their respective experimental data, marked as 3 for the fundamental shear mode. It is worth noting that the agreement is also good for the secondary modes, marked as 1, 2, 4, and 5. The calculated parameters show the validity of the reproduction of all details of the spectra.     worth noting that the agreement is also good for the secondary modes, marked as 1, 2, 4, and 5. The calculated parameters show the validity of the reproduction of all details of the spectra.     Figure 6. Fundamental radial extensional resonance mode of a thickness-poled and excited thin disk. The FEA-modeled (thick lines) impedance spectrum in an (R, G) plot of a disk with a diameter of 12 mm and a thickness of 1.15 mm. The experimental data (symbols) and the reconstructed spectrum after the calculation of the parameters in a 1D model (thin lines) are also shown for comparison. Black symbols and lines (G data) can be read in the left-y-axis, and the blue ones (R data) are in the right-y-axis. Figure 6. Fundamental radial extensional resonance mode of a thickness-poled and excited thin disk. The FEA-modeled (thick lines) impedance spectrum in an (R, G) plot of a disk with a diameter of 12 mm and a thickness of 1.15 mm. The experimental data (symbols) and the reconstructed spectrum after the calculation of the parameters in a 1D model (thin lines) are also shown for comparison. Black symbols and lines (G data) can be read in the left-y-axis, and the blue ones (R data) are in the right-y-axis.   are in the right-y-axis.   FEA-modeled impedance spectrum in an (R, G) plot of a cylinder with a diameter of 6 mm and a length of 15 mm (shown as thick lines). The experimental and reconstructed spectra after the calculation of the parameters (1D model) are also shown as symbols and thin lines, respectively, for comparison. Black symbols and lines (G data) can be read in the left-y-axis, and blue ones (R data) are in the right-y-axis.

Conclusions
A non-standard strategy for obtaining an uncoupled shear resonance mode of a thickness-poled and longitudinally excited shear plate of PIC700 eco-piezoceramic was described here. It allows for the use of an analytical expression of the resonance, which is strictly valid for a monomodal resonance, on samples with an aspect ratio as low as L/t = 4.52, whereas the standard measurement procedure requires aspect ratios higher than 20. This strategy led to an accurate determination of all related shear material complex coefficients from a manageable aspect ratio plate, along with the material losses.
All piezo-dielectric and elastic matrices of the PIC700 piezoceramic for the alternative forms of the sets of constitutive equations of piezoelectricity were obtained from additional measurements of the disks and cylinders. An analysis of the impedance curves was made by the automatic iterative method. A quantitative evaluation of the validity of the complex coefficients, obtained at each resonance using the regression factor of the recalculated spectra ( 2 ), was accomplished. The selection of spectra with a high 2 allows for the use of resonators of a manageable aspect ratio. The matrices of the coefficients are consistent and show the validity of PIC700 material modeling using an FEA in a frequency interval from 100 kHz to a few MHz.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.