PHYSICAL REVIEW B 108, 235171 (2023) Understanding the role of Hubbard corrections in the rhombohedral phase of BaTiO3 G. Gebreyesus ,1,* Lorenzo Bastonero ,2 Michele Kotiuga ,3 Nicola Marzari,2,3 and Iurii Timrov 3,† 1Department of Physics, School of Physical and Mathematical Sciences, College of Basic and Applied Sciences, University of Ghana, P. O. Box LG63, Legon - Accra, Ghana 2Bremen Center for Computational Materials Science, and MAPEX Center for Materials and Processes, University of Bremen, D-28359 Bremen, Germany 3Theory and Simulation of Materials (THEOS), and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland (Received 8 September 2023; revised 22 November 2023; accepted 22 November 2023; published 27 December 2023) We present a first-principles study of the low-temperature rhombohedral phase of BaTiO3 using Hubbard- corrected density-functional theory. By employing density-functional perturbation theory, we compute the onsite Hubbard U for Ti(3d) states and the intersite Hubbard V between Ti(3d) and O(2p) states. We show that applying the onsite Hubbard U correction alone to Ti(3d) states proves detrimental, as it suppresses the Ti(3d)–O(2p) hybridization and drives the system towards a cubic phase. Conversely, when both onsite U and intersite V are considered, the localized character of the Ti(3d) states is maintained, while also preserving the Ti(3d)–O(2p) hybridization, restoring the rhombohedral phase of BaTiO3. The generalized PBEsol+U+V functional yields good agreement with experimental results for the band gap and dielectric constant, while the optimized geometry is slightly less accurate compared to PBEsol. Zone-center phonon frequencies and Raman spectra are found to be significantly influenced by the underlying geometry. PBEsol and PBEsol+U+V provide satisfactory agreement with the experimental Raman spectrum when the PBEsol geometry is used, while PBEsol+U Raman spectrum diverges strongly from experimental data highlighting the adverse impact of the U correction alone in BaTiO3. Our findings underscore the promise of the extended Hubbard PBEsol+U+V functional with first-principles U and V for the investigation of other ferroelectric perovskites with mixed ionic-covalent interactions. DOI: 10.1103/PhysRevB.108.235171 I. INTRODUCTION BaTiO3 (BTO) holds a prominent position among the ex- tensively studied ABO3 perovskite materials due to its wide range of technological applications in electronics, electrome- chanical energy conversion, nonlinear optics, and nonvolatile data storage [1–7]. It possesses a paraelectric cubic perovskite structure with no net polarization at high temperatures and un- dergoes a series of three ferroelectric phase transitions as the temperature decreases [8,9]. Between 394 K and 278 K, BTO adopts a tetragonal structure with the polarization along the 〈100〉 crystal direction; as the temperature further decreases between 278 K and 183 K, it transforms into an orthorhombic structure with polarization along 〈110〉; finally, below 183 K, BTO adopts a rhombohedral structure with the R3m space group with the polarization along 〈111〉 [8–11]. These phase transitions have been the subject of extensive experimental and theoretical investigations inspiring both the displacive [9,12] and order-disorder [13,14] models of ferroelectric tran- sitions. Computational studies based on density-functional the- ory (DFT) [15,16] have been instrumental in exploring the * ghagoss@ug.edu.gh †Present address: Laboratory for Materials Simulations (LMS), Paul Scherrer Institut (PSI), CH-5232 Villigen PSI, Switzerland; iurii.timrov@psi.ch structural, electronic, optical, and vibrational properties of various phases of BTO. These studies have employed di- verse exchange-correlation functionals, such as local-density approximation (LDA), generalized-gradient approximation (GGA), hybrid, and meta-GGA functionals [17–33]. Despite many results in good agreement with experiments, none of the aforementioned functionals provide an overall quantitative description of BTO. Consequently, the pursuit of even more accurate functionals has persisted. For instance, the PBEsol functional [34] has shown high accuracy in predicting the structure of the rhombohedral phase of BTO; however, it un- derestimates the band gap [30,31]. On the other hand, hybrid functionals like PBE0 [35] and HSE06 [36,37] improve the band-gap prediction, but they overestimate the lattice constant and the atomic distortions associated with ferroelectricity [20,27]. As a result, any physical property dependent on atomic distortions (e.g., phonons) is significantly influenced by the choice of the functional [38]. In turn, this affects the accuracy of vibrational properties, leading to shifts in computed Raman spectral peaks or incorrect intensities as compared to experiments. Hence, there is a pressing need to search for novel functionals capable of providing an accurate characterization of the structural, electronic, and vibrational properties of BTO simultaneously. One of the widely used approaches to model transition- metal oxides is DFT + U [39–41], where the Hubbard U correction is applied to selected states (typically to partially filled d states) to alleviate self-interaction errors [42–44]. The 2469-9950/2023/108(23)/235171(15) 235171-1 ©2023 American Physical Society https://orcid.org/0000-0002-3413-2349 https://orcid.org/0000-0001-9374-1876 https://orcid.org/0000-0002-2592-7342 https://orcid.org/0000-0002-6531-9966 https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevB.108.235171&domain=pdf&date_stamp=2023-12-27 https://doi.org/10.1103/PhysRevB.108.235171 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) application of DFT + U to BTO has only gained prominence in the last decade [27,45–52]. The reluctance to apply this approach with the Hubbard U correction on Ti(3d) states in BTO in early DFT-based works was primarily due to the fact that Ti ions are in the 4+ oxidation state (OS), resulting in a d0 electronic configuration. Thus, in BTO the 3d states of Ti are nominally empty, suggesting that the application of Hubbard U would have a negligible effect. Conversely, in BTO with oxygen vacancies or in BTO oxyhydrides, some Ti ions undergo reduction to 3+ OS with a d1 configura- tion, which provides a stronger motivation for the use of Hubbard corrections [45,47,50]. However, earlier DFT studies of the pristine BTO employing LDA and GGA revealed a significant hybridization between Ti(3d) and O(2p) states, indicating that the four electrons of the Ti(3d) states are not entirely transferred to the neighboring O(2p) states (see, e.g., Ref. [23]). In fact, this hybridization plays a crucial role in the covalency of the Ti–O bonds, which is essential for the manifestation of ferroelectricity [17,23,53–55]. Therefore, the effect of the U correction on Ti(3d) states should not be presumed to be exactly zero. In fact, applying the Hubbard U correction to the Ti(3d) states localizes them more on Ti ions and erroneously suppresses hybridizations with the O(2p) states (see Fig. 6 in Ref. [51]). Consequently, the rhombohedral distortion disappears, and instead, the cubic phase of BTO becomes more stable (see Fig. S6 in Ref. [27]). This unequivocally demonstrates that the effect of onsite U corrections on Ti(3d) states is detrimental in pristine BTO and underscores the crucial importance of preserving hy- bridizations with O(2p) states. In this regard, the extended DFT + U + V formulation [56] proves highly attractive as it allows for the application of an onsite U correction to the Ti(3d) states to reduce self-interaction errors while en- suring the intersite Ti(3d)–O(2p) hybridizations thanks to the V corrections, thereby preserving the fundamental co- valency of the Ti–O bonds essential for ferroelectricity in BTO [53]. The main challenge of the DFT + U approach is the lack of a priori knowledge of the value for U . Even though there are various first-principles methods to compute it [57], U is still often determined empirically. Despite numerous DFT + U studies of BTO [27,45–52], no consensus has been reached regarding which states should undergo the Hubbard correction and which value of U should be used. Additionally, differ- ences between various types of Hubbard projector functions are often overlooked when comparing the U values from different DFT + U studies of various materials [58,59]. For example, in previous DFT + U investigations of BTO, Ti(3d) states were corrected using the Hubbard U value of 3 eV [50] or 4.49 eV [27,45,47], while in Refs. [46,49] the Hubbard U correction was applied to O(2p) states with the value of 8 eV. Furthermore, in Ref. [52] the U correction of 8 eV was applied to both Ti(3d) and O(2p) states. Additionally, the choice of Hubbard projector functions also varied in these works: Most studies employed projector-augmented wave (PAW) Hubbard projector functions, except for Ref. [45] where a different type of projector was utilized (not reported in that paper but either nonorthogonalized or orthogonalized atomic orbitals [60]). Moreover, the value of U was determined empirically in all the aforementioned DFT + U studies, with only one exception [45], where it was computed using linear-response theory [61]. As a result of this large variation in U values and the ambiguity surrounding the correction of Ti(3d) or O(2p) states or both, there exists a wide spread of results that of- ten contradict both each other and experimental observations. Furthermore, as of now, no DFT + U + V studies of BTO have been conducted to investigate the significance of intersite Hubbard V corrections between Ti(3d) and O(2p) states. In this study, we present a comprehensive first-principles investigation of the low-temperature rhombohedral phase of BTO, focusing on its structural, electronic, and vibrational properties using three functionals: PBEsol, PBEsol+U , and PBEsol+U+V . To determine the onsite U and intersite V Hubbard parameters, we employ a rigorous first-principles approach based on linear-response theory [61], recast in terms of density-functional perturbation theory (DFPT) [62,63] in a basis of Löwdin-orthogonalized atomic orbitals (Hub- bard projector functions). This approach eliminates any empirical input and potential ambiguities typically present in Hubbard-corrected DFT studies. Furthermore, the self- consistent procedure is employed for computing the Hubbard parameters [63] to ensure the mutual consistency of the crystal and electronic structures. In agreement with Ref. [27], we find that applying the onsite U correction solely to Ti(3d) states results in a pronounced suppression of hybridization with neighboring O(2p) states, driving the system into a cubic phase. However, our DFT + U calculations, utilizing first-principles U value, provide the cubic structure that still exhibits dynamical instability with imaginary phonon modes around the � point, in contrast to the results of Ref. [27]. Conversely, the introduction of intersite Hubbard V inter- actions between the Ti(3d) and O(2p) states dramatically alters the overall picture by restabilizing the rhombohedral phase thanks to the restored covalency of the Ti–O bonds. The optimized geometry obtained from PBEsol+U+V is found to be somewhat less accurate than that from PBEsol, whereas the projected density of states (PDOS) is qualitatively very similar in both cases. In addition, the Born effective charges (BEC) are slightly smaller for certain components within PBEsol+U+V as compared to PBEsol, while the di- electric constant and band gap from PBEsol+U+V exhibit good agreement with experimental data, surpassing the ac- curacy of PBEsol predictions by a significant margin. On the contrary, all predictions obtained from PBEsol+U are systematically less accurate compared to the PBEsol+U+V case. Last, the zone-center phonon frequencies and Raman spectra are found to be highly sensitive to the underlying geometry. The PBEsol and PBEsol+U+V Raman spectra are found to be in satisfactory agreement with experiments provided that the PBEsol geometry is used, while the Ra- man spectrum from PBEsol+U differs dramatically from the experimental one. The remainder of this paper is organized as follows: Sec- tion II presents the computational details; Sec. III contains an in-depth analysis of the results, encompassing the structural properties, PDOS, BEC, the dielectric tensor, phonon disper- sions, and Raman spectra; and, finally, Sec. IV contains the concluding remarks. 235171-2 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) FIG. 1. Crystal structure of BTO in the (a) cubic and (b) rhombo- hedral perovskite structures with the Ba, Ti, and O atoms represented in green, gray, and red, respectively. While the cubic structure is nonpolar, the rhombohedral phase exhibits a polarization along the 〈111〉 direction, denoted by the dark blue arrow. (c) Schematic illus- tration of the Ti and O atomic displacements in the rhombohedral structure relative to the cubic one projected onto the ab plane. The total displacements are depicted with yellow arrows, while the components relative to the cubic structure are indicated with black arrows: Dashed lines for �Ti and solid and dotted lines for the two O displacements �O and �O′ . For clarity, Ba atoms are omitted. The rotational symmetry around the 〈111〉 direction results in identical displacements when viewed in the ac and bc planes. II. COMPUTATIONAL DETAILS All calculations are performed using the plane-wave pseudopotential method as implemented in the QUANTUM ESPRESSO distribution [64–66]. We use the exchange- correlation functional constructed using GGA with the PBEsol parametrization [34], and we employ the following pseudopotentials from the Pslibrary v1.0.0 [67]: Ba.pbesol-spn-rrkjus_psl.1.0.0.UPF, Ti.pbesol-spn- rrkjus_psl.1.0.0.UPF, and O.pbesol-n-rrkjus_psl.1.0.0.UPF. The Kohn-Sham wave functions and potentials are expanded in plane waves up to a kinetic-energy cutoff of 80 and 800 Ry, respectively. The first Brillouin zone is sampled using uniform k-point meshes of size 8 × 8 × 8 and 18 × 18 × 18 centered at the � point for the ground state and PDOS calculations, respectively. PDOS is obtained using the tetrahedron method [68]. The low-temperature rhombohedral phase of BTO is de- scribed by its lattice constant, rhombohedral angle, and symmetry-preserving internal atomic displacements along the 〈111〉 direction. The atomic positions of the rhombohedral crystal structure can be expressed based on those in the cubic structure (in crystal coordinates) as follows: Ba : (0.0 + �Ba, 0.0 + �Ba, 0.0 + �Ba ), Ti : (0.5 + �Ti, 0.5 + �Ti, 0.5 + �Ti), O1 : (0.5 + �O, 0.5 + �O, 0.0 + �O′ ), O2 : (0.5 + �O, 0.0 + �O′ , 0.5 + �O), O3 : (0.0 + �O′ , 0.5 + �O, 0.5 + �O), where �Ba, �Ti, �O, and �O′ represent the displacements of the Ba, Ti, and O ions from their atomic positions in the cubic phase (see Fig. 1). We perform a full structural TABLE I. Comparison of the optimized lattice parameter a, the rhombohedral angle α, and the atomic displacements (in crystal coordinates) �Ti, �O, and �O′ (see Fig. 1) computed in this work using PBEsol, PBEsol+U , and PBEsol+U+V and as obtained in previous computational studies [20,31,100,101] and in experiments at 15 K [102]. �Ba = 0 in all cases. The Hubbard U correction is applied only to the Ti(3d) states while V is applied between Ti(3d) and O(2p). a (Å) α (deg) �Ti �O �O′ PBEsol 3.998 89.87 −0.012 0.011 0.018 PBEsol+U 3.990 90.00 0.000 0.000 0.000 PBEsol+U+V 4.017 89.77 −0.014 0.012 0.020 LDA [100] 3.931 89.92 −0.007 0.010 0.014 PBE [20] 4.073 89.71 −0.015 0.014 0.025 PBEsol [31] 3.998 −0.012 0.012 0.019 PBE0 [20] 4.029 89.73 −0.015 0.013 0.024 B1-WC [101] 3.991 89.86 Expt. [102] 4.004 89.84 −0.013 0.011 0.019 relaxation to calculate these displacements using three functionals (PBEsol, PBEsol+U , and PBEsol+U+V ). The Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [69] is employed for geometry optimization with convergence thresholds set to 10−8 Ry, 10−5 Ry/bohr, and 0.01 Kbar for the total energy, forces, and pressure, respectively. For the PBEsol+U and PBEsol+U+V calculations, we determine the onsite U and intersite V Hubbard parameters self-consistently using DFPT [62,63], implemented in the HP code [70] of QUANTUM ESPRESSO [64–66]. The reader is invited to check Refs. [62,63] for the detailed explanations on how the Hubbard parameters are defined and computed using DFPT. As Hubbard projector functions, we employ the Löwdin-orthogonalized atomic orbitals [60,71,72]. We use uniform �-centered k- and q-point meshes of size 6 × 6 × 6 and 4 × 4 × 4, respectively, in computing the Hubbard param- eters with an accuracy of ∼0.01 eV for U and ∼0.001 eV for V . Employing a self-consistent procedure [63], our cal- culations encompass cyclic iterations that involve structural optimizations and successive recalculations of Hubbard pa- rameters for each new geometry. By doing so, we obtain structural parameters that are optimized using PBEsol+U and PBEsol+U+V using the respective self-consistent Hubbard parameters and that are listed in Table I. This methodology has demonstrated its efficacy in yielding accurate results for various transition-metal compounds [73–80]. The computed values of the Hubbard parameters are as follows: U = 4.52 eV for Ti(3d) states within PBEsol+U and U = 5.21 eV for Ti(3d) states with intersite V values ranging from 1.17 to 1.37 eV (depending on the interatomic distance) between Ti(3d) and O(2p) states within PBEsol+U+V . These values are consistent with previous studies of SrTiO3 using the same approach [81]. It is important to note that Hubbard parameters are not universal [62,63,70]: Their values depend on the chemi- cal environment of transition-metal elements, oxidation state [43], spin configuration [74,76], and other factors. When it comes to the material phase’s influence on computed Hub- bard parameters, variations can arise, especially if distinct 235171-3 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) phases trigger notable changes in the chemical environment of transition-metal elements. For instance, in the case of the tetragonal phase of BTO, we obtain the following values: U = 4.50 eV for Ti(3d) states within PBEsol+U , and U = 5.25 eV for Ti(3d) states with intersite V values ranging from 1.29 to 1.32 eV between Ti(3d) and O(2p) states within PBEsol+U+V . As can be seen, the differences between the Hubbard parameters of tetragonal and rhombohedral phases are very small (< 1%), which is not surprising since the crystallographic differences between these two phases are small. Consequently, our computed Hubbard parameters for the rhombohedral phase are transferrable to other phases of BTO, provided identical computational setups, like pseudopo- tentials and Hubbard projector functions, are used. They can, at the very least, serve as very good starting point when refin- ing these parameters for other phases using DFPT. However, in general it is important to recompute Hubbard parameters for different phases of the same material and verify how significant are the changes. To compute the phonon frequencies, we utilize the frozen- phonon method implemented in the PHONOPY package [82], using atomic displacements of 0.01 Å. The 3 × 3 × 3 super- cell is employed to calculate the interatomic force constants. We verified that adopting larger supercells of size 4 × 4 × 4 does not significantly alter the phonon dispersions. The first Brillouin zone is sampled using a uniform 4 × 4 × 4 k-point mesh centered at the � point. To correct for the nonanalyt- icity of the dynamical matrix as q → 0 we include up to dipole-dipole interactions which depend on the BEC and the dielectric tensor [83]. We note that Royo et al. [32] investi- gated the effect of higher-order corrections, which we discuss but ultimately neglect in the results presented here. For the cal- culation of the BEC and the dielectric tensor, we employ two approaches: the finite electric field method [84–86] and DFPT [83,87–89]. The method of Refs. [84–86] for computing BEC, dielectric, Raman, and nonlinear optical susceptibility ten- sors is used as implemented in the aiida-vibroscopy [90] package, exploiting AiiDA [91,92]. The directional-sampling technique described in Ref. [90] is used to effectively sample the Brillouin zone with two sets of k-point distances, namely the “parallel” and the “orthogonal” (referring to the direction of the applied electric field). For the parallel and orthogonal distances, we set values of 0.1 Å and 0.3 Å, respectively (we note that 0.1 Å corresponds to the 16 × 16 × 16 uniform k-point mesh). The BEC and the dielectric tensor are defined in the real-space rhombohedral BTO reference system with the z axis parallel to the C3 rhombohedral axis (i.e., the 〈111〉 direction), while the x axis is perpendicular to it. The zone- center phonons and the phonon dispersions are computed for phonon wave vectors q defined in reciprocal space of the rhombohedral BTO reference system. To calculate the nonresonant Raman intensities, we adopt the Placzek approximation (first-order processes) [93]. In order to compare with experiments from Refs. [94,95], for Raman spectrum calculations we use the real-space tetragonal BTO reference system with the z axis being parallel to the c tetragonal axis, and the x axis is perpendicular to it. Within this framework, we consider the backscattering geometry with a parallel polarization configuration, namely the z(xx)z̄ setup in Porto notation [96]. This entails that the incident light’s propagation direction aligns with z, while the scattered light follows z̄. Both the incident and scattered light’s polarization direction align with x. The Raman scattering amplitude for a phonon mode ν is expressed as [97] Iν ∝ (ωL − ων ) nν + 1 ων |ei · ↔ Aν · es|2, (1) where nν represents the Bose-Einstein occupation, ↔ Aν denotes the Raman susceptibility tensor, ωL is the incident laser fre- quency, ων denotes the frequency of the phonon mode ν, and ei and es are the polarization vectors of the incident and scattered light, respectively. The temperature relevant for the Bose-Einstein occupation and the laser frequency are adjusted to match the experimental conditions. For the z(xx)z̄ setup, ei and es are parallel to each other, both lying along the x direction. While the propagation directions of the incident (ki) and scattered (ks) light do not appear explicitly in Eq. (1), their difference q = ks − ki is utilized to specify the direction for the limit q → 0 when computing the nonanalytic component of the dynamical matrix and the Raman susceptibility tensor [98]. To model the experimental backscattering geometry, q is chosen to be parallel to the z axis in reciprocal space of the tetragonal BTO reference system, as the photon propagates back and forth along this direction. To generate the Raman spectra, we employ a Lorentzian smearing function with a constant broadening parameter of 8 cm−1. The data used to generate the results presented in this paper are accessible in the Materials Cloud Archive [99]. III. RESULTS AND DISCUSSION A. Structural properties In this section, we present the results of structural op- timizations for the low-temperature rhombohedral phase of BTO using PBEsol, PBEsol+U , and PBEsol+U+V func- tionals. A summary of our findings and a comparison with previous computational studies based on DFT with (semi-)local and hybrid functionals [20,31,100,101], as well as experimental data [102], are presented in Table I. It is well known that LDA underestimates the lattice parameters [100], while PBE tends to overestimate them [20], which is also ob- served in the case of BTO. Our PBEsol results are consistent with previous PBEsol studies [31], showing very good agree- ment with the experimental lattice parameter a, rhombohedral angle α, and atomic displacements. On the other hand, previ- ous DFT studies using PBE0 [31] and B1-WC [101] hybrid functionals either significantly overestimate or underestimate a, respectively, while α is closer to the experimental value. As mentioned in Sec. I, the PBEsol+U and PBEsol+U+V predictions for structural properties are in stark contrast to each other: The former drives the structure towards the cu- bic phase, while the latter preserves the rhombohedral phase (see Table I). The PBEsol+U and PBEsol+U+V structural parameters are discussed in more detail in the following. Let us compare the accuracy of the structural parameters obtained using PBEsol and PBEsol+U+V . As can be seen in Table I, PBEsol underestimates a by only ∼0.2%, and α is overesti- mated by only ∼0.03%. These results are significantly better than those obtained with LDA and PBE [20,100]. On the other hand, PBEsol+U+V overestimates a by ∼0.3% and underes- 235171-4 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) FIG. 2. Optimized lattice parameter a (a) and the rhombohedral angle α (b) as a function of the Hubbard U parameter applied to Ti(3d) states using PBEsol+U . The ab initio value of U is 4.52 eV and it was computed using DFPT, as discussed in Sec. II. Experi- mental values of the structural parameters are also indicated: For the rhombohedral phase these are arh exp = 4.004 Å and αrh exp = 89.84◦ at 15 K [102], and for the cubic phase these are ac exp = 4.001 Å and αc exp = 90.00◦ at 400 K [103]. timates α by ∼0.08%. Such seemingly minor deterioration in the accuracy of structural predictions within PBEsol+U+V compared to PBEsol has a significant impact on the lattice vibrational properties, as evidenced in Secs. III D and III E. Motivated by Refs. [46,49], we examine the impact of applying the Hubbard U correction to O(2p) states. First, we determine U for the O(2p) states from first principles using DFPT in a self-consistent fashion, and we find the values of 8.6 and 8.7 eV within PBEsol+U and PBEsol+U+V , respectively. Applying the U correction to both O(2p) and Ti(3d) states within PBEsol+U and PBEsol+U+V results in the optimized geometry of rhombohedral BTO adopting a cubic structure. This occurs as U localizes O(2p) electrons, diminishes Ti-O hybridization, suppressing covalency, and ul- timately stabilizing the cubic phase. However, this contradicts experimental observations where the rhombohedral phase pre- vails. Consequently, we abstain from computing additional properties such as electronic structure and vibrational spectra that involve the U correction on O(2p). In the rest of the paper we neglect the effect of U on O(2p) states. The stabilization of the cubic phase of BTO within PBEsol+U requires further analysis. To investigate this, we performed structural optimizations starting from the experi- mental rhombohedral structure while varying the value of U in the range from 0 to 10 eV. Figure 2 illustrates the variation of lattice parameter a and rhombohedral angle α as a function of U . Increasing U from 0 to 4.5 eV leads to decreasing of a in a quasilinear fashion by ∼0.2% and a quasilinear increase in α by ∼0.1%, leading it towards 90◦. Further increasing U beyond 4.5 eV results in a larger a (also quasilinear change), while α remains stable at 90◦. Hence, a critical value of U is observed at approximately 4.5 eV, where the cubic phase becomes stable. Interestingly, our first-principles value of U is found to be 4.52 eV, which aligns well with this criti- cal value. It is important to note that these U values were determined using Löwdin-orthogonalized atomic orbitals as Hubbard projectors (see Sec. II), and hence they may not be directly applicable with other types of Hubbard projectors, necessitating reevaluation for such cases. Therefore, we find that by increasing U from 0 to 4.5 eV the unit cell volume is decreasing and the rhombohedral distortion smoothly van- ishes, while when further increasing U the cubic structure remains stable and its volume increases. Why is there such a nonmonotonic behavior of the cell volume as a function of U? The increase in U results in a more localized nature of Ti(3d) states and a reduction in their hybridization with the neighboring O(2p) states, consequently diminishing the cova- lent character of Ti(3d)–O(2p) bonds. This, in turn, causes the Ti and O ions to occupy high-symmetry positions of the cubic phase. As a result of such an interplay between structural and electronic degrees of freedom, overall the unit cell volume tends to decrease. However, when U is increased beyond the critical value, a typical behavior is observed: Larger U values lead to even stronger localization of Ti(3d) states and a more pronounced ionic character of interactions, which con- sequently expands the lattice (increase in a). As mentioned earlier, the inclusion of intersite Hubbard V interactions between Ti(3d) and O(2p) states plays a crucial role in preserving the hybridizations between these states and, consequently, in maintaining the rhombohedral distortion of the lattice. As indicated in Sec. II, within PBEsol+U+V , the U value on Ti(3d) states is found to be 5.21 eV, which is higher than the value obtained within PBEsol+U . This behavior is common and is related to changes in electronic screening when computing both U and V [63,74,76]. Despite the fact that the computed U value within PBEsol+U+V exceeds the critical U value obtained within PBEsol+U , the cubic phase does not appear. This observation underscores the strong impact of intersite V interactions in stabilizing the rhombohedral phase, even though the V values are much smaller than the U values. In the next section, we utilize the optimized structural parameters obtained for each of the considered functionals (if not otherwise stated) to analyze the effects of Hubbard corrections on the electronic and vibrational properties of BTO. This analysis provides further insights into the role of extended Hubbard functionals with both U and V interactions in accurately predicting various properties of the rhombohe- dral phase of BTO. B. Projected density of states and band gap Using the optimized structural parameters for each func- tional, here we analyze the respective PDOS and the band-gap values. Figure 3 illustrates the total DOS and PDOS for Ti(3d) and O(2p) states. It can be seen that the valence band maximum primarily originates from the O(2p) states, while the conduction band minimum is mainly associated with the Ti(3d) states. However, due to the hybridization between Ti(3d) and O(2p) states, there is some contribution from the Ti(3d) states in the valence band region and, vice versa, the O(2p) states contribute to the conduction band region. This hybridization effect is crucial for the covalency of the Ti–O bonds in BTO, as previously discussed in Sec. I. The non- vanishing contribution of Ti(3d) states in the valence region 235171-5 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) FIG. 3. Computed total density of states and PDOS for Ti(3d) and O(2p) states using three functionals: (a) PBEsol, (b) PBEsol+U , and (c) PBEsol+U+V . The zero of the energy is set at the top of the valence bands in all cases. supports our earlier assertion that applying the Hubbard U correction to these states will not result in a negligible effect. As can be seen in Figs. 3(b) and 3(c), the application of Hub- bard corrections introduces noticeable changes in the PDOS; the primary impact is on the band-gap value, and additionally there are some subtle adjustments in the shape of the PDOS. As discussed earlier, applying the Hubbard U correction solely to the Ti(3d) states leads to their increased localization and a reduction in their hybridization with the neighbor- ing O(2p) states. In the scenario where this hybridization is strongly suppressed, there would be vanishing energy overlap between Ti(3d) and O(2p) states. Consequently, the valence band region would solely consist of O(2p) states, while the Ti(3d) states would remain entirely unoccupied. Indeed, it can be noticed in Fig. 3(b) that the application of U does slightly reduce the Ti(3d) weight in the valence region and increase it in the conduction band region, and the opposite effect is ob- served for the O(2p) states. In fact, the integrated intensity of the Ti(3d) states in the valence region is decreased by ∼10% when applying our first-principles Hubbard U correction. Let us now analyze the impact of the Hubbard corrections on the band gap. It is important to recall that DFT is not a theory for spectral properties. Importantly, in Ref. [104] it is shown that DFT + U with the Hubbard U parameter determined using linear-response theory, can markedly en- hance agreement with experimental band gaps. This improve- ment occurs when the Hubbard correction is applied to edge states (in cases where the system is already insulating at the DFT level) or states near the Fermi level (in situations where the system exhibits unphysical metallic behavior in standard DFT calculations), resulting in a Koopmans-like linearization [105]. Table II summarizes our results and compares them with previous studies and experimental data. To the best of our knowledge, there is no direct experimental data available for the band gap of the rhombohedral BTO at low temperatures. However, the experimental band gap for the cubic phase is known to be 3.2 eV at 420 K, and it increases with decreasing temperature at a rate of 4.5 × 10−4 eV/K [106], leading to an extrapolated value of 3.4 eV at 0 K. We employ this extrapo- lated value for comparisons, noting that precise measurements of the rhombohedral BTO’s band gap at low temperatures are required. We note that this extrapolated value was also used for comparisons in other works, e.g., in Refs. [27,107]. Our PBEsol calculation yields a band gap that underestimates the extrapolated experimental value by approximately 33%. This is consistent with the trends observed in LDA [100] and PBE [20] calculations. Hybrid functionals, such as HSE06 [27] and PBE0 [20], overestimate the band gap by approximately 9% and 44%, respectively. Thus, HSE06 has been the most accurate so far in predicting the band gap of rhombohedral BTO. Surprisingly, the band gap obtained with PBEsol+U is very close to the PBEsol band gap. This result can be attributed to the cancellation of two effects: the decrease in the gap due to the vanishing rhombohedral distortion and the increase in the gap caused by the application of the U correc- tion which generally leads to larger band gaps. To elucidate this, we performed a PBEsol calculation using the PBEsol+U geometry (i.e., cubic structure), resulting in a band gap of 1.60 eV, while performing a PBEsol+U calculation using the PBEsol+U geometry yields a band gap of 2.26 eV. On the other hand, PBEsol+U+V provides a band gap in ex- cellent agreement with the extrapolated reference value, with a deviation of approximately 2%. However, it is important to remark that the experimental band gap from Ref. [106] was determined through an optical absorption spectroscopy experiment, introducing the potential relevance of excitonic effects. Indeed, Ref. [3] demonstrated that incorporating ex- citonic effects (via solving the Bethe-Salpeter equation on top of GW ) in the cubic and tetragonal phases of BTO in- duces a shift at the optical absorption onset, resulting in a band-gap reduction of approximately 0.3–0.5 eV compared to GW calculations. Consequently, a similar order of magni- tude adjustment in the band gap might be expected for the rhombohedral phase of BTO due to excitonic effects. All the theoretical band gaps in Table II exclude excitonic effects, im- plying an expected reduction on the order of 0.3–0.5 eV. This TABLE II. Comparison of the band gap (in eV) computed in this work using PBEsol, PBEsol+U , and PBEsol+U+V and in previous works using (semi-)local and hybrid functionals and as measured in experiments. The experimental band gap was estimated using the extrapolation technique as explained in the main text. This work Ref. [100] Ref. [20] Ref. [27] Ref. [20] Ref. [106] PBEsol +U +U+V LDA PBE HSE06 PBE0 Expt. 2.28 2.26 3.46 2.23 2.7 3.69 4.9 3.4 235171-6 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) TABLE III. Comparison of the computed BEC of Ba, Ti, and O using PBEsol, PBEsol+U , and PBEsol+U+V and those from previous studies [22,25,28]. The xx and yy components are equal for Z∗ Ba and Z∗ Ti, and hence we report only the former. The BEC values in Refs. [22] and [25] are calculated at experimental lattice parameters [102], with relaxed atomic positions. Functional Z∗xx Ba Z∗zz Ba Z∗xx Ti Z∗zz Ti Z∗ O|| Z∗ O⊥ PBEsol 2.79 2.75 6.61 5.72 −5.11 −1.99 PBEsol+U 2.77 2.77 6.57 6.56 −5.20 −2.07 PBEsol+U+V 2.80 2.74 6.18 5.08 −4.47 −1.95 LDA [22] 2.79 2.74 6.54 5.61 −5.05 −1.98 LDA [25] 2.78 2.74 6.61 5.77 −5.10 −1.99 B1-WC [28] 2.71 2.68 6.41 5.55 −4.94 −1.95 adjustment aligns HSE06 more closely with the experimental band gap, whereas the PBEsol+U+V band gap is likely to de- viate more significantly from the experimental value. Overall, while PBEsol+U+V does not exhibit significant improve- ments in structural properties compared to PBEsol (which already shows remarkable agreement with experiments), it significantly enhances the accuracy in predicting the band-gap value. C. Born effective charges and dielectric tensor In this section, we present and discuss the computed BEC and dielectric tensor using the three functionals considered in this work. Table III provides a comparison of our calcu- lations with previous computational studies [22,25,28]. In the rhombohedral crystal structure, the computed Z∗ Ti and Z∗ Ba are diagonal tensors with xx and yy components being equal due to symmetry, while the BEC of three O atoms have nonzero off-diagonal components. To simplify the analysis, we di- agonalize Z∗ O and denote its largest and doubly degenerate smallest eigenvalues as O|| and O⊥, respectively (which are the same for three O atoms). O|| represents displacements of O ions along the Ti–O bond, while O⊥ refers to displacements perpendicular to the bond. As shown in Table III, for all functionals, the values of Z∗ Ti diagonal components and Z∗ O|| are substantially larger than their nominal values (+4 and −2, respectively) [22,30]. This observation is consistent with other ABO3 perovskites and is a fingerprint of covalency of the Ti– O bonds [21,22,28,53,108]. In contrast, the computed Z∗ Ba and Z∗ O⊥ values are not significantly different from their nominal values (+2 for Ba), suggesting an ionic character of the Ba–O bonds [21,22,28]. Applying the onsite U correction alone does not lead to large changes in the BEC for all the elements as compared to PBEsol (except for Z∗zz Ti ). When considering PBEsol+U+V , however, we find that Z∗xx Ti , Z∗zz Ti , and Z∗ O|| are smaller by 7–13% compared to PBEsol. This reduction in the BEC due to +U + V corrections arises from the redistribution of electronic charge between Ti and O ions and sensitivity to changes in structural parameters (see Table I) [22]. Let us now discuss the dielectric constant of the low- temperature rhombohedral BTO. Again due to symmetry, the xx and yy components of the dielectric tensor are equal. In Table IV, we present a comparison of our computed values TABLE IV. Diagonal components of the dielectric tensor ε∞ xx and ε∞ zz (note that ε∞ yy = ε∞ xx ) and its average value [〈ε〉 = (ε∞ xx + ε∞ yy + ε∞ zz )/3] as computed in this work using PBEsol, PBEsol+U , and PBEsol+U+V and as obtained in previous works [20,28,30,109,110] and measured in experiments [111]. In Ref. [110] a scissor � correction was used. Functional ε∞ xx ε∞ zz 〈ε〉 PBEsol 6.14 5.69 5.99 PBEsol+U 5.88 5.88 5.88 PBEsol+U+V 5.54 5.00 5.36 LDA [109] 6.16 5.69 6.00 LDA [20] 6.05 5.82 5.97 LDA+� [110] 5.57 5.51 5.55 PBE [20] 5.42 4.80 5.21 PBEsol [30] 6.11 5.63 5.95 PBE0 [20] 4.62 4.13 4.46 B1-WC [28] 4.96 4.60 4.84 Expt. [111] 5.20 with those from previous studies [20,28,30,109,110] and ex- periments [111]. First, we want to remark that quite often [20,28,101,112] in the literature the theoretical components of the dielectric tensor of the rhombohedral BTO are compared with the following experimental values: ε∞ xx = 6.19 and ε∞ zz = 5.88 [113]. However, as was rightly pointed out in Ref. [26], the experimental data of Ref. [113] corresponds to the room- temperature tetragonal phase of BTO, and thus it should not be used for comparisons with the low-temperature rhombo- hedral phase. On the other hand, we follow Ref. [110] and use an average experimental dielectric constant 〈ε〉 = 5.20 which was obtained from the average refractive index of 2.28 [111]. However, it is crucial to emphasize that more accurate measurements of the dielectric tensor components for the low- temperature rhombohedral phase of BTO are needed. As can be seen in Table IV, our PBEsol dielectric tensor shows good agreement with previous LDA [22] and PBEsol [30] studies. The average dielectric constant obtained using PBEsol overes- timates the experimental value by ∼15%. When including the +U correction, 〈ε〉 slightly decreases, but it still overestimates the experimental value by ∼13%. However, it is essential to consider that the PBEsol+U geometry corresponds to the cu- bic phase at 0 K, making this comparison less straightforward. We note in passing that the experimental 〈ε〉 of the cubic phase of BTO is 5.40 [114]. The PBEsol+U+V functional yields 〈ε〉 that is in substantially closer agreement with the experimental value, overestimating it by only ∼3%. Hence, PBEsol+U+V provides the most accurate prediction for 〈ε〉 compared to PBEsol and PBEsol+U . Finally, we note that hybrid functionals PBE0 [20] and B1-WC [28] underestimate 〈ε〉 by ∼14% and ∼7%, respectively, while 〈ε〉 from PBE [20] is by chance in excellent agreement with the experimental value. D. Phonons Here we explore the lattice vibrational properties of BTO using the three considered functionals. A group-theory anal- ysis reveals that the zone-center phonon frequencies of 235171-7 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) TABLE V. Phonon frequencies of the rhombohedral BTO (in cm−1) at the Brillouin zone center computed using PBEsol on top of the PBEsol geometry (PBEsol@PBEsol), and PBEsol+U+V on top of the PBEsol geometry (PBEsol+U+V @PBEsol) and on top of the PBEsol+U+V geometry (PBEsol+U+V @PBEsol+U+V ). The results from previous computational studies are also shown. The symmetry of the phonon modes is specified in the first column, while the common labeling [25] of the TO and LO modes is specified in the second column. The label “q” denotes the direction of the phonon wave vector q, which is either along the x or z axis in reciprocal space of the rhombohedral BTO reference system (see Sec. II). The 3E (TO3) and 2E (LO2) modes are degenerate, while the A2 mode is silent. Mode Label q PBEsol PBEsol+U+V PBEsol+U+V Ref. [20] Ref. [25] Ref. [20] Ref. [20] Ref. [101] @PBEsol @PBEsol @PBEsol+U+V LDA LDA PBE PBE0 B1-WC TO 1E TO1 x, z 166 168 162 145 163 165 180 125 1A1 TO1 x 170 171 165 191 167 167 181 192 2E TO2 x, z 207 243 264 191 210 264 264 217 3E TO3 x, z 294 290 288 306 293 299 319 318 2A1 TO2 x 262 292 310 200 259 309 332 285 4E TO4 x, z 472 483 478 489 470 474 493 497 3A1 TO3 x 516 537 548 508 512 544 571 542 LO 1E LO1 x 176 180 177 190 174 176 190 193 1A1 LO1 z 180 186 183 193 178 182 197 199 2E LO2 x, z 294 290 288 306 293 299 319 318 3E LO3 x 441 443 437 460 441 441 468 472 2A1 LO2 z 460 464 463 471 461 467 498 491 3A1 LO3 z 681 705 701 707 676 689 739 740 4E LO4 x 691 714 711 713 687 705 750 730 A2 x, z 277 275 270 297 277 279 301 307 rhombohedral BTO can be decomposed as �opt = 3A1 + A2 + 4E [112], where both A1 and E modes are infrared and Raman active, while the A2 mode is silent. The zone-center phonon frequencies obtained from PBEsol and PBEsol+U+V are compared with previous computational studies [20,25,101] in Table V. In order to scrutinize the impact of lattice geome- try (a and α) on phonon frequencies, we present the outcomes of PBEsol+U+V phonon calculations for both the PBEsol+U+V optimized geometry (referred to as PBEsol+U+V @PBEsol+U+V ) and the PBEsol optimized geometry (PBEsol+U+V @PBEsol). It is worth noting that in both cases, atomic positions are optimized using PBEsol+U+V to ensure that forces acting on atoms are vanishing. Conversely, PBEsol+U is omitted from this com- parison due to its inability to sustain the rhombohedral phase of BTO (refer to Sec. III A). The nonanalytic com- ponent of the dynamical matrix up to the dipolar order is taken into account, which leads to the splitting between the longitudinal optical (LO) and transverse optical (TO) modes [83]. This nonanalytic component is derived using the BEC and dielectric tensor, as discussed in Sec. III C. In the case of PBEsol+U+V @PBEsol, the BEC and di- electric tensor are not recalculated; instead, those from the PBEsol+U+V @PBEsol+U+V case are used. The phonon mode data presented in Table V demonstrate distinctive behaviors among the considered functionals in this study and in previous works. We have not compared our zone-center phonons to the experimental data taken on single crystals in Refs. [94,95] as these samples are polydomain with the z axis oriented along the c axis of the tetragonal BTO reference system. This reorientation leads both to a mixing of the phonon modes listed in Table V (see Secs. II and III E for more details) and to potential differences in the resulting LO- TO splitting. Experimental data for rhombohedral BTO single crystals with a single ferroelectric domain, specifically for q directions along the x and z axes in reciprocal space of the rhombohedral BTO reference system (where, in real space, z is parallel to the C3 rhombohedral axis), would be highly desirable for evaluating the accuracy of computational data obtained from various functionals, as presented in Table V. Let us now shift to the phonon dispersion of the rhombohedral BTO lattice. In Fig. 4, we present a comprehensive comparison of the computed phonon frequencies as a function of the phonon wave vector q along high-symmetry directions in the Brillouin zone obtained from PBEsol@PBEsol, PBEsol+U+V @PBEsol+U+V , and PBEsol+U+V @PBEsol. It can be seen that there is a discontinuity in the phonon dispersion near � for certain optical phonon modes, due to the LO-TO splitting that depends on the direction of the phonon wave vector q in the long-wavelength limit (q → 0) [83]. The three acoustic branches display remarkable congruence across all cases, except in proximity to the � point. This correspondence stems from the fact that, up to around 150 cm−1, the primary contribution to the vibrational density of states (vDOS) can be attributed to the motion of Ba ions [100], which remains largely unaffected by the Hubbard corrections. However, above ∼150 cm−1, the vDOS is mainly influenced by both Ti and O atoms, leading to pronounced differences in the phonon bands due to the Hubbard corrections. Specifically, within PBEsol+U+V @PBEsol+U+V and PBEsol+U+V @PBEsol, the highest phonon bands are shifted to higher frequencies (compared 235171-8 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) FIG. 4. (a) Phonon dispersion of the rhombohedral BTO computed using PBEsol and PBEsol+U+V and using different optimized geometries. (a) PBEsol+U+V calculation using the PBEsol+U+V lattice geometry (PBEsol+U+V @PBEsol+U+V ) and (b) PBEsol+U+V calculation using the PBEsol lattice geometry (PBEsol+U+V @PBEsol). On both panels the PBEsol results are exactly the same and were obtained using the PBEsol lattice geometry. to PBEsol@PBEsol), consistently with similar obser- vations in other materials [115,116]. The effect of the lattice parameters on the PBEsol+U+V phonon calculations can be seen by comparing Figs. 4(a) and 4(b). Although PBEsol+U+V @PBEsol+U+V and PBEsol+U+V @PBEsol phonon dispersions exhibit qualitative similarity, notable quantitative disparities exist. It is important to pay special attention to the behavior of the acoustic phonon bands near the � point. As previously noted in Ref. [32], piezoelectric materials like BTO exhibit imaginary phonon frequencies around � when considering the nonanalytic part of the dynamical matrix only up to the dipolar order. To address this issue, higher multipolar orders of the nonanalytic correction to the dynamical matrix must be included [32], which are though not considered in our study. In our calculations, we find sizable imaginary phonon frequencies along the �-F high-symmetry direction when using the PBEsol functional. However, these imaginary frequencies are substantially reduced when using PBEsol+U+V , becoming imperceptible on the plot due to the large frequency scale. Such differences between the PBEsol and PBEsol+U+V results do not arise from the variations in the lattice geometry, because both PBEsol+U+V @PBEsol+U+V and PBEsol+U+V @PBEsol do not have a large negative bump along �-F like in the PBEsol@PBEsol case. Hence, such differences between PBEsol and PBEsol+U+V originate purely from changes in the electronic structure due to the U and V corrections. Furthermore, both PBEsol and PBEsol+U+V yield minimal imaginary phonon frequencies close to � along the �-L direction, although these frequencies are not readily visible on the large frequency scale shown in Fig. 4. These observations suggest that higher-order corrections to the dynamical matrix are essential in BTO and similar materials, regardless of the functional employed. As was shown in Ref. [117] for the cubic BTO, there is a giant LO-TO splitting. Conversely, in the rhombohedral BTO we find that there is a strong mixing of LO modes when incorporating the nonanalytic part to the dynamical matrix. Therefore, there is no one-to-one connection between the highest LO mode and one of the TO modes, and hence it is not straightforward how to determine the LO–TO splitting. Nevertheless, we believe that it is useful to point out what is the optical phonon bandwidth, i.e., the difference between the highest LO and lowest TO modes, which is 525 cm−1 in PBEsol and about 550 cm−1 in PBEsol+U+V . It is instructive to compare the phonon dispersion for the cubic phase of BTO using PBEsol and PBEsol+U . Figure 5 shows such a comparison when using the PBEsol+U lattice geometry for both functionals (i.e., PBEsol@PBEsol+U and PBEsol+U@PBEsol+U , respectively). In our PBEsol calcu- lations for the cubic BTO, we observe large imaginary phonon frequencies at various points in the Brillouin zone, including the high-symmetry points �, M, and X, which is consis- tent with previous studies [29,33,118,119]. When employing PBEsol+U with the first-principles U value, the imaginary phonon frequencies vanish everywhere in the Brillouin zone except in the vicinity of the � point. Since soft phonon modes measure the structural instability, this result shows that PBEsol+U tends to stabilize the cubic phase at 0 K. A similar observation was reported in Ref. [27] (see Fig. S6 in this reference). However, it is worth noting that no imaginary FIG. 5. Phonon dispersion of the cubic BTO computed using PBEsol and PBEsol+U , both using the PBEsol+U lattice geometry. 235171-9 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) phonon frequencies were found in that study, even around �. Despite the fact that the value of the Hubbard U parameter is extremely similar in both studies (4.52 eV in our case vs 4.49 eV in Ref. [27]), the remaining differences might be at- tributed to the following factors: (i) we use PBEsol+U in this study, while PBE + U was used in Ref. [27]; (ii) we employ Löwdin-orthogonalized atomic orbitals, while in Ref. [27], the PAW functions were used to build the Hubbard projectors; (iii) different pseudopotentials were used in both studies; and (iv) we use the PBEsol+U optimized lattice parameter of 3.990 Å (see Table I), whereas in Ref. [27], a lattice parameter of 4.077 Å was used. Although it is not clear which factor plays the dominant role, the conclusion remains that soft phonon modes around � are highly sensitive to numerical details of calculations. Additionally, we verified that using larger values of the U parameter (e.g., testing with U = 6 eV) within our computational setup completely eliminates the soft phonon modes around �. Finally, we address the comparison between the tendency of PBEsol+U to stabilize the cubic phase of BTO at 0 K and the effect of strain within PBEsol to achieve a similar stabilization. Previous studies have shown that by compress- ing the lattice (reducing the lattice parameter by ∼3%), the cubic phase of BTO can be stabilized at 0 K when using (semi-)local functionals since no imaginary phonon fre- quencies are present [28,33,118]. This phenomenon can be explained as follows: When BTO is compressed, the in- teratomic distances decrease, altering the energy landscape. Specifically, the minima of the symmetric double well of the energy potential move closer together, while the saddle point becomes lower, resulting in the stabilization of the high- symmetry cubic phase. On the other hand, when applying the first-principles U correction, the lattice tends to expand, as expected (see the discussion in Sec. III A). We found that the optimized lattice parameter for the cubic BTO using PBEsol+U is increased only by ∼0.3% compared to PBEsol [120]. Therefore, from a purely geometrical perspective, it is not evident why PBEsol+U would stabilize the cubic phase of BTO. As we discussed earlier, the reason lies purely in the electronic behavior. The localization of the Ti(3d) states suppresses the Ti(3d)–O(2p) hybridization that is crucial for the covalency of the Ti–O bonds and stabilization of the rhombohedral phase of BTO. Consequently, as the strength of the Hubbard U correction is increased, the dynamical sta- bility of the rhombohedral phase diminishes, while that of the cubic phase increases. Thus, this intricate electronic effect is responsible for the observed stabilization of the cubic phase by PBEsol+U . E. Raman spectra In this section, we present a comparison between the computed and experimental nonresonant Raman spectra of rhombohedral BTO. The spectra were experimentally mea- sured at 10 K [94] and 80 K [95,121]. Although all these experimental spectra display similar qualitative features, we choose to compare with the higher-resolution Raman spec- trum on single crystals of BTO from Ref. [95]. It is important to note that the Raman spectrum of single crystals in Ref. [95] closely resembles the spectrum presented in Ref. [94]. The FIG. 6. Raman spectra computed using PBEsol, PBEsol+U , and PBEsol+U+V on top of different optimized geometries. The ex- perimental Raman spectra were obtained for single crystals at 80 K [95]. All spectra are shifted vertically for the sake of clarity. Vertical dashed lines correspond to positions of peaks in the experimental Raman spectrum, and they are numbered for convenience. latter reference explicitly states that the measurements were conducted using the backscattering geometry with the z(xx)z̄ setup, where z is along the c axis of the tetragonal BTO phase. We use the same setup in our calculations, as de- tailed in Sec. II. On the other hand, we are aware of only two theoretical studies of the Raman spectra of the rhom- bohedral BTO, and they are based on DFPT using the LDA functional [25,98]. In Ref. [25] the authors considered the x(zz)y and z(xy)z̄ setups for single crystals with the z axis aligned along the C3 rhombohedral axis, while in Ref. [98] the averaged Raman spectra of polycrystalline rhombohedral BTO were computed. Hence, the computational predictions of these studies are not directly comparable with our Raman spectrum calculations for single crystals and the experiments of Refs. [94,95]. Figure 6 illustrates the computed Raman spectra using the three considered functionals along with the experimental spectrum from Ref. [95]. The Hubbard corrections exert a twofold influence on the Raman spectra: (i) via alterations in the optimized lattice geometry and (ii) via changes in the electronic structure. To disentangle these distinct effects, separate Raman spectrum computations were executed for PBEsol+U and PBEsol+U+V , each employing two dis- tinct lattice geometries: one optimized utilizing PBEsol+U and PBEsol+U+V , correspondingly, and the other using the PBEsol optimized lattice geometry. 235171-10 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) We begin our analysis of Fig. 6 by comparing the Raman spectrum obtained from PBEsol@PBEsol with the exper- imental one. The peak 3 is well reproduced (its position coincides very well with the experimental one), peak 2 is red-shifted by 19 cm−1, peak 1 is barely visible and it is blue-shifted by 5 cm−1, while peak 4 is not visible at all. Next, peaks 5, 6, and 7 are all blue-shifted by 15, 19, and 28 cm−1, respectively. We recall that here we use a constant broadening with the Lorenzian function, hence we do not account for the various lifetimes of all peaks, which might explain the differences in the widths of the peaks. The differences in the positions of peaks between PBEsol@PBEsol and experiments might be attributed to the following factors: (i) errors stem- ming from the approximations when using PBEsol@PBEsol and (ii) our calculations are performed assuming a single ferroelectric domain while measurements in Ref. [95] were carried out on a polydomain single crystal sample. Overall, our PBEsol@PBEsol spectrum looks qualitatively similar to the experimental spectrum, although some disparities in peak positions and relative intensities persist. Let us now delve into the analysis of the Raman spectra computed employing the Hubbard corrections. As depicted in Fig. 6, the Raman spectrum computed using PBEsol+U+V with the PBEsol+U+V optimized lattice geometry (PBEsol+U+V @PBEsol+U+V ) displays a no- table deterioration in comparison to PBEsol@PBEsol and the experiment. Namely, peaks 2 and 3 are red-shifted even more, peak 5 does not show significant changes, peak 6 is now red-shifted with respect to the experimen- tal peak, while peak 7 is now in much closer agreement with the experimental peak position. A substantial improve- ment emerges when the Raman spectrum is calculated using PBEsol+U+V with the lattice geometry optimized using PBEsol (PBEsol+U+V @PBEsol). More specifically, even though the peaks 2 and 3 are still more red-shifted com- pared to experiments than in the PBEsol@PBEsol case, but now the peaks 5, 6, and 7 are in remarkable agree- ment with the experimental peaks. Therefore, while some peaks are more accurately described using PBEsol@PBEsol, others are better described using PBEsol+U+V @PBEsol, and hence none of these approximations provides the best description of all peaks in the Raman spectrum simulta- neously. In contrast, the Raman spectrum computed using PBEsol+U on top of the PBEsol+U lattice geometry (PBEsol+U@PBEsol+U ) diverges dramatically from the experimental one. To ascertain the influence of the under- lying lattice geometry on the resulting Raman spectrum, we conducted a PBEsol+U Raman spectrum calculation employing the PBEsol lattice geometry. This resultant spec- trum, denoted as PBEsol+U@PBEsol, still exhibits poor alignment with the experimental data. This infers that the de- ficiency in the Raman spectrum of PBEsol+U@PBEsol+U is not predominantly attributed to the PBEsol+U lattice ge- ometry, which remains cubic instead of rhombohedral (as elaborated in Sec. III A). The primary rationale for the ineffectiveness of PBEsol+U resides in the fact that even with the adoption of the PBEsol lattice geometry, the atomic positions regress to the high-symmetry cubic orientations during the PBEsol+U relaxation, thereby influencing the electronic structure significantly due to the +U correction. Consequently, the covalency of the Ti(3d)–O(2p) bonding is suppressed, leading to a drastic deterioration of the Raman spectrum. Overall, we ascertain that PBEsol@PBEsol and PBEsol+U+V @PBEsol yield the most accurate description of the Raman spectrum for rhombohedral BTO. The inclusion of the +U + V corrections tends to red-shift low-wave-number Raman peaks but yields an excellent description of the high-wave-number Raman peaks, provided the PBEsol geometry is used. In stark contrast, the application of the +U correction alone proves to be highly detrimental, leading to a Raman spectrum that deviates dramatically from the experimental data. IV. CONCLUSIONS We have presented a detailed first-principles investigation of the low-temperature rhombohedral phase of BTO, focus- ing on its structural, electronic, and vibrational properties using PBEsol, PBEsol+U , and PBEsol+U+V . The onsite U and intersite V Hubbard parameters were computed using density-functional perturbation theory in a basis of Löwdin- orthogonalized atomic orbitals. Our findings reveal that the application of the Hubbard U correction to Ti(3d) states lo- calizes them and suppresses the Ti(3d)–O(2p) hybridizations, thus favoring the stabilization of the cubic phase in agree- ment with previous studies [27]. The inclusion of intersite Hubbard V interactions between the Ti(3d) and O(2p) states preserves covalency of the Ti–O bonds, ultimately leading to the stabilization of the rhombohedral phase in agreement with experimental observations. The optimized geometry in PBEsol+U+V is found to be slightly less accurate than in PBEsol, resulting in a detri- mental impact on the lattice vibrational properties. The Born effective charges for Z∗ Ti components and Z∗ O|| are found to be smaller by 7–13% compared to PBEsol results. The PDOS in PBEsol+U+V is qualitatively similar to PBEsol, while the PBEsol+U+V band gap and dielectric constant show good agreement with experiments by surpassing the accuracy of PBEsol predictions. The zone-center phonon frequencies, representing the positions of Raman peaks, are found to be very sensitive to the underlying geometry. The PBEsol and PBEsol+U+V functionals provide the Raman spectra in sat- isfactory agreement with the experimental one, provided that the PBEsol geometry is used. Conversely, PBEsol+U tends to stabilize the cubic BTO at 0 K due to the suppression of the Ti(3d)–O(2p) hy- bridizations, resulting in various properties that are in poorer agreement with experiments compared to PBEsol+U+V . Strikingly, the Raman spectrum computed using PBEsol+U differs dramatically from the experimental one. Hence, the application of the +U correction alone is found to be highly detrimental in the rhombohedral BTO. Therefore, our study has uncovered the crucial signif- icance of intersite Hubbard interactions in preserving the covalent features present in the rhombohedral phase of BTO. These findings could potentially extend beyond BTO and have broader implications for other ABO3 perovskites exhibiting similar properties. However, further investigations employing PBEsol+U+V for other materials are needed to fully explore 235171-11 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) the potential of this approach. Thus, our study contributes to laying the foundation for future research in this direction, offering new insights into tailoring the properties of these ma- terials. The predictive power of PBEsol+U+V holds promise for breakthroughs in perovskite materials, unlocking numer- ous opportunities for advanced technologies and applications. ACKNOWLEDGMENTS We thank Lorenzo Monacelli, Atsushi Togo, and Eu- gene Roginskii for fruitful discussions. This research was supported by the NCCR MARVEL, a National Centre of Competence in Research, funded by the Swiss Na- tional Science Foundation (Grant No. 205602). Computer time was provided by the Swiss National Supercomput- ing Centre (CSCS) under Project No. s1073, the Centre for High Performance Computing (CHPC), South Africa, and the supercomputer Lise and Emmy at NHR@ZIB and NHR@Göttingen as part of the NHR infrastructure under the projects hbc00053 and hbi00059. The authors gratefully acknowledge support from the Deutsche Forschungsgemein- schaft (DFG) under Germany’s Excellence Strategy (EXC 2077, Grant No. 390741603, University Allowance, Univer- sity of Bremen) and Lucio Colombi Ciacchi, the host of the “U Bremen Excellence Chair Program,” and the Abdus Salam International Centre for Theoretical Physics. [1] W. Peng, J. A. Zorn, J. Mun, M. Sheeraz, C. J. Roh, J. Pan, B. Wang, K. Guo, C. W. Ahn, Y. Zhang, K. Yao, J. S. Lee, J.-S. Chung, T. H. Kim, L.-Q. Chen, M. Kim, L. Wang, and T. W. Noh, Constructing polymorphic nanodomains in BaTiO3 films via epitaxial symmetry engineering, Adv. Funct. Mater. 30, 1910569 (2020). [2] F. Rubio-Marcos, D. A. Ochoa, A. Del Campo, M. A. García, G. R. Castro, J. F. Fernández, and J. E. García, Reversible optical control of macroscopic polarization in ferroelectrics, Nat. Photon. 12, 29 (2018). [3] S. Sanna, C. Thierfelder, S. Wippermann, T. P. Sinha, and W. G. Schmidt, Barium titanate ground- and excited-state properties from first-principles calculations, Phys. Rev. B 83, 054112 (2011). [4] Edited by R. Waser, Nanoelectronics and Information Technol- ogy: Advanced Electronic Materials and Novel Devices, 3rd ed. (Wiley-VCH, Weinheim, 2003). [5] U. De, K. R. Sahu, and A. De, Ferroelectric materials for high temperature piezoelectric applications, Solid State Phenom. 232, 235 (2015). [6] S. K. Mahadeva, K. Walus, and B. Stoeber, Piezoelectric paper for physical sensing applications, in Proceedings of the 28th IEEE International Conference on Micro Electro Mechanical Systems (MEMS ’15) (IEEE, Los Alamitos, CA, 2015), p. 861. [7] G. Vasta, T. J. Jackson, and E. Tarte, Electrical properties of BaTiO3 based ferroelectric capacitors grown on oxide sacrifi- cial layers for micro-cantilevers applications, Thin Solid Films 520, 3071 (2012). [8] H. F. Kay, Preparation and properties of crystals of barium titanate, BaTiO3, Acta Crystallogr. 1, 229 (1948). [9] W. J. Merz, The electric and optical behavior of BaTiO3 single- domain crystals, Phys. Rev. 76, 1221 (1949). [10] A. Devonshire, Theory of ferroelectrics, Adv. Phys. 3, 85 (1954). [11] S.-E. Park, S. Wada, L. E. Cross, and T. R. Shrout, Crystallographically engineered BaTiO3 single crystals for high-performance piezoelectrics, J. Appl. Phys. 86, 2746 (1999). [12] A. von Hippel, R. G. Breckenridge, F. G. Chesley, and L. Tisza, High dielectric constant ceramics, Ind. Eng. Chem. 38, 1097 (1946). [13] I. Bersuker, On the origin of ferroelectricity in perovskite-type crystals, Phys. Lett. 20, 589 (1966). [14] A. Chaves, F. C. S. Barreto, R. A. Nogueira, and B. Zeks, Thermodynamics of an eight-site order-disorder model for ferroelectrics, Phys. Rev. B 13, 207 (1976). [15] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964). [16] W. Kohn and L. J. Sham, Self-consistent equations includ- ing exchange and correlation Effects, Phys. Rev. 140, A1133 (1965). [17] R. E. Cohen and H. Krakauer, Lattice dynamics and origin of ferroelectricity in BaTiO3: Linearized-augmented-plane- wave total-energy calculations, Phys. Rev. B 42, 6416 (1990). [18] R. D. King-smith and D. Vanderbilt, A first-principles pseu- dopotential investigation of ferroelectricity in barium titanate, Ferroelectr. 136, 85 (1992). [19] J. J. Wang, P. P. Wu, X. Q. Ma, and L. Q. Chen, Temperature- pressure phase diagram and ferroelectric properties of BaTiO3 single crystal based on a modified Landau potential, J. Appl. Phys. 108, 114105 (2010). [20] R. A. Evarestov and A. V. Bandura, First-principles calcula- tions on the four phases of BaTiO3, J. Comput. Chem. 33, 1123 (2012). [21] E. Goh, L. Ong, T. Yoon, and K. Chew, Structural and response properties of all BaTiO3 phases from density functional theory using the projector-augmented-wave methods, Comput. Mater. Sci. 117, 306 (2016). [22] P. Ghosez, X. Gonze, P. Lambin, and J.-P. Michenaud, Born effective charges of barium titanate: Band-by-band decompo- sition and sensitivity to structural features, Phys. Rev. B 51, 6765(R) (1995). [23] P. Ghosez, X. Gonze, and J.-P. Michenaud, First-principles characterization of the four phases of barium titanate, Ferroelectr. 220, 1 (1999). [24] A. W. Hewat, Structure of rhombohedral ferroelectric barium titanate, Ferroelectr. 6, 215 (1973). [25] P. Hermet, M. Veithen, and P. Ghosez, Raman scattering intensities in BaTiO3 and PbTiO3 prototypical ferroelectrics from density functional theory, J. Phys.: Condens. Matter 21, 215901 (2009). [26] A. Mahmoud, A. Erba, K. E. El-Kelany, M. Rérat, and R. Orlando, Low-temperature phase of BaTiO3: Piezoelectric, dielectric, elastic, and photoelastic properties from ab initio simulations, Phys. Rev. B 89, 045103 (2014). 235171-12 https://doi.org/10.1002/adfm.201910569 https://doi.org/10.1038/s41566-017-0068-1 https://doi.org/10.1103/PhysRevB.83.054112 https://doi.org/10.4028/www.scientific.net/SSP.232.235 https://doi.org/10.1016/j.tsf.2011.11.032 https://doi.org/10.1107/S0365110X4800065X https://doi.org/10.1103/PhysRev.76.1221 https://doi.org/10.1080/00018735400101173 https://doi.org/10.1063/1.371120 https://doi.org/10.1021/ie50443a009 https://doi.org/10.1016/0031-9163(66)91127-9 https://doi.org/10.1103/PhysRevB.13.207 https://doi.org/10.1103/PhysRev.136.B864 https://doi.org/10.1103/PhysRev.140.A1133 https://doi.org/10.1103/PhysRevB.42.6416 https://doi.org/10.1080/00150199208016068 https://doi.org/10.1063/1.3504194 https://doi.org/10.1002/jcc.22942 https://doi.org/10.1016/j.commatsci.2016.01.037 https://doi.org/10.1103/PhysRevB.51.6765 https://doi.org/10.1080/00150199908007992 https://doi.org/10.1080/00150197408243970 https://doi.org/10.1088/0953-8984/21/21/215901 https://doi.org/10.1103/PhysRevB.89.045103 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) [27] N. Tsunoda, Y. Kumagai, and F. Oba, Stabilization of small polarons in BaTiO3 by local distortions, Phys. Rev. Mater. 3, 114602 (2019). [28] Y.-S. Seo and J. S. Ahn, Pressure dependence of the phonon spectrum in BaTiO3 polytypes studied by ab initio calcula- tions, Phys. Rev. B 88, 014114 (2013). [29] Y. Zhang, J. Sun, J. P. Perdew, and X. Wu, Comparative first- principles studies of prototypical ferroelectric materials by LDA, GGA, and SCAN meta-GGA, Phys. Rev. B 96, 035143 (2017). [30] S. Mellaerts, J. W. Seo, V. Afanas’ev, M. Houssa, and J.-P. Locquet, Origin of supertetragonality in BaTiO3, Phys. Rev. Mater. 6, 064410 (2022). [31] S. F. Yuk, K. C. Pitike, S. M. Nakhmanson, M. Eisenbach, Y. W. Li, and V. R. Cooper, Towards an accurate description of perovskite ferroelectrics: Exchange and correlation effects, Sci. Rep. 7, 43482 (2017). [32] M. Royo, K. R. Hahn, and M. Stengel, Using high multipolar orders to reconstruct the sound velocity in piezoelectrics from lattice dynamics, Phys. Rev. Lett. 125, 217602 (2020). [33] M. Kotiuga, S. Halilov, B. Kozinsky, M. Fornari, N. Marzari, and G. Pizzi, Microscopic picture of paraelectric perovskites from structural prototypes, Phys. Rev. Res. 4, L012042 (2022). [34] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and sur- faces, Phys. Rev. Lett. 100, 136406 (2008). [35] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys. 110, 6158 (1999). [36] J. Heyd, G. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003). [37] J. Heyd, G. Scuseria, and M. Ernzerhof, Erratum: “Hy- brid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)], J. Chem. Phys. 124, 219906 (2006). [38] K. Moseni, R. B. Wilson, and S. Coh, Electric field control of phonon angular momentum in perovskite BaTiO3, Phys. Rev. Mater. 6, 104410 (2022). [39] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I , Phys. Rev. B 44, 943 (1991). [40] V. Anisimov, F. Aryasetiawan, and A. Lichtenstein, First- principles calculations of the electronic structure and spectra of strongly correlated systems: The LDA + U method, J. Phys.: Condens. Matter 9, 767 (1997). [41] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the struc- tural stability of nickel oxide: An LSDA + U study, Phys. Rev. B 57, 1505 (1998). [42] H. J. Kulik, M. Cococcioni, D. A. Scherlis, and N. Marzari, Density functional theory in transition-metal chemistry: A self-consistent Hubbard U approach, Phys. Rev. Lett. 97, 103001 (2006). [43] H. Kulik and N. Marzari, A self-consistent Hubbard U density-functional theory approach to the addition-elimination reactions of hydrocarbons on bare FeO+, J. Chem. Phys. 129, 134314 (2008). [44] A. Cohen, P. Mori-Sánchez, and W. Yang, Insights into current limitations of density functional Theory, Science 321, 792 (2008). [45] J. Zhang, G. Gou, and B. Pan, Study of phase stabil- ity and hydride diffusion mechanism of BaTiO3 oxyhy- dride from first-principles, J. Phys. Chem. C 118, 17254 (2014). [46] P. Erhart, A. Klein, D. Aberg, and B. Sadigh, Efficacy of the DFT + U formalism for modeling hole polarons in perovskite oxides, Phys. Rev. B 90, 035204 (2014). [47] X. Liu, T. S. Bjorheim, and R. Haugsrud, Formation and mi- gration of hydride ions in BaTiO3−xHx oxyhydride, J. Mater. Chem. A 5, 1050 (2017). [48] H. Chen and A. Millis, Design of new Mott multiferroics via complete charge transfer: promising candidates for bulk photovoltaics, Sci. Rep. 7, 6142 (2017). [49] Y. Watanabe, Calculation of strained BaTiO3 with different exchange correlation functionals examined with criterion by Ginzburg-Landau theory, uncovering expressions by crystal- lographic parameters, J. Chem. Phys. 148, 194702 (2018). [50] S. Majumder, P. Basera, M. Tripathi, R. Choudhary, S. Bhattacharya, K. Bapna, and D. Phase, Elucidating the origin of magnetic ordering in ferroelectric BaTiO3−δ thin film via electronic structure modification, J. Phys.: Condens. Matter 31, 205001 (2019). [51] N. Din, T. Jiang, S. Gholam-Mirzaei, M. Chini, and V. Turkowski, Electron–electron correlations and structural, spectral and polarization properties of tetragonal BaTiO3, J. Phys.: Condens. Matter 32, 475601 (2020). [52] V. Dwij, B. De, G. Sharma, D. Shukla, M. Gupta, R. Mittal, and V. Sathe, Revisiting eigen displacements of tetragonal BaTiO3: Combined first principle and experimental investiga- tion, Phys. B: Condens. Matter 624, 413381 (2022). [53] R. E. Cohen, Origin of ferroelectricity in perovskite oxides, Nature (Lond.) 358, 136 (1992). [54] K. Tkacz-Śmiech, A. Koleżyński, and W. S. Ptak, Chem- ical bond in ferroelectric perovskites, Ferroelectr. 237, 57 (2000). [55] R. Resta, M. Posternak, and A. Baldereschi, Towards a quan- tum theory of polarization in ferroelectrics: The case of KNbO3, Phys. Rev. Lett. 70, 1010 (1993). [56] V. L. Campo Jr. and M. Cococcioni, Extended DFT + U + V method with on-site and inter-site electronic interactions, J. Phys.: Condens. Matter 22, 055602 (2010). [57] B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, Hubbard-corrected DFT energy functionals: The LDA + U description of correlated systems, Int. J. Quantum Chem. 114, 14 (2014). [58] Y.-C. Wang, Z.-H. Chen, and H. Jiang, The local projection in the density functional theory plus U approach: A critical assessment, J. Chem. Phys. 144, 144106 (2016). [59] M. Kick, K. Reuter, and H. Oberhofer, Intricacies of DFT + U , not only in a numeric atom centered orbital framework, J. Chem. Theory Comput. 15, 1705 (2019). [60] I. Timrov, F. Aquilante, L. Binci, M. Cococcioni, and N. Marzari, Pulay forces in density-functional theory with extended Hubbard functionals : From nonorthogonalized to orthogonalized manifolds, Phys. Rev. B 102, 235159 (2020). 235171-13 https://doi.org/10.1103/PhysRevMaterials.3.114602 https://doi.org/10.1103/PhysRevB.88.014114 https://doi.org/10.1103/PhysRevB.96.035143 https://doi.org/10.1103/PhysRevMaterials.6.064410 https://doi.org/10.1038/srep43482 https://doi.org/10.1103/PhysRevLett.125.217602 https://doi.org/10.1103/PhysRevResearch.4.L012042 https://doi.org/10.1103/PhysRevLett.100.136406 https://doi.org/10.1063/1.478522 https://doi.org/10.1063/1.1564060 https://doi.org/10.1063/1.2204597 https://doi.org/10.1103/PhysRevMaterials.6.104410 https://doi.org/10.1103/PhysRevB.44.943 https://doi.org/10.1088/0953-8984/9/4/002 https://doi.org/10.1103/PhysRevB.57.1505 https://doi.org/10.1103/PhysRevLett.97.103001 https://doi.org/10.1063/1.2987444 https://doi.org/10.1126/science.1158722 https://doi.org/10.1021/jp5030334 https://doi.org/10.1103/PhysRevB.90.035204 https://doi.org/10.1039/C6TA06611A https://doi.org/10.1038/s41598-017-06396-5 https://doi.org/10.1063/1.5022319 https://doi.org/10.1088/1361-648X/ab06d5 https://doi.org/10.1088/1361-648X/abaa81 https://doi.org/10.1016/j.physb.2021.413381 https://doi.org/10.1038/358136a0 https://doi.org/10.1080/00150190008216232 https://doi.org/10.1103/PhysRevLett.70.1010 https://doi.org/10.1088/0953-8984/22/5/055602 https://doi.org/10.1002/qua.24521 https://doi.org/10.1063/1.4945608 https://doi.org/10.1021/acs.jctc.8b01211 https://doi.org/10.1103/PhysRevB.102.235159 G. GEBREYESUS et al. PHYSICAL REVIEW B 108, 235171 (2023) [61] M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA + U method, Phys. Rev. B 71, 035105 (2005). [62] I. Timrov, N. Marzari, and M. Cococcioni, Hubbard parame- ters from density-functional perturbation theory, Phys. Rev. B 98, 085127 (2018). [63] I. Timrov, N. Marzari, and M. Cococcioni, Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formu- lations, Phys. Rev. B 103, 045141 (2021). [64] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos et al., Quantum ESPRESSO: A modular and open-source software project for quantum sim- ulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009). [65] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio Jr., A. Ferretti, A. Floris, G. Fratesi, G. Fugallo et al., Advanced ca- pabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017). [66] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum ESPRESSO toward the exascale, J. Chem. Phys. 152, 154105 (2020). [67] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci. 95, 337 (2014). [68] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Improved tetra- hedron method for Brillouin-zone integrations, Phys. Rev. B 49, 16223 (1994). [69] R. Fletcher, Practical Methods of Optimization, 2nd ed. (Wiley, Chichester, 1987). [70] I. Timrov, N. Marzari, and M. Cococcioni, HP—A code for the calculation of Hubbard parameters using density-functional perturbation theory, Comput. Phys. Commun. 279, 108455 (2022). [71] P.-O. Löwdin, On the non-orthogonality problem con- nected with the use of atomic wave functions in the theory of molecules and crystals, J. Chem. Phys. 18, 365 (1950). [72] I. Mayer, On Löwdin’s method of symmetric orthogonaliza- tion, Int. J. Quant. Chem. 90, 63 (2002). [73] I. Timrov, P. Agrawal, X. Zhang, S. Erat, R. Liu, A. Braun, M. Cococcioni, M. Calandra, N. Marzari, and D. Passerone, Electronic structure of Ni-substituted LaFeO3 from near edge x-ray absorption fine structure experiments and first-principles simulations, Phys. Rev. Res. 2, 033265 (2020). [74] R. Mahajan, I. Timrov, N. Marzari, and A. Kashyap, Im- portance of intersite Hubbard interactions in β-MnO2: A first-principles DFT + U + V study, Phys. Rev. Mater. 5, 104402 (2021). [75] I. Timrov, F. Aquilante, M. Cococcioni, and N. Marzari, Accurate electronic properties and intercalation voltages of olivine-type Li-ion cathode materials from extended Hubbard functionals, PRX Energy 1, 033003 (2022). [76] R. Mahajan, A. Kashyap, and I. Timrov, Pivotal role of in- tersite Hubbard interactions in Fe-doped α-MnO2, J. Phys. Chem. C 126, 14353 (2022). [77] I. Timrov, M. Kotiuga, and N. Marzari, Unraveling the effects of inter-site Hubbard interactions in spinel Li-ion cathode ma- terials, Phys. Chem. Chem. Phys. 25, 9061 (2023). [78] F. Haddadi, E. Linscott, I. Timrov, N. Marzari, and M. Gibertini, On-site andinter-site Hubbard corrections in magnetic monolayers: The case of FePS3 and CrI3, arXiv:2306.06286 [Phys. Rev. Mater. (to be published)]. [79] L. Binci, M. Kotiuga, I. Timrov, and N. Marzari, Hybridization driving distortions and multiferroicity in rare-earth nickelates, Phys. Rev. Res. 5, 033146 (2023). [80] P. Bonfà, I. Onuorah, F. Lang, I. Timrov, L. Monacelli, C. Wang, X. Sun, O. Petracic, G. Pizzi, N. Marzari, S. Blundell, and R. De Renzi, Magnetostriction-driven muon localisation in an antiferromagnetic oxide, arXiv:2305.12237 [Phys. Rev. Lett. (to be published)]. [81] C. Ricca, I. Timrov, M. Cococcioni, N. Marzari, and U. Aschauer, Self-consistent DFT + U + V study of oxy- gen vacancies in SrTiO3, Phys. Rev. Res. 2, 023313 (2020). [82] A. Togo and I. Tanaka, First principles phonon calculations in materials Science, Scr. Mater. 108, 1 (2015). [83] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density- functional perturbation Theory, Rev. Mod. Phys. 73, 515 (2001). [84] P. Umari and A. Pasquarello, Ab initio molecular dynamics in a finite homogeneous electric field, Phys. Rev. Lett. 89, 157602 (2002). [85] I. Souza, J. Iniguez, and D. Vanderbilt, First-principles ap- proach to insulators in finite electric fields, Phys. Rev. Lett. 89, 117602 (2002). [86] P. Umari and A. Pasquarello, Infrared and Raman spectra of disordered materials from first principles, Diam. Relat. Mater. 14, 1255 (2005). [87] P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Ab initio calculation of phonon dispersions in semiconductors, Phys. Rev. B 43, 7231 (1991). [88] X. Gonze and C. Lee, Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997). [89] J. Tóbik and A. Dal Corso, Electric fields with ultrasoft pseudo-potentials: Applications to benzene and anthracene, J. Chem. Phys. 120, 9934 (2004). [90] L. Bastonero and N. Marzari, Automated all-functionals in- frared and raman spectra, arXiv:2308.04308 (2023). [91] S. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. Yakutovich, C. Andersen, F. Ramirez, C. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari et al., AiiDA 1.0, a scalable computational in- frastructure for automated reproducible workflows and data provenance, Sci. Data 7, 300 (2020). [92] M. Uhrin, S. P. Huber, J. Yu, N. Marzari, and G. Pizzi, Work- flows in AiiDA: Engineering a high-throughput, event-based engine for robust and modular computational workflows, Comput. Mater. Sci. 187, 110086 (2021). 235171-14 https://doi.org/10.1103/PhysRevB.71.035105 https://doi.org/10.1103/PhysRevB.98.085127 https://doi.org/10.1103/PhysRevB.103.045141 https://doi.org/10.1088/0953-8984/21/39/395502 https://doi.org/10.1088/1361-648X/aa8f79 https://doi.org/10.1063/5.0005082 https://doi.org/10.1016/j.commatsci.2014.07.043 https://doi.org/10.1103/PhysRevB.49.16223 https://doi.org/10.1016/j.cpc.2022.108455 https://doi.org/10.1063/1.1747632 https://doi.org/10.1002/qua.981 https://doi.org/10.1103/PhysRevResearch.2.033265 https://doi.org/10.1103/PhysRevMaterials.5.104402 https://doi.org/10.1103/PRXEnergy.1.033003 https://doi.org/10.1021/acs.jpcc.2c04767 https://doi.org/10.1039/D3CP00419H https://arxiv.org/abs/2306.06286 https://doi.org/10.1103/PhysRevResearch.5.033146 https://arxiv.org/abs/2305.12237 https://doi.org/10.1103/PhysRevResearch.2.023313 https://doi.org/10.1016/j.scriptamat.2015.07.021 https://doi.org/10.1103/RevModPhys.73.515 https://doi.org/10.1103/PhysRevLett.89.157602 https://doi.org/10.1103/PhysRevLett.89.117602 https://doi.org/10.1016/j.diamond.2004.12.007 https://doi.org/10.1103/PhysRevB.43.7231 https://doi.org/10.1103/PhysRevB.55.10355 https://doi.org/10.1063/1.1729853 https://arxiv.org/abs/2308.04308 https://doi.org/10.1038/s41597-020-00638-4 https://doi.org/10.1016/j.commatsci.2020.110086 UNDERSTANDING THE ROLE OF HUBBARD CORRECTIONS … PHYSICAL REVIEW B 108, 235171 (2023) [93] S. Lee, Placzek-type polarizability tensors for Raman and resonance Raman scattering, J. Chem. Phys. 78, 723 (1983). [94] D. A. Tenne, X. X. Xi, Y. L. Li, L. Q. Chen, A. Soukiassian, M. H. Zhu, A. R. James, J. Lettieri, D. G. Schlom, W. Tian, and X. Q. Pan, Absence of low-temperature phase transitions in epitaxial BaTiO3 thin films, Phys. Rev. B 69, 174101 (2004). [95] O. Maslova, Y. Yuzyuk, M. Jain, and S. Barannikova, Lattice dynamics of barium titanate: Single crystal, ceramic, and poly- crystalline film, Phys. Status Solidi B 257, 1900762 (2020). [96] T. C. Damen, S. P. S. Porto, and B. Tell, Raman effect in zinc oxide, Phys. Rev. 142, 570 (1966). [97] P. Brüesch, Phonons: Theory and Experiments II (Springer, Berlin, 1986). [98] M. Popov, J. Spitaler, V. Veerapandiyan, E. Bousquet, J. Hlinka, and M. Deluca, Raman spectra of fine-grained materi- als from first principles, npj Comput. Mater. 6, 121 (2020). [99] G. Gebreyesus, L. Bastonero, M. Kotiuga, N. Marzari, and I. Timrov, Understanding the role of Hubbard corrections in the rhombohedral phase of BaTiO3, Materials Cloud Archive 2023.187 (2023), doi: 10.24435/materialscloud:vz-7q. [100] H.-Y. Zhang, Z.-Y. Zeng, Y.-Q. Zhao, Q. Lu, and Y. Cheng, First-principles study of lattice dynamics, structural phase transition, and thermodynamic properties of barium titanate, Z. Naturforsch. A 71, 759 (2016). [101] Y.-S. Seo and J. S. Ahn, First-principles investigations on polytypes of BaTiO3: Hybrid calculations and pressure depen- dences, J. Kor. Phys. Soc. 62, 1629 (2013). [102] G. H. Kwei, A. C. Lawson, S. J. L. Billinge, and S. W. Cheong, Structures of the ferroelectric phases of barium titanate, J. Phys. Chem. 97, 2368 (1993). [103] H. Megaw, Temperature changes in the crystal structure of barium titanium oxide, Proc. R. Soc. A: Math. Phys. Eng. Sci. 189, 261 (1947). [104] N. Kirchner-Hall, W. Zhao, Y. Xiong, I. Timrov, and I. Dabo, Extensive benchmarking of DFT + U calculations for predict- ing band gaps, Appl. Sci. 11, 2395 (2021). [105] N. L. Nguyen, N. Colonna, A. Ferretti, and N. Marzari, Koopmans-compliant spectral functionals for extended sys- tems, Phys. Rev. X 8, 021051 (2018). [106] S. H. Wemple, Polarization fluctuations and the optical- absorption edge in BaTiO3, Phys. Rev. B 2, 2679 (1970). [107] P. Erhart and K. Albe, Thermodynamics of mono- and di- vacancies in barium titanate, J. Appl. Phys. 102, 084111 (2007). [108] J. D. Axe, Apparent ionic charges and vibrational eigenmodes of BaTiO3 and other perovskites, Phys. Rev. 157, 429 (1967). [109] P. Ghosez, First-principles study of the dielectric and dynam- ical properties of Barium Titanate, Ph.D. thesis, Universite Catholique de Louvain, 1997. [110] T. Paoletta and A. A. Demkov, Pockels effect in low- temperature rhombohedral BaTiO3, Phys. Rev. B 103, 014303 (2021). [111] F. Eltes, G. E. Villarreal-Garcia, D. Caimi, H. Siegwart, A. A. Gentile, A. Hart, P. Stark, G. D. Marshall, M. G. Thompson, J. Barreto, J. Fompeyrine, and S. Abel, An integrated optical modulator operating at cryogenic temperatures, Nat. Mater. 19, 1164 (2020). [112] X. Wu, D. Vanderbilt, and D. R. Hamann, Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory, Phys. Rev. B 72, 035105 (2005). [113] B. Wang and C. Sun, Precise measurement of thermal-induced refractive-index change in BaTiO3 on the basis of anisotropic self-diffraction, Appl. Opt. 40, 672 (2001). [114] G. Burns and F. Dacol, Polarization in the cubic phase of BaTiO3, Solid State Commun. 42, 9 (1982). [115] A. Floris, S. de Gironcoli, E. K. U. Gross, and M. Cococcioni, Vibrational properties of MnO and NiO from DFT + U - based density functional perturbation theory, Phys. Rev. B 84, 161102(R) (2011). [116] A. Floris, I. Timrov, B. Himmetoglu, N. Marzari, S. de Gironcoli, and M. Cococcioni, Hubbard-corrected density functional perturbation theory with ultrasoft pseudopotentials, Phys. Rev. B 101, 064305 (2020). [117] W. Zhong, R. D. King-Smith, and D. Vanderbilt, Giant LO-TO splittings in perovskite ferroelectrics, Phys. Rev. Lett. 72, 3618 (1994). [118] Y. Xie, H.-T. Yu, G.-X. Zhang, and H.-G. Fu, Lattice dynamics investigation of different transition behaviors of cubic BaTiO3 and SrTiO3 by first-principles calculations, J. Phys.: Condens. Matter 20, 215215 (2008). [119] P. Ghosez, X. Gonze, and J.-P. Michenaud, Ab initio phonon dispersion curves and interatomic force constants of barium titanate, Ferroelectr. 206, 205 (1998). [120] The optimized lattice parameter for the cubic BTO at 0 K using PBEsol is found to be 3.979 Å. [121] V. Buscaglia, M. Buscaglia, M. Viviani, T. Ostapchuk, I. Gregora, J. Petzelt, L. Mitoseriu, P. Nanni, A. Testino, R. Calderone, C. Harnagea, Z. Zhao, and M. Nygren, Raman and AFM piezoresponse study of dense BaTiO3 nanocrystalline ceramics, J. Eur. Ceram. Soc. 25, 3059 (2005). 235171-15 https://doi.org/10.1063/1.444775 https://doi.org/10.1103/PhysRevB.69.174101 https://doi.org/10.1002/pssb.201900762 https://doi.org/10.1103/PhysRev.142.570 https://doi.org/10.1038/s41524-020-00395-3 https://doi.org/10.24435/materialscloud:vz-7q https://doi.org/10.1515/zna-2016-0149 https://doi.org/10.3938/jkps.62.1629 https://doi.org/10.1021/j100112a043 https://doi.org/10.1098/rspa.1947.0038 https://doi.org/10.3390/app11052395 https://doi.org/10.1103/PhysRevX.8.021051 https://doi.org/10.1103/PhysRevB.2.2679 https://doi.org/10.1063/1.2801011 https://doi.org/10.1103/PhysRev.157.429 https://doi.org/10.1103/PhysRevB.103.014303 https://doi.org/10.1038/s41563-020-0725-5 https://doi.org/10.1103/PhysRevB.72.035105 https://doi.org/10.1364/AO.40.000672 https://doi.org/10.1016/0038-1098(82)91018-3 https://doi.org/10.1103/PhysRevB.84.161102 https://doi.org/10.1103/PhysRevB.101.064305 https://doi.org/10.1103/PhysRevLett.72.3618 https://doi.org/10.1088/0953-8984/20/21/215215 https://doi.org/10.1080/00150199808009159 https://doi.org/10.1016/j.jeurceramsoc.2005.03.190