Demonstrating Quantum Computation for Quasiparticle Band Structures

Introduction.— Simulating solid-state materials accurately has been a long-standing issue. A widely used approach relies on first-principles calculations based on the density functional theory (DFT) [1, 2]. Despite its success, DFT fails to describe accurately the electron correlation, and band gaps for semiconductors and insulators are systematically underestimated. Many body perturbation theories such as the GW approximation overcomes this problem partially [3], while strongly correlated systems are still not properly described. Recently, wavefunction-based theories originally developed in the field of quantum chemistry have been successfully applied to solids [4, 5, 6, 7]. Moreover, the equation-of-motion coupled-cluster theory with single, double and triple excitations (EOM-CCSDT) is arguably more accurate than the GW approximation [8]. Nevertheless, the high computational costs of such calculations prohibit practical applications to solid-state materials.

Owing to the recent advances in quantum technologies, quantum computations offer a promising way to efficient and accurate simulations of quantum systems. The two leading quantum algorithms for solving electronic-structure problems are the quantum phase estimation (QPE) [9, 10, 11] and the variational quantum eigensolver (VQE) [12]. QPE is expected to provide an exponential speedup in exact full configuration interaction (FCI) calculations. However, QPE requires a quantum circuit with many gates, and thus it is sensitive to the noise on quantum devices, hindering its practical applications until the realization of large-scale fault-tolerant quantum computers. On the other hand, VQE can be implemented with a shallower quantum circuit, allowing its application on near-term noisy quantum devices as demonstrated for small molecules [12, 13, 14]. There is therefore a scenario in which the practical application of quantum computers might be realized even in the present noisy intermediate-scale quantum (NISQ) era [15].

Recently, VQE has been extended to periodic materials [16, 17, 18, 19, 20, 21, 22, 23]. Interestingly, the quasiparticle band structures can now be obtained using the quantum subspace expansion (QSE)-type formalisms [24, 19, 20], which are variants of EOM. So far, experimental demonstrations of band structures on quantum devices are limited to the tight-binding (one-body) Hamiltonian [17, 23]. However, the inclusion of two-body terms is crucial to treat electron correlations, and it is indispensable for practical first-principles calculations.

In this paper, we experimentally demonstrate the first-principles calculation of a quasiparticle band structure on noisy quantum devices. We treat an ab initio Hamiltonian including the two-body terms. Calculations are based on QSE in conjunction with VQE. To mitigate noise effects on calculation results, we use the qubit-reduction and error-mitigation techniques. Owing to this methodology, we successfully reproduce the ideal noise-free results. Our demonstration is an important step toward practical band-structure calculation on quantum computers.

Ab initio Hamiltonian.— We consider the second quantized representation of ab initio Hamiltonians of periodic systems which is given by

Here, ${\hat{c}}_{p{\mathbf{k}}}^{\dagger}$ (${\hat{c}}_{p{\mathbf{k}}}^{\vphantom{\dagger}}$) is the fermionic creation (annihilation) operator on the $p$th crystalline orbital (CO) with the crystal momentum ${\mathbf{k}}$ obtained through the Hartree-Fock (HF) calculation. The complex coefficients $t_{pq}^{\mathbf{k}}$ and $v_{pqrs}^{{\mathbf{k}}_{p}{\mathbf{k}}_{q}{\mathbf{k}}_{r}{\mathbf{k}}_{s}}$ are one-body and two-body integrals between COs, respectively. The prime on the summation symbol accounts for the momentum conservation law ${\mathbf{k}}_{p}+{\mathbf{k}}_{q}-{\mathbf{k}}_{r}-{\mathbf{k}}_{s}={\mathbf{G}}$, where ${\mathbf{G}}$ is a reciprocal lattice vector of the unit cell.

Hybrid quantum-classical algorithms.— To calculate quasiparticle band structures, we use the QSE method combined with VQE as reported in Refs. [19, 20]. In QSE, the ground state is used as a reference state $|\psi\rangle$. With an excitation operator ${\hat{O}}_{i}$, we construct a low-energy subspace Hamiltonian ${\mathbf{H}}^{\rm sub}$ which is given by $H_{ij}^{\rm sub}=\langle\psi|{\hat{O}}_{i}^{\dagger}{\hat{H}}{\hat{O}}_{j}^{\vphantom{\dagger}}|\psi\rangle$ in the matrix representation. Then we obtain the spectrum of energy eigenvalues (quasiparticle bands) by solving the generalized eigenvalue problem ${\mathbf{H}}^{\rm sub}{\mathbf{C}}={\mathbf{S}}^{\rm sub}{\mathbf{C}}{\mathbf{E}}$. Here, ${\mathbf{C}}$ is the matrix consisting of the eigenvectors and ${\mathbf{E}}$ is the diagonal matrix whose elements are eigenvalues. In addition, the matrix ${\mathbf{S}}^{\rm sub}$ is given by $S_{ij}^{\rm sub}=\langle\psi|{\hat{O}}_{i}^{\dagger}{\hat{O}}_{j}^{\vphantom{\dagger}}|\psi\rangle$. We evaluate the matrix elements $H_{ij}^{\rm sub}$ and $S_{ij}^{\rm sub}$ using quantum devices. The diagonalization for solving the generalized eigenvalue problem is performed on classical computers. As the excitation operator ${\hat{O}}_{l{\mathbf{k}}}^{\vphantom{\dagger}}$, we choose ${\hat{c}}_{l{\mathbf{k}}}^{\vphantom{\dagger}}$ for obtaining the valence bands at each ${\mathbf{k}}$, where $l$ runs over the occupied orbitals. Similarly, ${\hat{c}}_{l{\mathbf{k}}}^{\dagger}$ is used for the conduction bands, where $l$ runs over the unoccupied orbitals.

To prepare the ground state $|\psi\rangle$, we use the VQE method which is based on the variational method with a parameterized wavefunction ansatz. The wavefunction ansatz takes the form of $|\psi({\bm{\theta}})\rangle={\hat{U}}({\bm{\theta}})\ket{\psi_{0}}$, and the parameters ${\bm{\theta}}$ are variationally optimized to minimize the energy expectation value. Here, $\ket{\psi_{0}}$ is an input quantum state and ${\hat{U}}({\bm{\theta}})$ is a unitary operator implemented on a quantum circuit.

Settings.— We consider Si that has the diamond crystal structure with the experimental lattice constant of 5.43 Å, where two atoms are contained in the primitive cell. For the Brillouin zone sampling, a 1$\times$1$\times$1 k-point grid is centered at the target k-point on the band path. For each target k-point, the Hamiltonian is constructed separately. The GTH (Goedecker-Teter-Hutter) pseudpotential [25] is used together with the GTH-SZV (single-zeta valence) basis set. To reduce the number of required qubits, we use the active space approximation and take the highest occupied and lowest unoccupied COs into account for each k-point. In addition, we map the target fermionic Hamiltonians to the qubit ones with the technique reported in Ref. [26], which tapers off two qubits by exploiting $Z_{2}$ symmetries. Thus, the Si crystal is solved by a two-qubit system.

In Fig. 1, we show the quantum circuit used in the present study. It is motivated by the hardware-efficient ansatz [13], but we consider a simplified circuit with the minimal depth to reduce the computational cost and the effect of noise. The variational parameters are optimized by using the sequential minimal optimization [27] on a noiseless statevector simulator.

Quantum devices used in our demonstration are Rigetti’s Aspen-M processors[28]. In particular, we adopt qubits with high fidelities which are shown in Table 1 as examples. In the QSE calculations, we apply the readout-error mitigation (REM) [29, 30, 31, 32]. Moreover, to cancel out the bias coming from time-varying noise, we repeat the same calculation 40-70 times and average the outcomes (See also Supplemental Material [33]).

Results.— In Fig. 2, we show the quasiparticle band structure obtained by QSE. The result labeled as “VQE” represents the main one where the QSE calculations were performed on quantum devices. For comparisons, we include the results where the exact ground states were obtained using FCI or the complete active space configuration interaction (CASCI) with the same active space as VQE. For these two results, the QSE calculations were accurately performed using a statevector simulator. As seen in this figure, the VQE result agrees well with the corresponding FCI and CASCI ones, showing that the band structure is accurately obtained using actual quantum devices.

To apply REM, we constructed a calibration matrix before each calculation. More specifically, we prepared all the computational basis states and measured each probability distribution in the same basis (See also the Supplemental Material for the detailed definition of the calibration matrix). In Fig. 3 (a), we show all the records of the calibration matrices (diagonal elements) measured throughout the present study (December 20, 2022 to March 27, 2023), indicating that the noise is time-varying and its effect on outcomes also varies at every measurement. Figure 3 (b) shows a histogram of the calculated band energies at the $\Gamma$ point. Even though REM is applied, most of the results still deviate from the ideal noise-free ones obtained by the statevector simulator. Since the statistical error in each calculation was much smaller than the deviation (See also Supplemental Material [33]), the bias is mainly due to the effect of device noise. However, we found that averaging the outcomes at various times leads to error cancellations, which explains our well-reproduced ideal values.

Finally we address the accuracy of the VQE ground state. Figure 4 shows optimization processes of the VQE state $|\psi({\bm{\theta}})\rangle$ at the $\Gamma$ point. The energy optimized by the noiseless sampling simulator agrees with the exact ground-state energy obtained by CASCI within the chemical accuracy of $\pm 1$ kcal/mol (= $\pm 0.0434$ eV), showing that the ansatz we employ describes the ground state accurately. For this particular k-point, we also performed the optimization using actual quantum devices as a demonstration. As a result, we found that the optimization with REM is not sufficient to reach the exact ground-state energy. To overcome it, we also applied the zero-noise extrapolation (ZNE) method [34, 35, 36, 14] at the cost of introducing the redundant gates and measurements (See also Supplemental Material [33]). With this approach, we eventually observed a tendency of converging to the exact ground-sate energy.

Summary and Outlook.— To summarize, we have demonstrated $\it ab\ initio$ calculations of a quasiparticle band structure using actual quantum devices. The difficulty arises from the noise inherent to quantum devices. To circumvent it, a qubit-reduction strategy combining the pseudopotential, the active space approximation, and the symmetry-based qubit tapering technique were implemented. Owing to this strategy, we could solve the Si crystal with a two-qubit system. Moreover, the error-mitigation techniques were successfully applied and the ideal noise-free results were reproduced. Since the band structure is an essential physical characterization of the periodic materials, we believe that the present demonstration is a crucial step towards practical applications of quantum computers to solid-state materials.

Scaling up the quantum computations is certainly an indispensable subsequent step. The REM method used in the present study requires manipulations with probability vectors of size $2^{n}$, which limits the scalability ($n$ is the number of qubits). Problems on large systems might be tackled with scalable error-mitigation methods [31, 32, 37, 38]. In addition, a promising alternative approach is the application of recently developed hybrid quantum-classical algorithms which circumvents the VQE optimization difficulties. For example, a quantum algorithm combined with a quantum Monte Carlo was developed and demonstrated for the ground-state calculations of the diamond solid using 16-qubit systems [39]. In another example, a hybrid algorithm called quantum-selected configuration interaction was shown to calculate excited states as well as ground states [40]. These are interesting directions toward practical applications of quantum computers to larger systems.

Acknowledgements.— The authors thank Amazon Web Services for supporting this work through their Amazon Braket service. All computations in the present work are performed using the quantum-circuit simulation library Qulacs [41] with QURI Parts [42], a library for developing quantum algorithms. The Hamiltonians are constructed by the PySCF [43] and the OpenFermion [44] software.

Supplemental Material