Spontaneous polarization, n-type half-metallicity, low lattice thermal conductivity, and high structure stabilities in F@O-doped PbTiO₃

Abstract

Materials displaying ferroelectric and ferromagnetic properties offer a fascinating foundation for future spintronic innovations like multi-state random-access memory and data storage devices. Herein, numerous features of pristine and F@O-doped-PbO (PO)/TiO₂ (TO)-layer PbTiO₃ perovskite oxides are investigated via ab initio calculations. The calculated negative formation enthalpies and elastic coefficients verify the thermodynamic and mechanical stability of the structures, respectively. The pristine motif exhibits a giant spontaneous polarization (P) of 88 μC cm⁻², having a non-magnetic insulating state with an indirect energy gap (E_g) of 2.11 eV. It appears that the dopant reduces the structural distortions, lowering P to 42.56/42.78 μC cm⁻² in the F@O-doped PO/TO-layer-based structure. The most notable aspect of the present study is that F-doping in the PO layer induces an n-type half-metallic (HM) ferromagnetic (FM) behavior with a total magnetic moment (m_t) of 1.0 μ_B. It is found that Ti ions are the main contributor to the magnetism, which is confirmed by spin-magnetization density isosurface plots. Additionally, a large E_g of 2.49 eV in the spin-minority channel is found, which is large enough to ensure the HM state and avoid reverse leakage current. Conversely, the F@O-doped TO-layer-based system transforms to an FM semiconductor with an E_g of 0.23 eV. Interestingly, the highest figure of merit of 0.72/0.56/0.49 is predicted at 700/700/400 K for the pristine/F@O-doped PO/TO-layer-based structures with the inclusion of lattice thermal conductivity. Thus, due to high structural stability, high P, half-metallicity, and low thermal conductivity, the F@O motif emerges as a promising candidate for various potential applications in spintronics and energy conversion.

---

1 Introduction

Materials that demonstrate concurrent existence of multiple fundamental ferroic phases in one state, like magnetic and ferroelectric (FE)/ferroelastic, are referred to as multiferroics (MFs).^1–3 The interaction of various ferroic orders like magnetoelasticity,⁴ magnetocaloric effects,⁵ piezoelectricity,⁶ and the magnetoelectric effect (ME),^7,8 which allows remote polarization (P) modification with a magnetic field or remote magnetism with an electric field, adds further adaptability to the development of new multipurpose devices. Combined FE and ferromagnetic (FM) phases in a compound are becoming more desirable due to their effective contribution to future-oriented technologies, such as tunable and variable memory.^2,9 MFs are generally categorized into two types. In type I, ferroelectricity and magnetism originate from two different cations, such as in BiFeO₃ (BFO). In type II ferroelectrics, such as RMnO₃, electric P is induced by specific magnetic ordering, providing an incredible ME response.¹⁰ In particular, the enhanced interactions involving elasticity, spin, and charge in ferroic materials have attracted much curiosity owing to their fascinating physical aspects and numerous multi-faceted applications, like electric current sensors,¹¹ nonvolatile memory elements,¹² magnetic storage devices,¹³ AC/DC magnetic field sensors,¹⁴ solar cells,¹⁵ and MF magnetic recording read heads.¹⁶ The P of approximately 10 μC cm⁻², which is considerably lower than the P of roughly 100 μC cm⁻² for type I MFs, is a major issue, even though type II MFs have an extraordinarily strong ME coupling.¹⁷ The challenges of low P in type II and weak ME coupling in type I MFs are real obstacles in device applications.¹³ The earliest identified MF was nickel iodine boracite (Ni₃B₂O₁₃I),¹⁸ which demonstrated simultaneous occurrence of FM and FE behavior at about 64 K.¹⁹ However, the immense dimensions of the unit cell, asymmetric configuration, and absence of magnetoelectricity at lower temperatures (temp.) were impractical.²⁰ Hence, single-phase multifunctional materials are very rare, as long-range Coulomb interactions have been eliminated due to conduction of the electrons in the material, resulting in FE instability and distortion.²¹

In this vein, BFO is a well-known type I MF perovskite oxide (PO) that exhibits a MF phase at room temp., holding a substantial P of approximately 80 μC cm⁻² with a T_C (Curie temp.) of 830 °C and tilted G-type antiferromagnetic (AFM) character with a 370 °C Néel temp. (T_N).^24,25 Unfortunately, the weak ME coupling²⁶ and current-leakage issues²⁷ make it less attractive for technical applications. Similarly, BiMnO₃,²⁸ MnWO₄,²⁹ and DyMnO₃ (ref. 30) are among the recently discovered transition metal (TM) MF POs. Ferroelectricity triggered by magnetism has been observed recently in LnMn₂O₅,³¹ SmFeO₃,³² and DyFeO₃.³³ An array of mechanisms, such as exchange coupling in BaFe₂Se₃,³⁴ a spin-driven effect in BaYFeO₄,³⁵ and charge ordering in RMnO₃ (R = Er, Ho, In, Lu, Tm, Sc, and Y)³⁶ and RMn₂O₅ (R = rare earth elements),³⁷ have been explored, which trigger the MF phenomenon. Regarding POs, PbTiO₃ (PTO) is considered the best FE candidate due to its insulating and non-magnetic nature. It has an enormous P, ranging between 75 and 92 μC cm⁻² (ref. 38–40), and an energy gap (E_g) of 3.88 eV.⁴¹ The primary cause for this significant P is the shifting of the Ti⁴⁺ ion with respect to the oxygen lattice.^42,43 Similar to other FE materials like BaTiO₃ (BTO)⁴⁴ and KNbO₃,⁴⁵ PTO also suffers from current leakage, especially in thin films when the thickness decreases to micro- or nanometers.⁴⁶ This drawback is considered the main hindrance in technological applications such as FE memories.⁴⁷ In-spite of the toxic effects of lead (Pb), PTO provides potential for higher-temperature application as it exhibits the highest thermodynamic stability among all FE materials, having a Curie temperature close to 500 °C.⁴⁸ Furthermore, it has exceptional FE applications, including FE field-effect transistors⁴⁹ and nonvolatile FE random-access memory,⁵⁰ due to its exceptionally large P.⁴⁰ Additionally, owing to its large pyroelectric coefficient, it exhibits outstanding pyroelectric features.^51,52

It is believed that fully occupied d-orbitals are essential for FE materials, whereas partially filled unpaired d- and/or f-orbital electrons are responsible for magnetism in all magnetic materials. The fundamental condition to create a polar interaction, the d⁰ rule, requires unoccupied d-orbitals for P to arise. Therefore, it is necessary to violate the d⁰ rule to observe the FE state and magnetism in a single phase.^53,54 Consequently, combining these two states in the same phase proves extremely challenging. This leads to a range of approaches to induce the MF effect, including strain engineering (strain-induced multiferroicity in SrMnO₃ (ref. 55)), and doping (Mn-doped (dop.) KNbO₃,⁵⁶ and Fe-dop. BaTiO₃, which demonstrates the coexistence of FE and FM phases at room temp⁵⁷). Additionally, Bi³⁺ substitution in YFeO₃ (ref. 58) and Co/Nd doping at the Fe/Bi position in BFO⁵⁹ have been observed to strengthen MF effects. In a theoretical investigation, an enormous P of 70.32 μC cm⁻² and net magnetic moment (m_t) of 2.9 μ_B were found when PTO was dop. with Zn/Mn TM at the Ti position.⁶⁰ Likewise, Nd/Fe-doping at the Pb/Ti site,⁶¹ V-doping at the Ti site,⁶² and 3/4/5d-electron doping at the Ti site^60,63 in MF PTO systems have been previously investigated.^64,65

Motivated by the above discussion, the structural durability, FE, electronic, magnetic, and thermoelectric (TE) aspects of both the pristine (prst.) and F@O-dop. PTO systems are explored via density functional theory (DFT)-based computations. Initially, the structural durability of both structures is examined using enthalpy of formation (ΔH_f). Applying a stress–strain technique, the elastic constants, C_ij, are determined. The FE features are demonstrated by determining the values of P and the ∠O–Ti–O bond angles. The electronic structures are then illustrated by plotting the total and partial density of states (TDOS and PDOS) and band structures. Additionally, m_t and partial spin magnetic moments (m_s), as well as three-dimensional spin-magnetization density isosurface plots, have been utilized for analyzing magnetic traits. Ultimately, transport characteristics are disclosed on the basis of various thermoelectric (TE) parameters, which demonstrate the material’s energy conversion efficiency.

2 Computational and structural details

Spin-polarized (SP) DFT-based computations have been carried out employing the full-potential linearized augmented plane wave technique as employed in the Wien2k code.⁶⁶ An exchange-correlation functional, namely the generalized gradient approximation (GGA),⁶⁷ adopting the Hubbard parameter (U) of 4.4 eV,⁶⁸ has been employed on the Ti-3d states to analyze the substantial correlational impact of the Ti-3d states. However, hybrid exchange-correlation functionals like B3PW or B3LYP allow one to achieve an exceptional agreement with the experimental findings with respect to the band gaps of numerous perovskites, such as BTO, CaTiO₃, PbTiO₃, SrTiO₃, BaZrO₃, CaZrO₃, PbZrO₃, and SrZrO₃.⁶⁹ The ultimate value of l_max associated with wavefunction expansion within the atomic sphere is set to 12, and the plane-wave limit is assumed to be R_mt × K_max = 7 and G_max = 24. To ensure that the charge is not escaping from the core and the total energy (E_t) converges, we assign K_max to be the highest possible value. Fully relaxed crystal structures are achieved by lowering the force of the atoms to 5 mR per a. u. and an E_t convergence of 10⁻⁵ Ry is assumed to be self-consistent.

The prst. PTO, with space group number 99 (P4mm) and experimental lattice parameters of a = b = 3.905 Å and c = 4.156 Å, has been demonstrated as being a stable phase at room temp. The Pb, Ti, O₁, and O₂ are arranged at (0, 0, 0), (0.5, 0.5, 0.5 + Δ), (0.5, 0.5, Δ), and (0.5, 0, 0.5 + Δ), sites, respectively. In this case, Δ is a distortion along the c-axis and the experimentally reported values for Ti, O₁, and O₂ are 0.040, 0.112, and 0.112, respectively.⁷⁰ The crystalline makeup of the prst. cell consists of 5 atoms including one Pb, Ti and O₁, and two O₂, and it is replicated eight times in order to generate the F@O-dop. systems. This results in a 2 × 2 × 2 supercell comprising 40 atoms, with 8Pb, 8Ti, and 24O. To achieve a 12.5% total impurity concentration, one F atom is replaced at the O-site in both PbO and TiO₂ (PO and TO) layers, which is simple to carry out in experiments.⁶³ The completely relaxed crystalline configurations of the prst. and F@O-dop. (in the PO/TO layers) PTO structures in a 2 × 2 × 2 supercell configuration are displayed in Fig. 1(a), (b) and (c), respectively.

Fig. 1 Graphical diagrams of the PbTiO₃ structures in a 2 × 2 × 2 supercell for the (a) pristine, (b) F@O-doped PbO (PO) layer, and (c) F@O-doped TiO₂ (TO) layer.

3 Results and discussion

To explore the thermodynamic crystal stability of the prst. system, ΔH_f is initially determined by utilizing the competitive chemical potential (μ) components of PO and TO:

ΔH^PO_f ≥ Δμ_Pb + Δμ_O (1)

ΔH^TO_f ≥ Δμ_Ti + 2Δμ_O (2)

ΔH^PTO_f = Δμ_Pb + Δμ_Ti + 3Δμ_O (3)

Here ΔH^PO_f/ΔH^TO_f/ΔH^PTO_f denote the ΔH_f of the PO/TO/PTO, which can be calculated as:

To lessen the overestimation in the anticipated ΔH_f values, the changes suggested by Wang et al.⁷¹ have additionally been adopted. Thus, our estimated ΔH_f of −2.47 eV for the PO phase correlates significantly with the experimentally determined value of −2.27 eV (ref. 72) and the theoretical estimations of −2.59 and −2.437 eV.^73,74 In addition, the estimated ΔH_f of −10.22 eV for the TO phase is also analogous to the previously reported experimental values of −9.79/−9.75 eV (ref. 72) and the theoretically obtained −10.36/−9.783/−10.38 eV.^73–75 Ultimately, the predicted ΔH_f for the prst. PTO is −12.56 eV, which coincides well with the highly comparable structure of BTO (−17.20 eV)⁷² and the previously reported DFT results of −13.34 and −12.342 eV.^73,74 Consequently, PTO’s structural stability has been confirmed by all of the negative values of ΔH_f in the aforementioned circumstances, which further support eqn (1)–(3). Despite the negative ΔH_f, the thermodynamic durability with respect to the component phases, PO and TO, needs to be confirmed by determining the reaction energy (E^PTO_React.). PO + TO → PTO is a permissible chemical reaction, and the associated E^PTO_React. is calculated as:

E^PTO_React. = ΔH^PTO_f − (ΔH^PO_f + ΔH^TO_f) (7)

Utilizing the prior established ΔH_f of PO, TO, and PTO, E^PTO_React. is −0.41 eV per f. u., which indicates that the reaction is exothermic and the growth of the prst. structure is empirically stable relative to its constituents. Additionally, it has been found that E^PTO_React. is linked to the PTO decomposition energy (E^PTO_decomp.), which is obtained via the convex-hull approach:⁷⁶

The number of atoms contributing to the chemical reaction is denoted by n in this particular case. As a result, the predicted E^PTO_decomp. is 0.09 eV per atom, which validates the structural stability of the prst. motif regarding its constituents. Next, the energy of defect generation (Ω[F_O]), in eV per defect is evaluated to confirm the stable and practical substitution of F in PTO at the O-site:

Ω[F_O] = E^PTO_t − E^PTFO_t − μ_F + μ_O (9)

where μ_F/μ_O denote the chemical potential of the F/O atom and indicate the E_t of the prst./F-dop. structure. The current calculations satisfy the requirement for maximum variance of μ, where the μ for O is modified from O-rich (i.e., ) to O-poor (i.e., ). Additionally, (space group no. 15) is employed to represent the most persistent magnetizing phases. Consequently, for O-rich and O-poor circumstances, the calculated Ω[F_O] values are 1.3 and 2.1 eV per defect, respectively, for the F@O-dop.PO layer. Similarly, the Ω[F_O] values are 1.5 and 2.5 eV per defect for O-rich and O-poor conditions, respectively, in the F@O-dop.TO layer. Finally, the projected threshold for Ω of less than 3 eV per defect⁷² implies the feasible growth of F-dop. structures in the lab.

Additionally, by executing six finite lattice distortions, elastic tensors C_ij are generated from the strain–stress relation⁷⁷ for the prst. and F@O-dop.PO/TO-layer-based PTO structures. They reveal information regarding the mechanical durability of solids and dynamical features of a material, like the Gruneisen factor, specific heat, Debye temperature, and thermal expansion. They also provide crucial information about elastic parameters, including internal strain, load deflection, sound velocity, fracture toughness, and thermo-elastic stress.⁷⁸ To verify the mechanical strength of the prst. and F@O-dop.PO/TO-layer-based PTO structures, the Born criteria for a tetragonal phase are utilized:⁷⁹

The calculated elasticity coefficients are provided in Table 1, where the prst. and F@O-dop.PO/TO-layer-based PTO structures meet the aforementioned conditions for mechanical robustness and rigidity is assured. There are actually two Cauchy relationships for tetragonal crystals, C₁₂ = C₆₆ and C₃₃ = C₄₄. It is apparent from Table 1 that these two relationships are not followed by the determined elastic constants. These Cauchy relationships are only fulfilled if the atoms interact solely through the central forces, demonstrating that the non-central forces have significance in the prst. and F@O-dop.PO/TO-layer-based PTO structures because of their covalent nature.⁷⁹

Table 1 The calculated elastic constants (C_ij) in GPa for the pristine and F@O-doped PO/TO-layer-based PbTiO₃ structures

Elastic constants	C11	C12	C13	C33	C44	C66
PbTiO3 (our work)	309.18	120.28	96.33	132.87	84.06	121.29
PbTiO3 (theory)^22	313.1	110.71	92.40	197.0	82.72	85.41
PbTiO3 (expt.)^23	235.00	105.00	65.10	104.00	—	—
F@O-(PO layer)	306.01	126.10	95.26	132.40	88.47	112.68
F@O-(TO layer)	314.49	122.05	102.54	138.07	86.67	106.71

---

Various mechanical parameters, such as the Young modulus (Y), bulk modulus (B), shear modulus (G), Poisson ratio (ν), Pugh ratio (B/G), Frantsevich ratio (G/B), tetragonal shear modulus (C′), Kleinman parameter (ζ), machinability index , linear compressibility (β_a/β_c), volumetric compressibility (β_V), and Vicker hardness (H_V), have been calculated utilizing the above calculated elastic constants. These factors play a vital role in analyzing the strength, elasticity, plasticity, and rigidity of materials. B indicates the rigidity against volumetric deformation. Y determines the stiffness of a material under compression/tension⁸⁰ and also measures the resistance to thermal shocks. A higher Y value indicates a more covalent nature.⁸¹ G measures the capacity to resist the plastic deformation and strong directional bonding between atoms, reflected by a large value of the shear modulus.⁸² Based on the given Voigt, Reuss, and Hill approximation equations,⁸³ the calculated values of B/G/Y of the prst./F@O-dop.PO/TO-layer-based structures are 128.14/128.37/133.13, 81.88/81.04/79.99 and 202.50/200.85/199.93 GPa, respectively, which indicate the hardness and stiffness of our reported materials as mentioned in Table 1S of the SI. Our calculated results are comparable to the other theoretically reported B/G/Y values of 152.67/84.7/213.7 for the prst. PTO structure.⁸⁴ The calculated values of G/B/Y for the prst. and F@O-dop.PO/TO-layer-based systems obey the sequence of Y > B > G, indicating that these materials possess high stiffness and less elasticity and offer less resistance to shape deformation compared to volumetric deformation. Xifan et al. ⁸⁵ reveals that a material exhibits either ductile or brittle character if the B to G ratio exceeds or is less than 1.75. Thus, B to G ratios of 1.56/1.58/1.66 suggest a brittle nature of the prst. and F@O-dop.PO/TO-layer-based structures (see Table 2S of SI).

9B_V = C₁₁ + C₂₂ + C₃₃ + 2(C₁₂ + C₁₃ + C₂₃) (10)

15G_V = (C₁₁ + C₂₂ + C₃₃) − (C₁₂ + C₁₃ + C₂₃) + 3(C₄₄ + C₅₅ + C₆₆) (12)

15 G_R = 4(S₁₁ + S₂₂ + S₃₃) − 4(S₁₂ + S₁₃ + S₂₃) + 3(S₄₄ + S₅₅ + S₆₆) (13)

Additionally, ν generally varies between −1 and 0.5, and materials possessing a value less than 0.26 are brittle; those having a value greater than 0.26 are ductile. Hence, computed values of 0.24/0.24/0.25 for the prst. and F@O-dop.PO/TO-layer-based structures further confirm the stability of the brittle nature (see Table 2S of the SI). It also provides information about the bonding nature, as ν = 0.25 indicates ionic, lower than 0.25 covalent, and higher than 0.25 a metallic nature.⁸⁶ The covalent bonding character is demonstrated by the calculated values of 0.24/0.24/0.25 for the prst. and F@O-dop.PO/TO-layer-based structures (see Table 2S of the SI). The Frantsevich ratio (G/B), additionally, provides details regarding a material’s ductile/brittle nature. If G/B is lower than 0.57, the material is considered ductile, while G/B greater than 0.57 indicates a brittle nature.⁷⁸ Hence, computed G/B values of 0.64/0.63/0.60 for the prst. and F@O-dop.PO/TO-layer-based structures also further confirm the brittle nature (see Table 2S of the SI). To evaluate the dynamical stability of a crystalline solid, tetragonal shear modulus is a key parameter.^87,88 It also indicates the crystal stiffness, has an association with slowed transversal sound waves and is essential for modification of a system.⁸⁹ The crystal is dynamically stable when it shows a positive value, while it is dynamically unstable if it possesses a negative value. The calculated positive values of 94.45/89.96/96.22 GPa for the prst. and F@O-dop.PO/TO-layer-based structures suggest dynamic stability of the structures (see Table 2S of the SI). The internal-strain variable, frequently referred to as the Kleinman parameter (ζ),⁹⁰ determines a material’s stability toward stretching and bending. This dimensionless parameter usually exhibits values between 0 and 1. It defines how easy bond bending is compared to bond stretching. Significant contributions from bond stretching and bond bending are demonstrated by this variable’s lower and upper ranges, respectively. Consequently, a ζ close to zero shows that bond stretching dominates under stress, while bond bending is dominant at values close to 1. The values of 0.53/0.55/0.53 for the prst. and F@O-dop.PO/TO-layer-based structures, calculated by using eqn (17), are listed in Table 2S of the SI, suggesting that the rigidity of the materials is more affected by bond bending than bond stretching or contraction.

The machinability, , of materials has gained a lot of attention in today's industries, as it relates to several factors, such as the intrinsic features of the material being worked, the nature of the cutting-tool material, the geometry of the tool material, the conditions and type of cutting, and the rigidity and capacity of the machine tool. It is calculated using eqn (18),⁹¹ which may also be utilized to assess the plasticity of a material⁹² and lubrication capabilities. Materials having a small C₄₄ and a large B/C₄₄ value are associated with excellent dry lubrication qualities, smaller input forces, less friction, and large values of plastic strain. The calculated values of 1.52/1.45/1.53 for the prst. and F@O-dop.PO/TO-layer-based structures suggest a high degree of machineability (see Table 2S of the SI). The response of a material or structure against pressure, either negative (tensile) or positive (compressive), is referred to as compressibility. Linear compressibility defines the variation along a specific axis and volumetric compressibility describes the change in volume against applied pressure.⁹³ The values of linear compressibility along the a and c axes⁹⁴ and volumetric compressibility for the prst. and F@O-dop.PO/TO-layer-based structures, calculated by using eqn (19) and (20), are listed in Table 2S of the SI. The compressibility along the a-axis is smaller compared to that along the c-axis in the prst. and F@O-dop.PO/TO-layer-based structures and the F@O-doped PO is more compressible in comparison to the other materials.

In terms of applications, hardness is a basic mechanical parameter. It is characterized as the material’s ability to resist deformation. We explore the hardness of the prst. and F@O-dop.PO/TO-layer-based structures by employing the semi-empirical equations of hardness suggested by Chen et al.⁹⁵ (eqn (21)) and Tian et al.⁹⁶ (eqn (22)), where k represents G/B, and G and B represent the shear and bulk moduli, respectively. The obtained results are provided in Table 3S of the SI, where the F@O-dop.TO-layer-based system has the lowest values and the prst. PTO structure possesses the highest one, which indicates that the prst. material exhibits higher hardness as compared to the dop. one.

H_V = 2(k²G)^0.585 − 3 (21)

H_V = 0.92k^1.137G^0.708 (22)

The directional dependency of a material’s mechanical properties is often characterized by elastic anisotropy. Anisotropic elastic properties are crucial, as they are related to multiple physical features, including the development of plastic deformations, the motion of cracks, and the production of microcracks in solids. The universal anisotropy parameter (A^U) indicates if the material’s physical properties change depending on the direction. Materials with A^U = 0 are referred to as isotropic. A material’s anisotropy increases as the value of A^U varies from zero. The calculated A^U values from eqn (23) for the prst. and F@O-dop.PO/TO-layer-based structures are 0.86 and 0.79/0.77, respectively (see Table 3S of the SI). In tetragonal phases, a comparatively small elastic anisotropy behavior is observed as these values deviate slightly from 0.⁹⁷ Furthermore, we employ the bulk and shear moduli to figure out the anisotropic nature in terms of anisotropy indices. A_B and A_G equal to 0 denote absolute elastic isotropy, and A_B and A_G equal to 1 represent the optimal anisotropic state. The A_B/A_G values for the tetragonal prst., F@O-dop.PO/TO-layer-based structures are 0.030/0.073, 0.036/0.067, and 0.026/0.067, respectively, which are quite far from 1, further showing slight shear and compression anisotropy for the phases (see Table 3S of the SI). The Zener anisotropy factor (A^eq)⁹⁸ possesses a value of 1 for an isotropic crystal. The calculated values of 2.28 and 2.21/2.19 for the prst. and F@O-doped PO/TO structures indicate that these materials have anisotropic characteristics (see Table 3S of the SI).

Additionally, to investigate the anisotropy, three-dimensional illustrations of the elastic features, including the Young modulus, linear compressibility, shear modulus, and Poisson ratio, are generated by utilizing the ELATE code.⁹⁹ Comprehensive knowledge about the material’s mechanical properties in multiple directions is obtained using these 3D models. Isotropy is represented by a spherical shape, whereas deviation from a sphere indicates anisotropic phenomena in these models. Along with the anisotropy factor (A^U), the 3D visualizations provide the highest and lowest values of the elastic features, which are provided in Table 4S of the SI and shown in Fig. 2. Utilizing the data from Table 4S of the SI and the 3D images, it is concluded that the Young modulus, linear compressibility, shear modulus, and Poisson ratio demonstrate anisotropic phenomena.

Fig. 2 Three-dimensional representations of the various anisotropic elastic traits of the pristine and F@O-doped-PbO (PO)/TiO₂ (TO)-layer-based PbTiO₃ structures.

Sound velocity in a material is another crucial parameter that affects its electrical and thermal conductivity. In the fields of physics, materials science, geology, seismology, musical instrument design, and medicine, the acoustic behavior of materials has received significant attention. The longitudinal, transverse, and average sound velocities (v_l/v_t/v_avg.) are calculated using mechanical parameters like B and G:¹⁰⁰

where ρ indicates the mass density. The calculated values of v_l, v_t, and v_avg. for the prst./F@O-dop.PO/TO-layer-based structures are listed in Table 5S, 6S, and 7S of the SI, respectively, in the Voigt, Reuss, and Hill approximations. A material either possesses identical or a different acoustic impedance in comparison to its surroundings, which is important for transducer design, reducing noise in aircraft engines, manufacturing operations, and underwater acoustic applications. The difference in acoustic impedance between two materials at the interface determines the amount of sound energy that is reflected and transmitted while sound moves between them. Consequently, the majority of the sound is transmitted if the two impedances are roughly comparable, while the majority of the sound is reflected if they differ greatly. The acoustic impedance (Z) is calculated as:

It is clear that materials having large values of density and shear modulus possess large Z. The calculated values of Z listed in Tables 5S/6S/7S of the SI for the structures show that the prst. structure has the highest Z value as compared to the dop. structures. The radiation factor (I) can be calculated as:

The computed values of I listed in Tables 5S/6S/7S of the SI for the structures indicate that the prst. structure has the highest I value as compared to the dop. structures. The Debye temp. (Θ_D), the most significant solid thermophysical parameter, is a crucial measure for investigating the vibrational characteristics of a material and is strongly correlated with physical features, including the melting point, specific heat, elastic constants, thermal conductivity, and lattice vibration. Usually, it directly correlates with strong interatomic bonds, smaller average atomic masses, greater melting temp., higher hardness, and greater mechanical wave velocities. Therefore, the Θ_D is determined employing the v_avg. from Table 5S/6S/7S of the SI in the Voigt, Reuss, and Hill approximations:¹⁰⁰

where ρ indicates the mass density, n is the number of atoms in the unit cell, N_A is the Avogadro constant, h and k_B are the Planck and Boltzmann constants, and M is the molecular mass. The calculated 456 and 451/454 K values of Θ_D for the prst. and F@O-dop.PO/TO-layer-based structures are listed in Table 8S of the SI. The prst. system exhibits the highest value, which indicates that the prst. material is harder than the dop. materials, which is in excellent agreement with the calculated hardness parameters listed in Table 3S of the SI. Another crucial dimensionless thermophysical variable that evaluates lattice anharmonicity is the Grüneisen parameter (γ). Higher values of γ represent a greater level of anharmonicity. It contains specific details regarding the temp. influence on phonon frequencies and damping. Additionally, it facilitates investigation of thermal expansion impacts in a crystal. The calculated ν values listed in Table 2S of the SI are employed to determine γ:¹⁰¹

The computed values of γ, 1.44 and 1.45/1.50 for the prst. and F@O-dop.PO/TO-layer-based structures, are listed in Table 8S of the SI. This indicates that the F@O-dop.TO-layer-based structure exhibits the highest degree of anharmonicity compared to the prst. and F@O-dop.PO-layer-based structure. A higher anharmonicity value improves the phonon–phonon interactions and ultimately reduces the lattice thermal conductivity (κ_L) for the F@O-dop.TO-layer-based structure in comparison to the prst. and F@O-dop.PO-layer-based structures. These findings coincide with our calculated κ_L (discussed later). Investigating the melting temp. of novel materials that can be utilized at various temps. is a fascinating and crucial field of study in theoretical research. Greater bonding energy and reduced thermal expansion are characteristics of materials having higher melting temp. We calculate the melting temp. (T_m) as:¹⁰²

The expected melting temp. values for the prst. and F@O-dop.PO/TO-layer-based structures are 1480 and 1470/1504 K, respectively, demonstrating that these are excellent candidate materials for applications at higher temp. (see Table 8S of the SI). Several physical characteristics, such as specific heat, thermal conductivity (κ), E_g variation with temp., and electron effective mass are related to thermal expansion of materials. Therefore, the thermal expansion coefficient (α) is a crucial factor for the epitaxial development of crystals, which is calculated as:

The determined α values of 1.95 × 10⁻⁵ and 1.97 × 10⁻⁵/2 × 10⁻⁵ K⁻¹ for the prst. and F@O-dop.PO/TO-layer-based structures show that the prst. structure has the lowest value, which indicates its dimensional stability with respect to temp. as compared to the other dop. materials (see Table 8S of the SI). Another crucial thermodynamic variable is heat capacity (C_P), which is calculated as:

where the number of atoms in a unit volume is symbolized by 1/Ω. It generally takes more energy for the system to attain an equilibrium state when the C_P is greater. The calculated values of ρC_P for the prst. and F@O-dop.PO/TO-layer-based structures are listed in Table 8S of the SI. Phonons are crucial for a number of physical characteristics, like electrical conductivity (σ), κ, thermopower, and C_P. The precise position of the peak in the phonon distribution curve is identified as the dominant phonon wavelength (λ_dom.). A material having a higher v_avg., smaller ρ, and larger G exhibits a longer λ_dom.. The λ_dom.¹⁰³ is calculated as:

The calculated values of λ_dom. at 300 K, 149.63 and 148.01/149.01 for the prst. and F@O-dop.PO/TO-layer-based PbTiO₃ structures, respectively, show that the prst. structure has the highest λ_dom. value.

Subsequently, the FE features of the prst. and F@O-dop. PTO systems are examined on the basis of the magnitudes of P. This refers to the amount of dipole moment in a unit volume along a particular axis, which arises from broken symmetry occurring due to the corresponding shifts of atoms from their original places. It is calculated as:

The estimated O–Ti–O bond angles and their associated P in the prst. as well as F@O-dop.PO/TO-layer-based structures are shown in Fig. 3. Based on our findings, it is revealed that P significantly depends upon the O–Ti–O bond angles, which further define the structural distortions. The determined O–Ti–O bond angle of 160.9° that is illustrated in Fig. 3(a) for the prst. structure leads to a significant P of 88 μC cm⁻², as shown in Fig. 3(b), which closely matches the experimentally obtained values ranging from 75 to 92 μC cm⁻² (ref. 38, 39 and 104). This is also consistent with the theoretical results based on the polar phonon mode, the Berry phase theory for the PAW (projector-augmented wave), USPP (ultra-soft pseudo-potentials), and LDA (local density approximation), which are 75, 78.6, 81.2 and 85 μC cm⁻², respectively.^38,39

Fig. 3 Calculated (a) ∠O–Ti–O average bond angle and (b) spontaneous polarization (P) in the pristine and F@O-doped-PbO (PO)/TiO₂ (TO)-layer-based PbTiO₃ structures.

For the F@O-dop. structures, P is considerably smaller as compared to that of the prist. one because of a significant reduction in structure deformities. Fig. 3(b) illustrates that the F@O-dop. system for the PO layer keeps a P of 42.56 μC cm⁻² and Fig. 3(a) shows a ∠O–Ti–O bond angle of 161.5°. Likewise, the F@O-dop.TO layer exhibits a ∠O–Ti–O bond angle of 161.3° along with a P of 42.78 μC cm⁻². Hence, it is concluded that as the ∠O–Ti–O bond angle value increases (moving toward 180°), structural distortions decrease, and hence so does P. Previously, it was revealed that the tetragonal behavior of the PTO structure decreases when F is substituted at the O site as the ratio of the dop. structure declines to 1.056, in comparison to 1.064 (prst. one).¹⁰⁵ The higher electro-negativity of F in comparison to O explains the decline in tetragonal behavior as it promotes the establishment of ionic interaction between Pb/Ti and F, with a corresponding decrease in covalent interactions between the cations and O.¹⁰⁶ Ferroelectricity and polarized deformation in FE perovskites depend on the hybridization of the cations at the B positions with anions O.¹⁰⁷ It also reduces the distortion of the F-connected Ti ions and lessens the hybridization between Ti and O,¹⁰⁵ which negatively influence the P magnitude.

Now, to analyze the electronic aspects of the systems, we display the calculated SP TDOS and band structures. In Fig. 4(a), the prst. one demonstrates insulating behavior with an E_g of 2.11 eV between the valence band (VB) and conduction band (CB) edges in both the spin-majority/spin-minority channels (N^↑/N^↓), which is an underestimate compared to the experimentally reported value of 3.88 eV.¹⁰⁸ This is a well-known problem of exchange-correlation functionals. However, our determined E_g is close to another theoretically calculated E_g of 2.02 eV for the prst. PTO system, determined via the PBE-GGA method employing CASTEP.¹⁰⁹ Furthermore, the insulating character of the prst. motif is validated by visualizing the SP band structures in Fig. 1S of the SI, which display an indirect E_g of 2.11 eV between the X and Γ symmetry points. In addition, TDOS degeneracy verifies the non-magnetic behavior of the system.

Fig. 4 GGA + U calculated spin-polarized total density of states (TDOS) of the (a) pristine and (b and c) F@O-doped-PbO (PO)/TiO₂ (TO)-layer-based PbTiO₃ structures.

Next, to investigate the F@O-doping impact on the electronic structure of the PTO, we plot the SP TDOS for both doped PO/TO-layer-based systems in Fig. 4(b) and (c). For the PO layer, it is found that certain states are observed to be crossing the Fermi level (E_F) from the CB to VB in N^↑ and there is a significant E_g of 2.49 eV within the VB, confirming the n-type conducting character of the corresponding channel, and there is a large E_g of 2.49 eV in N^↓ (see Fig. 4(b)). Overall, the F@O-dop.PO-layer-based system is termed an n-type HM FM. The large E_g in the non-metallic channel (N^↓) reduces reverse leakage current and ensures 100% SP at E_F in the conducting channel, making it an extremely desirable material for spintronics.^110,111 Furthermore, the calculated TDOS for the F@O-dop.TO-layer-based system is shown in Fig. 4(c), wherein a narrower/wider E_g of 0.23/2.41 eV exists in the N^↓/N^↑ channel, further supporting the SC character of the structure. F-substitution at the oxygen site provides one extra electron to the system because it has an electronic configuration of 2s²2p⁵ and oxygen has 2s²2p⁴. This extra electron resides at the Ti site, which changes its valence from +4 to +3 and the system exhibits metallicity as well magnetism due to partially filled Ti 3d states. In earlier studies, it was demonstrated that the substitution of oxygen by F transforms Ti⁴⁺ to Ti³⁺ (3d¹).^112,113 Those unpaired 3d electrons (3d¹) are responsible for the observed ferromagnetism in F@O-dop. PTO structures. In the same way, F@O substitution yields FM in the BTO structure.¹¹⁴ Likewise, a Ti⁴⁺ to Ti³⁺ transformation is observed in response to F@O substitution in PTO¹⁰⁵ and SrTiO₃.¹¹⁵ Hence, a similar phenomenon is observed in the present case.

To better understand the role of atomic orbitals in the HM electronic state of the F@O-dop.PO-layer-based system, we present the 3d orbital resolved PDOS in Fig. 5. It has been found that each Ti atom shows a contribution to the electronic state around the E_F, except for Ti5 and Ti6 (located at the CB edge (CBE)), and the Ti4 atom exhibits a significant contribution around E_F (see Fig. 5(a)). Therefore, we plotted the d_z²/d_x²−y²/d_yz+xz/d_xy sub-states of the Ti4 ion. It is evident that the d_z²/d_x²−y² states reside at the CBE and do not contribute to the metallicity (see Fig. 5(b)). However, d_yz+xz/d_xy states significantly grow at E_F, resulting in the conducting behavior in N^↑. Similarly, to get improved knowledge about the role of Ti ions in the SC electronic state of the F@O-dop.TO-layer-based system, we displayed the 3d orbital resolved PDOS in Fig. 2S of the SI. It has been determined that Ti3-4 states in the CB grow significantly toward E_F and defect states are induced by doping at the VBE due to Ti1-2 and Ti3-4 atoms in N^↓, resulting in E_g reduction and leading to the system displaying a SC nature (see Fig. 2S(a) of the SI). Consequently, we additionally displayed the d_z²/d_x²−y²/d_xy/d_xz/d_yz sub-states of the Ti ion. It is apparent that the d_xz states mainly contribute at the VBE and expand toward E_F in the CB (see Fig. 2S(b) of the SI). However, the d_z²/d_x²−y²/d_xy/d_yz states reside behind and do not contribute to the SC behavior in N^↓.

Fig. 5 GGA + U calculated spin-polarized partial density of states (PDOS) for the (a) 3d states of the Ti-ions and (b) orbital resolved Ti4-3d sub-states in the F@O-dop.PO-layer-based PbTiO₃ structure.

In addition, to confirm the n-type HM/SC FM state for the F@O-dop.PO/TO-layer-based system, we present the SP band structures in Fig. 6. From Fig. 6(a), specific states in N^↑ are found to be crossing the E_F from the CB to VB in the F@O-dop.PO-layer-based system, and a significant E_g within the VB verifies the n-type metallic character of the channel. However, no states are passing through the E_F in N^↓ and there is a large E_g of 2.49 eV, as presented Fig. 6(a′), which coincides with the TDOS that appears in Fig. 4(b). For the F@O-dop.TO-layer-based system, a large E_g of 2.41 eV appears in N^↑ (see Fig. 6(b)), and a smaller one of 0.23 eV lies in N^↓, as shown in Fig. 6(b′). This confirms the SC character of the motif, which is identical to the TDOS presented in Fig. 4(c). Hence, the n-type HM/SC FM state of the F@O-dop.PO/TO-layer-based system is verified via the SP TDOS and band-structure investigations.

Fig. 6 GGA + U calculated spin-polarized band structures of the spin-majority/spin-minority channel for the (a/a′) F@O-doped PbO(PO)-layer and (b/b′) F@O-doped TiO₂(TO)-layer-based PbTiO₃ structures.

Furthermore, non-degeneracy in the TDOS demonstrates the magnetic character for both doped systems. The calculated m_t in the F@O-dop.PO-layer-based motif is 1.00 μ_B, where a substantial contribution comes from the interstitial sites and the Ti atoms have m_s values of 0.35 and 0.16 μ_B, respectively. Hence, the integral m_t of the system ultimately indicates an HM FM character. Similarly, the m_t of the F@O-dop.TO-layer-based system is 1.00 μ_B, where the interstitial sites and Ti atoms play a substantial role with m_s values of 0.31 and 0.36 μ_B, respectively. Even though the system has an integral m_t, however, it exhibits a SC character. Additionally, Fig. 7 shows 3D spin magnetization density iso-surfaces for the highly prominent F@O-dop.PO-layer-based system, having an iso-value of ±0.001 e Å⁻³. It is apparent that significant density magnitude appears only around the Ti ions, which further confirms that these ions contribute the most to the m_t, and a minor proportion comes from the oxygen ions. Furthermore, it qualitatively validates the m_t related to the Ti ions.

Fig. 7 GGA + U computed three-dimensional spin-magnetization iso-surfaces for the F@O-dop.PO-layer-based PbTiO₃ structure.

Now, we discuss the TE behavior of the prst. and F@O-dop.PO/TO-layer-based structures. Based on their prolonged functional lifetime, energy conservation, thermal recycling, and environment-friendly characteristics, TE materials have received a lot of attention.^116,117 Without harming the environment, these materials directly convert around 60% of waste heat into usable energy.^118,119 Reasonable size and lack of mechanical moving parts make them appropriate for converting heat into electrical energy in remote locations.¹²⁰ In this respect, POs are regarded as excellent TE materials due to having high chemical stability,¹²¹ σ, and Seebeck coefficients (S), along with low electronic thermal conductivity (κ_e),¹²² which are crucial for device applications such as power generation for deep-space detectors,¹²³ TE coolers,¹²⁴ and laser diodes.¹²⁵ To boost TE efficiency, materials possessing higher S values and lower κ are highly preferred.¹²⁶ Therefore, to assess the transformation of heat into electrical energy, we estimated multiple TE variables in the prst. and F@O-dop.PO/TO-layer-based systems within a constant relaxation time (τ) of 10⁻¹⁴ s, employing the BoltzTrap computational code.¹²⁷ It is generally believed that S, σ, and κ_e, which are additionally associated with the carrier concentration, are the main factors that determine a material’s TE responses.¹²⁸ The computed essential TE variables, including S, , κ, power factor (PF), and dimensionless ZT with respect to temp. within the 200–700 K range are plotted in Fig. 8.

Fig. 8 Multiple thermoelectric variables against temperature for the pure and F@O-doped-PbO/TiO₂ (PO/TO)-layer-based PbTiO₃ structures: (a) σ/τ, (b) S, (c) σS²/τ, (d) κ_e/τ, (e) κ_L and (f) ZT.

The fluctuation in with temp. is displayed in Fig. 8(a). At room temp., the prst./F@O-dop.PO/TO-layer-based systems achieve values of 0.9/46.5/1.1 × 10¹⁸ (Ω ms)⁻¹, which increase for the prst./F@O-dop.TO-layer-based system as the temp. rises and reach 3.2/2.3 × 10¹⁸ (Ω ms)⁻¹ at 700 K. This could occur as a result of bond breaking and a rise in electron motion,¹²⁶ enhancing the material’s propensity for transferring charge carriers across the VB and CB. It might be argued that higher temp. promotes further conductive movements of charge carriers from the VB into the CB, which enhances σ. For the F@O-dop.PO-layer-based system, it decreases gradually and reaches 40.84 × 10¹⁸ (Ω ms)⁻¹ at 700 K. Another important parameter, , commonly referred to as thermopower, couples the variation in potential ΔV established across the terminals of the material with the temp. difference ΔT (ref. 80) and is plotted in Fig. 8(b). It provides insights on the material’s charge carriers and reveals whether they are mostly holes (positive S) or electrons (negative S).¹²² At room temp., the prst./F@O-dop.TO-layer based system exhibits values of 257/197 × 10⁻⁶ V K⁻¹, which gradually decrease as the temp. goes higher and reach 243/88 × 10⁻⁶ V K⁻¹ at 700 K. In comparison, the F@O-dop.PO-layer-based material demonstrates negative values of S from 400–700 K, indicating that electrons are the majority charge carriers.¹²⁹ Additionally, S has a value of 6.06 × 10⁻⁶ V K⁻¹ at room temp. and becomes −21.18 × 10⁻⁶ V K⁻¹ with the rise in temp. to 400 K, then smoothly decreases to −105.85 × 10⁻⁶ V K⁻¹ at 700 K. More negative S values at elevated temp. compared to ambient temp. show that the carrier concentration becomes higher when the temp. increases. Consequently, even at higher temp., the system shows an improved TE conversion ratio.

Next, , which indicates how effectively a substance transforms thermal waste into energy,¹³⁰ is shown in Fig. 8(c). The values of the PF are 5.9/0.17/4.4 × 10¹⁰ W K⁻² (ms)⁻¹ for the prst./F@O-dop.PO/TO-layer-based systems at room temp. These directly correlate with the rise in temp. and climb to 18.6/45.76 × 10¹⁰ W K⁻² (ms)⁻¹ at 700 K in the prst./F@O-dop.PO-layer-based systems. In contrast, the PF slightly rises to 4.4 × 10¹⁰ W K⁻² (ms)⁻¹ at 400 K then declines until it reaches 1.8 × 10¹⁰ W K⁻² (ms)⁻¹ at 700 K in the F@O-dop.TO-layer-based system. Materials with substantial PF values are potentially desired in TE industries, since they are capable of greater effective conversion of thermal waste into electrical power.¹²⁶ Next, κ, characterized as κ_e + κ_L, where the first term is associated with carrier transportation while the next describes the phonon-driven temp. impact,¹³¹ is illustrated in Fig. 8(d and e). The prst./F@O-dop.PO/TO-layer-based systems attain values of 0.22/2.76/0.62 × 10¹⁴ W K⁻¹ at room temp. As the temp. rises, increases and climbs to 1.72/5.72/0.59 × 10¹⁴ W K⁻¹ for the prst./F@O-dop.PO/TO-layer-based systems (see Fig. 8(d)). This pattern shows that carrier mobility improves as temp. increases. The charge carriers gain more energy and participate actively in conduction as the temp. rises. The primary source of heat conduction in insulators and SCs is phonon-induced heat conduction. Although phonon-induced lattice vibrations and spin-polarized electrons both contribute to heat conduction in HMs, the phonon contribution generally dominates.¹³² Therefore, it is crucial to investigate the temporally dependent modifications in κ_L for the prst. and F@O-dop.PO/TO-layer-based PTO systems, which are displayed in Fig. 8(e) utilizing the Slack model:⁸⁰

In this case, γ represents the Grüneisen parameter, M_avg. is the average molar mass per atom, V is volume per atom, T is temp. in Kelvin, and A is a constant that depends on γ, which further depends on ν:

At room temp., the prst./F@O-dop.PO/TO-layer-based systems have κ_L values of 0.23/0.056/0.053 W m⁻¹ K⁻¹, which continuously decrease to 0.10/0.024/0.02 W m⁻¹ K⁻¹ in the prst./F@O-dop.PO/TO-layer-based systems as the temp. increases to 700 K. Improved anharmonic phonon scattering, which is also known as the phonon rattling effect,¹³³ is the main cause of the κ_L decline at higher temp. Compared to the prst. structure, the F@O substitution produces more ionic and shorter Ti–O/F bonds, which strengthens the anharmonic scattering phenomenon.¹⁰⁵ In addition, F (1.33 Å) and O (1.40 Å) have distinct ionic radii, which distorts the Ti–O and Ti–F bonds in the F@O-dop. structures.¹¹⁴ These distortions enhance phonon scattering possibilities with the transmission of heat. Moreover, F@O substitution adds more electrons to the dop. structures, which may enhance electron–phonon interactions, reduce the mean free path and improve phonon scattering.¹¹⁵ Along with this, the atomic mass of O (16 amu) is close to that of F (19 amu); this slight mass variation causes a mass imbalance, which also facilitates the mass fluctuation scattering mechanism¹³⁴ and enhances phonon scattering. Also, it is found that the κ_L of the prst. and F@O-dop.PO/TO-layer-based structures reduces gradually with rising temp., as improved multifaceted scattering phenomena based upon point defects, crystal deformations, and mass variation mainly contribute with the increase in temp.¹³⁴

Finally, ,¹²³ presented in Fig. 8(f), immediately correlates the product of the PF and estimates the material energy conversion efficiency. Materials with greater PF and smaller κ exhibit greater ZT, which ultimately improves the energy utilization effectiveness of the material. The prst./F@O-dop.PO/TO-layer-based systems possess ZT values of 0.15/0.026/0.39 at 200 K. For the prst./F@O-dop.PO material, it increases with the rise in temp. and attains values of 0.72/0.56 at 700 K. For the F@O-dop.TO one, it initially increases to 0.49 at 400 K and then gradually decreases and reaches 0.21 at 700 K. This suggests that these are appropriate materials for TE activity at elevated temp.

4 Conclusion

In conclusion, first-principles calculations are employed to investigate the thermodynamic, mechanical, ferroelectric, electronic, magnetic, and thermoelectric (TE) properties of the pristine (prst.) and F@O-doped (dop.)-PbO (PO)/TiO₂ (TO)-layer-based PbTiO₃ systems. The negative formation enthalpy of the structures confirms the thermodynamic stability and exothermic reaction. The calculated elastic constants of the structures verify the mechanical stability, brittle nature, and the existence of covalent bonding between atoms. Various thermophysical parameters, such as the Debye temperature, Grüneisen parameter, thermal expansion coefficient, melting temperature, and dominant phonon mode also indicate the hardness of the structures. The prst. system displays a non-magnetic insulating nature with an energy-gap (E_g) of 2.11 eV, having spontaneous polarization (P) of 88 μC cm⁻². Surprisingly, the F@O-dop.PO material displays a ferromagnetic (FM) half-metallic behavior with a P of 42.56 μC cm⁻², and the F@O-dop.TO material is a FM semiconductor having an E_g = 0.23 and P = 42.78 μC cm⁻². The reduced P value in the doped structure is due to the reduction of the structural distortions. Remarkably, the prst. and F@O-dop.PO/TO-layer-based structures demonstrate the highest figures of merit of 0.72 and 0.56/0.49 at 700 and 700/400 K, respectively, including lattice-thermal conductivity. This indicates that these are effective candidates for TE applications at higher temperatures.

Author contributions

S. Haider: writing – original draft, investigations, formal analysis, data curation. Bassem F. Felemban: validation, visualization, formal analysis. Hafiz Tauqeer Ali: writing, investigations, and formal analysis. S. Nazir: writing – review and editing, validation, supervision, project administration, conceptualization.

Conflicts of interest

The authors declare no competing interests.

Data availability

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ra05969k.

Acknowledgements

The computational work was supported by the University of Sargodha, Sargodha, Pakistan. The authors extend their appreciation to Taif University, Saudi Arabia, for supporting this work through project number (TU-DSPP-2024-133).

Notes and references

1. K. F. Wang, J. M. Liu and Z. F. Ren, Adv. Phys., 2009, 58, 321–448.
2. W. Eerenstein, N. Mathur and J. F. Scott, Nature, 2006, 442, 759–765.
3. S. Picozzi and C. Ederer, J. Phys.: Condens. Matter, 2009, 21, 303201.
4. D. J. Lovinger, E. Zoghlin, P. Kissin, G. Ahn, K. Ahadi, P. Kim, M. Poore, S. Stemmer, S. J. Moon, S. D. Wilson and R. D. Averitt, Phys. Rev. B, 2020, 102, 085138.
5. R. S. Silva Jr, J. E. Rodrigues, J. Gainza, F. Serrano-Sánchez, L. Martínez, Y. Huttel, J. L. Martínez and J. A. Alonso, Inorg. Chem., 2024, 63, 7007–7018.
6. C. R. Bowen, H. A. Kim, P. M. Weaver and S. Dunn, Energy Environ. Sci., 2014, 7, 25–44.
7. C. W. Nan and J. M. Liu, Natl. Sci. Rev., 2019, 6, 620.
8. P. Bongurala and V. Gorige, J. Magn. Magn. Mater., 2019, 477, 350–355.
9. R. Ramesh and N. A. Spaldin, Nat. Mater., 2007, 6, 21–29.
10. C. Lu, M. Wu, L. Lin and J.-M. Liu, Natl. Sci. Rev., 2019, 6, 653–668.
11. H. Palneedi, V. Annapureddy, S. Priya and J. Ryu, Actuators, 2016, 5, 9.
12. N. Ortega, A. Kumar, J. F. Scott and R. S. Katiyar, J. Phys.: Condens. Matter, 2015, 27, 504002.
13. M. Kumar, S. Shankar, A. Kumar, A. Anshul, M. Jayasimhadri and O. Thakur, J. Mater. Sci.: Mater. Electron., 2020, 31, 19487–19510.
14. M. M. Vopson, Crit. Rev. Solid State Mater. Sci., 2015, 40, 223–250.
15. S. Sharma, M. Kumar, A. Laref, J. M. Siqueiros and O. R. Herrera, Mater. Lett., 2023, 331, 133490.
16. Y. Zhang, Z. Li, C. Deng, J. Ma, Y. Lin and C.-W. Nan, Appl. Phys. Lett., 2008, 92, 152510.
17. D. Khomskii, Physics, 2009, 2, 20.
18. W. Prellier, M. Singh and P. Murugavel, J. Phys.: Condens. Matter, 2005, 17, R803.
19. E. Ascher, H. Rieder, H. Schmid and H. Stössel, J. Appl. Phys., 1966, 37, 1404–1405.
20. S. NV, K. Vinayakumar and K. Nagaraja, Coatings, 2020, 10, 1221.
21. Y. J. Pu, Z. C. Huang, H. C. Xu, D. F. Xu, Q. Song, C. H. P. Wen, R. Peng and D. L. Feng, Phys. Rev. B, 2016, 94, 115146.
22. S. Kuma and M. M. Woldemariam, Adv. Condens. Matter Phys., 2019, 2019, 3176148.
23. A. Hachemi, H. Hachemi, A. Ferhat-Hamida and L. Louail, Phys. Scr., 2010, 82, 025602.
24. J.-G. Park, M. D. Le, J. Jeong and S. Lee, J. Phys.: Condens. Matter, 2014, 26, 433202.
25. S. M. Selbach, M.-A. Einarsrud and T. Grande, J. Mater. Chem., 2009, 21, 169–173.
26. H. Liu and X. Yang, Ferroelectrics, 2017, 507, 69–85.
27. J. Wang, J. Neaton, H. Zheng, V. Nagarajan, S. Ogale, B. Liu, D. Viehland, V. Vaithyanathan, D. Schlom and U. Waghmare, et al., Science, 2003, 299, 1719–1722.
28. T. Shishidou, N. Mikamo, Y. Uratani, F. Ishii and T. Oguchi, J. Phys.: Condens. Matter, 2004, 16, S5677.
29. O. Heyer, N. Hollmann, I. Klassen, S. Jodlauk, L. Bohaty, P. Becker, J. A. Mydosh, T. Lorenz and D. Khomskii, J. Phys.: Condens. Matter, 2006, 18, L471.
30. C. Lu and J.-M. Liu, J. Materiomics, 2016, 2, 213–224.
31. C. Rao, A. Sundaresan and R. Saha, J. Phys. Chem. Lett., 2012, 3, 2237–2246.
32. S. Sahoo, P. Mahapatra and R. Choudhary, J. Phys. D: Appl. Phys., 2015, 49, 035302.
33. A. Panchwanee, V. R. Reddy, A. Gupta and V. Sathe, Mater. Chem. Phys., 2017, 196, 205–212.
34. S. Dong, J.-M. Liu and E. Dagotto, Phys. Rev. Lett., 2014, 113, 187204.
35. J. Z. Cong, S. P. Shen, Y. S. Chai, L. Q. Yan, D. S. Shang, S. G. Wang and Y. Sun, J. Appl. Phys., 2015, 117, 174102.
36. J. M. Wesselinowa and S. Kovachev, J. Appl. Phys., 2007, 102, 043911.
37. M. Balli, B. Roberge, P. Fournier and S. Jandl, Crystals, 2017, 7, 44.
38. T. Abe, S. Kim, C. Moriyoshi, Y. Kitanaka, Y. Noguchi, H. Tanaka and Y. Kuroiwa, Appl. Phys. Lett., 2020, 117, 252905.
39. T. Shimada, Y. Umeno and T. Kitamura, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 77, 094105.
40. H. S. Bhatti, S. T. Hussain, F. A. Khan and S. Hussain, Appl. Surf. Sci., 2016, 367, 291–306.
41. V. Zelezny, D. Chvostová, D. Šimek, F. Máca, J. Mašek, N. Setter and Y. H. Huang, J. Phys.: Condens. Matter, 2015, 28, 025501.
42. T. Shimada, S. Tomoda and T. Kitamura, J. Phys.: Condens. Matter, 2010, 22, 355901.
43. R. E. Cohen and H. Krakauer, Ferroelectrics, 1992, 136, 65–83.
44. J.-M. Oh, H.-J. Kim and S.-M. Nam, J. Korean Phys. Soc., 2010, 57, 1096–1101.
45. S. Hazra, T. Schwaigert, A. Ross, H. Lu, U. Saha, V. Trinquet, B. Akkopru-Akgun, B. Z. Gregory, A. Mangu and S. Sarker, et al., Adv. Mater., 2024, 36, 2408664.
46. S. Yan-Chao, Z. Wei-Xi, W. Ying, L. Kai-Wu and J. Xiao-Wei, Acta Phys. Sin., 2018, 67, 187701.
47. L. Delimova, E. Guschina, N. Zaitseva, S. Pavlov, D. Seregin, K. Vorotilov and A. Sigov, J. Mater. Res., 2017, 32, 1618–1627.
48. V. Bhide, K. Deshmukh and M. Hegde, Physica, 1962, 28, 871–876.
49. H. Jiao, X. Wang, S. Wu, Y. Chen, J. Chu and J. Wang, Appl. Phys. Rev., 2023, 10, 011310.
50. H. Kohlstedt, Y. Mustafa, A. Gerber, A. Petraru, M. Fitsilis, R. Meyer, U. Böttger and R. Waser, Microelectron. Eng., 2005, 80, 296–304.
51. B. Jiang, J. Peng, L. Bursill and W. Zhong, J. Appl. Phys., 2000, 87, 3462–3467.
52. A. Sobhani-Nasab, M. Rangraz-Jeddy, A. Avanes and M. Salavati-Niasari, J. Mater. Sci., 2015, 26, 9552–9560.
53. B. Behera, B. C. Sutar and N. R. Pradhan, Emergent Mater., 2021, 4, 847–863.
54. W. Zhou and A. Ariando, Jpn. J. Appl. Phys., 2020, 59, SI0802.
55. J. H. Lee and K. M. Rabe, Phys. Rev. Lett., 2010, 104, 207204.
56. M. Manikandan, K. Saravana Kumar and C. Venkateswaran, J. Appl. Phys., 2015, 118, 234105.
57. P. Pal, K. Rudrapal, S. Mahana, S. Yadav, T. Paramanik, S. Mishra, K. Singh, G. Sheet, D. Topwal and A. R. Chaudhuri, et al., Phys. Rev. B, 2020, 101, 064409.
58. O. R. González, F. S. D. Jesús, F. P. García, C. A. C. Escobedo, M. R. Cardona and A. M. B. Miró, Materials, 2019, 12, 2054.
59. I. S. Banu and S. D. Lakshmi, J. Mater. Sci.: Mater. Electron., 2017, 28, 16044–16052.
60. S. Gilani, S. Nazir, M. Zulfiqar and Y. Zhao, Phys. Scr., 2022, 97, 075816.
61. F. Craciun, F. Cordero, M. Cernea, V. Fruth, I. Atkinson, N. Stanica, B. Vasile, R. Trusca, A. Iuga and P. Galizia, et al., Ceram. Int., 2019, 45, 9390–9396.
62. F. Ricci and V. Fiorentini, J. Phys.: Conf. Ser., 2013, 470, 012013.
63. S. Nazir and A. Tariq, J. Appl. Phys., 2019, 125, 093905.
64. S. M. Selbach, T. Tybell, M.-A. Einarsrud and T. Grande, Chem. Mater., 2007, 19, 6478–6484.
65. L. Ju, T. Sabergharesou, K. G. Stamplecoskie, M. Hegde, T. Wang, N. A. Combe, H. Wu and P. V. Radovanovic, J. Am. Chem. Soc., 2012, 134, 1136–1146.
66. P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen and L. D. Mark, J. Chem. Phys., 2020, 152, 074101.
67. A. Dal Corso, A. Pasquarello, A. Baldereschi and R. Car, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 53, 1180.
68. R. Tesch and P. M. Kowalski, Phys. Rev. B, 2022, 105, 195153.
69. R. I. Eglitis and R. Jia, Materials, 2023, 16, 7623.
70. A. Glazer and S. Mabud, Acta Crystallogr., Sect. B, 1978, 34, 1065–1070.
71. L. Wang, T. Maxisch and G. Ceder, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 73, 195107.
72. H. Javed and S. Nazir, Mater. Chem. Phys., 2022, 289, 126469.
73. Z. Alahmed and H. Fu, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 76, 224101.
74. F.-F. Ge, W.-D. Wu, L.-H. Cao, X.-M. Wang, H.-P. Wang, Y. Dai, H.-B. Wang and J. Shen, Ferroelectrics, 2010, 401, 154–160.
75. H.-S. Lee, T. Mizoguchi, T. Yamamoto, S.-J. L. Kang and Y. Ikuhara, Acta Mater., 2007, 55, 6535–6540.
76. S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl and C. Wolverton, npj Comput. Mater., 2015, 1, 1–15.
77. M. Ziati and H. Ez-Zahraouy, J. Phys. Opt. Sci., 2021, 3, 2–8.
78. A. Bakar, A. O. Alrashdi, M. M. Fadhali, A. Afaq, H. Yakout and M. Asif, J. Mater. Res. Technol., 2022, 19, 4233–4241.
79. E. Karaca, E. Arslan, H. M. Tütüncü and G. P. Srivastava, Philos. Mag., 2017, 97, 1866–1883.
80. M. Y. Sofi, M. S. Khan and M. A. Khan, Mater. Adv., 2024, 5, 4913–4931.
81. M. Rajagopalan, S. P. Kumar and R. Anuthama, Phys. B, 2010, 405, 1817–1820.
82. S. Mitro, M. Hadi, F. Parvin, R. Majumder, S. Naqib and A. Islama, J. Mater. Res. Technol., 2021, 11, 1969–1981.
83. W. Xiao, Y. Yang, Z. Pi and F. Zhang, Coatings, 2022, 12, 73.
84. S. Gadisa and M. Woldemariam, Adv. Condens. Matter Phys., 2019, 2019, 1–12.
85. W. Xifan, D. Vanderbilt and D. Hamann, Phys. Rev. B:Condens. Matter Mater. Phys., 2005, 72, 035105.
86. M. Benkabou, H. Rached, A. Abdellaoui, D. Rached, R. Khenata, M. H. Elahmar, B. Abidri, N. Benkhettou and S. Bin-Omran, J. Alloys Compd., 2015, 647, 276–286.
87. M. Naher, M. Mahamudujjaman, A. Tasnim, R. Islam and S. Naqib, Solid State Sci., 2022, 131, 106947.
88. M. Ali, M. Hossain, M. Uddin, A. Islam, D. Jana and S. Naqib, J. Alloys Compd., 2021, 860, 158408.
89. J. Worgull, E. Petti and J. Trivisonno, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 15695.
90. L. Kleinman, Phys. Rev., 1962, 128, 2614.
91. Z. Sun, D. Music, R. Ahuja and J. M. Schneider, Phys. Rev. B:Condens. Matter Mater. Phys., 2005, 71, 193402.
92. M. Phasha, P. Ngoepe, H. Chauke, D. Pettifor and D. Nguyen-Mann, Intermetallics, 2010, 18, 2083–2089.
93. W. Miller, K. E. Evans and A. Marmier, Appl. Phys. Lett., 2015, 106, 231903.
94. M. Naher and S. Naqib, Sci. Rep., 2021, 11, 5592.
95. X.-Q. Chen, H. Niu, D. Li and Y. Li, Intermetallics, 2011, 19, 1275–1281.
96. Y. Tian, B. Xu and Z. Zhao, Int. J. Refract. Met. Hard Mater., 2012, 33, 93–106.
97. H. Dong, C. Chen, S. Wang, W. Duan and J. Li, Appl. Phys. Lett., 2013, 102, 182905.
98. M. K. Alam, M. N. M. Nobin and M. L. Ali, RSC Adv., 2024, 14, 1267–1283.
99. R. Gaillac, P. Pullumbi and F. X. Coudert, J. Phys.:Condens. Matter, 2016, 28, 275201.
100. S.-B. Chen, S.-D. Guo, B. Lv and Y. S. Ang, Phys. B, 2024, 689, 416201.
101. Y. Jain and R. Kurchania, Mater. Sci. Semicond. Process., 2025, 186, 109030.
102. M. Fine, L. Brown and H. Marcus, Scr. Metall., 1984, 18, 951–956.
103. D. R. Clarke, Surf. Coat. Technol., 2003, 163, 67–74.
104. J. A. Brehm, H. Takenaka, C.-W. Lee, I. Grinberg, J. W. Bennett, M. R. Schoenberg and A. M. Rappe, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 89, 195202.
105. Z. Pan, Y.-W. Fang, S. A. Nikolaev, L. Wu, J. Zhang, M. Ye, J. Liu, X. Ye, X. Wang and T. Nishikubo, et al., Mater. Horiz., 2025, 12, 6804–6811.
106. H. Pritchard and H. Skinner, Chem. Rev., 1955, 55, 745–786.
107. L. Xu, Z. Wang, B. Su, C. Wang, X. Yang, R. Su, X. Long and C. He, Crystals, 2020, 10, 434.
108. S. Palagummi and F. Yuan, Sens. Actuators, A, 2018, 279, 743–752.
109. T. Lu, B. Huang, Z. Wen, L. Huang, X. Wei and Z. Zou, Dig. J. Nanomater. Biostruct., 2022, 17, 715–730.
110. Q. Sun, J. Li, Y. Li, Z. Yang and R. Wu, Appl. Phys. Lett., 2021, 119, 062404.
111. X. Wang, Z. Cheng, J. Wang, L. Wang, Z. Yu, C. Fang, J. Yang and G. Liu, RSC Adv., 2016, 6, 57041–57047.
112. S. Yoon, E. H. Otal, A. E. Maegli, L. Karvonen, S. K. Matam, S. G. Ebbinghaus, B. Walfort, H. Hagemann, S. Pokrant and A. Weidenkaff, J. Alloys Compd., 2014, 613, 338–343.
113. S. Yoon, A. E. Maegli, L. Karvonen, S. K. Matam, A. Shkabko, S. Riegg, T. Großmann, S. G. Ebbinghaus, S. Pokrant and A. Weidenkaff, J. Solid State Chem., 2013, 206, 226–232.
114. S. Yoon, K. Son, S. G. Ebbinghaus, M. Widenmeyer and A. Weidenkaff, J. Alloys Compd., 2018, 749, 628–633.
115. H. W. Kang and S. B. Park, Chem. Eng. Sci., 2013, 100, 384–391.
116. L. E. Bell, Science, 2008, 321, 1457–1461.
117. P. Wu, Z. He, M. Yang, J. Xu, N. Li, Z. Wang, J. Li, T. Ma, X. Lu and H. Zhang, et al., Int. J. Thermophys., 2021, 42, 111.
118. S. R. Brown, S. M. Kauzlarich, F. Gascoin and G. J. Snyder, Chem. Mater., 2006, 18, 1873–1877.
119. J. R. Sootsman, D. Y. Chung and M. G. Kanatzidis, Angew. Chem., Int. Ed., 2009, 48, 8616–8639.
120. T. Sugahara, M. Ohtaki and T. Souma, J. Ceram. Soc. Jpn., 2008, 116, 1278–1282.
121. C. Soumya and P. Pradyumnan, Mater. Today Commun., 2023, 35, 106197.
122. A. Ayyaz, G. Murtaza, M. Naeem, A. Usman, S. M. Ramay, M. Irfan and H. Irfan, J. Phys. Chem. Solids, 2024, 188, 111936.
123. Y. Liu, Y. Tang, L. Shen, W. Ge, P. Yang and S. Deng, Ceram. Int., 2023, 49, 27772–27779.
124. R. Funahashi and S. Urata, Int. J. Appl. Ceram. Technol., 2007, 4, 297–307.
125. G. Kieslich, G. Cerretti, I. Veremchuk, R. P. Hermann, M. Panthöfer, J. Grin and W. Tremel, Phys. Status Solidi A, 2016, 213, 808–823.
126. D. Abdullah and D. C. Gupta, Sci. Rep., 2024, 14, 26168.
127. A. Hanif, S. A. Aldaghfag, A. Aziz, M. Yaseen and A. Murtaza, Int. J. Energy Res., 2022, 46, 10633–10643.
128. R. Chasmar and R. Stratton, Int. J. Electron., 1959, 7, 52–72.
129. R. Meenakshi, R. A. S. Srinivasan, A. Amudhavalli, R. R. Palanichamy and K. Iyakutti, Phys. B, 2024, 678, 415741.
130. M. Wolf, R. Hinterding and A. Feldhoff, Entropy, 2019, 21, 1058.
131. M. Yaseen, H. Ambreen, R. Mehmood, M. Iqbal, J. Iqbal, T. Alshahrani, S. Noreen and A. Laref, Phys. B, 2021, 615, 412857.
132. J. W. Vandersande and C. Wood, Contemp. Phys., 1986, 27, 117–144.
133. Y. Ouyang, Z. Zhang, D. Li, J. Chen and G. Zhang, Ann. Phys., 2019, 531, 1800437.
134. H. Yu, S. Gao, H. Ma and X. Jia, Ceram. Int., 2025, 51, 46135–46141.

---