Quantum chemical insights into hexaboride electronic structures: correlations within the boron $p$-orbital subsystem

The notion of strong electronic correlations arose in the context of $d$-metal oxides such as NiO but can be exemplified on systems as simple as the H$_2$ molecule. Here we shed light on correlation effects on B$_6^{2-}$ clusters as found in $M$B$_6$ hexaborides and show that the B 2$p$ valence electrons are fairly correlated. B$_6$-octahedron excitation energies computed for CaB$_6$ and YbB$_6$ agree with peak positions found by resonant inelastic x-ray scattering, providing a compelling picture for the latter. Our findings characterize these materials as very peculiar $p$-electron correlated systems and call for more involved many-body investigations within the whole hexaboride family, both alkaline- and rare-earth compounds, not only for $N$- but also ($N\pm1$)-states defining e. g. band gaps.


Introduction
The physical properties of M B 6 hexaborides are remarkably diverse and have received persistent interest [1,2], as a thorough understanding of this peculiar family of materials could not be achieved yet. Some of its members possess properties relevant for applications, generally determined by specificities of their electronic structures: CaB 6 is a semiconductor [3], YB 6 turns superconducting for temperatures lower than 8 K [4], LaB 6 is widely used as thermionic electron emitter [5], EuB 6 displays colossal magneto-resistance [6], the heavy-fermion system Ce 1−x La x B 6 has been intensively studied in the context of multipolar phases and magnetically hidden orders [7,8], while SmB 6 is believed to host correlated topological states [9]. However, even for the seemingly simplest member of the family, CaB 6 , there are important electronic-structure features that are not at all clear, for example, the size and nature of its band gap: densityfunctional computations yield a vanishing band gap [10][11][12][13][14] while experiment indicates a gap of ∼1 eV [3]. The band structures of SmB 6 and YbB 6 are also a topic of active debate, especially in relation to their surface states [9,[15][16][17]; for SmB 6 , even predictions made for the symmetry of the Sm 3+ -ion ground-state term are recently contradicted by experiment [18,19].
In this context, we argue that one aspect constantly neglected in the study of these compounds is correlations within the B 2p electronic subsystem. To address such physics, we here employ ab initio multi-configurational wavefunction theory. For CaB 6 and YbB 6 , representative divalent hexaborides, the ground-state wavefunctions turn out strongly multi-configurational, with large admixture of higher-lying configurations. Trying to understand the nature of B 2p N -particle states, a rich Nparticle excitation spectrum is evidenced in the ab initio calculations. These computed p-p excitation energies compare well with peak positions in B K-edge RIXS spectra reported by Denlinger et al. [3,20]. Analysis of the many-body wavefunctions in terms of site-centered valence orbitals suggests that the ab initio data can be mapped onto a six-site, half-filled effective Hubbard model. Unfolded to the extended lattice, this might provide a convenient frame for further insights into the electronic properties of the B 2p subsystem. The computationally more expensive option is using quantum chemical methods to determine the nature of (N±1) quasi-particle states [21][22][23] defining e. g. band gaps.

Correlated electronic structure of CaB6
The crystal structure of cubic M B 6 hexaborides is CsCl-like, with P m3m symmetry [1,2]. A given B 6 octahedron is surrounded by eight M nearest-neighbor sites defining a cube (see Fig. 1). M can be an alkalineearth atom, Y, or a rare-earth ion (e.g., La, Ce, Sm, Yb). Interestingly, even for CaB 6 , a conventional closedshell p-electron compound at first sight, standard band-  [10][11][12][13][14] while a gap of ∼1 eV is found experimentally [3]. One obvious question is then what kind of correlations are missed in (LDA-based) density functional theory in this class of materials but the answer is not evident. From post-LDA GW calculations, both finite [24] and vanishing [13,14] band gaps were reported. These contradictory sets of data suggest that the results are sensitive to specific details in the numerical implementation. Moreover, the GW scheme is not designed for dealing with shortrange correlations; the latter turn out to be massive, as discussed below.
For an isolated B atom, the Aufbau principle implies an [1s 2 ]2s 2 2p 1 valence configuration; on a B 2− 6 octahedral unit, there are 6 × 3 + 2 = 20 electrons associated with the B 2s and 2p shells. But those 20 electrons can be now regarded as occupying symmetry-adapted orbitals distributed over several B sites of a given octahedron. In cubic geometry (i. e., regular octahedron), there are a 1g , e g , and t 1u symmetry-adapted functions associated with the six s atomic orbitals and a 1g , e g , t 1u , and t 2g linear combinations related to the set of 18 p atomic orbitals.
To gain insight into the B 2s/2p electronic structure, we here rely on ab initio quantum chemical methods [25] as implemented in the Orca program package v5.0 [26]. The material model is a point-charge embedded B 6 2cluster, for both CaB 6 and YbB 6 (see Fig. 1). Details on the technicalities of this embedding process are given in the Methods section. The quantum chemical investigation was initiated as a preliminary Hartree-Fock (HF) calculation that provides starting orbitals for a much more sophisticated post-HF treatment: complete-active-space self-consistent-field (CASSCF) computations and subsequent multi-reference configuration interaction (MRCI) with single and double excitations [25]. To make the highly-demanding MRCI calculations feasible, we enabled the use of symmetry, according to the D 2h subgroup; B 1s orbitals were not correlated.
Starting from the closed-shell HF solution, basic information on the nature of the low-lying excited states and on the nature of most important correlation effects was gathered from a series of CASSCF calculations in which different numbers of orbitals were considered as active. In quantum chemical terminology, "active" is the set of orbitals within which all possible determinants are generated on the basis of a given number of electrons. These are chemically or physically relevant orbitals, to e.g. bonding or spectroscopy, respectively; less relevant orbitals are assigned either double (core/semi-core levels) or zero (virtual levels) occupancy. By analyzing different active orbital spaces, the six molecular orbitals shown in Fig. 2 the A 1g , E g , and T 1u irreducible representations. We found that for the aspects discussed here a six-orbital active space (a 1g + e g + t 1u ) accommodating six electrons (for lower-energy orbitals on the B 6 2octahedron, double occupancy is imposed in CASSCF) already provides a balanced description; we refer to this active space as CAS (6,6).
Using this active space, a total of one septet, five quintets, 15 triplets and eight singlet states were considered in a state-averaged [25] CASSCF calculation. Results computed this way for CaB 6 are listed in Table I. Electronic configurations are here expressed in terms of symmetry-adapted molecular-like orbitals (a 1g , e g , t 1u ); this is common practice when discussing the electronic structure of octahedral atomic units and also the usual way output numerical data are printed by quantum chemical programs. A first observation is the pronounced multi-configurational structure of the computed wavefunctions. The leading ground-state configuration, for example, has a weight of only 61% in the many-body ground-state wavefunction ( 1 A 1g state) and a similar situation is observed for many of the excited states; the remaining weight has to do in each case with other electron configurations. The multi-configurational character of most of the states in Table I is also reflected in the natural-orbital occupation numbers, obtained through diagonalization of the density matrix [25]: for the 1 A 1g ground state, for instance, the a 1g , e g , and t 1u occupation numbers are 1.85, 1.70 (×2), and 0.25 (×3), respectively; at the HF level the occupation of the antibonding (see Fig. 2(a) ) t 1u levels is 0 since the HF ground-state wavefunction is 100% a 2 1g e 4 g . Also illustrative for the effect of electronic correlations: the (restricted open-shell) HF 1 A 1g -5 T 2g splitting, for instance, is less than 1 eV, much smaller than the CASSCF/MRCI values provided in Table I.
An instructive exercise is re-expressing the CAS(6,6) ground-state wavefunction in terms of site-centered orbitals. To achieve this, the symmetry-adapted a 1g , e g , and t 1u composites were transformed to a set of six symmetry-equivalent functions using the orbital "localization" module in Orca [26]. Two of these symmetry-equivalent site-centered orbitals are displayed in Fig. 2(b). Their asymmetric shape clearly indicates p-s mixing: dominant p but also s character, with substantial weight at the opposite B site. In this basis, the ground-state wavefunction amounts to 54% |↑↓↑↓↑↓ , i. e., 54% determinants with single orbital occupation, and 36% |↑↓↑↓ 0 2 + ... , 7% |↑↓ 0 2 0 2 + ... , and 3% |0 2 0 2 0 2 + ... configurations, where ↑, ↓, 2, and 0 stand for spin-up electron, spin-down electron, zero orbital occupation, and double occupation, respectively. As such, the wavefunction can in principle be mapped onto an effective six-site Hubbard-type model (t-U or t-U -V ). Estimating hopping (t) and Coulomb repulsion (U , V ) integrals [27][28][29] defining such effective models is beyond the scope of this letter. But the quantum chemical wavefunctions and excitation energies reported here provide useful ab initio benchmarks for tuning the relevant parameters. Remarkably, p-electron Hubbard-U values as large as 9 eV were estimated for e. g. carbon and phosphorus allotropes [30,31]. Good insight into correlation effects on the B 2− 6 unit is attained by inspecting the wavefunction describing two opposite B sites (the hopping integral is larger for opposite-corner sp hybrids). To this end the following numerical test can be made: allow all possible occupations for two sp hybrids on opposite corners of the octahedron (c.f. those pictured in Fig. 2(b)), but restrict the occupation of each of the other four sp hybrids to 1 (without re-optimizing the orbitals). The weight of con-figurations with 0, 2 and 2, 0 orbital occupation for the opposite-corner orbitals is only 19%, indicative of sizable correlations. For comparison, in the H 2 molecule at equilibrium H-H interatomic distance the weight of those ionic configurations in the HF wavefunction (uncorrelated limit) is 50%, being reduced to about 44% in the exact wavefunction [32]. Support to our computational results is provided by the good agreement between calculated N -particle excitation energies and excitation energies extracted from B K-edge RIXS experiments [3]. This is apparent by comparing numerical (MRCI level of theory) and experimental values in the last two columns of Table I and Fig. 3(a). In other words, the numerically obtained excitation energies fit the measured excited-state relative energies; this good agreement allows for an one-to-one assignment of the main features for less than ≈4 eV energy loss in the RIXS spectra. But the interpretation of the electronic ground state itself and of the lower N -particle charge excitations resulting from this assignment is peculiar, with no correspondence in the extensive literature available on M B 6 hexaborides.
That computations beyond the CAS(6,6) level are necessary is visible from the large corrections brought through the MRCI treatment to some of the excitation energies, in particular, to the relative energies of singlet and triplet states at 2.4-3.5 eV and of the 7 A 1u septet. These corrections amount to up to 1.5 eV, see Table I, and point to configuration-interaction (CI) effects that are not yet captured with a six-orbital activespace CASSCF. This is evident, e. g., from the substantial reduction of the weight of the leading configuration in the CI wavefunction of each of those states. Moreover, we found that more approximate post-CASSCF CI schemes such as spectroscopy-oriented CI (SORCI) or second-order perturbation theory (PT2) treatments such as N -electron valence PT2 correction (NEVPT2) yield inaccurate excitation energies (see Supplementary Table  1). Relative shifts of this magnitude when improving the description of electronic correlations (by appropriately enlarging the active orbital space in CASSCF and/or post-CASSCF correlation treatment) are typically found in d-electron oxides, e.g., for d 6 states in LaCoO 3 [33].

Correlated electronic structure of YbB6
For comparison, similar quantum chemical computations were carried out for the N -particle B 2− 6 electronic structure in the rare-earth hexaboride YbB 6 . The relative energies of the first twelve excited states and leading electronic configurations are listed in Table II. In agreement with the computational results obtained for CaB 6 , the ground-state wavefunction is of strong multiconfigurational character, with its leading configuration amounting to only 60%. Also in this case, the calculated MRCI energies can be compared with B K-edge X-ray spectra [20]; MRCI N -particle excitation energies (a) CaB 6 (b) YbB 6 are overlaid for this purpose onto the experimental plots of Denlinger et al. [20] in Fig. 3. It is seen that the relative energy of the well-defined first peak is rather well reproduced in both compounds. At higher energies, the level structure is denser; qualitatively, there is reasonably good correspondence between the peak found experimentally in each of the two materials at 2-3 eV and the excitations computed in that energy range. The B p-p excitation energies computed for YbB 6 are slightly lower as compared to CaB 6 . The likely explanation is that the larger effective charge of Yb ions in YbB 6 (2 or slightly larger than 2 according to experimental estimates [17,[34][35][36], a representation that is also adopted for the surroundings of the B 6 octahedron in our calculations) stabilizes the t 1u B-cluster orbitals as compared to the e g (also a 1g ) components and reduces this way the e g → t 1u excitation energies -the t 1u composites have stronger antibonding character, are more extended, and therefore sense more effectively the nearest-neighbor cations.

Conclusions
While strong correlation effects within a whole manifold of electron configurations were already pointed out for the case of d-ion clusters [37][38][39][40][41], we here describe similar physics occurring on p-electron octahedral clusters as found in M B 6 hexaborides. In particular, we highlight the role of correlation effects on B 6 octahedral units. For divalent CaB 6 and YbB 6 , our quantum-chemically computed valence-shell excitation energies agree well with peak positions in earlier B K-edge RIXS measurements  [3,20], providing an interpretation of the RIXS spectrum. Shifts to lower energies of the N -particle excited states in YbB 6 are ascribed to a relative stabilization of the antibonding t 1u cluster orbitals as a result of larger cation charges in the Yb compound. Given (i) the limitations of density functional theory in describing the electronic structure of hexaboride compounds, documented in the case of CaB 6 in literature extending over a few decades [3,[10][11][12][13]24], and (ii) more recent controversies on the electronic structures of e. g. SmB 6 [9,18,19] and YbB 6 [15][16][17], an extension of the present investigation of Nparticle excitations to quasi-particle band structures is highly desirable. The foundations of quasi-particle bandstructure calculations with ab initio wavefunction-based methods were laid [22,23,42] and such an approach provides better control over the underlying approximations as compared to main stream computational schemes based on density functional theory. The failures of the latter in hexaborides appear to be related to the correlated nature of the B 2− 6 electronic wavefunctions.

Methods
As material model, a B 6 2cluster was embedded in a field of M (with M = Ca, Yb) and B point charges (see Fig. 1 and Supplementary Table 2 for the cluster coordinates). The point charge field was created using the Ewald program [43,44]. For the lattice parameters, experimental data were used: lattice vectors a = 4.1537Å and an x-axis B atom position x B = 0.201 for CaB 6 [45][46][47] and parameters of a = 4.1792Å and x B = 0.202 for YbB 6 [17,48]. To account for covalency effects, an initial charge value of +1.8 e was chosen for Ca ions in the embedding, which corresponds to an ionic charge of −0.3 e for B. In contrast, in YbB 6 initial effective charges of +2.0 e for Yb [17,[34][35][36] and −0.33 e for B were assumed. For both embeddings, these charges were slightly adjusted by Ewald to ensure convergence of the Madelung sum. The final obtained point charge values do not deviate more than 0.01 e from the initial values. Charge neutrality was ensured in the overall system by slightly adjusting a set of about 150 B point charges at the boundary of the whole array of approximately 40000 point charges. To eight M nearest neighbors and to six adjacent B 6 octahedra around the reference B 6 2quantum cluster (see Fig. 1) we assigned capped effective core potentials (cECPs); the pseudopotentials of Kaupp et al. [49], Dolg et al. [50] and Bergner et al. [51] were employed for Ca, Yb and B, respectively. Their main purpose is to prevent electron density from the quantum cluster being artificially polarized by nearby point charges. Consequently, for the quantum cluster, an allelectron correlation-consistent polarized valence quadruple zeta (cc-pVQZ) basis set [52] was used in conjunction with the CASSCF method described above.

Data Availability Statement
Supplementary data (excitation energies of CaB 6 using further methods, cluster and point charge geometries) is attached and Orca calculation outputs are deposited at the ioChem-BD repository under the DOI 10.19061/iochem-bd-6-146. The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.