English

• Piezoelectric materials have garnered extensively focus owing to their remarkable ability to facilitate bidirectional conversion between mechanical and electrical energy rooted in their non-centrosymmetric structure [1-3]. This dynamic interplay of energy conversion makes them indispensable in a diverse range of modern industrial applications [4-8]. Ferroelectrics, especially for perovskite-based ceramics and crystals, have attracted considerable interest for their splendid piezoelectric properties [9-12]. Through the formation of phase boundaries and introduction of relaxor behavior, further enhancement of piezoelectric response would be realized [13-15]. Owing to the existence of several equivalent polarization axis, polydomain state can universally be found in ferroelectric materials, which cannot be completely textured into a single domain in single crystals through complex poling procedure. Therefore, extrinsic contribution originated from micro-scale domain wall motion plays an ignorable role for the piezoelectricity, especially for relaxor ferroelectrics in which high-density low-symmetry domain walls separate the polar nanoregions, this increases difficulty for the deconvolution of intrinsic piezoelectric response from each functional unit and atom [16-20]. Moreover, the depolarization and polymorphic phase transition on heating not only limits the application temperature range (usually <400 ℃ for perovskites) far lower than their melting point but also further complicates the contributions and mechanism for piezoelectricity and thermal stability [21,22].

  The urgent need for high-temperature piezoelectric materials has spurred extensive research. Fig. 1a presents several prominent candidates along with their approximate maximum operational temperatures. Below 1000 ℃, ceramics are predominantly favored due to their cost-effectiveness in preparation. Among these, Aurivillius-phase ferroelectric ceramics, which comprise a combination of perovskite and bismuth-oxygen layers, are particularly dominant owing to the tunability of their components and layer structures (e.g., though the high-entropy design, the superior performance acquired) [23]. However, at temperatures exceeding 1000 ℃, with the exception of a few systems, such as the perovskite-like layer-structured (PLS) [24], phase transitions of depolarization and a reduction in resistivity (falling below 10⁶ Ω cm at 700 ℃ in the Aurivillius structure) occur, rendering these materials less viable. As a result, non-ferroelectric materials have become more suitable for high-temperature applications. Currently the mainstream non-ferroelectric piezoelectrics could be divided into several systems, that Gallium orthophosphate (GaPO₄) [25,26], Ca₂Al₂SiO₇ (CAS) [27], aluminum nitride (AlN) [25,28], Tourmaline [29], Ba₂TiSi₂O₈ (BTS) [30,31], langasite family crystals [25,32-34] and rare-earth calcium oxyborate crystals (ReCOB) [25,35-38]. Among them, ReCa₄O(BO₃)₃ (ReCOB, with Re: Y, Er, Nd, Gd, Pr and La in this work) crystals possess crucial status due to them not only remain the onefold steady structure and monotonic variation of piezoelectricity before melting point (exceeding ~1400 ℃) but also equip remarkable piezoelectric performance combined with high electrical resistivity (10⁷–10⁹ Ω cm at 700 ℃), making them irreplaceable in the field of ultra-high temperature vibration sensors [39]. Most importantly, ReCOB crystals are non-ferroelectric materials, the absence of multiple domain structure in ferroelectrics with different possible symmetries largely limits the application of design strategies, such as phase boundary design and domain engineering, for improving piezoelectricity. Every coin has two sides, the lack of external influence in the structure provides the foundation for quantitatively calculating piezoelectric contributions and facilitates the direct comparison between theoretical models and experimental results. In contrast to ferroelectrics, where the polarization direction can change under the influence of an external electric field even in polycrystalline ceramic, the single crystals or texture ceramics are needed for non-ferroelectric materials to maintain the macroscopic non-centrosymmetric symmetry. The high cost for crystal preparation and absence of guiding theory makes it very difficult to develop a new generation of non-ferroelectric crystals. Therefore, new theoretical analytical methods need to be developed to understand the intrinsic piezoelectric origin.

  Figure 1

  

  Figure 1.  (a) Relationship between piezoelectric coefficient (d₃₃ or d₂₆) and operating temperature for various ceramics and single crystals. (b) The structure of ReCOB and its three octahedral piezoelectric subunits, the orange arrows reveal the polarization vectors. (c) The iDPC-STEM image of NdCOB perpendicular to z zone axis.

  DownLoad: Full-Size Img PowerPoint

  Recent advancements in deep-learning have ushered in a transformative paradigm for simulating and understanding the dynamic behavior of materials [40,41]. The integration of deep-learning potentials (DP) with molecular dynamics (MD) simulations represents a cutting-edge methodology that enables the precise prediction of atomic-level interactions [42]. By training these models using atomic coordinates derived from Density Functional Theory (DFT) calculations, the neural network directly targets the corresponding system energy, providing a highly efficient and accurate approach to energy computation [43]. The training process of these DP models virtually mirrors the establishment of the effective Hamiltonian and empirical potentials commonly employed in traditional MD simulations, but from the deep-learning perspective that accounts for the inherent non-linearity and high-dimensionality of the interatomic interactions rather than relying on certain combinations of empirical formulas. Besides, one of the critical strengths of DP-based MD simulations is their ability to handle perturbations in the large-scale system, enabling the examination of localized structural variations and their subsequent impact on macroscopic material properties, thereby particularly powerful in investigating phenomena such as phase transitions, ferroelectric/piezoelectric behaviors and topological phases [44-46].

  Even though algorithms have developed rapidly in recent years, it still lacks of a reliable model to calculated the macroscopic piezo-response directly from chemical composition. On the other hand, ferroelectrics with polydomain and ceramics with polycrystal exhibits large extrinsic contributions, which affects the piezo-response at microscale making them to be the inappropriate candidates before building a mature model. Therefore, based on the above discussion, it can be found that the ReCOB crystal is an ideal system to make a pioneering exploration using an advanced algorithm for the sub-cell quantitative decomposition for both piezoelectricity and its thermal stability. As a result, in this work, an algorithm model based on DP is built to excavate the universal internal mechanisms driving the piezoelectricity, which is further proved to be reliable by comparing with experimental results. The new method would not only guidance the discovery of new ReCOB ultra-high temperature crystal with higher piezoelectricity to meet the requirements of next-generation device applications but also be used to give a deeper understand for the intrinsic contributions in designing higher piezoelectricity materials.

  The monoclinic ReCOB framework (Cm space group) consists of alternating lanthanide (Re) and calcium (Ca) octahedra embedded within a borate-oxygen network (Fig. 1b). Each Re³⁺ ion occupies a highly asymmetric octahedral site, coordinated by six oxygen atoms. The Ca²⁺ ions are classified into two categories: One category is similar to the Re³⁺ ions that possesses enormous distortion, while the other resides in a comparatively regular octahedral configuration. To describe the octahedral geometry more effectively, displacement vectors are used to quantify the degree of distortion. These vectors are defined as the difference between the central cation coordinate and the average position of its adjacent coordinating oxygen atoms:

  Δr=rcation−1n∑i=1nrOi

  (1)

  Table 1 presents the corresponding displacement vectors for various piezo-units in the NdCOB crystal, aligning with the intuitionistic octahedral morphology illustrated in Fig. 1b. Moreover, variations in the lanthanide elements directly influence the nominal radii of the substituted Re³⁺ cations, in Table S2, it demonstrates a high positive degree of correlation between them. Atomic-resolution aberration-corrected integrated differential phase contrast (iDPC) scanning transmission electron microscopy (STEM) was utilized to examine the localized structure, as depicted in Fig. 1c, for NdCOB crystals along the [010] zone axis. The consistency between the ideal structure and the iDPC-STEM mapping ensures the authenticity of the results, with the lateral view providing a clear depiction of the ideal displacement distribution in the two-dimensional plane.

  Table 1

  

  Table 1.  The piezoelectric subunit information of NdCOB. For each element, ionic displacement is defined as the difference between the center cationic coordinate and the average value of nearest neighbor coordination oxygen atoms.

  DownLoad: CSV

  Cation	Coordination number	X axis (Å)	Y axis (Å)	Z axis (Å)	Average ionic displacement (Å)
  Nd	6	−0.3102	0.0000	−0.1401	0.3403
  Ca1	6	0.0562	0.1058	0.0063	0.1200
  Ca1	6	0.0562	−0.1058	0.0063	0.1200
  Ca2	6	−0.3244	−0.1785	−0.2240	0.4328
  Ca2	6	−0.3244	0.1785	−0.2240	0.4328
  B	3	0.0067	0.0000	0.0002	0.0067
  B	3	0.0018	−0.0013	0.0063	0.0066
  B	3	0.0018	0.0013	0.0063	0.0066

  The piezoelectric coefficient d_ij, a fundamental parameter in piezoelectric materials quantifies the relationship between the generated charge induced in the material versus the corresponding applied perturbations, which could be explained as the variation in polarization to external stress [1-3]. The formula can be written as

  dαj=∂Pα∂σj|E=ΔPαΔσj|E

  (2)

  In which ΔPα present the variation of polarization along α direction (α=1,2,3), which be defined as the ionic displacement Δrα multiples the Born effective charge Zαα* along certain direction, in formula as

  Pα=Zαα*·Δrα

  (3)

  and Δσj is the variation of stress (j=1,…,6). For instance, the shear piezoelectric coefficient d₂₆ is specifically defined as the slope of the function that links the polarization along the y-axis to the shear stress applied in the xy-plane, which can be displayed as follow:

  d26=ΔPyΔτxy

  (4)

  Previous experimental studies have demonstrated that the piezoelectric coefficient d₂₆ in ReCOB crystals attains the maximum value within its piezoelectric tensor [25,35-38,46]. To validate this observation, we performed first-principles calculations leveraging the Density Functional Perturbation Theory (DFPT) method to obtain the piezoelectric matrix. Table 2 presents the static piezoelectric tensor for NdCOB derived from DFPT calculations, where d₂₆ indeed exhibits the highest value (For further details, refer to Table S4 in Supporting information). However, ReCOB crystals exhibit not only different piezoelectric coefficient with various lanthanides but also significant dissimilarity of temperature sensibility, remaining the urgency to clarify the deeper mechanism to guide the material design and hindering the industrial application due to the demand of stability. With the training-sets constituted from DFT calculations, the Deep-learning Potential (DP) was constructed to depict the interatomic dynamic phenomenon with the aid of molecular dynamics (MD) simulation method. The performance of the model reveals the remarkable consistency between DP forecast and DFT result in energy and force as shown in Fig. S2 (Supporting information), moreover taking NdCOB as the example (Fig. S3 in Supporting information), the phonon spectrum and density of state further confirm the ability of DP models to describe dynamic characteristics.

  Table 2

  

  Table 2.  The static piezoelectric matrix from DFPT calculation of NdCOB crystal. The zero values are due to the symmetric restraint (unit: pC/N).

  DownLoad: CSV

  d 11: 2.2836	d 12: 3.4207	d 13: −3.7334	d 14: 0	d 15: 7.1512	d 16: 0
  d 21: 0	d 22: 0	d 23: 0	d 24: 14.2390	d 25: 0	d 26: 18.9848
  d 31: 3.1742	d 32: 3.6269	d 33: −4.0093	d 34: 0	d 35: −1.2673	d 36: 0

  Leveraging the shear-strain method (Fig. S4 in Supporting information), the variable process above could be described in MD simulation. Figs. 2a-c present the relationship of three piezo-units (Re³⁺, Ca1²⁺ and Ca2²⁺) between stress and ionic displacement in ErCOB, GdCOB and NdCOB (for the others refer to Figs. S5-S10 in Supporting information) with the same strain alteration. Hence, the partial d₂₆ of each piezo-unit could be excluded as Figs. 2d-f, where the lanthanides reveal disparate tendency both in temperature and types while remain tow calcium site nearly unchanged. The practical d₂₆ of ReCOB crystals is the summation of each piezo-unit in formular as

  d26=∑ici*d26i

  (5)

  while ci and d26i denote the number of i cation and its partial d₂₆ within primitive cell (with the order being 2 for Re, 4 for Ca1, and 4 for Ca2, respectively). Figs. 2g and h present a comparison between the DP forecast and experimental measurements of the piezoelectric properties in ideal ReCOB crystals. The qualitative agreement, even within the same order of magnitude, underscores the potential of deep learning techniques in materials science, demonstrating their capacity to accurately predict material properties. However, it is important to note that some degree of deviation exists between the predicted and experimental values. This discrepancy may be partially attributed to the similarity in ionic radii between Re³⁺ and Ca²⁺, which can lead to issues related to disorder occupation in the crystal lattice, a factor not considered in this study but deemed essential for more precise predictions [47-49]. Nevertheless, leveraging DP method, the application of the shear-strain method enables the contribution extraction for individual cations, it may provide a neoteric perspective compared to traditional experimental approaches.

  Figure 2

  

  Figure 2.  (a-c) The raw-data of the relationship between the ionic displacement along y zone axis and the shear stress applied in the xy-plane at different temperatures within ErCOB, GdCOB and NdCOB, where the slope is in direct proportion to the shear piezoelectric coefficient d₂₆. (d) The contribution of lanthanide site to d₂₆ calculated from (a-c) at various temperature, revealing its temperature sensibility. (e, f) The contribution of two calcium sites to d₂₆ remaining temperature independency opposite to the lanthanide. (g) The total d₂₆ (the summation of lanthanide and calcium) obtained from DP calculation, exhibiting qualitative consistency with the trend observed in the thermal alterations of the experimental results. (h) The experimental measurements.

  DownLoad: Full-Size Img PowerPoint

  Phonon behavior is tightly tied to their dynamic properties, such as heat conduction, electron-phonon coupling in superconductors, and phonon softening effects in ferro/piezoelectric materials, therefore, it would offer valuable insights for the thermal stability of piezoelectricity [50-54]. In Fig. 3a, the PDOS derived from DPMD simulations for NdCOB, reveals that the phonon density associated with the lanthanide element is primarily concentrated in the low-frequency range of 0–15 THz. Utilizing Fourier Transform to Velocity Autocorrelation Function (FTVAF) method, the PDOS spectra for ReCOB with various lanthanide substitutions can be plotted across a temperature range from 100 K to 900 K, as shown in Figs. 3b and c for Er³⁺ and Nd³⁺ as examples (further details refer to Figs. S11 and S12 in Supporting information) [55,56]. This approach accounts for anharmonic phonon contributions, meaning that the shape of the PDOS is predominantly influenced by anharmonic phonon behavior. Given the definition of weighted average ω¯ and weighted variance (ω−ω¯)2‾ for PDOS to further aid in understanding these effects, that

  ω¯=∫ω·PDOS(ω)·dω∫PDOS(ω)·dω

  (6)

  (ω−ω¯)2‾=∫(ω−ω¯)2·PDOS(ω)·dω∫PDOS(ω)·dω

  (7)

  Figure 3

  

  Figure 3.  Phonon dynamic. (a) The phonon density of states (PDOS) distribution of NdCOB acquired from DPMD simulation, where the lanthanides almost distribute at the region of low-frequency. (b, c) The PDOS of Er³⁺ and Nd³⁺ via various temperatures obtained by the Fourier transform to velocity autocorrelation function method with DPMD. (d) The average weighted phonon vibration frequency ω¯ of PDOS. (e) The corresponding average weighted variance (ω−ω¯)2‾ of PDOS, which reflects the degree of localization of phonons.

  DownLoad: Full-Size Img PowerPoint

  The morphology of PDOS curve can be described quantitatively as shown in Figs. 3d and e, where the alteration in ω¯ reflects the influence of anharmonic interactions across various temperatures, meanwhile the variance (ω−ω¯)2‾ provides insights into the degree of localization or delocalization of the phonon frequency distributions. If the system is consisted solely of harmonic phonons, the PDOS would remain temperature-independent, while the morphology of PDOS spectrum will remain unchanged. In the meantime, as shown in Fig. 3d, the ω¯ for ReCOB with different lanthanides emerges a clear redshift via temperature increases, indicating the softening behavior in phonons due to anharmonic interactions. Additionally, as shown in Fig. 3e, the changes in (ω−ω¯)2‾ correlate strongly with the variations in the shear piezoelectric coefficient contribution from the Re³⁺ subcell (Fig. 2d), implying possibly physical dependency in some extent. Specifically, the variance decreases on heating for Nd³⁺, suggesting a reduced sensitivity of the phonon frequency distribution to different wave vectors, leading to more localized phonon behavior. In contrast, for Er³⁺, the trend is reversed, with phonon states becoming more delocalized.

  Deep-learning potentials are employed to describe the dynamics of ReCa₄O(BO₃)₃, where the lanthanide ions are found to primarily influence the material’s piezoelectric properties and their thermal evolution. However, to fully understand the underlying mechanisms, it is crucial to explore the contribution to piezoelectricity from a localized atomic perspective. To address this, density functional theory (DFT) calculations are utilized to analyze the electronic distribution and chemical bonding interactions at the atomic level.

  Figs. 4a and b present the Born effective charge (BEC) and Bader charge along y axis for each cation [57,58]. Overall, doping with different lanthanide elements induces significant anisotropy in the charge distribution of Re³⁺ ions, particularly when compared to the two different occupancy states of the Ca²⁺ sites. Notably, both charge analysis methods consistently reveal that the variation in Re³⁺ ions follows a trend closely mirroring the experimental piezoelectric coefficient measured in room temperature [23-25]. This alignment suggests a strong correlation between the charge dynamics of Re³⁺ ions and the material’s piezoelectric behavior.

  Figure 4

  

  Figure 4.  Electronic behavior analysis. (a) The Born effective charge and (b) Bader charge in various piezo-unit of ReCOB, where the lanthanide site reveals charge sensitivity opposite to calcium site. The electronic density of states in (c) YCOB and (d) ErCOB, with orbital labels for each element indicated near the peaks. The fonts within the box represent the vanishing orbital states due to electronic transfer and bond formation. Near the Fermi level (located at 0 eV), the regions below and above respectively represent the valence band and conduction band (denoted by "*"), where the occurrence of 4f state shrinks the band gap.

  DownLoad: Full-Size Img PowerPoint

  For BEC, it serves as an indicator of the polarization response of ions in a crystal under external perturbations. A larger BEC suggests that the ions are more likely to undergo substantial displacements in response to the external field, which is of particular significance in ferroelectric and piezoelectric materials. This can be described by the relationship between the atomic force and the electric field as

  Zij*=−eδFiδEj

  (8)

  where Fi represents the force along the i-direction on an ion and Ej corresponds to the component of the electric field along the j-direction. In a given system, a larger BEC leads to a stronger force on the atoms under the same electric field. Physically, both external electric fields and mechanical stress exert similar effects, inducing forces that cause ions to adjust to a new equilibrium. This means that a higher BEC results in larger atomic displacements under a fixed external stimulation, which in turn reflects on the enhancement of the piezoelectric coefficient.

  Bader charge analysis is used to examine the gain or loss of charge by ions, based on the zero-flux surface [58], helping to understand charge transfer within the system. Among the lanthanides in Fig. 4b, Pr³⁺ has the highest valence of +2.132, suggesting a stronger ionic character due to electron loss, bringing it closer to its nominal charge of +3. In contrast, Er³⁺ demonstrates a stronger covalent character. Covalent bonds, being more directional, possess stronger bond configuration and insensitive to external perturbations from traditional knowledge. Another most obvious manifestation is from the perspective of stress, as shown in Figs. 2b-d, ErCOB requires significantly higher stress (~1.2 GPa) compared to NdCOB (~0.9 GPa) under the same maximum strain, leading to a reduction in piezoelectric response owing to both the smaller ionic displacements and enhancement of stress.

  Figs. 4c and d present the electronic density of states (EDOS) for YCOB and ErCOB (see Fig. S18 in Supporting information for others). Generally, the electrons near the Fermi level are involved in the formation of chemical bonds. Due to the high electronegativity of oxygen, electron density flows from adjacent metallic elements toward oxygen. For the elements Y, Er, and Ca, the covalent electrons are Y: 5s, Y: 4d, Er: 6s, Er: 5d, and Ca: 4s, respectively. These orbitals exhibit significant depletion near the Fermi surface, leaving only the O 2p orbitals, which indicates the formation of chemical bonds through electron transfer and ionicity. However, the 4f electrons of Er³⁺ near the Fermi level disrupt this condition (Fig. 4d). This deviation is attributed to the localization of the 4f electrons (Fig. S19 in Supporting information), which, despite their energy proximity to the Fermi surface, do not directly participate in bond formation. Additionally, the Er: 4f* conduction bands above the Fermi level result in a narrowing of the band gap, thereby enhancing the intrinsic conductivity of the material. Our DFT calculations demonstrate that the band shifts in response to varying lanthanide doping align well with the tendency of experimental observations, as shown in Fig. S18 (ignoring the contribution of oxygen vacancies’ migration [59]). Although the 4f electrons do not engage directly in chemical bonding, they exert a significant indirect influence on the energy levels of other orbitals due to the shielding effect they impose on the atomic nucleus, as evidenced in Table S6 (Supporting information). For instance, in Fig. 4d, the localized deep orbital Er: 5s breaks the symmetry of the system due to the unequal distribution of the spin-up (5 electrons) and spin-down (7 electrons) 4f states. Furthermore, the localization of the 4f electrons induces asymmetry in EDOS, characterized by a non-uniform distribution between spin-up and spin-down states, which in turn contributes to the material’s unique magnetic and optical properties [60,61]. Thus, while the 4f electrons do not directly participate in chemical bonding, they modulate the energy levels of surrounding orbitals and influnce the ionic radius.

  Valence bond analysis may provide other perspective for understanding the mechanism. Crystal orbital Hamilton population (COHP) analysis, a quantitative method for analyzing atomic interactions in a crystal, further supports these observations [62]. The integrated COHP below the Fermi level (ICOHP) provides a measure of bond strength. As shown in Fig. 5a, the average ICOHP for Er³⁺ with its coordinated oxygen atoms is higher (greater absolute value), suggesting stronger chemical bonding than in Pr³⁺, being consistent with covalent situation of Er³⁺. Additionally, the presence of anti-bonding states below the Fermi level in Pr³⁺, indicates weaker bond stability and stronger ionicity in Pr³⁺ in comparison to Er³⁺ as well. Furthermore, Crystal Orbital Bond Indices (COBI) analysis could measure the property to a single bond [63]. A value closer to zero indicates a more ionic bond, while a value approaching one reflects a covalent nature. In Fig. 5c, we display the local bonding features of each piezoelectric unit in NdCOB, with orange fonts representing bond COBI and black fonts showing bond lengths. Overall, Ca²⁺ demonstrates a strong ionic character, while the bond between B³⁺ and oxygen is almost purely covalent. Meanwhile, the bond between Nd³⁺ and oxygen has a hybrid character, resembling that of Ti⁴⁺ in ferroelectric BaTiO₃, sitting in an intermediate state. This hybrid nature is sensitive to changes in bond length, meaning the physical properties of the Re³⁺ unit can fluctuate with temperature due to the thermal expansion and the non-homophonic phonons, while the ionic calcium remains largely unaffected by such changes. For Nd³⁺, Fig. 2d highlights notable directional anisotropy in the COBI values. Along the y-axis, the Nd-O5 and Nd-O6 bonds exhibit a COBI value of 0.3533 (bond length of 2.476 Å). In contrast, within the xz-plane, two strong covalent bonds appear: Nd-O1 (COBI: 0.5642) and Nd-O3 (COBI: 0.5788). The strength of these covalent bonds restricts the movement of Nd³⁺ within the xz-plane, resulting in lower piezoelectric coefficients d_aj along the x and z axes when a = 1 and 3 (Table 2). The Differential Charge Density of NdCOB within two specific planes is illustrated in Fig. 5b, where the lighter regions represent electron accumulation opposite to the electron depletion in darker regions. The prominent lighter region along the line connecting boron and oxygen highlights significant electronic hybridization, which is indicative of a covalent bonding nature. Additionally, the variation in the distribution of light and dark regions around the Nd-O1 and Nd-O2 bonds aligns with the observations in Fig. 5c, that Nd-O1 exhibit more covalent than Nd-O2. Regarding calcium, the clear gap between the calcium cation and the adjacent oxygens suggests a predominantly ionic bonding character.

  Figure 5

  

  Figure 5.  Chemical bond analysis. (a) Average integrated crystal orbital Hamilton population (iCOHP) Re site with its coordination oxygen, the value directly related to the bond strength. (b) The charge density difference mapping of NdCOB in various plane respectively containing Nd and Ca2 cations. (c) Crystal Orbital Bond Indices (COBI) to each cation of NdCOB with its adjacent oxygen, the black and organ fonts represent bond length and COBI respectively, the COBI value evolves from 0 to 1 reflecting a transition from ionicity to covalency.

  DownLoad: Full-Size Img PowerPoint

  In summary, through multi-perspective conjoint analysis, we elucidate the mechanism and contribution of piezoelectric properties of ReCOB systems from each sub-unit. Leveraging the DP guidance for the dynamic calculation, the shear piezoelectricity d₂₆ was obtained by the stress-strain method aligning the experimental result, that highlights the great capacity of artificial neural networks in materials science. The lanthanide site in ReCOB contribute the discrepancy of piezoelectricity and its thermal altering, while the contribution from calcium remains unchanged. The localization of phonons is closely linked to the temperature-dependent drift of piezoelectric performance, suggesting that the distribution of chemical bonds, modulated by anharmonic phonons, plays a pivotal role in driving the temperature stability. For static calculations, charge and chemical bonds analysis reveals that the strength of the chemical bonds of Re³⁺ cations exhibits a clear inverse correlation with piezoelectric performance, suggesting that bonds with a totally more ionic character are more favorable for enhancing piezoelectricity. The 4f orbital does not participate in the formation of chemical bonding, but by altering the energy-level and the ionic radius by the shielding effect on atomic nuclei. This work presents a distinctive approach to understanding the mechanisms in non-ferroelectric materials and provides essential insights and guidance for the design of new materials.

  CRediT authorship contribution statement

  Xingshuai Ma: Writing – original draft, Visualization, Software, Methodology, Investigation, Data curation, Conceptualization. Chuanrui Huo: Visualization, Formal analysis. Xiaoniu Tu: Supervision, Investigation, Data curation. Hua Tan: Supervision, Methodology, Formal analysis, Data curation, Conceptualization. He Qi: Writing – review & editing, Writing – original draft, Supervision, Funding acquisition, Conceptualization. Shi Liu: Writing – review & editing, Investigation, Formal analysis, Conceptualization. Jun Chen: Writing – review & editing, Supervision, Funding acquisition.

  Declaration of competing interest

  The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

  Acknowledgments

  This work was supported by the National Key R&D Program of China (No. 2022YFB3204000) and the Beijing Outstanding Young Scientist Program (No. JWZQ20240101015).

  Supplementary materials

  Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cclet.2025.111895.

  
  
• 1. [1]

     H. Jaffe, J. Am. Ceram. Soc. 41 (1958) 494–498. doi: 10.1111/j.1151-2916.1958.tb12903.x

  2. [2]

     G.H. Haertling, J. Am. Ceram. Soc. 82 (1999) 797–818. doi: 10.1111/j.1151-2916.1999.tb01840.x

  3. [3]

     W. Jaffe, H.Jaffe Cook, Piezoelectric Ceramics, Academic Press, New York, 1971.

  4. [4]

     G.T. Hwang, H. Park, J.H. Lee, et al., Adv. Mater. 26 (2014) 4880–4887. doi: 10.1002/adma.201400562

  5. [5]

     L. Chen, H. Liu, H. Qi, J. Chen, Prog. Mater. Sci. 127 (2022) 100944. doi: 10.1016/j.pmatsci.2022.100944

  6. [6]

     H.D. Li, C. Tian, Z.D. Deng, Appl. Phys. Rev. 1 (2014) 041301. doi: 10.1063/1.4900845

  7. [7]

     C. Bowen, H.A. Kim, P. Weaver, S. Dunn, Energy Environ. Sci. 7 (2014) 25–44. doi: 10.1039/C3EE42454E

  8. [8]

     C. Qiu, B. Wang, N. Zhang, et al., Nature 577 (2020) 350–354. doi: 10.1038/s41586-019-1891-y

  9. [9]

     Y. Saito, H. Takao, T. Tani, et al., Nature 432 (2004) 84–87. doi: 10.1038/nature03028

  10. [10]

      W. Liu, X.B. Ren, Phys. Rev. Lett. 103 (2009) 257602. doi: 10.1103/PhysRevLett.103.257602

  11. [11]

      J. Scott, Science 315 (2007) 954–959. doi: 10.1126/science.1129564

  12. [12]

      K. Xu, J. Li, X. Lv, et al., Adv. Mater. 28 (2016) 8519–8523. doi: 10.1002/adma.201601859

  13. [13]

      F. Li, M.J. Cabral, B. Xu, et al., Science 364 (2019) 264–268. doi: 10.1126/science.aaw2781

  14. [14]

      S. Wada, S. Suzuki, T. Noma, et al., Jpn. J. Appl. Phys. 38 (1999) 5505. doi: 10.1143/JJAP.38.5505

  15. [15]

      X. Lv, J. Zhu, D. Xiao, X.X. Zhang, J. Wu, Chem. Soc. Rev. 49 (2020) 671–707. doi: 10.1039/c9cs00432g

  16. [16]

      S.J. Zhang, F. Li, J. Appl. Phys. 111 (2012) 031301. doi: 10.1063/1.3679521

  17. [17]

      C.T. Nelson, P. Gao, J.R. Jokisaari, et al., Science 334 (2011) 968–971. doi: 10.1126/science.1206980

  18. [18]

      L.Z. Li, L. Xie, X.Q. Pan, Rep. Prog. Phys. 82 (2019) 126502. doi: 10.1088/1361-6633/ab28de

  19. [19]

      J. Li, R. Rogan, E. Üstündag, K. Bhattacharya, Nat. Mater. 4 (2005) 776–781. doi: 10.1038/nmat1485

  20. [20]

      J. Fu, A.W. Xie, T.Y. Li, R.Z. Zuo, Acta Mater. 230 (2022) 117862. doi: 10.1016/j.actamat.2022.117862

  21. [21]

      T. Zheng, J. Wu, D. Xiao, J. Zhu, Prog. Mater Sci. 98 (2018) 552–624. doi: 10.1016/j.pmatsci.2018.06.002

  22. [22]

      M. Waqar, H.J. Wu, J.S. Chen, K. Yao, J. Wang, Adv. Mater. 34 (2022) 2106845. doi: 10.1002/adma.202106845

  23. [23]

      J. Wu, X.S. Ma, D.H. Zhou, et al., Adv. Mater. 37 (2025) 2419134. doi: 10.1002/adma.202419134

  24. [24]

      H. Yan, H. Ning, Y. Kan, P. Wang, M.J. Reece, J. Am. Ceram. Soc. 92 (2009) 2270–2275. doi: 10.1111/j.1551-2916.2009.03209.x

  25. [25]

      H.F. Zu, H.Y. Wu, Q.M. Wang, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63 (2016) 486–505. doi: 10.1109/TUFFC.2016.2527599

  26. [26]

      D. Balitsky, E. Philippot, V. Balitsky, et al., Molecules 25 (2020) 4518. doi: 10.3390/molecules25194518

  27. [27]

      H. Takeda, M. Hagiwara, H. Noguchi, et al., Appl. Phys. Lett. 102 (2013) 242907. doi: 10.1063/1.4811163

  28. [28]

      T. Kim, J. Kim, R. Dalmau, et al., IEEE Trans. Ultrason. Ferroelectr. Freq. Control 62 (2015) 1880–1887. doi: 10.1109/TUFFC.2015.007252

  29. [29]

      C.S. Pandey, J. Schreuer, J. Appl. Phys. 111 (2012) 013516. doi: 10.1063/1.3673820

  30. [30]

      Y.N. Li, C. Jiang, G.L. Wang, et al., J. Alloys Compd. 1019 (2025) 179263. doi: 10.1016/j.jallcom.2025.179263

  31. [31]

      C.Y. Shen, H.J. Zhang, H.J. Cong, et al., J. Appl. Phys. 116 (2014) 044106. doi: 10.1063/1.4891827

  32. [32]

      S.J. Zhang, Y.Q. Zheng, H.K. Kong, et al., J. Appl. Phys. 105 (2009) 114107. doi: 10.1063/1.3142429

  33. [33]

      T. Iwataki, H. Ohsato, K. Tanaka, et al., J. Eur. Ceram. Soc. 21 (2001) 1409–1412. doi: 10.1016/S0955-2219(01)00029-2

  34. [34]

      H. Fritze, J. Electroceram. 26 (2011) 122–161. doi: 10.1007/s10832-011-9639-6

  35. [35]

      F. Yu, S. Hou, X. Zhao, S. Zhang, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61 (2014) 1344–1356. doi: 10.1109/TUFFC.2014.3042

  36. [36]

      S. Zhang, F. Yu, J. Am. Ceram. Soc. 94 (2011) 3153–3170. doi: 10.1111/j.1551-2916.2011.04792.x

  37. [37]

      F.P. Yu, S.J. Zhang, X. Zhao, et al., IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58 (2011) 868–873. doi: 10.1109/TUFFC.2011.1881

  38. [38]

      F. Yu, G. Yao, G. Wang, X. Liu, X. Zhao, J. Chin. Ceram. Soc. 50 (2022) 563–574.

  39. [39]

      S.J. Zhang, Y.T. Fei, E. Frantz, et al., IEEE Trans. Ultrason. Ferroelectr. Freq. Control 55 (2008) 2703–2708. doi: 10.1109/TUFFC.2008.985

  40. [40]

      Y. LeCun, Y. Bengio, G. Hinton, Nature 521 (2015) 436–444. doi: 10.1038/nature14539

  41. [41]

      K. Choudhary, B. DeCost, C. Chen, et al., npj Comput. Mater. 8 (2022) 59. doi: 10.1038/s41524-022-00734-6

  42. [42]

      L.F. Zhang, J.Q. Han, H. Wang, et al., Phys. Rev. Lett. 120 (2018) 143001. doi: 10.1103/PhysRevLett.120.143001

  43. [43]

      H. Wang, L. Zhang, J. Han, Comput. Phys. Commun. 228 (2018) 178–184. doi: 10.1016/j.cpc.2018.03.016

  44. [44]

      J. Yang, S. Liu, Phys. Rev. B 110 (2024) 214112. doi: 10.1103/PhysRevB.110.214112

  45. [45]

      Y.H. Hu, J.Y. Yang, S. Liu, Phys. Rev. Lett. 133 (2024) 046802. doi: 10.1103/PhysRevLett.133.046802

  46. [46]

      J. Wu, J. Yang, Y.J. Liu, et al., Phys. Rev. B 108 (2023) L180104. doi: 10.1103/PhysRevB.108.L180104

  47. [47]

      L.Y. Bai, Z.J. Liu, Q.S. Bao, et al., J. Mater. Chem. C 13 (2025) 12793–12800. doi: 10.1039/d5tc00878f

  48. [48]

      F.P. Yu, S.J. Zhang, X. Zhao, et al., J. Phys. D: Appl. Phys. 44 (2011) 135405. doi: 10.1088/0022-3727/44/13/135405

  49. [49]

      X.Y. Lu, L.L. Li, S.W. Tian, et al., J. Mater. Chem. C 8 (2020) 10109–10120. doi: 10.1039/d0tc01073a

  50. [50]

      G. Chen, Nat. Rev. Phys. 3 (2021) 555–569. doi: 10.1038/s42254-021-00334-1

  51. [51]

      M.N. Luckyanova, J. Garg, K. Esfarjani, et al., Science 338 (2012) 936–939. doi: 10.1126/science.1225549

  52. [52]

      F. Giustino, Rev. Mod. Phys. 89 (2017) 015003. doi: 10.1103/RevModPhys.89.015003

  53. [53]

      A. Lanzara, P. Bogdanov, X. Zhou, et al., Nature 412 (2001) 510–514. doi: 10.1038/35087518

  54. [54]

      P.C. Xie, Y.X. Chen, W.N. E, R. Car, Phys. Rev. B 111 (2025) 094113. doi: 10.1103/PhysRevB.111.094113

  55. [55]

      J.M. Dickey, A. Paskin, Phys. Rev. 188 (1969) 1407–1418. doi: 10.1103/PhysRev.188.1407

  56. [56]

      B. Fultz, Prog. Mater Sci. 55 (2010) 247–352. doi: 10.1016/j.pmatsci.2009.05.002

  57. [57]

      X. Gonze, C. Lee, Phys. Rev. B 55 (1997) 10355. doi: 10.1103/PhysRevB.55.10355

  58. [58]

      M. Yu, D.R. Trinkle, J. Chem. Phys. 134 (2011) 064111. doi: 10.1063/1.3553716

  59. [59]

      S. Tian, L. Li, X. Lu, et al., Acta Mater. 183 (2020) 165–171. doi: 10.1016/j.actamat.2019.11.013

  60. [60]

      N.D. Kelly, S.E. Dutton, Inorg. Chem. 59 (2020) 9188–9195. doi: 10.1021/acs.inorgchem.0c01098

  61. [61]

      Y.Q. Liu, Z.P. Wang, F.P. Yu, et al., Opt. Express 25 (2017) 11867–11893. doi: 10.1364/OE.25.011867

  62. [62]

      R. Dronskowski, P.E. Bloechl, J. Phys. Chem. 97 (1993) 8617–8624. doi: 10.1021/j100135a014

  63. [63]

      P.C. Müller, C. Ertural, J. Hempelmann, R. Dronskowski, J. Phys. Chem. C 125 (2021) 7959–7970. doi: 10.1021/acs.jpcc.1c00718

• 

  Figure 1  (a) Relationship between piezoelectric coefficient (d₃₃ or d₂₆) and operating temperature for various ceramics and single crystals. (b) The structure of ReCOB and its three octahedral piezoelectric subunits, the orange arrows reveal the polarization vectors. (c) The iDPC-STEM image of NdCOB perpendicular to z zone axis.

  下载: 全尺寸图片 幻灯片
  

  Figure 2  (a-c) The raw-data of the relationship between the ionic displacement along y zone axis and the shear stress applied in the xy-plane at different temperatures within ErCOB, GdCOB and NdCOB, where the slope is in direct proportion to the shear piezoelectric coefficient d₂₆. (d) The contribution of lanthanide site to d₂₆ calculated from (a-c) at various temperature, revealing its temperature sensibility. (e, f) The contribution of two calcium sites to d₂₆ remaining temperature independency opposite to the lanthanide. (g) The total d₂₆ (the summation of lanthanide and calcium) obtained from DP calculation, exhibiting qualitative consistency with the trend observed in the thermal alterations of the experimental results. (h) The experimental measurements.

  下载: 全尺寸图片 幻灯片
  

  Figure 3  Phonon dynamic. (a) The phonon density of states (PDOS) distribution of NdCOB acquired from DPMD simulation, where the lanthanides almost distribute at the region of low-frequency. (b, c) The PDOS of Er³⁺ and Nd³⁺ via various temperatures obtained by the Fourier transform to velocity autocorrelation function method with DPMD. (d) The average weighted phonon vibration frequency ω¯ of PDOS. (e) The corresponding average weighted variance (ω−ω¯)2‾ of PDOS, which reflects the degree of localization of phonons.

  下载: 全尺寸图片 幻灯片
  

  Figure 4  Electronic behavior analysis. (a) The Born effective charge and (b) Bader charge in various piezo-unit of ReCOB, where the lanthanide site reveals charge sensitivity opposite to calcium site. The electronic density of states in (c) YCOB and (d) ErCOB, with orbital labels for each element indicated near the peaks. The fonts within the box represent the vanishing orbital states due to electronic transfer and bond formation. Near the Fermi level (located at 0 eV), the regions below and above respectively represent the valence band and conduction band (denoted by "*"), where the occurrence of 4f state shrinks the band gap.

  下载: 全尺寸图片 幻灯片
  

  Figure 5  Chemical bond analysis. (a) Average integrated crystal orbital Hamilton population (iCOHP) Re site with its coordination oxygen, the value directly related to the bond strength. (b) The charge density difference mapping of NdCOB in various plane respectively containing Nd and Ca2 cations. (c) Crystal Orbital Bond Indices (COBI) to each cation of NdCOB with its adjacent oxygen, the black and organ fonts represent bond length and COBI respectively, the COBI value evolves from 0 to 1 reflecting a transition from ionicity to covalency.

  下载: 全尺寸图片 幻灯片

  Table 1.  The piezoelectric subunit information of NdCOB. For each element, ionic displacement is defined as the difference between the center cationic coordinate and the average value of nearest neighbor coordination oxygen atoms.

  Cation	Coordination number	X axis (Å)	Y axis (Å)	Z axis (Å)	Average ionic displacement (Å)
  Nd	6	−0.3102	0.0000	−0.1401	0.3403
  Ca1	6	0.0562	0.1058	0.0063	0.1200
  Ca1	6	0.0562	−0.1058	0.0063	0.1200
  Ca2	6	−0.3244	−0.1785	−0.2240	0.4328
  Ca2	6	−0.3244	0.1785	−0.2240	0.4328
  B	3	0.0067	0.0000	0.0002	0.0067
  B	3	0.0018	−0.0013	0.0063	0.0066
  B	3	0.0018	0.0013	0.0063	0.0066

  下载: 导出CSV

  Table 2.  The static piezoelectric matrix from DFPT calculation of NdCOB crystal. The zero values are due to the symmetric restraint (unit: pC/N).

  d 11: 2.2836	d 12: 3.4207	d 13: −3.7334	d 14: 0	d 15: 7.1512	d 16: 0
  d 21: 0	d 22: 0	d 23: 0	d 24: 14.2390	d 25: 0	d 26: 18.9848
  d 31: 3.1742	d 32: 3.6269	d 33: −4.0093	d 34: 0	d 35: −1.2673	d 36: 0

  下载: 导出CSV
•