Quantum simulations of materials on near-term quantum computers

Quantum computers hold promise to enable efficient simulations of the properties of molecules and materials; however, at present they only permit ab initio calculations of a few atoms, due to a limited number of qubits. In order to harness the power of near-term quantum computers for simulations of larger systems, it is desirable to develop hybrid quantum-classical methods where the quantum computation is restricted to a small portion of the system. This is of particular relevance for molecules and solids where an active region requires a higher level of theoretical accuracy than its environment. Here we present a quantum embedding theory for the calculation of strongly-correlated electronic states of active regions, with the rest of the system described within density functional theory. We demonstrate the accuracy and effectiveness of the approach by investigating several defect quantum bits in semiconductors that are of great interest for quantum information technologies. We perform calculations on quantum computers and show that they yield results in agreement with those obtained with exact diagonalization on classical architectures, paving the way to simulations of realistic materials on near-term quantum computers.


Introduction
In the last three decades, atomistic simulations based on the solution of the basic equation of quantum mechanics have played an increasingly important role in predicting the properties of functional materials, encompassing catalysts and energy storage systems for energy applications, and materials for quantum information science. Especially in the case of complex, heterogeneous materials, the great majority of first-principles simulations are conducted using density functional theory (DFT), which is in principle exact but in practice requires approximations to enable calculations. Within its various approximations, DFT has been extremely successful in predicting numerous properties of solids, liquids and molecules, and in providing key interpretations to a variety of experimental results; however it is often inadequate to describe so-called strongly-correlated electronic states [1,2]. We will use here the intuitive notion of strong correlation as pertaining to electronic states that cannot be described by mean-field theories.
Several theoretical and computational methods have been developed over the years to treat systems exhibiting strongly-correlated electronic states, including dynamical mean-field theory [3,4] and quantum Monte-Carlo [5,6]; in addition, ab initio quantum chemistry methods, traditionally developed for molecules, have been recently applied to solid state problems as well [7]. Unfortunately, these approaches are computationally demanding and it is still challenging to apply them to complex materials containing defects and interfaces, even using high-performance computing architectures.
Quantum computers hold promise to enable efficient quantum mechanical simulations of weakly and strongly-correlated molecules and materials alike [8,9,10,11,12,13,14,15]; in particular when using quantum computers, one is able to efficiently simulate systems of interacting electrons, with an exponential speedup over exact diagonalization on classical computers (full configuration interaction, FCI, calculations). Thanks to decades of successful experimental efforts, we are now entering the noisy intermediate-scale quantum (NISQ) era [16], with quantum computers expected to have on the order of 100 quantum bits (qubits); unfortunately this limited number of qubits still prevents straightforward quantum simulations of realistic molecules and materials, whose description requires hundreds of atoms and thousands to millions of degrees of freedom to represent the electronic wavefunctions. An important requirement to tackle complex chemistry and material science problems using NISQ computers is the reduction of the number of electrons treated explicitly at the highest level of accuracy [17,18]. For instance, building on the idea underpinning dynamical mean field theory [3,4], one may simplify complex molecular and material science problems by defining active regions (or building blocks) with stronglycorrelated electronic states, embedded in an environment that may be described within mean-field theory [19,20,21].
In this work, we present a quantum embedding theory built on DFT, which is scalable to large systems and which includes the effect of exchange-correlation interactions of the environment on active regions, thus going beyond commonly adopted approximations. In order to demonstrate the effectiveness and accuracy of the theory, we compute ground and excited state properties of several spin-defects in solids including the negatively charged nitrogen-vacancy (NV) center [22,23,24,25,26,27,28], the neutral silicon-vacancy (SiV) center [29,30,31,32,33,34] in diamond, and the Cr impurity (4+) in 4H-SiC [35,36,37]. These spin-defects are promising platforms for solid-state quantum information technologies, and they exhibit strongly-correlated electronic states that are critical for the initialization and read-out of their spin states [38,39,40,41,42,43]. Our quantum embedding theory yields results in good agreement with existing measurements. In addition, we present the first theoretical predictions for the position and ordering of the singlet states of SiV and of Cr, and we provide an interpretation of experiments which have so far remained unexplained.
Importantly, we report calculations of spin-defects using a quantum computer [44,45]. Based on the effective Hamiltonian derived from the quantum embedding theory, we investigated the stronglycorrelated electronic states of the NV center in diamond using quantum phase estimation algorithm (PEA) [46,8] and variational quantum eigensolvers (VQE) [47,48,49], and we show that quantum simulations yield results in agreement with those obtained with classical FCI calculations. Our findings pave the way to the use of near term quantum computers to investigate the properties of realistic heterogeneous materials with first-principles theories.

General strategy
We summarize our strategy in Fig. 1. Starting from an atomistic structural model of materials (e.g. obtained from DFT calculations or molecular dynamics simulations), we identify active regions with strongly-correlated electrons, which we describe with an effective Hamiltonian that includes the effect of the environment on the active region. This effective Hamiltonian is constructed using the quantum embedding theory described below, and its eigenvalues can be obtained by either classical algorithms such as exact diagonalization (FCI) or quantum algorithms.

Embedding Theory
A number of interesting quantum embedding theories have been proposed over the past decades [50]. For instance, density functional embedding theory has been developed to improve the accuracy and scalability of DFT calculations [51,52,53,54,55]. Density matrix embedding theory (DMET) [56,57,58] and various Green's function based approaches [59,60], e.g. dynamical mean field theory (DMFT), have been developed to describe systems with strongly-correlated electronic states. At present, ab initio calculations of materials using DMET and DMFT have been limited to relatively small unit cells (a few tens of atoms) of pristine crystals, due to their high computational cost [61,62]. In this work, we present a quantum embedding theory that is applicable to strongly-correlated electronic states in realistic heterogeneous materials and we apply it to systems with hundreds of atoms. The theory, inspired by the constrained random phase approximation (cRPA) approach [63,64,65], does not require the explicit evaluation of virtual electronic states [66,67], thus making the method scalable to materials containing thousands of electrons. Furthermore, cRPA approaches contain a specific approximation (RPA) to the screened Coulomb interaction, which neglects exchange-correlation effects and may lead to inaccuracies in the description of dielectric screening. Our embedding theory goes beyond the RPA by explicitly including exchange-correlation effects, which are evaluated with a recently developed finite-field algorithm [68,69].
The embedding theory developed here aims at constructing an effective Hamiltonian operating on an active space (A), defined as a subspace of the single-particle Hilbert space: Here t eff and V eff are one-body and two-body interaction terms that take into account the effect of all the electrons that are part of the environment (E) in a mean-field fashion, at the DFT level. An active space can be defined, for example, by solving the Kohn-Sham equations of the full system and selecting a subset of eigenstates among which electronic excitations of interest take place (e.g. defect states within the gap of a semiconductor or insulator). To derive an expression for V eff that properly accounts for all effects of the environment including exchange and correlation interactions, we define the environment density response function (reducible polarizability) difference between the polarizability of the Kohn-Sham system χ 0 and its projection onto the active space χ A 0 (see Supplementary Information (SI)). χ E thus represents the density response outside the active space. The term f = V + f xc is often called the Hartree-exchange-correlation kernel, where V is the Coulomb interaction and the exchange-correlation kernel f xc is defined as the derivative of the exchange-correlation potential with respect to the electron density. We define the effective interactions between electrons in A as given by the sum of the bare Coulomb potential and a polarization term arising from the density response in the environment E. When the RPA is adopted, the exchange-correlation kernel f xc is neglected in Eq. 2 and the expression derived here reduces to that used within cRPA. We represent χ E and f on a compact basis obtained from a low-rank decomposition of the dielectric matrix [66,67] that allows us to avoid the evaluation and summation over virtual electronic states. Once V eff is defined, the one-body term t eff can be computed by subtracting from the Kohn-Sham Hamiltonian a term that accounts for Hartree and exchange-correlation effects in the active space (see SI).

Embedding theory applied to spin-defects
The embedding theory presented above is general and can be applied to a variety of systems for which active regions, or building blocks, with strongly-correlated electronic states may be identified: for example active sites in inorganic catalysts or organic molecules or defects in solids and liquids (e.g. solvated ions in water). Here we apply the theory to spin-defects including NV and SiV in diamond and Cr in 4H-SiC. Most of these defects' excited states are strongly-correlated (they cannot be represented by a single Slater determinant of single-particle orbitals), as shown e.g. for the NV center in diamond by Bockstedte et al. [70]. We demonstrate that our embedding theory can successfully describe the manybody electronic structure of different types of defects including transition metal atoms; our results not only confirm existing experimental observations but also provide a detailed description of the electronic structure of defects not presented before, which sheds light into their optical cycles.
We first performed spin-restricted DFT calculations using hybrid functionals [71] to obtain a meanfield description of the defects and of the whole host solid. The spin restriction ensures that both spin channels are treated on an equal footing and that there is no spin-contamination when building effective Hamiltonians. Based on our DFT results, we then selected active spaces that include single-particle defect wavefunctions and relevant resonant and band edge states. We verified that the size of the chosen active spaces yields converged excitation energies (see SI). We then constructed effective Hamiltonians (Eq. 1-2) by taking into account exchange-correlation effects, and we obtained many-body ground and excited states using classical (FCI) and, for selected cases, quantum algorithms (PEA, VQE). All calculations were performed at the spin triplet ground state geometries obtained by spin-unrestricted DFT calculations, thus obtaining vertical excitation energies (equal to the sums of zero phonon line (ZPL) and Stokes energies). It is straightforward to extend the current approach to compute potential energy surfaces at additional geometries, so as to include relaxations and Jahn-Teller effects [70,34]. In Fig. 2 we present atomistic structures, single-particle defect levels and the many-body electronic structure of three spin-defects. Several relevant vertical excitation energies are reported in Table 1, and additional ones are given in the SI. In the following discussion, lower-case symbols represent single-particle states obtained from DFT and upper-case symbols represent many-body states.
For the NV in diamond, we constructed effective Hamiltonians (Eq. 1) by using an active space that includes a 1 and e single-particle defect levels in the band gap and states near the valence band maximum (VBM). Our FCI calculations correctly yield the symmetry and ordering of the low-lying The vertical excitation energies reported in Table 1 show that including exchange-correlation effects yields results in better agreement with experiments than those obtained within the RPA.
In the case of the SiV in diamond, we built effective Hamiltonians using an active space with the e u and e g defect levels and states near the VBM, including resonant e u and e g states. Effective Hamiltonians including or neglecting exchange-correlation effects yield similar results, with the excitation energies obtained beyond RPA being slightly higher. We predicted the first optically-allowed excited state to be a 3 E u state with vertical excitation energy of 1.59 eV, in good agreement with the sum of 1.31 eV ZPL measured experimentally [29] and 0.258 eV Stokes shift estimated using an electron-phonon model [34]. Remarkably, our calculations predicted a 3 A 2u state 11 meV below the 3 E u state, supporting a recent experimental observation by Green et al. [33], which proposed a 3 A 2u -3 E u manifold with 7 meV separation in energy. The small difference in energy splitting between our results and experiment is likely due to geometry relaxation effects not yet taken into account in our study. In addition to states of u symmetry generated by e u → e g excitations, we observed a number of optically dark states of g symmetry (grey levels in Fig. 2b) originating from the excitation from the e g level and the VBM states to the e g level.
Despite significant efforts [31,32,33,34], several important questions on the singlet states of SiV remain open. These states are crucial for a complete understanding of the optical cycle of the SiV center. Our predicted ordering of singlet states of SiV is shown in Fig. 2b. We find the vertical excitation energies of the 1 A 1u state to be slightly higher than that of the 3 A 2u -3 E u triplet manifold, suggesting that the intersystem crossing (ISC) from 3 A 2u or 3 E u to singlet states may be energetically unfavorable (first-order ISC to lower 1 E g and 1 A 1g states are forbidden). We note that the 1 E u and 1 A 2u states are much higher in energy than 1 A 1u and are not expected to play a significant role in the optical cycle. In addition the first-order ISC process from the lowest energy singlet state 1 E g to the 3 A 2g ground state is forbidden by symmetry. Overall our results indicate that the 3 A 2g state is populated through higher-order processes and therefore the spin-selectivity of the full optical cycle is expected to be low.
While more detailed studies including spin-orbit coupling are required for definitive conclusions, our predictions shed light on the strongly-correlated singlet states of SiV and provide a possible explanation for the experimental difficulties in measuring optically-detected magnetic resonance (ODMR) of SiV.
We now turn to Cr in 4H-SiC, where we considered the hexagonal configuration. We constructed is expected to be small given the large Debye-Waller factor [37]. There is currently no experimental report for the triplet excitation energies of Cr in 4H-SiC, but our results are in good agreement with existing experimental measurements for Cr in GaN, a host material with a crystal field strength similar to that of 4H-SiC [36]. We predict the existence of a 3 E + 3 A 1 manifold at 1.4 eV and a 3 E + 3 A 2 manifold at 1.7 eV above the ground state (Fig. 2c), resembling the 3 T 2 manifold (1.2 eV) and 3 T 1 manifold (1.6 eV) for Cr in GaN observed experimentally [72]. We note that in many cases it is challenging to study materials containing transition metal elements with DFT [73]. The agreement between FCI results and experimental measurements clearly demonstrates that the embedding theory developed here can effectively describe the strongly-correlated part of the system, while yielding at the same time a quantitatively correct description of the environment.

Quantum simulations of spin-defects
The results presented in the previous section were obtained using classical algorithms. We now turn to the use of quantum algorithms. To perform quantum simulations with PEA and VQE, we constructed a minimum model of an NV center including only a 1 and e orbitals in the band gap. This model (4 electrons in 6 spin orbitals) yields excitation energies within 0.2 eV of the converged results using a larger active space. In Fig. 3 we show the results of quantum simulations.
We first performed PEA simulations with a quantum simulator (without noise) [44] to compute the energy of 3 A 2 , 3 E, 1 E and 1 A 1 states. We used molecular orbital approximations of these states derived from group theory [24] as initial states for PEA, which are single Slater determinant for 3 A 2 (M S = 1) and 3 E (M S = 1) states, and superpositions of two Slater determinants for 1 E and 1 A 1 states. As shown in Fig. 3a, PEA results converge to classical FCI results with an increasing number of ancilla qubits.
We then performed VQE simulations with a quantum simulator and with the IBM Q 5 Yorktown quantum computer [45]. We estimated the energy of the 3 A 2 ground state manifold by performing state, so we started with a higher energy reference state |ae xēx e y in the 3 E manifold. We used unitary coupled-cluster single and double (UCCSD) ansatzs [47] to represent the trial wavefunctions. Fig. 3b and Fig. 3c show the estimated ground state energy as a function of the number of VQE iterations, where VQE calculations performed with the quantum simulator correctly converges to the ground state energy in both the M S = 1 and M S = 0 case. Despite the presence of noise, whose characterization and study will be critical to improve the use of quantum algorithms [74], the results obtained with the quantum computer converge to the ground state energy within a 0.2 eV error. Calculations of excited states with quantum algorithms will be the focus of future works.

Discussion
With the goal of providing a strategy to solve complex materials problems on NISQ computers, we proposed a first-principles quantum embedding theory where appropriate active regions of a material and their environment are described with different levels of accuracy, and the whole system is treated quantum mechanically. In particular, we used hybrid density functional theory for the environment, and we built a many-body Hamiltonian for the active space with effective electron-electron interactions that include dielectric screening and exchange-correlation effects from the environment. Our method overcomes the commonly used random phase approximation, which neglects exchange-correlation effects; importantly it is applicable to heterogeneous materials and scalable to large systems, thanks to the novel algorithms used here to compute response functions [68,69]. We emphasize that the embedding theory presented here provides a flexible framework where multiple effects of the environment may be easily incorporated. For instance, dynamical screening effects can be included by considering a frequency-dependent screened Coulomb interaction, evaluated using the same procedure as the one outlined here for static screening; electron-phonon coupling effects can be incorporated by including phonon contributions in the screened Coulomb interactions. Furthermore, for systems where the electronic structure of the active region is expected to influence that of the host material, a self-consistent cycle in the calculation of the screened Coulomb interaction of the environment can be easily added to the approach.
We presented results for spin-defects in semiconductors obtained with both classical and quantum algorithms, and we showed excellent agreement between the two sets of techniques. Importantly, for selected cases we showed results obtained using a quantum simulator and a quantum computer, which agree within a relatively small error, in spite of the presence of noise in the quantum hardware. We predicted several excited states not observed before, in particular our results for the SiV in diamond and Cr in SiC represent the very first theoretical predictions of their singlet states, and provide important insights into their full optical cycle. We also demonstrated that a treatment of the dielectric screening beyond the random phase approximation leads to an accurate prediction of excitation energies.
The method proposed in our work enables calculations of realistic, heterogeneous materials using the resources of NISQ computers. We demonstrated quantum simulations of strongly-correlated electronic states in considerably larger systems (with hundreds of atoms) than previous studies combining quantum simulation and quantum embedding [17,18,19,20,21]. We have studied solids with defects, not just pristine materials, which are of great interest for quantum technologies. The strategy adopted here is general and may be applied to a variety of problems, including the simulation of active regions in molecules and materials for the understanding and discovery of catalysts and new drugs, and of aqueous solutions containing complex dissolved species. We finally note that our approach is not restricted to strongly-correlated active regions and will be useful also in the case of weakly correlated systems, where different regions of a material may be treated with varying levels of accuracy. Hence we expect the strategy presented here to be widely applicable to carry out quantum simulations of materials on near-term quantum computers.

Density Functional Theory
All ground state DFT calculations are performed with the Quantum Espresso code [75] using the planewave pseudopotential formalism. Electron-ion interactions are modeled with norm-conserving pseudopotentials from the SG15 library [76]. A kinetic energy cutoff of 50 Ry is used. All geometries are relaxed with spin-unrestricted DFT calculations using the PerdewBurke-Ernzerhof (PBE) functional [77] until forces acting on atoms are smaller than 0.013 eV/Å. NV and SiV in diamond are modeled with 216-atom supercells; Cr in 4H-SiC is modeled with a 128-atom supercell. The Brillouin zone is sampled with the Γ point.

Construction of effective Hamiltonians
Construction of effective Hamiltonians is performed with the WEST code [67], starting from wavefunctions of spin-restricted DFT calculations. For this step, we remark that the use of hybrid functional is important for an accurate mean-field description of defect levels, even though the geometry of defects are well represented at the PBE level. We used a dielectric dependent hybrid (DDH) functional [71] which self-consistently determines the fraction of exact exchange based on the dielectric constant of the host material. The DDH functional was shown to yield accurate band gaps of diamond and silicon carbide [71], as well as optical properties of defects [39,40,78,79,80]. FCI calculations [82] on the effective Hamiltonian are carried out using the PySCF [7] code.

Quantum simulations
In order to carry out quantum simulations, a minimum model of the NV center is constructed by applying the embedding theory with a 1 and e orbitals beyond the RPA.
In PEA simulations, the Jordan-Wigner transformation [83] is used to map the fermionic effective Hamiltonian to a qubit Hamiltonian, and Pauli operators with prefactors smaller than 10 −6 a.u. are neglected to reduce the circuit depth, which results in less than 10 −4 a.u. (0.003 eV) change in eigenvalues. In order to achieve optimal precision, the Hamiltonian is scaled such that 0 and 2.5 eV are mapped to phases φ = 0 and φ = 1 of the ancilla qubits, respectively. We used the first-order Trotter formula to split time evolution operators into 4 time slices.
In VQE simulations, the parity transformation [9] is adopted. For the simulation of the M S = 1 state, the resulting qubit Hamiltonian acts on 4 qubits and there are 2 variational parameters in the UCCSD ansatz. For the simulation of the M S = 0 state, we fixed the occupation of the a orbital and the resulting qubit Hamiltonian acts on 2 qubits. We replicated the exponential excitation operator twice, with parameters in both replicas variationally optimized. Such a choice results in 6 variational parameters, providing a sufficient number of degrees of freedom for an accurate representation of the strongly-correlated M S = 0 state. Parameters in the ansatz are optimized with the COBYLA algorithm [84].
Quantum simulations are performed with the QASM simulator and the IBM Q 5 Yorktown quantum computer using the IBM Qiskit package [44]. Each quantum circuit is executed 8192 times to obtain statistically reliable sampling of the measurement results.    [22]. † Ref [23]. ‡ Estimated by Ref [28] with a model for intersystem crossing. § Computed with Stokes energy from Ref [34]. || Ref [33]. ¶ Ref [35].