## 2 Materials and Methods  
Ceramics were prepared from autocompustion sol-gel nanopowders synthesized at 500 \({}^{\circ}\)C, from acetate and nitrates, as explained elsewhere [18]. Pressed pellets were sintered in two steps using moderate maximum temperatures (1050 \({}^{\circ}\)C for 1-2 h) to avoid loss of volatiles. Thin sintered ceramic disks and plates were electroded at the major faces with silver paste, sintered at 400 \({}^{\circ}\)C, and poled in silicon oil bath from 120 \({}^{\circ}\)C to 50 \({}^{\circ}\)C, cooling with an applied field of 30 KV\(\cdot\)cm\({}^{-1}\), for 30 min. Shear plates were re-electroded for the electrical measurements. Measurements of complex impedance of disks and plates were made using precision impedance analyzers (Hewlett-Packard, Palo Alto, CA, USA). Model HP4194A was used for the measurements as a function of the temperature, until the inhibition of the piezoelectricity, in situ at an ad hoc sample-holder with a heating module. The sample holder and the temperature control were specially designed for this application. The 16 \({}^{\circ}\)C corresponds to room temperature, after that a 5 \({}^{\circ}\)C step was used starting from 35 \({}^{\circ}\)C up to 105 \({}^{\circ}\)C, since in the first interval (16-35 \({}^{\circ}\)C) the variation is irrelevant to the aim of the work. Model HP4192A-LF, connected to a computer with the appropriate software implemented [19], was used for the in situ iterative analysis of resonance curves. To obtain the complete parameters set, a minimization of the difference between the numerical and the experimental impedance curve is performed. The literature does not provide a complete matrix characterization of these ceramics that could be used as initial conditions for the numerical calculation of the impedance spectra of a thin disk. In order to determine such an initial condition, comprising 10 material coefficients including all losses--piezoelectric, dielectric, and elastic--the iterative analysis of complex impedance curves [12; 13; 19] of thin disks and shear plates, thickness-poled, was used. To determine d\({}_{33}\), a quasi-static measurement was carried out on a thin disk at 100 Hz using a Bernicourt piezometer (Channel products Inc., Chesterland, OH, USA). After this first condition, the optimization methodology presented in [6; 7] was used. Although there may be more efficient solutions, the methodology based in the Nelder-Mead algorithm [20] is generally robust and easy to implement, which are reasons why it was selected here. This methodology can determine all the complex coefficients in the linear piezoelectric model using only one sample. FEM simulations are performed using square with four nodes using complex numbers for all parameters in the model. The geometry is assumed axisymmetric and the material is assumed homogenous and isotropic. After performing a convergence analysis, the element size is fixed in 10 \(\upmu\)m, giving 40 elements for one wave length in the thickness mode frequency. Additional details about the implementation of the methodology and the FEM code are in the work of [PERSON] et al. [6; 21].This is suitable for this case when we are looking for the simultaneous evolution of the parameters as a function of the temperature including the coupled modes. There are different ways to express the constitutive equations, here we use the method adopted by the IEEE standard [5]. \[T_{p}=c_{pq}^{E}S_{q}-e_{kp}E_{k} \tag{1}\]  
\[D_{i}=e_{iq}S_{q}+\varepsilon_{ik}^{S}E_{k} \tag{2}\]  
Due to symmetry, only 10 independent parameters are used to represent a piezoelectric material with axial 6 mm macroscopic symmetry. All parameters are complex numbers in order to take in to account the energy losses. The determination involves five elastic constants (\(c_{11}^{E}\), \(c_{12}^{E}\), \(c_{13}^{E}\), \(c_{44}^{E}\), \(c_{33}^{E}\)), three piezoelectric (\(e_{31}\), \(e_{33}\), \(e_{15}\)), and two dielectric (\(\varepsilon_{11}^{S}\), \(\varepsilon_{33}^{S}\)).