Spin-flip Bethe-Salpeter equation approach for ground and excited states of open-shell molecules and defects in solids

Open-shell systems such as magnetic molecules or defects with a triplet ground state are challenging to describe in electronic structure methods, but are of great interest for quantum information applications. We demonstrate a spin-flip approach within the Bethe-Salpeter equation to calculate ground and excited states of open-shell molecules and defected solids. The approach works in periodic boundary conditions without any need for embedding or selection of a subspace. Our benchmark results for the torsion potential-energy surface of ethylene and the optical excitations of the diamond NV$^{-}$ center show excellent agreement with the literature, and a low or moderate level of spin contamination.

Introduction.-Defects in solids have been a focus of great recent interest for quantum information. A canonical example is the diamond NV − center, consisting of substitutional N atom adjacent to a vacancy, and having an extra negative charge. This system has wideranging applications, from qubits [1][2][3][4][5][6] to nanosensing [7][8][9][10], and its optical properties are well studied both experimentally [11][12][13][14][15] and theoretically [16][17][18][19][20]. catalysts, and single-photon emitters There has been a search for other defects with appropriate electronic structure, in bulk as well as two-dimensional materials recently, especially on the computational side [21][22][23]. However, computing their electronic structure poses challenges due to their open-shell electronic configurations, meaning electrons singly occupy one or more energy levels.The complications of such systems are well appreciated in the quantum-chemistry community [24,25] but only recently have gathered attention in condensed-matter physics.
Density-functional theory (DFT) is the standard workhorse method for electronic structure, but it can only control for the total magnetic quantum number M S and not the total spin quantum number S, meaning that solutions are not necessarily eigenstates ofŜ 2 and therefore do not belong to a proper spin manifold (spin contamination). Further, DFT is single-reference, with the non-interacting Kohn-Sham (KS) system implying a single Slater determinant for the many-body wavefunction. Many investigations of open-shell defects rely on the simple but crude approximation of constrained DFT [19,21,26,27], in which different ground or excited states are described by simply altering the occupations of single-particle KS states. There is not even a distinction between a singlet and triplet of M S = 0, a significant limitation in the ability to describe spin physics.
Multi-reference methods from quantum chemistry that explicitly account for the multiconfigurational character of the many-body wavefunction are applicable for molecular systems but are quite challenging to apply directly to solid-state systems [28] including defects. Instead, many works have instead used embedding schemes [29], in which a finite subspace of in-gap defect states is selected, and then treated with configuration interaction or other correlated methods, in the presence of the solid environment [19,[30][31][32][33]. Difficult questions of how to correct double-counting of correlation arise in such theories [34,35].
The spin-flip method [36] allows for the calculation of ground-and excited-state energies for systems with multiconfigurational character, such as open-shell systems or molecules breaking bonds, from a single-reference wavefunction. The method computes excitation energies from a high-spin reference state, where S = M S , to a set of lower-spin target states with M S − 1, as diagrammed in Fig. S1 [37]. While DFT's Hohenberg-Kohn Theorem [38] requires a non-degenerate ground state, DFT is able to compute the lowest-energy state for a given symmetry [20,39] The choice S = M S > 0 for the high-spin reference state is made throughout this work. The target states with one spin flipped provide a basis of transitions. Some target states may include effective double excitations with respect to the ground state, allowing for the description of conical intersections as per Brillouin's Theorem [40].
The spin-flip method applied to TDDFT (SF-TDDFT) has demonstrated success for the description of antiferromagnetic molecules [41,42], conical intersections [43], and ordering of near-degenerate ground-and excitedstates. However, SF-TDDFT relies on the kernel f xc to describe the interactions between electrons and holes. Conventional options give incorrect long-range behavior [44,45]. Moreover, there is a zoo of possible options and unclear which to choose [46,47]. Collinear semilocal kernels, as would be typically used in TDDFT, do not mix different transitions, and so to be useful SF-TDDFT requires hybrid functionals with an unusually large fraction of exact exchange [46] or noncollinear functionals [48].
Instead, in this Letter we draw on the many-body perturbation theory approach of GW -BSE [44,[49][50][51][52]  interaction W . This allows us to describe excitonic effects [52], valence excitations, Rydberg excitations, and charge-transfer excitations [53]. Our method works naturally in a periodic cell, like conventional GW /BSE, without the need for any selection of a subspace, embedding technique, or double-counting corrections.
Independent from the present effort [54], Monino and Loos [47] recently considered an SF-BSE approach, applied to atoms and molecules with a detailed comparison to different quantum chemistry and TDDFT approaches. Our work marks a major advance in applying SF-BSE to defects in solids, demonstrating high accuracy for potential-energy surfaces, and identifying the problematic nature of open-shell GW for use in BSE calculations.
The target states are the collection of the spin-flip excited states which preserve particle number, |N, I . As single excitations, these can be described with the Bethe-Salpeter equation, which can be written in the usual Tamm-Dancoff approximation [52] for the case of a spinflip from α to β: (1) with quasiparticle energies E QP for conduction states c and valence v, an excitation energy Ω, and an eigenvector A. The BSE kernel K has direct and exchange terms, K D and K X , respectively. As with the triplet excitations in usual BSE [52], the K X matrix elements must be zero for spin-flips due to the orthogonality of the α and β spinors.
Single-particle excitations of the high-spin reference state, in general, have multiple quasiparticle energies E QP for a given electron addition/removal [20]. This issue has been neglected in some previous calculations [60]. It remains an open question how such quasiparticle multiplets (no longer indexable by α, β) should be correctly used in BSE calculations, an issue which is not addressed in the SF-BSE work of Monino and Loos [47]. Errors on the order of the multiplet splittings (∼ 1 − 2 eV in molecules and defects [20]) could result. Using quasiparticle energies from the Generalized Plasmon Pole (GPP) approximation [51] sidesteps this difficulty but at the cost of ignoring the correct effects of the multiplet splittings [20]. Instead we take the pragmatic approach of simply using the KS energy eigenvalues E KS instead of E QP . (By contrast, neglecting the kernel produces completely wrong results [37].) The ground-state energy is the sum E DFT tot (Ref)+Ω Imin , using the spin-flip state with the lowest energy, possibly negative. As we will see, the gap correction due to GW is mostly just a constant offset [37], which does not affect differences between various values of Ω I .
Ethylene.-Ethylene is a well-studied system, both for vertical transition energies from its ground-state and for its torsion potential, and is a classic test for spin-flip methods [46]. The torsion potential twists the π bond, making the π and π * orbitals become degenerate at 90 • [37], making it an open-shell system with a triplet ground state.
We obtained the M S = 1 reference state by constraining occupations. We performed extensive tests of convergence, finding similar characteristics to ordinary BSE, and demonstrating how SF-BSE does not require a priori knowledge to select an appropriate subspace.
The N, T, V, and Z potential energies for ethylene under torsion, listed with their highest contributing Slater determinants at zero torsion angle, are shown in Fig. 1. The excitation energies computed within GW -BSE required 500 empty states for the BSE kernel to achieve convergence, while SF-BSE converges with about 55 empty states [37]. We note that the SF-BSE Hamiltonian was not constructed using quasiparticle energies, while displaying remarkably high agreement with the (conventional) GW /BSE calculations (Table I). This is an important demonstration of consistency with the usual GW /BSE. Examination of the relation between KS eigenvalues of the high-spin reference and the closed-shell ground state shows that they actually differ by the quasiparticle corrections [37], emulating GW . This is suggestive of studies such as [61] in which a local functional, used in a particular way, is able to provide a derivative discontinuity to correct the gap. We note that by contrast neglecting quasiparticle energies with conventional BSE gives completely wrong excitation values of -3.00 and 1.27 eV for N → T and N → V, respectively. Our results also are in good agreement with experiment and with quantum Monte Carlo. Using GPP at 0 • also provides similar results [37].
Our N state can be compared with the usual DFT ground state at 0 • , and is slightly lower which seems to be due to improving static correlation from mixing of a doubly-excited configuration. This energy difference is comparable conceptually and numerically to the differences in estimating atom energies by GW as HOMO of a neutral atom or LUMO of a cation [62]. We find that restricted calculations for ethylene are generally worse [37], as for molecules in [47].  The torsion potential is in excellent agreement with the "gold standard" Two-Configuration Self-Consistent Field-CISD [46], and far closer than SF-TDDFT methods (Fig. 2). The barrier, which is the energy difference between the singlet N state at 90 • torsion and no torsion, is computed in SF-BSE as 3.25 eV and is 3.27 eV in TCSCF-CISD [46]. We note that using GPP for the potential-energy surfaces produces significantly worse results [37].
NV − center in diamond.-To demonstrate the applicability to defects, we consider the NV − center. There is no need for constraining occupations, as the DFT and true ground states are triplets. We calculate zero-phonon line (ZPL) transition energies from the total-energy differences E I − E G.S. in SF-BSE. A comparison between the in-gap defect states computed with restricted and unrestricted KS orbitals is in Fig. 3. Results using KS eigenvalues have an appreciable difference only for the excitation energy to 3 E, 0.30 eV. Unrestricted KS and GPP results (with 300 empty states in the Coulomb-hole sum [67]) are fairly close.
A comparison between SF-BSE and other methods is presented in Fig. 4. The 1 A 1 excitation energy shows appreciable difference from experiment in SF-BSE and other methods that cannot include a contribution from the determinant in which the up-and down-spin v orbitals are empty and the e x and e y orbitals are occupied, The ethylene torsion potential, comparing SF-BSE to SF-TDDFT and TCSCF-CISD results from [46]. The dashed curve is the difference between the SF-BSE and the TCSCF-CISD results, multiplied by ten to be seen more clearly.
FIG. 3. Energies of the in-gap many-body defect states for the NV − defect in diamond, computed with SF-BSE, with (a) restricted KS eigenvalues, (b) unrestricted KS eigenvalues, and (c) unrestricted GPP quasiparticle energies [51]. The red arrow is the experimental value [14] for the vertical optical transition, and the green double-headed arrow is the experimental value for the singlet splitting [15].
Spin contamination.-Excited-state theories, especially spin-flip methods, encounter the problem of spin contamination, which can be a diagnostic of poor energies or wavefunctions, or also simply cause difficulty in identifying the spin character of the calculated states. The eigenvector describing the excited state yields an expectation value ofŜ 2 that may differ from proper values of S (S + 1) such 0 (singlet) or 2 (triplet). Sources of spin contamination include the difference in spatial wavefunctions for the different spin channels for spin-polarized DFT [68], or the absence of particular transitions from the set of target states that are required to form an eigenvector ofŜ 2 [69]. Projection methods can be used to FIG. 4. Energies of the in-gap many-body defect states for the NV − defect in diamond, computed with SF-BSE compared to various other results in the literature: (b) [30], (c) [60], (d) and (e) [19], (f) and (g) [32], and (h) [33]. The red arrow is the experimental value [14] for the vertical optical transition, and the green double-headed arrow is the experimental value for the singlet splitting [15]. correct spin contamination and energies [70,71].
If a many-body wavefunction is available, the the Löwdin formula [37,72] can be used to compute N, I|Ŝ 2 |N, I for state I. Löwdin results using the KS Slater determinant have been found to be similar to those from density functional approximations to Ŝ 2 [73], and we find results very close to 2 for our high-spin reference states [37]. Within an excited-state theory, it is natural to evaluate N, I|Ŝ 2 |N, I via the difference ∆ Ŝ 2 I from the ground state value, which can be computed from only the excitation eigenvector. The result for ∆ Ŝ 2 I with a spin-flip (different from the spin-conserving expression) is given in [74] in the context of SF-TDDFT, and detailed in our Supplementary Material [37]. The related spin-conserving expression of [75] could be used to for GW -BSE in the presence of spin-orbit coupling [76]. We note that in TDDFT such expressions are an approximation, based on treating the KS orbitals as a Slater determinant [75]. By contrast, in BSE, we can rigorously calculate ∆ Ŝ 2 I consistently with the BSE derivation, as noted in [77], since only the exciton wavefunctions are needed [50,52].
We show results of Ŝ 2 I in Fig. 5. We find a small spin contamination comparable to TDDFT [69] that still allows identifying singlet and triplet states. Spin contamination is substantially larger in NV − [37], as others have found for solids [71]. States where the needed transitions are available still have values of Ŝ 2 I that enable identification of singlet and triplets.
Conclusion.-We have presented the SF-BSE method to calculate electronic structure for open-shell molecules and defects in solids. The ethylene torsion potential shows high agreement with a high-level multi-reference method, considerably improved from SF-TDDFT. Due to the problematic nature of GW multiplets for BSE, we use KS energies as input. The optical transition energies in the equilibrium geometry show remarkable agreement with conventional GW -BSE and experiment, while requiring roughly one-tenth of the empty orbitals and no single-particle quasiparticle energies. This surprisingly successful approach may be promising for optical calculations of closed-shell systems. The in-gap defect states for the NV − in diamond show good agreement with experiment and multi-reference methods, apart from the toolow 1 A 1 state. Spin contamination is small for ethylene and manageable for NV − . This method will be available in a future release of BerkeleyGW.
Our results show SF-BSE is promising for potentialenergy surfaces in defects [19]. This method, unlike embedding schemes, can easily provide forces through BSE excited-state forces schemes [78][79][80] for efficient exploration and investigation of Stokes shifts [81], Jahn-Teller distortions [26,82], and internal conversion [19]. Other future directions for the SF-BSE formalism includes use of quasiparticle multiplets [20] or spinors [76], multiplet reference states [69,74] or other reference states to access different excited states.
The authors thank Hrant Hratchian and Abdulrahman Zamani for useful discussions and John Parkhill for originally suggesting the concept This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, CTC and CPIMS Programs,