First-principles method justifying the Dieke diagram and beyond

We present a method to determine the model Hamiltonians to treat rare-earth multiplets in solids from the results of the quasiparticle self-consistent \textit{GW} (QSGW) method. We apply the method to trivalent Eu compounds EuCl$_3$, EuN, and Eu-doped GaN after examining free rare-earth ions. We solve the model Hamiltonian by the exact diagonalization. Our results justify applying the Dieke diagram to ions in solid, while its limitation is clarified. In particular, we show that the crystal fields cause sizable breaking of the Russell-Saunders coupling.


I. INTRODUCTION
The photoluminescence of rare-earth atoms (PLR) embedded in solids is one of the essential physical phenomena, which has a wide range of possible applications for light emitters. For example, Eu-doped GaN is a candidate for trichromatic LEDs, thanks to red-light emission by PLR of Eu [1][2][3]. It is favorable to assist the investigation of such PLR with the computations, that is, the computational materials design (CMD). However, CMD for PLR is not so easy because we do not have a reliable theoretical basis to perform computations for PLR even today. Historically, the model-Hamiltonian methods based on the atomic multiplet theory had been developed to explain complex spectra of PLR, where the model Hamiltonian was parametrized by the Slater integrals [4], which are treated as material-dependent parameters. The atomic multiplet theory to determine the model Hamiltonian was given by Racah [5], followed by further developments [6][7][8][9][10]. Particularly, a comprehensive study by Dieke and Crosswhite on doubly and triply ionized rare earth (RE) is well known. The energy level diagram for all lanthanides in the study called the Dieke diagram is taken as a standard for analyzing experimental results even now [11]. The material-dependent parameters to give the Dieke diagram are determined experimentally so as to reproduce the optical spectra of PLR. Since we had no means to determine the parameters by computation, the model-Hamiltonian methods were quite limited from the view of CMD.
One of the recent theoretical approaches to treating the multiplets was given by Ungur and Chibotaru, who successfully applied the complete active space self-consistent-field method (CASSCF) [12] to lanthanide complexes [13,14]. However, applying CASSCF to handle RE in solids such as semiconductors should be difficult. This is because we have to consider the screening effects of the Coulomb interaction between 4 electrons caused by the polarization of host semiconductors. Furthermore, 4 orbitals can hybridize well with atomic orbitals in solids: the 4 orbitals can be itinerant or localized depending on their environments [15][16][17][18]. We have to use a method reproducing the properties of both RE and solids on the same footing. In CASSCF, it is hopeless from the view of high computational demand to use sufficiently large enough active space to reproduce the properties of solids. CASSCF is applicable to only the cases where 4 electrons are very localized.
We have other first-principles-based methods to handle the multiplets until now. Most of these methods have two key steps. One is how to determine the model Hamiltonian, and the other is how to solve the model, namely the solver. For the solver, we have extensive developments until now. Assuming the LDA+ model Hamiltonian, the Hubbard-I approximation [19,20] (HIA) and the dynamical mean-field theory (DMFT) have been developed to predict multiplet properties and peak structures of the excitation spectra [21][22][23][24].
A serious problem is in the first step. The reliability of widely used LDA+ is quite questionable. We usually have to give the size of by hand. LDA+ assumes a very simple double-counting term. There are no parameters to control the center of 4 bands relative to the anion bands, that is, we simply assume that the center is determined in LDA. Considering these facts, we guess that LDA+ is often abused with little theoretical justifications as was analyzed by Lee, Kotani, and Ke [25]. Moreover, the multiplet excitation of RE is not controlled by but by the 2nd Slater integral 2 as we will explain later on.
In this paper, we will give a new method for the first step, deducing the parameters in the model Hamiltonian from the results of the quasiparticle self-consistent GW (QSGW) method [26]. QSGW is a reliable mean-field approximation in the sense that the one-particle Hamiltonian determined in QSGW gives a good independent-particle picture for a wide range of materials [26,27] including not only semiconductors but also 4 systems [28]. QSGW is roughly identified to be a "screened" Hartree-Fock method where the screening effect is internally determined self-consistently. Both excitation energies and quantities such as spin fluctuations can be reproduced well based on the one-particle Hamiltonian in QSGW [29]. Let us explain our core idea about how to deduce the parameters in the model Hamiltonian. Our core idea is that we require "QSGW applied to the model Hamiltonian" should reproduce the (model part of) one-body Hamiltonian given in QSGW. In contrast to LDA, we have no theoretical problem applying QSGW to the model Hamiltonian.
In this paper, we apply our new method to the multiplets of 4 orbitals in RE compounds, EuCl 3 , EuN, and Eu-doped GaN, after examining trivalent RE ions (RE 3+ ) in the supercell. The parameters of the model Hamiltonian for describing 4 electrons are obtained based on our core idea. With the exact diagonalization applied to the model Hamiltonian, we show the eigenvalues of multiplets of RE 4 orbitals and discuss the relation with experiments.

II. METHOD
We assume a model Hamiltonian of the multiplets with a fixed number of electrons of 4 orbitals. Here we neglect the hybridization of 4 orbitals with the other orbitals. Thus, the dimension of the model Hamiltonian is the number of the electronic configurations of 4 electrons, i.e., 14 where is the number of 4 electrons. Then, the model Hamiltonian is written as H 0 is a constant matrix to give the base level of 4 . The level is irrelevant to energy spectra since we do not consider the hybridizations. The non-spherical part of the one-body potential is given by the crystal field term H CF . In addition, we have the spin-orbit coupling (SOC) term H SOC , and the effective Coulomb interaction term H C . Surrounding atoms of the RE atom affect not only H CF but also H C through the size of interaction, and H SOC as well via the shape of 4 orbitals. H SOC and H C are given as Here indices , , ( = 1, 2, 3, 4) are for the magnetic quantum number, and for spins, andˆis the electronannihilation operator. H SOC is made of the strength of SOC and the angular-momentum and spin matrices , , · · · . We assume the effective Coulomb interactions 1 2 3 4 are given as where we have the Gaunt coefficients ( , ) [30] and the Slater-Condon parameters 2 for 4 orbitals [4,31]. We use the scaled-Slater-Condon parameters 0 , 0 and 2 to represent 2 as =7361.64 × 0.151 2 (5) in the manner of Ref. [32] for analyzing the RE 3+ elements. Here we fix the ratios 4 / 2 and 6 / 2 as in Ref. [33] respecting the case of Hydrogen. The spin-dependent term 0 is introduced to make a compromise in our fitting procedure to the QSGW results with keeping the spinspace symmetry of the model Hamiltonian.
We assume a general crystal field that corresponds to the symmetry of each material. Using Steven's operators , and so on [34], H CF is given as where 0 4 , 0 6 , 0 2 , 6 6 are the parameters to specify crystal fields. Thus, the model Hamiltonian H of Eq. (1) is specified by parameters , 0 , 0 , 2 , and . To apply our core idea shown in the introduction, we should apply QSGW to the model Hamiltonian. Here we neglect the correlation part since we expect little screening effects of 4 orbitals by themselves. That is, we apply the Hartree-Fock approximation to H instead, resulting in the Hartree-Fock model Hamiltonian (HFMH) H HF as where H HF C is the mean-field approximation to H C . denotes the opposite spin to , and . . . means the expectation values for the ground state. Based on our core idea, we compare H HF with the 4 part of the one-body Hamiltonian H QSGW 4 given in QSGW, so as to determine parameters in H .
We use the ecalj package[35] to perform QSGW calculations. In practice, we use 20% LDA mixing (QSGW80) so as to reduce too-large exchange effects. We can take QSGW80 as a quick remedy to include the effects of the vertex correction, which enhances the screening effect by ∼ 20% [36,37]. In fact, QSGW80 can reproduce the experimental band gap very well [27,38]. We obtain H QSGW 4 based on localized Wannier functions (it was implemented in ecalj.) as |H QSGW 4 | . Here and are the Wannier functions generated by the one-shot projection of 4 -type atomic-like seed orbitals [35,39].
In the QSGW calculations, we add the SOC term H QSGW SOC in the -only approximation. We obtain by an average Other parameters 0 , 0 , 2 , and of H are determined to minimize the difference of eigenvalues between H QSGW 4 and H HF . Then, we can obtain eigenvalues of H by the exact diagonalization [40].
In a summary, we mainly added two approximations to our core idea. One is H HF instead of applying QSGW to the model Hamiltonian. The other is the QSGW80 instead of QSGW. In addition, we fix the ratios 4 / 2 and 6 / 2 . We utilize the fixed ratio to enhance numerical stability and simple interpretation. In principle, our core idea is rather general for extracting an essential degree of freedom from the firstprinciples calculations. One of the advantages of our core idea is that we do not calculate effective interaction directly. The calculation of effective interaction is somehow complicated as discussed in Refs. 36 and 37, especially when we like to include vertex corrections.

III. RESULTS
Prior to discussion of the RE in compounds, we have examined free +3 in QSGW to confirm the performance of our method. For the calculation of the free +3 , we use a fcc supercell placing +3 at their centers, where +3 are separated by 7.07 Å. We assume a homogeneous background charge to keep charge neutrality.
In Figs. 1 (a) and 1 (b), we show the calculated band structure of Eu 3+ in QSGW. We superpose that of H QSGW 4 with green flat lines (seven lines per spin). PDOS of 4 orbitals is shown in Figs. 1 (c) and (d). Up to ∼ 18 eV, the calculated bands corresponding to 4 , 5 , and 6 are almost flat, indicating that our supercell is large enough. On the other hand, 6 bands around ∼ 20 eV have band width ∼ 1 eV because of the finite size of our supercell. The eigenvalue at Γ point around ∼ 21.5 eV is identified to be the vacuum level, the bottom of the scattering states. The band gap (LUMO-HOMO gap) is 13.7 eV, much larger than the LDA value of 4.85 eV. The majority bands of 4 orbitals except for = 3 are occupied, that is, the ground state is = + = −3 + 3 = 0, corresponding to 7 F 0 . Thus, our results are consistent with Hund's rule. Generally speaking, a good mean-field approach should satisfy Hund's rule since the ground states of atoms should be essentially described by the electronic configuration of a single Slater determinant.
We determine 0 , 0 , and 2 in H HF so as to reproduce eigenvalues of H QSGW 4 . We see the eigenvalues match well as shown in Fig. 1 (h). Figures. 1 (e)-(h) show our analysis of how the parameters affect the eigenvalues of H HF . Figure 1 (e) shows splits into seven states with degeneracy between and − . Figure (f) shows that 0 makes the difference between occupied and unoccupied states. Figure (g) shows that 0 with 2 is still not enough to reproduce H QSGW 4 . The parameter 0 is for describing larger effective interaction between occupied orbitals than that between occupied and unoccupied orbitals. This is reasonable because occupied 4 orbitals are localized more than unoccupied orbitals.
In Table I, we summarize obtained , 0 , 0 , and 2 for RE +3 in our method. The QSGW results are in the Supplemental Material [41]. We skipped Ce since QSGW did not give localized 4 eigenfunctions. Table I shows that and 2 in our method give good agreements with those of empirical studies: their differences are only ∼10%. This is a justification for our method. 0 , corresponding to of LDA+ , is changing systematically along the atomic number. Except for Gd, 0 is the largest in the middle of the first half and the latter half of the 4 series, namely at Sm when filling majority electrons, and at Er when minority. The latter half of the species shows a little larger values than the first half. 2 systematically becomes larger along the atomic number. For example, the 2 values of Tm and Pr are 0.080 eV and 0.043 eV, and the difference is 0.037 eV in our method, while the empirical study shows these 2 values are 0.056 eV and 0.038 eV, and the difference is as 0.018 eV.
In Fig. 2, we compare the excitation energies of free RE 3+ by the parameters of our method and by those of empirical studies in Table I. The latter exactly reproduced the original Dieke diagram [11]. In free RE 3+ , the excitation energies only depend on and 2 . Although the lowest excitation energies show good agreements (except for Tm with a little large error of ∼30%), we see disagreements overall. This is because the excitation energies are somehow sensitive to the small differences of and 2 in Table I.   TABLE I. Parameters in eV determined in our method for free trivalent ions. The empirical studies in the right column correspond to the empirical values of Dieke's review paper [11]. These values are derived from many experimental and theoretical data [6,7,[42][43][44][45][46][47][48][49].
Our  [50,51]. We show the QSGW results for EuCl 3 and EuN in Supplemental Material [41]. In Table II, we show , 0 , 0 , 2 , and . We see 0 and 0 are strongly reduced from those of the free Eu ion. However, 2 is not so different from that of free Eu ion. This finding suggests that 2 given by Dieke is reasonable probably even in other solids. However, we see some complicated behavior of 2 in detail: 2 =0.058 eV for the free Eu ion is reduced to be 2 =0.053 eV for EuN while a little enhanced to be 2 =0.075 eV for EuCl 3 . On the other hand, =0.195 eV   Table I for free Eu ion is very different from =0.093 eV for Eu-doped GaN. Considering the expression of SOC, this is due to the delocalization of 4 orbitals in Eu-doped GaN, suggesting hy-bridization with other surrounding orbitals. The parameters for CF should be material dependent. In particular, we see that 0 2 for the -axis anisotropy is rather large (=2.217 eV) for EuCl 3 , while 0 2 =0.346 eV for Eu-doped GaN. Let's focus on Eu-doped GaN. As shown in Figs. 3 (a) and (b), QSGW80 gives the band gap 3.78 (3.69) eV for the majority (minority) spin, corresponding to the experimental value of GaN 3.4 eV. Just above the bottom of the conduction band, we have an unoccupied band of = ±3 for the majority spin. In contrast to the unoccupied band, the occupied 4 orbitals are hybridized well with valence bands of GaN. Such hybridization is also seen in the calculations of the HSE functional [54]. In Fig. 3 (c), we compare the eigenvalues of H HF (blue lines) and those of H QSGW 4 (red lines). A little matching error indicates room for improving our method, while the error may not change our conclusions.
We show the eigenvalues of H in Fig. 3 (g) obtained by the exact diagonalization. Here we classify the eigenvalues by the expectation values of total angular momentum . In Fig. 3 (f), we show eigenvalues when we neglect H CF in H . Without H CF , each eigenstate is represented by Russell-Saunders states. In Fig 3 (f) without CF, there is a large gap between 7 F and 5 D. The size of this gap is ∼ 2 eV. This corresponds to the observed red emission which is experimentally identified to be the transition from 5 D 0 to 7 F 2 [2]. Comparison of Figs. 3 (f) and (g) shows that H CF is large enough to mix low-and higheigenstates, resulting in a middle size of . That is, H CF causes a sizable breaking of the Russell-Saunders coupling. However, H CF is not large enough to alter the overall structure. That is, we see remnants of 5 D and 7 F in Fig. 3 (g) while the excitation gap between them is preserved. This justifies the use of the Dieke diagram for the analysis of PLR.
In summary, we presented a new method to determine the model Hamiltonian from the QSGW calculations. By the exact diagonalization of the model Hamiltonian, we can justify the applicability of the Dieke diagram to RE in solids and its limitations. Along the line of our method, it is possible to include the hybridization of 4 electrons with others. In fact, our QSGW calculation shows we have large hybridization of occupied 4 orbitals with valence bands. We now have an extension of our method applied to multiplets of 3 orbitals working well [55].