Combining Localized Orbital Scaling Correction and Bethe-Salpeter Equation for Accurate Excitation Energies

We applied localized orbital scaling correction (LOSC) in Bethe-Salpeter equation (BSE) to predict accurate excitation energies for molecules. LOSC systematically eliminates the delocalization error in the density functional approximation and is capable of approximating quasiparticle (QP) energies with accuracy similar or better than the $GW$ Green's function approach and with much less computational cost. The QP energies from LOSC instead of commonly used $G_{0}W_{0}$ and ev$GW$ are directly used in BSE. We show that the BSE/LOSC approach greatly outperforms the commonly used BSE/$G_{0}W_{0}$ approach for predicting excitations with different characters. For the calculations for Truhlar-Gagliardi test set containing valence, charge transfer (CT) and Rydberg excitations, BSE/LOSC with the Tamm-Dancoff approximation provides a comparable accuracy to time-dependent density functional theory (TDDFT) and BSE/ev$GW$. For the calculations of Stein CT test set and Rydberg excitations of atoms, BSE/LOSC considerably outperforms both BSE/$G_{0}W_{0}$ and TDDFT approaches with a reduced starting point dependence. BSE/LOSC is thus a promising and efficient approach to calculate excitation energies for molecular systems.


I. INTRODUCTION
Electronic excitation energy is one of the most important quantities for describing the electronic excited states. It can be compared with the 0-0 energy, which is measured by the optical spectroscopy. Computationally, determining the excitation energy from the electronic structure theory plays a critical role for obtaining insights about various phenomena and processes in chemistry, biochemistry and material science, such as molecular drug delivery 1,2 and solar cells [3][4][5] . In past decades, many efforts have been devoted to develop accurate and efficient theoretical approaches to compute excitation energies. One of the most popular approaches is time-dependent density functional theory (TDDFT) [6][7][8] . TDDFT has been widely implemented in modern quantum chemistry packages to calculate energies, structures and other properties of excited states for molecular and periodic systems [9][10][11][12] . The success of TDDFT stems from the good compromise between the accuracy and the computational cost. However, TDDFT still suffers from several problems.
For example, it is known that TDDFT with commonly used density functional approximations (DFAs) fails to describe Rydberg excitations and charge transfer (CT) excitations 12,13 . This issue can be attributed to the incorrect description for the long-range behavior for the potential energy surface 14,15 . Efforts for correcting the long-range behavior includes using range-separated or Coulomb-attenuated functionals [16][17][18] and mixing the Hartree-Fock (HF) 19,20 exchange in DFAs 21 .
Besides the failures for describing CT and Rydberg excitations, TDDFT has an undesired dependence on the exchange-correlation (XC) kernel stemming from different choices of DFAs. The difference originating from using different DFAs can exceed 1.0 eV for valence excitation energies and even exceed 2.0 eV for Rydberg excitation energies 12,13 .
Recently, Bethe-Salpeter equation (BSE) [22][23][24][25][26][27] in the Green's function many-body perturbation theory 26,28 has gained increasing attention to compute the optical spectroscopy for molecular systems. The BSE approach commonly takes the energy levels computed from the GW approximation as the input 26,[28][29][30] and this approach is denoted as the BSE/GW approach. In the BSE/GW approach, the screened interaction is used instead of the bare Coulomb interaction to describe the electron-hole interaction. The screened interaction is formulated with the quasiparticle (QP) energies from the GW calculation. It is known that the QP energies from GW are more accurate for predicting HOMO-LUMO gap (fundamental gap) than the conventional Kohn-Sham (KS) orbital energies. Besides the improved gap prediction, energies also have more clear physical meanings that they are interpreted as the charged excitation energies, or excitation energies for electron removal and addition. It has been shown that the GW approximation substantially improves the accuracy of predicting energy levels over the KS density functional theory (KS-DFT) approach for both occupied and unoccupied states 26,31,32 , which are the key quantities to calculate the excitation energies. Because the correct long-range behavior in the BSE/GW approach and the importance of the dynamical screening in real systems, BSE/GW has been applied to calculate excitation energies for systems of different sizes [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] . However, the BSE/GW approach still has several challenges. First, although BSE has the same scaling as TDDFT 50,51 , which is O(N 4 ) (N is the size of the system), the preceding GW calculation is computationally expensive.
In the fully analytical treatment of GW , the scaling of solving the random phase approximation (RPA) equation is O(N 6 ) and the scaling of evaluating the self-energy is O(N 5 ) 29,52 . Thus, the computational-demanding GW calculation is the bottleneck of the BSE/GW approach. To reduce the cost of GW calculations, different techniques can be used. For example, the cost of formulating the response function can be reduced to O(N 4 ) by using the plasmon-pole models 53 or the analytic continuation 29,31,54 . The cost of evaluating the self-energy can be reduced to O(N 5 ) in the contour deformation approach 55 and O(N 4 ) in the analytical continuation approach 56 . Recently the cubic scaling implementations 57,58 of GW calculations have also gained increasing attention. Second, the performance the most used BSE/G 0 W 0 approach strongly depends on the choice of the DFA.
Because of the perturbative nature of the one-shot G 0 W 0 method, the accuracy of G 0 W 0 strongly depends on the starting point [59][60][61] . This undesired dependence is inherited in the BSE/G 0 W 0 approach. It has been shown that the accuracy of BSE/G 0 W 0 for predicting excitation energies of molecular systems is largely affected by the starting point 40 . G 0 W 0 based on range-separated functionals and tuned hybrid functionals provides more accurate QP energies 55,[62][63][64][65] , which lead to better excitation energies in BSE/G 0 W 0 . It has been shown that BSE/G 0 W 0 with optimally tuned hybrid functionals and range-separated functionals predicts accurate core electron excitation energies and low-lying excitation energies 66,67 . One path to improve the accuracy is to introduce the self-consistency into the GW calculations. BSE combining with the eigenvalue self-consistent GW (evGW ) approach 32 is shown to predict accurate excitation energies for organic molecules and CT systems with the reduced starting point dependence 36,38,40,68 . It has been shown that BSE/evGW provides comparable accuracy to TDDFT with hybrid functionals and the difference of excitation energies originating from using different DFAs in BSE/evGW is only around 0.1 to 0.3 eV 36,37,68 .
In practice, few iterations are necessary for evGW calculations to reach the convergence 32 . The 3 additional computational cost in evGW is only prohibited for large systems.
While Green's function theory provides QP energies by construction, there are parallel development in DFT. Within DFT, the physical meaning of the one-electron orbital energies of the frontier orbitals, namely HOMO and LUMO, has been established based on the three key theoretical results. First, the Janak theorem links Kohn-Sham orbital energies to the derivatives of the total energy with respect to the orbital occupation numbers, which are, however, not related to any physical observables directly 69 . Second, the derivatives of the total energy with respect to the total electron number, the chemical potentials, are respectively the negative of the first ionization potential (IP) and the first electron affinity (EA) for the exact functional based on the linear condition on the behavior of energy for fractional charges [70][71][72] . Third, the chemical potentials were established to be equal to the derivatives of the total energy with respect to the HOMO/LUMO orbital occupation numbers in the Kohn-Sham calculation with XC energy being functionals of the density, or the generalized Kohn-Sham calculation with XC energy being functionals of the noninteracting one-electron density matrix 73 . Combining these theoretical results, the HOMO and LUMO energies are thus the approximation to the negative of the first IP and the first EA, respectively, as first established by Cohen et. al 73 . This interpretation of the frontier orbital energies has been further extended for other orbitals: Kohn-Sham or generalized Kohn-Sham orbital energies are corresponding QP energies, from the DFA used 74 . This extension was based on extensive observation for a large set of molecules that DFAs, which were designed with minimal delocalization error and provide accurate prediction of IP and EA from the HOMO and LUMO energies, also predict other QP energies from the corresponding generalized Kohn-Sham orbital energies with similar excellent accuracy as HOMO and LUMO orbitals. Therefore, accurate QP energies can be provided from the ground state calculations of DFT 74 .
We want to leverage these recent developments of DFT within the BSE formalism for electronic excitations to achieve and improve the accuracy of the BSE/G 0 W 0 approach at more affordable computational cost, based on the localized orbital scaling correction (LOSC) 75 . Let us briefly review the development of LOSC and other DFAs. Over past decades, DFT [76][77][78] has become the most popular tool in the electronic structure theory. In DFT, the complicated electron correlation effect can be properly and efficiently described by DFAs, including local density approximations (LDAs) 79,80 , generalized gradient approximations (GGAs) [81][82][83] and hybrid functionals 81,84 . How-ever, the predictive power of DFT is impaired by some intrinsic deficiencies. It has been shown that the delocalization error 85,86 is responsible for many failures in DFT, such as the band-gap prediction 86 . The delocalization error in mainstream DFAs is manifested in small molecules as the violation of the Perdew-Parr-Levy-Balduz (PPLB) condition [70][71][72] showing that the total energy of a system as a function of the electron number should be piecewise linear between energies at integer points. In 2011, global scaling correction (GSC) 87 was developed to impose the PPLB condition by using canonical occupations and curvatures that are constructed from canonical orbitals. It has shown that GSC largely restores the linearity behavior and predicts accurate band gaps for systems of all sizes. However, GSC offers zero corrections to total energies at integers, which implies the correction is not size-consistent. To provide a size-consistent correction, local scaling correction (LSC) 88,89 was developed by using local fractional information. But LSC has numerical difficulties for capturing tiny fractions 75,88 . To combine merits of describing global fractions in GSC and local fractions in LSC, LOSC was developed to systematically eliminate the delocalization error in a size-consistent manner by utilizing orbitalets 75 . Orbitalets are defined as a set of orbitals localized in both physical and energy spaces that are obtained by the restrained Boys localization 75 . By using orbitalets, the LOSC correction can be applied to the global and local regions of the system in a dynamical way. It has been shown that LOSC successfully describes dissociation of cationic species, band gaps and photoemission spectrum 74,75 . As shown in the recent work from our group, accurate QP energies can be approximated from LOSC 74 . LOSC provides very similar or better photoemission spectrums and fundamental gaps to those obtained from the fully self-consistent GW (scGW ) approach 74 . As shown in Section.5 in the Supporting Information, G 0 W 0 with conventional DFAs gives underestimated fundamental gaps compared with evGW . LOSC provides larger fundamental gaps, which are similar to the evGW level. In addition, Fundamental gaps obtained from LOSC with hybrid functionals are closer to those obtained from evGW than LOSC with GGA functionals. Furthermore, there are further development of LOSC approach 90 and a open-source software 91

II. THEORY
A. Localized orbital scaling correction The LOSC correction to the total energy is expressed as 75 where λ is the local occupation matrix and κ is the curvature matrix. In Eq.1 and following equations, we use i, j for occupied orbitals, a, b for virtual orbitals and p, q for general orbitals.
The local occupation matrix λ in Eq.1 is defined as 75 where ρ s is the KS density matrix and {φ p } is the set of orbitalets. The diagonal elements of λ contain the information of the fractional electron distribution and the off-diagonal elements bring corrections to the unphysical interaction between the local fractions centered at different positions 75 .
Orbitalets {φ p } used in λ and κ are obtained the restrained Boys localization, which minimizes the following spread function 75 where w is the penalty function defined in Ref. 75 and U is the unitary transformation matrix from canonical orbitals to orbitalets. This restrained Boys localization mixes both occupied and virtual orbitals, which is different from the original Boys localization that only mixes occupied orbitals 93 .
In Eq.4, the first term that has the similar form as the spread function in the original Boys localization ensures orbitals localized in the physical space, and the second term forbids mixing between canonical orbitals that are far apart in the energy space by using the penalty function w.
More recently, a modified localization function was developed which leads to similar accuracy, but preserving the degeneracy of orbitals 90 . We have used the localization function of Eq.4 in this work.
With the LOSC correction defined in Eq.1 the total energy is expressed as And the LOSC correction to the Hamiltonian is where The LOSC orbital energies are obtained by diagonalizing H in Eq.6.
The LOSC correction can be applied in either the self-consistent manner (SCF LOSC) or the post-SCF manner (post-SCF LOSC). As shown in Ref. 75, SCF LOSC provides improved orbitals and corrects the wrong electron density caused by the delocalization error. As shown in Table   S18 in Section.6 in the Supporting Information, dipole moments of tested molecules obtained from SCF LOSC are similar to those obtained from quasiparticle-self-consistent GW (qsGW ) and scGW . This indicates that the electron density from LOSC is similar to those from qsGW and scGW , which are also more localized compared with the electron density from KS-DFT 32,94 . SCF LOSC with hybrid functionals provides better agreement with GW results, which is similar to fundamental gaps. However, as discussed in Ref. 32 and shown in numerical results in Table.S2 and Table.S7 in the Supporting Information, orbital updates have a minor effect on QP energies for the systems studied in this work. Therefore, we do not expect SCF to change the results.
In addition to the minor effects from using SCF LOSC orbitals, as shown in numerical results of where ω is the excitation energies. The A, B matrices in Eq.8 are defined as where {ε p } are input orbital energies. In Eq.8, v is the Coulomb interaction defined as where {ϕ p } is the set of input orbitals. And W is the screened interaction defined as where the dielectric function is calculated by the static response function χ 50,51 ε pq,rs = δ pr δ qs − v pq,rs χ rs,rs .
The Tamm-Dancoff approximation (TDA) in BSE is obtained by neglecting the B matrix in Eq.8, which is denoted as BSE-TDA. TDA has been used in TDHF and TDDFT to overcome the triplet instability problem 95 Table.I, using TDA leads to increased excitation energies and improves the accuracy for triplet excitations in BSE/LOSC and BSE/G 0 W 0 , which largely underestimate triplet excitation energies. However, one should be careful when using TDA in the BSE/GW approach. As shown in the present work and Ref.

A. Truhlar-Gagliardi test set
37, using TDA in BSE/evGW leads to similar or worse results and provides minor effects for singlet calculations. In addition, as shown in recent studies using TDA in BSE/GW can lead to blue-shifts in nanosized systems [115][116][117] and worse estimations for singlet-triplet energy gaps in organic molecules 36