A systematic variational approach to band theory in a quantum computer

Quantum computers promise to revolutionize our ability to simulate molecules, and cloud-based hardware is becoming increasingly accessible to a wide body of researchers. Algorithms such as Quantum Phase Estimation and the Variational Quantum Eigensolver are being actively developed and demonstrated in small systems. However, extremely limited qubit count and low fidelity seriously limit useful applications, especially in the crystalline phase, where compact orbital bases are difficult to develop. To address this difficulty, we present a hybrid quantum-classical algorithm to solve the band structure of any periodic system described by an adequate tight-binding model. We showcase our algorithm by computing the band structure of a simple-cubic crystal with one s and three p orbitals per site (a simple model for polonium) using simulators with increasingly realistic levels of noise and culminating with calculations on IBM quantum computers. Our results show that the algorithm is reliable in a low-noise device, functional with low precision on present-day noisy quantum computers, and displays a complexity that scales as Ω(M3) with the number M of tight-binding orbitals per unit-cell, similarly to its classical counterparts. Our simulations offer a new insight into the “quantum” mindset and demonstrate how the algorithms under active development today can be optimized in special cases, such as band structure calculations.


I. INTRODUCTION
The primary objective of computational chemistry is to calculate the eigenstates and eigenenergies of a chemical system.Increasing the accuracy of these calculations enables more accurate characterization of a wide range of chemical properties (ionization potential, equilibrium constants, IR spectra, etc.).However, the computational complexity of the most accurate algorithms scales exponentially with the number of basis orbitals, when using a classical computer.Therefore, large systems are typically treated approximately with a single-electron approximation, in which each electron independently interacts with an effective potential produced by all other electrons and atomic centers.Spin exchange forces can be treated through the Hartree-Fock self-consistent field method, and multi-body correlations can be treated through ad hoc correction terms as in density functional theory, but these methods fail when applied to highly-correlated systems.Therefore, computational chemists have begun to turn to quantum computers to reach higher levels of accuracy.
Quantum computers surpass Hartree-Fock by imposing fermionic statistics onto qubits and including electron correlation terms directly in the system Hamiltonian.Quantum Phase Estimation (QPE) is a quantum algorithm which extracts eigenstate energies from a simulated system on a noise-resilient quantum computer at a low cost of classical computational resources. 1,2However, quantum circuits designed for the QPE algorithm tend to require longer coherence times than are available in the era of Noisy Intermediate-Scale Quantum (NISQ) devices, so recent efforts focus on hybrid algorithms which balance quantum and classical resource costs.4][5][6] Many variants of VQE have arisen in recent literature, including methods such as Variational Quantum Deflation (VQD) capable of exploring excited states. 7However, resource costs for interesting systems still tend to exceed those currently available, and there is still some doubt whether error in hybrid algorithms can be made to converge sub-exponentially. 8,9riodic systems, such as the crystalline phase of a molecule or protein, are an especially interesting arena in the development of quantum algorithms.On the one hand, translational symmetry over the entire crystal seems to offer extremely powerful ways to reduce the computational complexity. 10On the other, periodic systems are typically considered infinite in extent and thus appear to require an exceedingly large number of resources to adequately approximate.For example, one may simulate a periodic system of N unit cells, each consisting of M orbitals (see for example Cade et al. 11 ), where M is typically comparable to the number of orbitals considered in a single molecular simulation, but N is large enough to approximate infinity.Alternatively, one may adopt a plane-wave basis, for which a quantum circuit is available to efficiently diagonalize the kinetic and potential operators directly. 12In either approach, the size of the basis, and therefore the number of qubits, must be very large to accurately represent the periodic system.The quantum resources required to simulate a crystal will thus tend to be many times larger than those required for a solitary molecule, and generally larger than the size of quantum computers available today.
In this work, we offer an easier transition to adapt quantum algorithms for periodic systems, by implementing correlation-free band structure calculations on NISQera quantum computers.Band structures are the funda-mental toolbox of materials scientists in the characterization and discovery of the electronic properties of crystalline solids.Band theory adopts the single-electron approximation (as in Hartree-Fock), for which a periodic Hamiltonian becomes separable in reciprocal space, reducing the system at any particular momentum k to the complexity of a single unit cell.In this way, the eigenstates of an electron with momentum k can be efficiently calculated in a classical computer; the energies of each eigenstate along a path varying k through reciprocal space form the band structure of the crystal.Integrating the band structure provides early insight into structural, electronic, optical, and thermal properties of the crystal. 13uantum algorithms to obtain classically verifiable results provide an invaluable tool in establishing foundational building blocks such as efficient mesaurement protocols and error mitigation, 14 as well as being generally easier to understand and replicated by researchers just breaking into quantum computation.With this motivation, we show how the single-electron approximation accommodates a systematic approach to efficiently apply VQE to solve the band structure of any periodic system.We will illustrate our procedure using an empirical tight-binding model for a simple cubic lattice, but our procedure is easily applied to any tight-binding Hamiltonian.For example, recent work in our group has demonstrated that "exact" tight-binding Hamiltonians can be derived from accurate electronic structure calculations using a procedure based on the projection of electronic wavefunctions on localized orbital bases (PAOFLOW). 15PAOFLOW is an advanced software tool that constructs tight-binding Hamiltonians from self-consistent electronic wavefunctions projected onto a set of atomic orbitals and provides numerous materials properties that otherwise would have to be calculated via approximate model approaches.Thus, an efficient quantum algorithm for the solution of the tightbinding Hamiltonian would provide an avenue for the calculation of advanced crystalline properties on a quantum computer.We do not expect our approach in this paper to offer immediate quantum advantage; rather, our purpose is to help chemists think the quantum way, motivating new, resource-efficient approaches to studying highly correlated systems.By considering the simplest available model, we can provide lower bounds on resource complexity and give insight into the practical difficulties chemists may expect when implementing quantum algorithms.
We have previously considered this topic, employing a VQE-based algorithm to iteratively calculate the band energies of a tight-binding silicon model. 16We now apply recent developments in the literature 7,17,18 to extend and improve upon our previous work.In Section II, we briefly outline the essential ideas we have taken from the literature.In Section III, we present our robust procedure for accurately calculating the band structure of any periodic system with a quantum computer.In Section IV, we demonstrate our procedure applied to a simple-cubic lattice, presenting data from a quantum simulator and preliminary results from IBM's open-access ibmq_athens and ibmq_santiago cloud devices.In Section V, we discuss the algorithmic complexity of our procedure and highlight the steps which may or may not be improved in later work.

II. BACKGROUND
In this section we briefly outline some essential techniques actively studied in the quantum computing literature.In particular, while QPE provides a robust strategy for measuring eigenenergies with minimal classical resources, its performance suffers greatly from the imperfect fidelity of NISQ devices.As such, we will focus mostly on VQE and its close cousin VQD when measuring eigenspectra, applying QPE when available as an optional refinement.

A. Quantum Phase Estimation
In the QPE algorithm, 1,2 a set of qubits (the "state register") are first prepared into an eigenstate |ψ of a unitary operator Û , such that Û |ψ = e 2πiφ |ψ .The unitary operator U is then repeatedly applied as a controlled quantum circuit so that the phase shift φ is encoded into another set of qubits (the "readout register").Measurements on the readout register give the binary expansion of φ.In molecular simulations one selects Û ≡ exp i Ĥτ , the operator which evolves a system with Hamiltonian Ĥ by a unit time τ .If one first transforms Ĥ to guarantee that all possible energies E fall within the interval [0, 2π/τ ), the measured phases φ map directly onto the eigenenergies of Ĥ.The algorithm is deterministic when the state register is prepared exactly into an eigenstate |ψ and its eigenvalue φ has an exact binary expansion, but it retains some probability of success when both conditions are relaxed.Thus, QPE can be adapted to discover eigenstates and eigenvalues a priori, at the cost of additional rounds of measurement.
Generally, Ĥ is given as a weighted sum of noncommuting Pauli words (Section III A).An exact circuit for Û ≡ exp i Ĥτ is not readily available, but it can be closely approximated by Suzuki-Trotter expansion, which factors Û into many small-time slices. 19The number of time slices scales polynomially with the accuracy required, and the depth of each time slice depends on the number of commuting groups in Ĥ.For this reason, QPE is extremely susceptible to errors arising from the low gate fidelity and short coherence times which plague NISQ devices.Alternative approaches scale more favorably with error at the cost of additional ancilla qubits. 20,21. Variational Quantum Eigensolver In the VQE algorithm, 3,4 one begins with a Hamiltonian Ĥ, represented as a weighted sum of non-commuting Pauli words, and a parameterized quantum circuit V (θ), the "ansatz", which prepares a set of qubits into an arbitrary state (Section III B).The ansatz is applied on an ensemble of states such that qubit measurements give the expectation value of each Pauli word in H.The weighted sum of expectation values gives the energy E(θ) of the arbitrary state prepared -this procedure is called "operator estimation". 4,22The parameters θ are varied by a classical optimization routine until E is minimized.According to the variational principle, this minimum is exactly the ground-state energy of the system when the ansatz V (θ) is robust (ie. it spans the full Hilbert space of the system), and if the classical optimization succeeds in producing the global minimum.
The algorithmic complexity of VQE depends on several factors.Measuring the expectation value of a Pauli word is a stochastic process, requiring a large number of measurements on the order of O( −2 ) for an acceptable sampling noise .These ensembles are usually measured for each Pauli word in Ĥ, although recent advances reduce the size of the ensemble by simultaneously measuring each commuting group of Pauli words 23 or by "classically shadowing" 24 the quantum circuit to require only a logarithmic number of measurements.Like in QPE, the efficiency of the algorithm is determined by the complexity of Ĥ, and every element of the ensemble requires a unique application of the ansatz, meaning that circuit depth and gate count should be kept minimal.The dimension of the ansatz also determines the efficiency and efficacy of the classical optimization.6][27] Because circuit depth and gate count are kept low, VQE is wellsuited to NISQ devices.

C. Variational Quantum Deflation
Variational Quantum Deflation (VQD) 7 is one approach for extending the VQE algorithm to explore excited states in addition to the ground-state.VQD begins as a typical VQE run to locate the ground state, and the ground-state parametrization θ 0 is recorded.The variational process is then repeated with an additional term in the optimization routine's cost function, which gives the overlap between the current ansatz and the ground-state, weighted by a factor β. States similar to the ground-state will be shifted into a higher effective energy, so that the optimization routine considers them unfavorable.Meanwhile, higher-energy eigenstates must be orthogonal to the ground-state, so their overlap contribution will be zero.Therefore, the next lowest energy that can be found is the first excited state.This process is repeated for each energy level, adding a new overlap term for each eigenenergy already found.Each overlap can be evaluated as the expectation value of a single commuting group of Pauli words in the Hamiltonian, so that the total number of additional measurements after finding M eigenvalues is Θ(M 2 ).If Ĥ consists of Ω(M ) commuting groups, measured for each of M energy levels, then the additional cost of the overlap circuits is negligible.

III. METHOD
Our objective is to calculate the band structure of a periodic system, as described by a tight-binding Hamiltonian Ĥ of the form: Each c † j and c j represent a creation and annihilation (ladder) operator on an atomic orbital φ j , centered on a coordinate r j in the crystal.The hopping parameters t αβ denote the energy cost of an electron transitioning from orbital φ β to orbital φ α .They are calculated from the overlap integrals between each pair of orbitals φ α and φ β , or they are selected to fit empirical observations.A general tight-binding Hamiltonian may also include multielectron correlation such as t αβγδ c † α c † β c γ c δ , but we neglect these terms in this work.
Our strategy is to transform Eq. 1 into reciprocal space and to apply VQD to solve for each eigenenergy at each momentum k along the desired path through reciprocal space.When sufficient quantum resources are available, we refine each band energy with QPE.Our procedure for mapping a single-electron periodic system onto a set of qubits is derived in Section III A. The variational ansatz we have selected, suitable for any band structure calculation, is described in Section III B. Details of implementing the quantum algorithm are presented in Section III C. Finally, we provide a step-by-step schematic of our algorithm and its relation to VQE in Fig. 2.

A. Qubit Mapping
The Hamiltonian in Eq. 1 consists of ladder operators acting on atomic orbitals.The Hamiltonians appearing in the quantum algorithms of Section II consist of Pauli words acting on qubits.We define a "Pauli word" Pi as an operator acting independently on each qubit with either the identity Î or one of the Pauli spin matrices X, Ŷ , Ẑ. Pauli words are a natural choice for representing physical operators in a quantum computer because their expectation values can be readily measured 16 and their unitary time evolution exp i Pi t can be readily implemented as a quantum circuit. 28Our goal in this section is to map our atomic orbitals onto a qubit basis, and our Hamiltonian to a weighted sum of Pauli words: The simplest conceivable mapping between orbitals and qubits is to identify each orbital with its own qubit.The qubit state |1 represents an occupied orbital, while |0 is empty.There are however an infinite number of orbitals in an infinite crystal, and quantum computers with an infinite number of qubits are beyond our engineering capabilities.We therefore reinterpret Ĥ as the Hamiltonian of an arbitrarily large supercell with periodic boundary conditions, consisting of N unit cells, each with M orbitals.
Hopping parameters are now dependent on the orbitals α, β and the displacement vector δ ≡ r ν α − r νβ between their atoms.As δ increases, t αβ tends to vanish, permitting a nearest-neighbor approximation in which one considers only a few of the smallest δ.
Eq. 3, when supplemented with two-electron correlation terms, is the form typically considered when applying quantum algorithms to periodic systems, requiring a total of M N qubits.In the single-electron approximation, however, we can reduce the size of the system to only M qubits by transforming into reciprocal space.Reciprocal space orbitals are characterized by their own ladder operators c † kj and ckj , related to c † νj and c νj by Fourier transform: Substituting Eqs. 4 into Eq.3, we obtain We simplify this sum by recalling r ν α = r νβ + δ.Then r ν α becomes a common factor of each k in the exponential, and we may exploit the orthogonality relation Each momentum k contributes an independent subystem with only M orbitals, whose eigenenergies may be solved independently.Classically, the values H αβ (k) in Eq. 7 form an M × M Hermitian matrix whose eigenvalues can be efficiently calculated with standard linear algebraic techniques in Θ(M 3 ) time.This work instead considers how to calculate these eigenvalues the "quantum" way.We focus on a specific Ĥk for the remainder of this section, with the understanding that our procedure must be repeated for each momentum k along the path of interest in reciprocal space.Eq. 6 has a form very similar to Eq. 1, except that it acts on reciprocal-space orbitals rather than atomic orbitals.We therefore adopt a "reciprocal-orbital" basis, in which each reciprocal-space orbital is identified with its own qubit.

Hamiltonian Mapping
Having transformed our Hamiltonian into reciprocal space (Eq.6), we must now consider mapping each ladder operator to a set of Pauli words.The ladder operators must satisfy the following: Meanwhile, the Pauli spin operators X, Ŷ , Ẑ act on a qubit's basis states in the following way: It is easy to verify that the following mapping suffices for a single qubit: In multi-electron systems, one typically adopts the Jordan-Wigner transformation, which retains the form of Eqs. 10 but appends a Ẑ operation on Θ(M ) other qubits to enforce fermionic antisymmetry.Alternatively, one may adopt the Bravyi-Kitaev transformation, which requires operations on only Θ(log M ) qubits, but uses a non-intuitive basis and involves non-adjacent interactions more difficult to simulate on certain qubit architectures.
We refer the reader to Seeley et al. 28 for an excellent introduction to both transforms.However, because we are considering single-electron systems, there are no other fermions to exchange with, and we may use Eqs. 10 directly, so that each ladder operator acts on only Θ(1) qubits.This significantly reduces the complexity of our Hamiltonian, as we shall presently see.
We may rewrite Eq. 6 to exploit the Hermiticity of Ĥk .
since the transpose term Eq. 12 provides the weighted sum of Pauli words required in the quantum algorithms of Section II.Eq. 12 consists of Θ(M 2 ) Pauli words.The complexity of each algorithm in Section II is determined in part by the number of commuting groups in Ĥ.In Eq. 12, all terms of the form Ẑα , Xα Xβ , and Ŷα Ŷβ each form commutative groups.Therefore, when Ĥk has no imaginary part, the energy can be determined with just 3 rounds of measurement.When Ĥk does have an imaginary part, we note that for fixed α, Ŷα Xβ>α and Xα Ŷβ>α each form commutative groups, so in general we have Θ(M ) commuting groups.Finally, we note that each of these commuting groups are qubit-wise commutative, meaning that each index of all Pauli words in the set has either the same spin operator or the identity.This accommodates a particularly simple procedure for measuring expectation values of each set simultaneously, requiring no additional overhead in the measurement circuit.

B. Ansatz
The VQE and VQD algorithms require an ansatza parameterized quantum circuit V (θ) preparing a trial state Ψ(θ) ≡ V (θ) |0 for energy measurements.A quantum circuit to span the full Hilbert space of M qubits requires 2(2 M − 1) parameters, and it will not generally have an efficient decomposition into one-and two-qubit gates.However, most applications to molecular simulation consider a system with fixed number of electrons.In the orbital basis, or in our reciprocal orbital basis, one need only consider that subset of Hilbert space spanned by the basis states whose Hamming weights match the number of electrons in our system.For example, in band structure calculations we consider just one electron, so we need only consider the space spanned by |10... , |010... , etc.. Gard et al. 17 provide a procedure for generating variational ansatz which conserve particle number, which is particularly simple when particle number is 1.We begin with M qubits labeled 0 through M − 1 in the state |0 .First, we apply an X gate to qubit 0, to set our ansatz with a single filled orbital.Then we apply the entangling parameterized A gate 17 such that each qubit is entangled directly or indirectly with qubit 0 (see Fig. 1).This ansatz requires M − 1 A gates, each contributing two independent parameters, for a total of Θ(M ) gates and parameters.The circuit is compatible with any quantum architecture exhibiting linear qubit connectivity and has a depth of Θ(M ).Alternatively, in a fully connected device, the A gates could be applied with a "divide-and-conquer" strategy, reducing the circuit depth to Ω(log M ).
Rather than assigning each orbital to its own qubit, we could assign each orbital to an individual basis state, requiring only Θ(log M ) qubits total.This "compact basis" is the approach of our previous work. 16While this is more efficient in the number of qubits, it must explore states with an arbitrary Hamming weight.The number of parameters required to span the space of interest is unchanged, and a suitable ansatz must be developed ad hoc.Additionally, the Hamiltonian for the compact basis consists of global operators acting on all qubits at once, and it would generally form the maximum number 3 log 2 M = M log 2 3 of commuting sets, requiring a less efficient measurement protocol.Reliance on a random ansatz and a global cost function made our previous procedure vulnerable to exponentially difficult optimizations induced by barren plateaus. 8,29Finally, the Hamiltonian for the compact basis is less-structured and more difficult to reduce based on symmetries in the Hamiltonian (for example our observation in Section III A that a real Ĥk results in Θ(1) rounds of measurement).

C. Band Structure Calculation
With our ansatz V (θ) (Fig. 1) and qubit Hamiltonian Ĥk (Eqs.7 and 12) prepared, we are ready to implement VQD for each momentum k along a path through reciprocal space.This path is usually constructed from highsymmetry segments in the crystal's First Brillouin Zone, because this proves sufficient to calculate many properties of interest.As briefly described in Section II, the idea is to vary the trial state prepared by our ansatz until the energy E ≡ Ĥ is minimized.We repeat the optimization for each band energy, adding additional terms to the cost function proportional to the overlap between the trial state and each previously found eigenstate, weighted by the constant factor β.
The expectation values Ĥ of a generic observable cannot be directly measured in the quantum computer.
The particle-number preserving A(θ, φ) gate from Gard et al. 17 αβ as input and outputs each band energy E l (θ).Optionally, each band energy may be refined with QPE.The operator Ω0 is the sum of all Pauli words spelled with letters Î and Ẑ.The operator V is the quantum circuit presented in Fig. 1.
Rather, the expectation value of each Pauli word Pi are measured independently, and the energy is evaluated from the weighted sum Ĥ = i a i Pi , with weights a i taken from Eq. 12. Obtaining the Pauli expectation values Pi is also somewhat indirect.First, the Pauli word Pi should be transformed so that it contains only letters Î or Ẑ -let us refer to the modified Pauli word as Qi .In practice, the transformation is easily accomplished by applying a "basis rotation" gate to each qubit before measurement.Next, each qubit is measured to be in one of the two computational basis states |0 or |1 .The bitstring obtained from concatenating the state of each qubit is itself an eigenstate of Qi , with eigenvalue +1 or -1.This procedure is applied to a large ensemble of qubits, each prepared independently with the ansatz and basis rotation gates.The expectation value Pi is the average of all the eigenvalues of Qi measured across the ensemble.
The ensemble necessarily has a finite size S, introducing an energy variance on the order of 2 ∈ O(1/S).In practice, the ensemble is usually prepared in sequence, resetting a single register of qubits after each round of measurement, relegating the sampling error a parameter in the time complexity of any VQE-based algorithm.Fortunately, the same ensemble may be used to calculate the expectation values of any Pauli word which is qubit-wise commutative with Pi .For simplicity, we assign S = 8096 for each commuting group in this work, although advanced methods exist which optimally distribute measurements to minimize the sampling error . 22any popular optimization routines (eg.SLSQP, BFGS) are gradient-based, and they have difficulty converging to the correct value in the presence of sampling noise.Therefore, we use COBYLA, a simplex-based algorithm implemented in the scipy python package, which we have empirically noted to give good results.We randomly generate our initial guess for the parameters θ, and we use the default tolerance parameters implemented by scipy.These choices are by far the simplest, but they are by no means optimal, and our results may be improved greatly by a more careful choice of optimization routine. 30efore we can implement the deflation procedure, we must select the constant β suitable for "deflating" each band energy.We do this with a systematic procedure, first maximizing the energy of our system to find the highest possible energy E max .We then minimize the energy to find E 0 and ∆ ≡ E max − E 0 .In theory, β = ∆ is a sufficiently high number to guarantee each eigenstate is projected sufficiently out of the optimization in later steps.In practice, we take β = 2∆ to insure against errors in the sampling and optimization process.
Higgott et al. 7 offer several strategies for computing the overlap, offering robustness against error at the cost of ancilla qubits or additional optimization steps.In this work, we choose the simplest, evaluating the overlap with an eigenstate Ψ(θ l ) by preparing the trial state Ψ(θ) and applying the adjoint circuit V † l ≡ V † (θ l ).The probability of measuring the bitstring 0..0 gives the overlap | Ψ(θ)|Ψ(θ l ) | 2 .In practice, the probability of measuring bitstring 0..0 is equivalent to the expectation value of an operator Ω 0 ≡ i Qi , the sum of all unique Pauli words spelled with the letters Î and Ẑ (eg.Î Î Î, Î Î Ẑ, ... Ẑ Ẑ Ẑ).All such operators are qubit-wise commutative and can be estimated with a single round of measurements.Therefore, we can implement the deflation procedure conveniently in the qiskit Python package provided by IBM, by solving for each band energy and then adding to our Hamiltonian the deflation operator β Vl Ω0 V † l .Initializing the Hamiltonian Ĥ0 ≡ Ĥk , our procedure can be formally summarized as follows: Each E l we find is recorded as the energy of the lth band at momentum k, and we repeat the procedure for each k in our selected path.Optimization routines do not always converge to the true minimum, and errors incurred early in the deflation procedure can propagate unfavorably to higher bands.Therefore, we include an optional QPE refinement to our algorithm, which applies QPE to each optimized state Ψ l ≡ Vl |0 .QPE has the effect of selecting the dominant eigenstate of Ψ l and giving the corresponding eigenenergy with high precision.Thus, as long as the optimization procedure is "good enough", we may update our energy calculations with the result of the QPE experiment.We have used the iterative version of QPE implemented in qiskit.Details of the algorithm can be found in Dobsícek et al. 2

IV. RESULTS
To demonstrate our procedure, we consider a basic model for a crystal in a simple cubic lattic structure (see Fig. 3).Each atom has s, p x , p y , and p z orbitals (M = 4), with a large energy gap between the s and p orbitals, and comparatively small hopping parameters between neighboring s and p orbitals.This is the simplest threedimensional periodic system accommodating multiple orbitals per atom, and the required number of qubits (one qubit per orbital, plus one ancilla qubit to implement QPE) fits nicely onto IBM's open-access five-qubit machines.It may also be considered a rough model for elemental Polonium, although more accurate models should take into account relativistic effects and Coulomb interaction between orbitals located on the same atom. 31,32he exact eigenenergies of our model at specific k-points along a high-symmetry path are calculated using standard linear algebraic techniques to diagonalize the matrix elements in Eq. 7. We compare this band structure to the results from the quantum algorithm presented in Section III in four different levels of simulation: 1. Statevector -quantum operations are simulated with unitary matrices, and expectation values are calculated exactly.
2. Sampling -expectation values are now calculated by sampling from a probability distribution.
3. Noisy -quantum operations and measurements are now applied with an error rate drawn from real quantum devices.
4. Calibrated -the same noisy simulator is used, but classical post-processing steps are applied to mitigate the error.
We also present preliminary results from IBM quantum devices.

A. Statevector Simulator
To validate our algorithm's capability of producing the correct band structure, we model the state of an n qubit system as a complex statevector and quantum operations as unitary matrices acting on the Hilbert space spanned by the 2 n dimensional basis vectors.Expectation values are evaluated analytically.Such a simulation gives the ideal behavior of a quantum computer, with perfect qubit fidelity and no sampling variance.Fig. 3 summarizes the results of over 32 randomly-seeded optimization runs, marking the median optimization with a diamond and the interquartile range with a bar.For a few momenta where the third and fourth bands are very close together, optimization tends to locate the wrong eigenstate, giving a small variance in results.For every other point, the diamond coincides perfectly with the classical solution and the bar is absent, demonstrating that our ansatz is robust, the deflation procedure is mathematically sound, and that our choice of optimization routine (COBYLA) is generally consistent in converging to the correct values on a smooth surface.

B. Sampling Simulator
We now consider long-term viability of our procedure by retaining perfect qubit fidelity but simulating realistic measurement.The same unitary matrices as in the statevector simulator are applied to an ensemble of states, which are "measured" by sampling from the resulting probability distribution a finite number of times.While mathematically equivalent, the sampling noise resulting from the stochastic measurement process can make the energy surface bumpy, which can have a detrimental effect on the optimization step.We have selected the COBYLA optimization algorithm because it is resistant to these bumps; nevertheless, the anomalous variance observed at nearly degenerate points in the statevector simulator is now commonplace.Fig. 4 shows our results on the noiseless qubit simulator over 32 randomly-seeded runs, clearly marking the median (asterisk) and mean (X) for both optimization (left) and QPE refinement (right).FIG.4: Sampling Simulator -Our method applied in the presence of sampling noise (high-fidelity qubits).The left column shows raw optimization results; the right column shows the energy obtained by QPE refinement.Gray dots denote the results from each of 32 trials, with each band given its own row.The asterisk and X denote the median and mean, respectively.
The smaller dots denote the results of individual trials.
The optimization results are extremely accurate and precise on the high-symmetry momenta but deviate slightly on the intermediate points.In fact, the highsymmetry points in this particular model each happen to have matrices H αβ (k) (Eq.7) which are already diagonalized, and the resulting cost function yields a wellbehaved surface which is reliably optimized, even in the presence of noise.Averaged results on the intermediate points still tend to be quite good, but individual trials can exhibit a large variance.However, the optimization does succeed in finding a point close enough to a correct eigenstate that the QPE refinement consistently extracts the dominant eigenvalue with high precision.The median QPE results prove to be as accurate as is permitted by the finite binary expansion calculated by the algorithm.

C. Noisy Simulator
We now consider the realistic application of our procedure on present-day quantum computers, which suffer from relatively short coherence times and are vulnerable to a number of error sources.This makes practical computations extremely difficult, even in systems requiring relatively few qubits.We model environmental effects by simulating a dephasing channel, randomly introducing a bit-flip and/or phase-flip on each qubit after each unitary operation 33 .We also model errors in the measurement process by randomly introducing a bit-flip with low probability prior to sampling each probability distribution.Fig. 5 shows our results on a simulator emulating the error rates characteristic of IBM's ibmq_athens quantum computer.Qubit noise has a clearly negative impact on the quality of results.Lowest-band optimization results tend to suffer a large systematic shift, characteristic of coherent noise in a quantum computer.Additionally, while average QPE refinement often improves energy es-FIG.5: Noisy Simulator -Our method applied while simulating low-fidelity qubits, without calibration.The left column shows raw optimization results; the right column shows the energy obtained by QPE refinement.Gray dots denote the results from each of 32 trials, with each band given its own row.The asterisk and X denote the median and mean, respectively.
timates, its results are now clustered with some variance around each nearby band, and on occasion (eg. the third band between M and Γ) the mean optimization result is more accurate.This is symptomatic of the long circuit requirements for QPE, and supports the widely-believed notion that variational algorithms are better suited to NISQ devices.

D. Calibrated Simulator
Although automated error correction procedures, based on redundant qubit registers and applied during calculation, are the most promising path toward practical quantum computation, several classical post-processing methods have already proven successful in mitigating error.Errors modeled by bit-flips in the measurement process ("readout error") can be mitigated in part or entirely, at the cost of additional calibration circuits 34 .Errors modeled by a dephasing channel as each unitary oper-ation is applied ("gate error") tend to result in systematic distortions of the energy surface, as we have seen in the optimization results of Fig. 5.These distortions can be mitigated by applying Zero-Noise Extrapolation (ZNE) 18 in which the same measurements are repeated several times, each time modifying the circuit to incur more noise.These measurements can then be extrapolated to a hypothetical circuit with zero noise, using Richardson extrapolation or a similar method.Fig. 6 shows our results on a noisy simulator, applying readout calibration and ZNE for each energy evaluation during the optimization.ZNE offers noticeable improvement in the highest and lowest bands (calculated independently), but appears less impactful on the intermediate bands (calculated after deflation), perhaps even increasing variance in the third band.This may be explained by noting that ZNE is designed to assuage systematic error, and this is what we tend to observe when we can rely on the variatonal principle, where energies cannot in principle be measured below the ground-state FIG.6: Calibrated Simulator -Our method applied while simulating low-fidelity qubits, along with rudimentary calibration.The left column shows raw optimization results; the right column shows the energy obtained by QPE refinement.Gray dots denote the results from each of 32 trials, with each band given its own row.The asterisk and X denote the median and mean, respectively.The squares and diamonds on the left denote the energies measured on quantum devices ibmq_santiago and ibmq_athens respectively, using the least-error optimization results obtained with the calibrated simulation data.Device architecture constrains the length of quantum circuits, and QPE results for either device could not be obtained.
energy.This is not always true because our energy estimates are linear combinations of stochastically evaluated Pauli expectation values, and on occasion we do observe trials which appear above the highest band, but these points are relatively rare, and the average values on the highest and lowest bands are shifted inwards.However, the deflation circuits V0 are somewhat different for each trial, depending on exactly what eigenstate was selected for the lowest band, and this, coupled with the coherent qubit error, has the effect of inducing a random noise on the intermediate bands.This explains why the intermediate bands seem to suffer a larger variance but reasonable average values.We note that many other error mitigation techniques besides readout calibration and ZNE have been proposed in the literature, and our results can likely be improved greatly by implementing more of them.Nevertheless, the best solution to combat random error remains averaging over more and more trials.
In addition to statistics from a calibrated simulator, Fig. 6 also shows data from the IBM devices ibmq_athens and ibmq_santiago.These are calibrated energy measurements of the eigenstates given by the least-error optimization runs on the (calibrated) noisy simulator.Results are generally consistent with the simulator, but our error mitigation is evidently even less effective on the real devices, and we note that ibmq_santiago is especially ill-behaved at certain points.This may be due to less effective thermal isolation from its environment at the hardware level.Finally, implementing the controlled-unitary operations necessary for the QPE procedure on a linear architecture introduces an overwhelming amount of overhead in the form of additional SWAP gates, making the QPE refinement part of our algorithm completely intractable on these devices.

V. DISCUSSION
We have presented an application of VQD to calculate the band structure of a periodic system.This algorithm is hypothetically successful in producing accurate results on a device with low noise and is functional to a limited extent on current NISQ devices.In this section, we carefully analyze the complexity of the algorithm.The classical approach to band structure includes up to Eq 7, at which point the calculated values H αβ (k) are arranged into a Hermitian matrix.This matrix can be diagonalized using row-reduction or a similar technique in Θ(M 3 ) steps, where M is the number of atomic orbitals per unit cell.This is the standard against which we must compare our quantum algorithm.
Quantum resources are employed in the VQD phase of our algorithm during the operator estimation procedure, for every evaluation of the energy E ≡ Ĥk .Each application of the ansatz from Fig. 1 requires M qubits, Θ(M ) entangling gates, and has a depth between Θ(log M ) and Θ(M ) layers, depending on qubit architecture.The Hamiltonian in Eq. 12 has Θ(M ) commuting groups, even including additional terms from the deflation procedure (Eq.16).Our implementation requires an ensemble size of O( −2 ) for each commuting group in Ĥ to obtain an expectation value accurate within , but since does not scale with M , we omit it in the present analysis.The ensemble states may be prepared sequentially, for a worst-case (linear architecture) execution time on the order of Θ(M 2 ).Alternatively, the ensemble states may be prepared in parallel, decreasing execution time at the cost of additional qubits.In the best case, implementing "classical shadowing" 24 reduces the number of required measurements to Θ(log M ), and a fully-connected architecture permits a circuit depth as low as Θ(log M ), bringing our algorithm into a sub-polynomial quantum resource requirement.However, the operator estimation procedure is still bounded by the number of Pauli words Θ(M 2 ) when measurement results are assembled into the energy E Operator estimation is repeated for each function evaluation in the optimization procedure.The number of function evaluations required depends on the optimization routine selected and the shape of the energy surface, so it is difficult to estimate.Quantum variational algorithms are notoriously vulnerable to so-called barren plateaus, regions in the parameter surface with a vanishing gradient expected to result in an exponential number of function evaluations for successful convergence. 8,29However, research into barren plateaus has focused mainly on densely-packed ansatze which alternate between parameterized rotation gates and entanglement among each qubit.The ansatz we have presented has a more constrained structure which rotates and entangles only two qubits at a time.Additionally, Cerezo et al. 8 found that local cost functions such as the Hamiltonian we have employed are much more resistant to barren plateaus compared to global cost functions.Therefore, we conjecture that the number of function evaluations required in our algorithm may be expected to scale polynomially with the number of ansatz parameters, in our case Θ(M ), provided sufficient error mitigation to suppress noise-induced barren plateaus, which are independent of the ansatz 35 .Thus, we include a factor of Θ(M c ), where c ≥ 1 depends on the optimization.The optimization is repeated for each of M energy levels; therefore, the VQD phase of our algorithm has a total run-time on the order of Θ(M 3+c ).
An optional QPE phase may be implemented to estimate the eigenvalue to an arbitrary binary precision t. 2 The implementation of QPE we have used requires M + 1 qubits and Ω(t) rounds of measurement (see Dobšíček et al. 2 for a tighter bound).Each round applies a quantum circuit approximating a unitary operator Ûj = exp i Ĥτ j .Each time slice in the Suzuki-Trotter expansion of Ûj on a linear architecture requires an entangling gate count of Θ(M 3 ) and a circuit depth of Θ(M 2 ). 36QPE is repeated for each of M energy levels, setting the best case run-time of the QPE phase of our algorithm on the order of Θ(M 3 ).Note also that the simulation time τ j scales exponentially with the accuracy of the phase estimation procedure, and the number of time-slices must scale accordingly to maintain an accurate Ûj .Thus, QPE tends to incur too much overhead for practical application on present-day NISQ devices.
Altogether, evaluating the band energies for each momentum k requires Ω(M 3 ) time steps, comparable to the classical approach.Even with a "perfect" optimizer in which the optimal parameters θ l are produced instantly (c = 0), the complexities of operator estimation and QPE alone exhibit the same scaling as classical diagonalization and incur significantly greater overhead from the finite accuracy .While in this form band structure calculations are not a strong candidate for quantum advantage, quantum computers are expected to provide a superior edge when including electron correlation terms such as t αβγδ c † α c † β c γ c δ in the Hamiltonian, which introduce factors of exponential complexity in the classical approach.However, such terms also appear to force us to abandon several simplifications we have made.First, transforming into reciprocal space no longer enables H to be separated into subsytems of size M , meaning many more qubits are required to accurately simulate a periodic system.Second, considering multiple electrons forces us to adopt a qubit mapping which enforces fermionic antisymmetry, greatly increasing the number of commuting groups in the Hamiltonian.Finally, our ansatz dimension, entangling gates, and circuit depth can no longer remain linear in the number of qubits while simultaneously remaining robust.Our hope is that this work will inspire similar simplifications to those we have made here, while remaining applicable to highly-correlated systems.

VI. CONCLUSION
In this work, we have presented a systematic algorithm for evaluating band structures on a quantum computer.We have demonstrated the viability of implementing this algorithm in noiseless qubit systems, and we have demonstrated several of the difficulties faced when implementing it on present-day NISQ devices.Given the analogy to the classical band structure problem, our algorithm evidently generalizes to solving the eigenvalues of any Hermitian matrix.Finally we have demonstrated how state-of-the-art quantum algorithms can be applied with drastically lower resource requirements to correlation-free crystalline systems and motivated similar approaches for highly-correlated systems less accessible to classical com-puting.

CONFLICTS OF INTEREST
There are no conflicts to declare.

FIG. 1 :
FIG.1:The ansatz V (θ) suitable for any band structure calculation.Each qubit is initialized in the |0 state; the output is an arbitrary superposition of states with a single qubit in the |1 state.

FIG. 3 :
FIG. 3: Statevector Simulator -The band structure of a simple cubic lattice with s and p orbitals (right inset) along the high-symmetry path XMΓ through the lattice's First Brillouin Zone (left inset).Solid curves denote classical (exact) diagonalization.Diamonds denote the median optimization result from applying our method on a noiseless statevector simulator 32 times with a different random seed.Bars (only visible between the third and fourth bands at nearly-degenerate momenta) denote interquartile ranges.Hopping parameters are 2 eV between adjacent s and p orbitals and 2 eV between colinear p orbitals.Each s orbital has a self-energy of -14 eV to generate a large band gap.