Optical properties of BaTiO₃ at room temperature: DFT modelling

Abstract

The optical properties of the tetragonal phase of BaTiO₃ have been studied using the density functional method. In the study of a static lattice, we employed the generalized gradient approximation functional combined with on-site Hubbard correlation (GGA + U) and a hybrid functional. To account for the thermal motion of the atoms, we performed ab initio molecular dynamics calculations using the GGA + U method. We calculated the optical absorption spectra both for the static lattice and along the molecular dynamics trajectory. The results show that considering the motion of atoms leads to a significant decrease in the calculated value of the threshold energy for optical absorption. This effect occurs for two main reasons. First, changes in the atomic configuration due to thermal motion make electronic transitions that were previously dark for a static lattice become bright. Second, the optical absorption threshold decreases due to fluctuations in the energy of electronic transitions caused by the motion of atoms. The calculations were performed separately for different k-points of the Brillouin zone. The dispersion of electron energy in the reciprocal space may explain some features observed in the photoluminescence spectra.

---

1 Introduction

Barium titanate (BaTiO₃) has been thoroughly studied as a ferroelectric,^1–3 high-k dielectric,^4,5 and photorefractive material^6,7 among the perovskite oxides, both in experimental and theoretical research. BaTiO₃ shows promise for various technological applications, including gate dielectrics,⁸ holographic storage,⁹ photocatalysts,^10,11 and multi-layer ceramic capacitors,¹² due to its unique properties.

In order of increasing temperature, BaTiO₃ has five distinct structures: rhombohedral, orthorhombic, tetragonal, cubic, and hexagonal. The ferroelectric and optical properties of the tetragonal phase at 290–390 K render it significant.

The PL response can be significantly enhanced by external influences, such as temperature,¹³ applied electric field,^14–16 or mechanical stress,¹⁷ especially when the host matrix is a ferroelectric material that is highly susceptible to these excitations. Ferroelectrics are materials that are inherently characterized by a spontaneous polarization that can be switched in response to an electric field. They are highly responsive to any structural changes induced by electrical, thermal, or mechanical means, and polarization is also dependent on the temperature and strain.

A broad PL spectrum is observed in BaTiO₃, which is highly crystalline and exhibits weak PL properties at room temperature.^18,19 The PL spectra are enhanced by structural deficiencies and reduced band gaps.^18,20,21 BaTiO₃ that is non-crystalline and/or disordered exhibits PL at approximately 550 nm when excited at 488 nm.

It was found experimentally that highly crystalline BaTiO₃ exhibits a small broad peak at approximately 500 nm, acute peaks at 540 and 610 nm, and a large broad peak above 800 nm in its PL spectra.¹⁹ Emission peaks were observed in the visible-light region, particularly at 300 K, and they dissipated as the temperature increased. The quantum yields for the emission peaks are less than 1% and demonstrated a significant decrease as the temperature increased. At room temperature, as the wavelength of the exciting radiation decreases, the intensities of the radiation peaks at 540 and 610 nm increase.

In this paper, we simulate the optical absorption and photoluminescence spectra at room temperature, taking into account the contribution of each k-point of the Brillouin zone. The thermal motion of atoms significantly influences the optical absorption of BaTiO₃. The displacement of atoms from equilibrium positions results in the brightening of previously dark optical transitions in the static lattice, and the motion of atoms causes fluctuations in the transition energies. By figuring out the quantum yield of photoluminescence, it was found that BaTiO₃ doesn't have many luminescent properties. It is likely that structure defects are responsible for the photoluminescence observed in the experiment.

2 Methods

2.1 Theoretical approaches

Computational characterization of materials' properties is accomplished through five distinct methods: ground-state electronic structure, linear optical absorption, thermal motion of nuclei, nonradiative transitions, and photoluminescence.

2.1.1 Ground-state electronic structure. Atomistic models are represented by set of nuclear coordinates {R_I}. For selected atomistic model, for a set of momentum values k one solves a set of self-consistent equations is solved iteratively for total density, exchange–correlation potential, Kohn–Sham orbitals ψ^KS_j(k) (bands), and their energies ε^KS_j(k). These self-consistent calculations are performed for both optimized geometry and along time dependent trajectory {R_I(t)}. The evaluation of optical properties, nonradiative energy dissipation, and photoluminescence is performed in the basis of ψ^KS_j(k).

2.1.2 Linear optical absorption. The linear optical properties were analyzed based on the transition-dipole-moment matrix elements at different k-points in the reciprocal space:

D_ij(k) = e〈ψ^KS_i(k)|r|ψ^KS_j(k)〉 (1)

for transitions between the initial state (i) and final state (j) calculated in the basis of Kohn–Sham (KS) orbitals ψ^KS_i(k) where k is a wave vector in reciprocal space, e is an elementary charge. The transition dipole moment was used for calculating the oscillator strength:where m and ℏ are the electron mass and Planck constant, respectively, and ν_ij is the frequency of transition between the i^th and j^th states. Using the oscillator strengths and assuming a lack of spin–orbit coupling, the absorption spectra can then be determined as

2.1.3 Thermal motion of nuclei. Heating and adiabatic molecular dynamics (MD) simulations are used to model the movement of atoms. In each time step, the atomistic system interacts with a thermostat set at a specific temperature, which then rescales the ionic momenta (P_I). Newton's equation of motion is used to calculate the MD trajectory of the adiabatic ground state as a function of time once the system reaches thermal stability. The nuclear momenta are used as input. The coordinates of the atoms are collected at each time step of the MD trajectory to enable an electronic structure analysis. For adiabatic MD calculations, the geometry optimized models were heated to 300 K by a Nosé–Hoover thermostat with repeated velocity rescaling. Then, microcanonical trajectories were generated using the Verlet algorithm.²²

2.1.4 Nonradiative dynamics. The time evolution of the electronic degrees of freedom is calculated by solving the equation of motion,Here, the F_ij correspond to the Fock matrix elements F_ij = δ_ijε_j (ε_j is the energy of j^th KS orbital) of adiabatic contribution and represents the electronic dissipative transitions caused by nonadiabatic electron–ion interaction. Nonadiabatic couplings (NAC) can be computed according to the on-the-fly procedure along the nuclear trajectory as a measure of violation of KS orbital orthogonality due to nuclear motionwhere Δt, r, and R_I(t) are the MD timestep, electronic, and ionic coordinates, respectively. The autocorrelation function M_ijkl(τ) is processed by averaging along the duration of the trajectory T_MD as A Fourier transform of the autocorrelation function provides elements of Redfield tensor, which control the dissipative dynamics of the density matrix

The nonequilibrium distribution of charge as a function of time and energy reads where (a,b) is corresponding to the initial photoexcitation from state a to b. The change in population from the equilibrium distribution is

Δn^(a,b)(ε,t) = n^(a,b)(ε,t) − n^eq(ε,t) (6)

where Δn > 0 in the equation describes the population gain and Δn < 0 the population loss.

2.1.5 Photoluminescence. One calculates the time-resolved and time-integrated PL as follows.

An emission event between a pair of orbitals may occur if inverse population conditions are satisfied where ρ_jj > ρ_ii, ε_j > ε_i. The technique for calculating nonequilibrium electron dynamics and photoluminescence spectra is described in more detail in our previous publications.^23–27

2.2 Computational details

The electronic structures, which encompass Coulomb, correlation, and exchange electron–electron interactions, as well as interactions between electrons and ions, were determined using density functional theory (DFT) through the self-consistent solution of the Kohn–Sham equations, as implemented in the Vienna ab initio simulation package (VASP).^28,29 The hybrid functional Heyd–Scuseria–Ernzerhof (HSE06)³⁰ and the Perdew–Burke–Ernzerhof (PBE)³¹ exchange–correlation functional under the generalized gradient approximation were employed. The precise Hartree–Fock exchange was determined by the mixing parameter of 0.25 for HSE06 functional. The Dudarev parametrization³² was used to account for the on-site Coulomb correlation of Ti-3d electrons, with a U–J value of 2.6 eV in the case of PBE functional.³³ The projector augmented-wave formalism,³⁴ which was devised using the concept of pseudopotentials in the plane-wave basis, was employed to perform the calculations. All MD simulations were performed with PBE + U functional at T = 300 K. The temperature was controlled using the Nosé–Hoover thermostat.³⁵

The modeling conducted in this study employed the tetragonal BaTiO₃ phase (P4mm, #99), which is not energetically favorable at zero and exists at ambient temperature. The Ba atom has a coordination number of six with oxygen anions, while the Ti atom has a coordination number of twelve. The TiO₆ octahedra form a skeletal structure by sharing their corners, as illustrated in Fig. 1. Additionally, Ti cation is located at the center of this structure.

Fig. 1 Crystal structure of tetragonal BaTiO₃.

The initial crystal structure was obtained from the Materials Project database.³⁶ Used was a 3 × 3 × 3 supercell of the primitive cell, which contained 135 atoms. The calculations were conducted using a 500 eV basis and a 2 × 2 × 2 k-point mesh centered on the Γ-point. The analysis of the nature of chemical bonds was carried out using the electron localization function (ELF).^37,38 We used the VASPKIT code to postprocess the calculated data from VASP.³⁹

The original methodology^24,40,41 was developed for monitoring excited state dynamics in nanostructures with surfaces exposed to (disordered) liquid solvent. While applying this methodology to a periodic crystal one may find interesting observations related to dynamic change of the symmetry. One may introduce tolerance limit for coordinates retaining certain symmetry within allowed error bar. However such exploration would go beyond the scope of this work.

3 Results and discussion

3.1 Electronic properties of static lattice

The results presented here were obtained using PBE + U calculations of the optimized BaTiO₃ supercell geometry. Table 1 displays the band gap values (E_g) at all irreducible points in the Brillouin zone. We found these numbers in two ways: by calculating the difference between the energies of the conduction band minimum (CBM) and valence band maximum (VBM) (E^CBM–VBM = E(CBM) − E(VBM)), and by calculating the smallest value of the optical transition energy with a nonzero oscillator strength (E^optic). We performed the calculations using the PBE + U and HSE06 functionals.

Table 1 Calculated values of E^CBM–VBM and E^optic energies at different points of the Brillouin zone using PBE and HSE06 functionals

k-Point	E^CBM–VBM, eV		E^optic, eV
	PBE + U	HSE06	PBE + U	HSE06
Γ = (0 0 0)	2.40	3.56	2.58	3.70
Z = (0 0 ½)	2.46	3.62	2.67	3.79
X = (½ 0 0)	2.69	3.89	3.23	4.33
M = (½ ½ 0)	2.51	3.73	3.11	4.34
R = (0 ½ ½)	2.77	3.97	3.38	4.47
A = (½ ½ ½)	2.72	3.95	3.42	4.67

---

The valence band maximum of BaTiO₃ is located at the M point, whereas the conduction band minimum is found at the Γ point. Consequently, BaTiO₃ is classified as a semiconductor with an indirect bandgap, in line with results reported in the literature based on PBE calculations.⁴² The calculated PBE + U indirect bandgap for BaTiO₃ is 1.94 eV, while the direct bandgap is 2.48 eV.

In Section 3.4, we reexamine the character of Ti–O bonding and find that it is rather ionic and indirect optical transitions are not expected. The electronic structure of the ground state of BaTiO₃ was simulated while considering the effects of spin–orbit coupling. The calculations indicated a slight reduction in the direct bandgap value, which was found to be 2.44 eV. The value of the indirect bandgap remained unchanged. Due to its minimal influence, spin–orbit coupling was not considered further.

The E_g values calculated using the PBE + U functional are significantly smaller than the corresponding values obtained using the HSE06 hybrid functional. The results for E^CBM–VBM calculated at the Γ-point using the HSE06 functional are close to the analogous value of 3.66 eV obtained using the B3LYP functional¹⁸ and differ greatly from the results of PBE0 calculations, which give 4.1 eV.⁴³ At the same time, the E^CBM–VBM and E^optic values differ from each other. The minimum E^optic value significantly exceeds the experimental value of the optical absorption threshold, equal to 3.2 eV.^18,44 The reasons for such a large difference between the theoretical and experimental data are discussed below in Section 3.2.

The PBE + U calculated element-projected band structure of BaTiO₃ is presented in Fig. 2. Near the Fermi level, the valence band (VB) is formed by the O-2p states, while the conduction band (CB) is formed by the Ti-3d states. The top of the VB is at the A-point and is very close to the M-point. The bottom of CB is at the Γ-point and is close to the Z-point. In this paper, we consider only direct optical transitions, although the possibility of indirect optical transitions was predicted theoretically due to the presence of a covalent Ti–O bond in the compound.⁴⁵

Fig. 2 Element-projected band structure of BaTiO₃ supercell. Red and blue symbols correspond to O and Ti atoms, respectively.

The PBE + U calculations of the electron density of states (DOS) at all the irreducible points in Table 1's Brillouin zone are shown in Fig. 3(a). At all considered k-points in VB and CB, there are sub-bands in which electron states are absent. This leads to energy intervals in the optical absorption spectrum in which electrons with a given value of the wave vector do not absorb electromagnetic radiation (Fig. 3(b), dotted lines). Considering all k values, the total optical absorption indicates that BaTiO₃ absorbs in the entire energy range above E^optic (Fig. 4(a)). Similar data for the static lattice calculated using the HSE06 functional are presented in Fig. 4(b).

Fig. 3 (a) DOS (b) absorption static (dotted lines) and at T = 300 K (solid lines) calculated at different k-points.

Fig. 4 (a) PBE + U calculated total linear absorption spectra for static BaTiO₃ lattice and at T = 300 K. (b) HSE06 calculated linear absorption spectra for static BaTiO₃ lattice: total absorption and contribution of each irreducible k-point.

3.2 Optical properties at T = 300 K

The above theoretical values of the band gap, calculated using hybrid functionals, differ significantly from the corresponding experimental values. To explain this contradiction, it is necessary to remember that the theoretical data indicated were obtained for static lattices without considering the motion of atoms at a finite temperature.

The oscillator strengths and energies of electron transitions along the MD trajectories were calculated at T = 300 K for all k-points in the Brillouin zone to consider the effect of atomic motion. The comparison of the results of static calculations using the optimized structure of the model supercell and the results of calculations along the MD trajectory at 300 K is carried out using the PBE + U functional since calculations using the hybrid HSE06 functional along the MD trajectory are too expensive. The optical absorption values averaged along the MD trajectory are shown in Fig. 3(b) compared to similar data for the static lattice. At all k-points, a significant (up to 0.8 eV in X-point) red shift in the optical absorption energy threshold value is observed.

To explain this effect, Fig. 5 shows the calculated values of the oscillator strengths and energies of electron transitions at the Γ- and Z-points. For the remaining points in k-space, the data are qualitatively the same. The selected k-points correspond to the lowest-energy electron transitions for the static lattice. The results shown in Fig. 5 represent data for two types of electron transitions – bright and dark transitions when considering the static crystal lattice. Thus, a change in the atomic configuration at a finite temperature enhances optical absorption by involving a more significant number of states in the process. Due to the thermal motion of atoms, some transitions that are dark in a static lattice become bright. While transitions that are bright in a static lattice do not turn dark, their oscillator strengths may decrease. For comparison, Fig. 5(a and c) displays the oscillator strengths for a static lattice as horizontal red dashed lines.

Fig. 5 Oscillator strengths (left axes) and electronic transition energies (right axes) are calculated along the MD trajectory at T = 300 K for the Γ- and Z-points. Panels (a) and (c) correspond to electronic transitions that are bright for static lattice calculations, while panels (b) and (d) present data for transitions that are dark for static lattice. Horizontal black dashed lines correspond to transition energy calculated for static lattice. Horizontal red dashed lines in panels (a) and (c) depict oscillator strengths for static lattice.

The transition of optical characteristics from dark to bright occurs due to changes in lattice symmetry caused by the thermal motion of atoms, where the interatomic distances experience contractions and elongations at the frequency of activated normal modes. The analysis reveals that throughout all structures along the MD trajectory, the lattice symmetry is P1. This change in symmetry is a result of a decrease in point symmetry due to atomic movement. One should note that the NVE molecular dynamics starts with equipartitioned activation of all nuclear degrees of freedom (normal modes). There is an expectation for independent contraction/elongation of all internuclear distances, bonds, dihedrals as function of time. The amplitude of these contractions/elongations is small, but phase is random. Strictly speaking, it will distort any absolute symmetry; there will be no rotation, or inversion that will reproduce the same structure, except translation by the lattice vector. Distortion of symmetry during ab initio MD is referred to as dynamical disorder. For analysis of keeping or breaking crystalline structure one may explore RDF as function of time, catching transitions from crystalline to amorphous phase or keeping “approximate” symmetry. Alternatively, one may “cool down” the model after the dynamics (optimization) and then analyze new equilibrium geometry, in a spirit of simulated annealing, for retention or change of original symmetry.

The reduction of symmetry leads to change in selection rules for matrix elements, eqn (1), and affects the absorption spectrum. Superimposed oscillations of all normal modes make these trajectories quasiperiodic; in a large enough model, the trajectories would gain a random nature.

Some of the optical transitions that are dark for the static lattice become bright as the system evolves along the MD trajectory. Such changes concern transitions with energies both below and above the value of E^optic. For transitions with energies below E^optic, the oscillator strengths decrease as the E^CBM–VBM transition energy is approached. The oscillator strengths for E^CBM–VBM transitions are less than 10⁻⁴. The second reason for the decrease in the optical absorption energy is thermal fluctuations in the energies of electronic transitions.

It is necessary to figure out the oscillator strengths and energies of electron transitions along the MD trajectory using the hybrid HSE06 functional to get a good idea of the optical absorption threshold. This task requires lengthy calculations. However, based on PBE + U calculations, one can assume that since the lowest energy of the optical absorption threshold corresponds to absorption at the Γ-point of the Brillouin zone. The shift in the onset of optical absorption due to the thermal motion of atoms is approximately 0.5 eV. After applying this value to the HSE06 results for a static lattice, we get 3.2 eV as the expected value of the optical absorption threshold. This value is consistent with the results of the experiments.

3.3 Nonradiative charge relaxation and photoluminescence

Fig. 6 illustrates the supercell calculated excited-state dynamics of electron–hole pairs. The results were obtained by processing NACs from adiabatic MD trajectories using spin-restricted DFT for the in the different k-points. As mentioned in eqn (6), the blue and yellow colors represent the gain in population of electrons and holes, respectively. The green color indicates no change, which corresponds to the ground-state charge density. The energy expectation values for electrons and holes in the excited states are shown by solid and dashed lines, respectively.

Fig. 6 Supercell calculated charge carrier dynamics for states in (a) Γ- and (b) Z-points in Brillouin zone. The arrows show the intermediate states where excitation is temporarily delayed during the relaxation process.

Each panel shows how electronic populations are distributed based on energy and time. When each panel starts, there is an immediate excitation that moves electrons to the CBM_K + y (y-th orbital above the CBM) in the conduction band and introduces holes to the VBM_K − x (x-th orbital below the VBM) orbital in the valence band. The index K denotes a point in the Brillouin zone. These excitations, VBM_K − x → CBM_K + y, indicate optical transitions with very high oscillator strengths. After the excitation, a nonradiative process called vibrational relaxation converts the excited electrons to the CBM_K and the holes to the VBM_K. It's important to note that the electron–hole pairs stay in the frontier orbitals for a specific period of time before returning to the ground electronic state. Arrows in Fig. 6 depict intermediate states that trap excited electrons during relaxation. As a result, a radiative transition from these intermediate states is also possible, leading to the appearance of peaks in the PL spectrum. Upon including noncollinear spin basis and spin–orbit coupling effect into account one expects some modification of transition energies, relaxation rates, recombination rates, and intensity of optical transitions. Such comparison was reported for lead halide perovskites.⁴⁶

The most probable transitions can be identified by analyzing the time-integrated PL spectra of excited-state trajectories in Fig. 7. Nevertheless, this approach uses transition energies obtained from HSE06 calculations, which results in the generation of PL features with narrow line widths. In order to achieve a more accurate match with experimental PL spectra, it is necessary to consider the impact of thermal broadening. This has been reported through the method known as molecular dynamics photoluminescence (MDPL).²⁵

Fig. 7 Time-integrated PL spectra computed along the excited-state dynamics trajectories for all irreducible k-points in the reciprocal space.

The MDPL spectrum, which considers the contribution from the states at each k-point, is displayed in Fig. 8. The PL spectra for each k-point are within the energy range ≤ E^optic. The minimum energy is determined by the fluctuations of the KS orbitals along the MD trajectory. The PL spectrum exhibits a wide energy extension, with several PL peaks present at all k-points. This is due to the fact that electronic excitations relax to states at the edges of the bandgap, and excited states occupy orbitals with transition energies below E^optic. Fluctuating oscillator strengths of optical transitions are attributed to instantaneous atomic configurations at each point of the MD trajectory. As excited electrons and holes approach the CBM and VBM, the oscillator strengths decrease in absolute value. Despite this it remains a relatively high probability of a radiative transition occurring before the excitations have reached the band gap edge. This effect results in a wide PL range and multiple emission peaks.

Fig. 8 MDPL spectra computed at different k-points of Brillouin zone. Total spectrum is a sum of contribution from each k-point.

The calculated PL spectrum qualitatively matches the measured spectrum.¹⁹ This range spans from 1.8 to 2.8 eV (440–690 nm) and shows the existence of multiple peaks, though not as clearly as observed in the experiment. However, the theoretical PL quantum yield (PLQY) is ∼10⁻⁸, while the experimental PLQY value is ∼10⁻³. This difference can be explained by the fact that experimentally it was shown that the PL peak decreases as oxygen-deficient BaTiO₃ samples are annealed in an oxygen flow.¹⁸ This finding is supported by theoretical modeling.¹⁸ Since our simulations were carried out using a perfect lattice, the observed difference between calculated and experimental data suggests that even a small concentration of defects greatly improves the PL properties of BaTiO₃ compared to the ideal lattice. Therefore, the complete PL behavior in BaTiO₃ is influenced by both defects and fluctuations in electron transitions near the band gap edge.

This significant difference suggests that even a small concentration of defects greatly improves the PL properties of BaTiO₃ compared to the ideal lattice. Therefore, the complete PL behavior in BaTiO₃ is influenced by both defects and fluctuations in electron transitions near the band gap edge.

3.4 Chemical bonding in BaTiO₃

In Gao et al.,⁴⁵ the analysis of chemical bonds between atoms in the BaTiO₃ lattice was based on the examination of charge density in the interatomic space. However, to fully understand the nature of these chemical bonds, we must consider the formation of electron pairs, which we can achieve through Electron Localization Function (ELF) analysis.

The ELF value ranges from 0 to 1, with a value of 0.5 indicating a homogeneous electron gas. An ELF of 1.0 signifies perfect localization of electrons, while a value of 0.5 denotes perfect delocalization. Additionally, ELF readings below 0.5 suggest areas with minimal electron pair density.

The calculated ELF projections for the Ti–O and Ba–O bonds are shown in Fig. 9. A narrow region of minimum ELF separating the valence shells of Ti and O indicates the presence of ionic bonding. The minimal ELF values in the regions between Ti and O ions range from 0.12 to 0.22. In contrast, the minimal ELF value between Ba and O ions is 0.08. This analysis suggests that BaTiO₃ behaves primarily as a purely ionic crystal, with no covalent bonds present.

Fig. 9 Projection of ELF on planes passing through Ti–O (left panel) and Ba–O (right panel) bonds in the lattice of BaTiO₃. Color code blue (ELF = 0) and red (ELF = 1.0) represents full absence and full presence of electron pair in actual point of space.

4 Conclusions

We have employed the density functional method to investigate the optical properties of the tetragonal phase of BaTiO₃. We conducted the calculations for both a static lattice and one that takes into account the thermal motion of atoms at ambient temperature. The results show a significant reduction in the calculated threshold energy for optical absorption when accounting for atom motion. There are two primary causes for this phenomenon. Initially, electronic transitions that were previously dark for a static lattice become luminous as a result of changes in the atomic configuration caused by thermal motion. Secondly, fluctuations in the energy of electronic transitions induced by atoms' motion lead to a decrease in the optical absorption threshold. The dispersion of electron energy in the reciprocal space may account for certain characteristics observed in the photoluminescence spectra.

Data availability

Atomistic models of the considered materials are available upon request.

Author contributions

Conceptualization, D. S. K. and T. M. I.; methodology, D. R. G. and T. M. I.; software, T. M. I. and D. S. K.; validation, D. S. K. and T. M. I.; formal analysis, D. R. G., A. U. A. and F. U. A.; investigation, A. U. A., D. R. G., F. U. A., T. M. I. and D. S. K.; resources, T. M. I. and D. S. K.; writing—original draft preparation, T. M. I.; writing—review and editing, T. M. I., D. R. G. and D. S. K.; visualization, A. U. A.; supervision, T. M. I.; project administration, F. U. A.; funding acquisition, T. M. I. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was carried out with the financial support of the Ministry of Science and Higher Education of the Republic of Kazakhstan: AP14869492 “Development of nanocrystalline metal oxide catalysts for hydrogen production”. The work by T. M. I. was performed under the state assignment of GEOKHI RAS. D. R. G. acknowledges support of NSF 1944921 for developing methods of quantum dynamics. D. S. K. acknowledges partial support from the U.S. Department of Energy DE-SC0022239 for computational exploration of photoluminescence of earth abundant metal compounds.

Notes and references

1. K. J. Choi, M. Biegalski, Y. L. Li, A. Sharan, J. Schubert, R. Uecker, P. Reiche, Y. B. Chen, X. Q. Pan, V. Gopalan, L.-Q. Chen, D. G. Schlom and C. B. Eom, Science, 2004, 306, 1005–1009.
2. V. Buscaglia, M. T. Buscaglia and G. Canu, in Encyclopedia of Materials: Technical Ceramics and Glasses, ed. M. Pomeroy, Elsevier, Oxford, 2021, pp. 311–344,  DOI:10.1016/B978-0-12-803581-8.12132-0.
3. M. E. Lines and A. M. Glass, Principles and Applications of Ferroelectrics and Related Materials, Oxford University Press, 2001.
4. M. T. Buscaglia, M. Viviani, V. Buscaglia, L. Mitoseriu, A. Testino, P. Nanni, Z. Zhao, M. Nygren, C. Harnagea, D. Piazza and C. Galassi, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 73, 064114.
5. S. S. Parizi, D. Caruntu, A. Rotaru and G. Caruntu, Mater. Adv., 2022, 3, 6474–6484.
6. D. Kip, Appl. Phys. B:Lasers Opt., 1998, 67, 131–150.
7. Z. Hasan, M. A. Rahman, D. K. Das and H. K. Rouf, Sci. Rep., 2023, 13, 10487.
8. H. Lee, J. F. McGlone, S. I. Rahman, C. Chae, C. Joishi, J. Hwang and S. Rajan, Phys. Status Solidi RRL, 2024, 18, 2400042.
9. W. R. Lopes, H. F. A. Medeiros, G. S. Santos, T. C. Araujo, J. F. Carvalho, P. V. dos Santos and M. T. de Araujo, J. Opt. Soc. Am. A, 2018, 35, 1919–1928.
10. S. Chandrappa, S. N. Myakala, N. A. Koshi, S. J. Galbao, S.-C. Lee, S. Bhattacharjee, D. Eder, A. Cherevan and D. H. K. Murthy, ACS Appl. Mater. Interfaces, 2024, 16, 8763–8771.
11. T. M. Inerbaev, A. U. Abuova, Z. Y. Zakiyeva, F. U. Abuova, Y. A. Mastrikov, M. Sokolov, D. Gryaznov and E. A. Kotomin, Molecules, 2024, 29, 2707.
12. Y. Y. Huang, Y. L. Yang, L. Y. Zhang, V. Laletin, V. Shur, R. Y. Jing and L. Jin, J. Alloys Compd., 2024, 1006.
13. M. K. Mahata, T. Koppe, K. Kumar, H. Hofsäss and U. Vetter, Sci. Rep., 2020, 10, 8775.
14. J. Hao, Y. Zhang and X. Wei, Angew. Chem., Int. Ed., 2011, 50, 6876–6880.
15. X. Tian, Z. Wu, Y. Jia, J. Chen, R. K. Zheng, Y. Zhang and H. Luo, Appl. Phys. Lett., 2013, 102, 042907.
16. H. Sun, X. Wu, D. F. Peng and K. W. Kwok, ACS Appl. Mater. Interfaces, 2017, 9, 34042–34049.
17. Y. Zhou, R. Xiong, P. Wang, X. Wu, B. Sa, C. Lin, M. Gao, T. Lin and C. Zhao, J Materiomics, 2022, 8, 586–595.
18. E. Orhan, J. A. Varela, A. Zenatti, M. F. C. Gurgel, F. M. Pontes, E. R. Leite, E. Longo, P. S. Pizani, A. Beltràn and J. Andrès, Phys. Rev. B:Condens. Matter Mater. Phys., 2005, 71, 085113.
19. T. Watanabe, D. Hoshi, M. Ishida and T. Ohba, Langmuir, 2018, 34, 10250–10253.
20. M.-S. Zhang, Z. Yin, Q. Chen, W. Zhang and W. Chen, Solid State Commun., 2001, 119, 659–663.
21. M. L. Moreira, M. F. C. Gurgel, G. P. Mambrini, E. R. Leite, P. S. Pizani, J. A. Varela and E. Longo, J. Phys. Chem. A, 2008, 112, 8938–8942.
22. L. Verlet, Phys. Rev., 1967, 159, 98–103.
23. T. M. Inerbaev, Y. Han, T. B. Bekker and D. S. Kilin, J. Phys. Chem. C, 2022, 126, 6119–6128.
24. T. M. Inerbaev, J. D. Hoefelmeyer and D. S. Kilin, J. Phys. Chem. C, 2013, 117, 9673–9692.
25. D. J. Vogel and D. S. Kilin, J. Phys. Chem. C, 2015, 119, 27954–27964.
26. Fatima, Y. Han, D. J. Vogel, T. M. Inerbaev, N. Oncel, E. K. Hobbie and D. S. Kilin, J. Phys. Chem. C, 2019, 123, 7457–7466.
27. D. R. Graupner and D. S. Kilin, J. Phys. Chem. C, 2024, 128, 3935–3944.
28. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
29. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
30. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
32. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B, 1998, 57, 1505–1509.
33. F. Maldonado, S. Jácome and A. Stashans, Comput. Condens. Matter, 2017, 13, 49–54.
34. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
35. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
36. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
37. A. Savin, R. Nesper, S. Wengert and T. F. Fässler, Angew Chem. Int. Ed. Engl., 1997, 36, 1808–1832.
38. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397–5403.
39. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033.
40. D. S. Kilin and D. A. Micha, J. Phys. Chem. Lett., 2010, 1, 1073–1077.
41. S. Huang and D. S. Kilin, J. Chem. Theory Comput., 2014, 10, 3996–4005.
42. V. B. Parmar, D. Raval, S. K. Gupta, P. N. Gajjar and A. M. Vora, Mater. Today: Proc., 2023 DOI:10.1016/j.matpr.2023.01.410.
43. R. A. Evarestov and A. V. Bandura, J. Comput. Chem., 2012, 33, 1123–1130.
44. S. H. Wemple, Phys. Rev. B, 1970, 2, 2679–2689.
45. H. Gao, J. Cao, L. Liu and Y. Yang, J. Mol. Struct., 2011, 1003, 75–81.
46. A. Forde, T. Inerbaev, E. K. Hobbie and D. S. Kilin, J. Am. Chem. Soc., 2019, 141, 4388–4397.

---