Quantum algorithm for band structures with local tight-binding orbitals

While the main thrust of quantum computing research in materials science is to accurately measure the classically intractable electron correlation effects due to Coulomb repulsion, designing optimal quantum algorithms for simpler problems with well-understood solutions is a useful tactic to advance our quantum “toolbox”. With this in mind, we consider the quantum calculation of a periodic system’s single-electron band structure over a path through reciprocal space. Previous efforts have used the Variational Quantum Eigensolver algorithm to solve the energy of each band, which involves numerically optimizing the parameters of a variational quantum circuit to minimize a cost function, constructed as the expectation value of a Hamiltonian operator. Traditionally, a unique Hamiltonian operator is constructed for each k-point, so that many cost functions, each with their own parameter space, must be optimized to generate a single band. Similarly, calculating higher bands than the ﬁrst has traditionally involved modifying the cost function with additional overlap terms to ensure higher-energy eigenstates are orthogonal to those of lower bands. In this paper, we adopt a direct space approach, using a novel hybrid ﬁrst/second-quantized qubit mapping which allows us to construct a single Hamiltonian, and a single cost-function, suitable for solving the entire band-structure. In contrast to previous approaches, the k-point and the band index are selected by additional parameters in our quantum circuit, rather than through modiﬁcations to the cost function. The result is a technically and conceptually simpler approach to band structure calculations on a quantum computer. Moreover, we expect that the tools developed herein will motivate new strategies for tackling highly-correlated materials beyond the grasp of classical computing.


Introduction
Band structures plot the energy eigenstates of an electron in the presence of a periodic potential as a function of momentum. Integrating or differentiating a band structure yields many useful properties of solid-state materials. Band structures are easily calculated classically under the single-electron approximation, but extending the approach to highly correlated materials has been challenging. This obstacle has motivated many materials scientists to consider quantum computation, which enables extensions to the the band structure approach which more accurately account for electron-electron interactions. The literature is full of distinct algorithms which treat periodic systems, including hybrid quantum-classical DMFT, 1-4 cyclic qubit arrays, 5 and "holographic" VQE. 6 Other methods developed for quantum chemistry are also readily applied to periodic systems expanded in a basis of plane waves 7 or Bloch atomic orbitals. [8][9][10] We refer the reader to Bauer et al. 11 for a thorough review of the most popular approaches.
Several recent papers 9,12,13 have demonstrated band structure calculations using the Variational Quantum Eigensolver (VQE) algorithm. 14-16 VQE is a popular quantum algorithm in the Noisy Intermediate-Scale Quantum (NISQ) era which consists of a low-depth parameterized quantum circuit preparing a variational ansatz |Ψ(θ )⟩, and a Hamiltonian operator H whose expectation value can be measured in the quantum computer. The variational parameters are adjusted until the energy E ≡ ⟨Ψ|Ĥ|Ψ⟩ is minimized. This energy is a good approximation of the ground state of the system described byĤ, provided that the ansatz is sufficiently expressive. Many ansatze have been developed for electronic structure calculations, balancing expressivity with efficient implementation on available quantum hardware. 10,[17][18][19][20][21] In band structure calculations, the ansatz |Ψ(θ )⟩ is constrained to a single-electron wavefunction, and the ground state is the energy of the lowest band at a particular momentum k. Higher bands are found by repeating the optimization with additional constraints to ensure the ansatz is orthogonal to previously located eigenstates.
In order to estimate the expectation value of the HamiltonianĤ in a quantum computer,Ĥ should be expressed as a linear combination of "Pauli words": Each Pauli wordP j = Î ,X,Ŷ ,Ẑ ⊗n is an n-qubit operator with a Pauli spin matrix associated with each qubit. The precise choice of c j depends on the choice of orbital basis, and how this basis is mapped onto the available logical qubits. Previous results for solving band structures each used somewhat different strategies for the qubit mapping, but all three adopted a Bloch atomic orbital basis. 9,12,13 Each Bloch atomic orbital |α k ⟩ is related to the local atomic orbitals {|α r ⟩} at each lattice point r in the crystal by a Fourier transformation: Bloch atomic orbitals form a valuable basis for band structure calculations because the basis states are intrinsically periodic, and any wavefunction constructed in this basis automatically satisfies the periodicity of the system. Furthermore, the single-electron Hamiltonian is separable in k, meaning that each set of orbitals |α k ⟩ for fixed k can be solved independently, reducing the effective size of the eigenvalue problem. However, this comes with the consequence that a new set of c j must be obtained for every k. Similarly, while previous results used different strategies for exploring higher-level bands, they each enforced orthogonality with previously located eigenstates by including additional terms in the cost function. This effectively generates a new set of c j for each band, and often increases the number of quantum measurements required to obtain each E.
In this paper, we present an algorithm which simplifies band structure calculations in a quantum computer by setting the cost function just once, for all k and for all bands. We accomplish this by adopting the basis of local atomic orbitals, and by enforcing the periodicity of our eigenstates directly through the ansatz. In this framework, k enters as an additional (fixed) parameter to our variational quantum circuit. Additionally, we present a novel procedure for exploring higher bands without the use of additional terms in the cost function. Instead, orthogonality with previously located eigenstates is enforced by the quantum circuit itself, incurring negligible overhead in circuit complexity and no overhead in the required number of measurements.
The Qubit Mapping, Quantum Circuit, and Cost Function sections present the technical details of our algorithm. The Examples section briefly presents examples using simulated results to validate our method. In the Conclusion, we discuss the advantages and disadvantages of our approach, and we suggest how this method could be adapted to efficiently treat highly-correlated systems.

Qubit Mapping
In this paper, we adopt a local atomic orbital basis. We associate with each unit cell a lattice coordinate r and a set of M orbitals {|α⟩}, such as the hydrogen-like orbitals |s⟩, |p x ⟩, etc. centered on each atom. The set {|α r ⟩} of all orbitals over all unit cells forms the basis for the crystal. In practice, we restrict ourselves to a supercell consisting of N total unit cells. The value of N determines the resolution in k we can obtain in our band structure, where the limit N → ∞ corresponds to a continuous k space. We index each unit cell with an integer coordinate ν counting from 0 to N, thereby adopting the finite basis {|α ν ⟩} with MN elements. We use the bold font to indicate that in a d-dimensional crystal, ν takes the form of a d-tuple.
We map our supercell onto a set of qubits using a novel hybrid first/second-quantized approach, factoring the basis orbital |α ν ⟩ into two sub-states: The |α⟩ state is mapped onto an "orbital register" M consisting of M qubits encoded with second quantization. Each local atomic orbital α is associated with a specific qubit q α ∈ M ; the state |α⟩ corresponds to the computational basis state in which q α is |1⟩ and all other qubits are |0⟩. The |ν⟩ state is mapped onto a "site register" N consisting of log N qubits encoded with first quantization. Each basis state of N is the binary representation of the integer coordinate ν; a d-dimensional crystal will have d sub-registers in N . The total number of qubits required by this mapping is M + log N. Fig. 1 shows a schematic of our qubit mapping.

Quantum Circuit
A general ansatz expanded in the local atomic orbital basis {|α r ⟩} may be written as: Bloch's theorem guarantees that single-electron eigenstates of a periodic system must themselves be periodic, with a phase factor determined by the momentum of the electron k.
φ αr = e ik·r φ α0 (5) Figure 1. A schematic illustrating the hybrid qubit mapping used in this paper. Each unit cell in a supercell is indexed with the "site register" using first quantization. Meanwhile, each qubit in the "orbital register" is associated with a single orbital in each unit cell using second quantization. For example, in this schematic, the qubit labelled q p x holds the amplitude for every p x orbital in the supercell, while the site register contributes the phase for each site. Figure 2. The quantum circuit applied on our site register N . The sub-circuit P loads the computational basis state |p⟩. The sub-circuit QFT implements the Quantum Fourier Transform on a linear architecture.
For ease of notation, let us assume a one-dimensional crystal with lattice constant a, so that each lattice coordinate r → aν for an integer coordinate ν. Imposing periodic boundary conditions on a supercell of size N requires φ α,aN = φ α,0 , implying for a momentum quantum number p. Under the qubit mapping described above, we can combine Eqs. 3, 4, 5, and 6: The first factor describes the amplitudes of each orbital on the principal unit cell and is expressed in the orbital register M , while the second factor describes the phases accumulated on each site and is expressed in the site register N . We have absorbed a factor of √ N from the site register into the amplitudes φ α ≡ √ Nφ α,0 so that both registers are normalized. The action on both registers is factorized and thus can be expressed with two independent quantum circuits.

Site Register
We wish to prepare the following state in the site register N , from the starting state |0⟩: This state is determined exactly by the momentum quantum number p, requiring no other variational parameters. Eq. 8 has the form of a discrete Fourier transform, and we can make use of the well-known Quantum Fourier Transform (QFT), an efficient implementation of the discrete Fourier transform as a quantum circuit. 22 We first prepare the state |p⟩ with a single layer of Pauli X gates applied to each bit for which the binary representation of p is |1⟩. We then apply QFT, preparing the state in Eq. 8 exactly. Fig. 2 presents an implementation of QFT optimized for linear qubit architectures, requiring Θ((log N) 2 ) entangling gates and Θ(log N) depth.

Orbital Register
We wish to prepare the following state in the orbital register M , from the starting state |0⟩: The amplitudes φ α are a priori unknown, and we will treat them as variational parameters. Our objective is to design a parameterized quantum circuit capable of expressing any arbitrary superposition of the basis states |α⟩, which are those M-qubit computational basis states with a Hamming weight of 1. The circuit V presented in our previously-published work 13 and re-presented in Fig. 3a is suitable for locating the lowest band. It consists of M − 1 so-called A gates (Fig. 3b), first presented in Gard et al., 21 each of which accepts a polar angle θ and an azimuthal angle ϕ, for a total of 2(M − 1) independent parameters. If desired, the ansatz may be constrained to real values by setting all azimuthal angles to 0.
The variational circuit applied on our orbital register M . We now wish to consider how to adapt V for higher bands, in a way which ensures orthogonality with previously located eigenstates. Let us define V 0 as that circuit V with parameters optimized to prepare the eigenstate |ψ 0 ⟩ of the lowest band, so that V 0 |10..0⟩ = |ψ 0 ⟩. If V 0 is applied to any other computational basis state, the result must be orthogonal to |ψ 0 ⟩. Therefore, if we prepare an arbitrary single-electron state |Ψ⟩ over all qubits but the first (ie. |Ψ⟩ = |0⟩ ⊗ |Ψ + ⟩ for an arbitrary single-electron state |Ψ + ⟩ with M − 1 orbitals), the state V 0 |Ψ⟩ will be orthogonal to the ground-state |ψ 0 ⟩ (Fig. 4). Thus, we may use the same variational form V for each band l, except that we omit qubit q l after each iteration, and we apply each previously optimized circuit V l .
The entire circuit for the orbital register can be compactly represented as in  Figure 4. The variational circuit applied on our orbital register M to locate the second-lowest band of any periodic system. The V 0 sub-circuit is as shown in Fig. 3a. The V 1 circuit is analogously constructed, save with one fewer qubit. Figure 5. The generalized variational quantum circuit applied on our orbital register M to calculate any band of a periodic system. Each diagonal corresponds to a particular energy level. Parameters for bands lower than the target energy level are fixed at their optimal values, while parameter for bands higher than the target energy level are fixed at zero. Thus, only the parameters for the target band must be optimized.

Cost Function
Our objective in this section is to present a cost-function suitable for band structure calculations over all k and all bands. We first write the tight-binding HamiltonianĤ in the basis of local atomic orbitals |α ν ⟩ described above, and map it onto a set of Pauli words so thatĤ → ∑ j c jPj . Ostensibly, the cost-function is E = ⟨Ĥ⟩ = ∑ j c j ⟨P⟩, where each ⟨P j ⟩ can be estimated independently in a quantum computer. 15 However, we will also present a measurement strategy which ensures that the number of circuit measurements required to evaluate E scales minimally with the number of commuting groups inĤ.

Hamiltonian Mapping
The single-electron tight-binding HamiltonianĤ is expanded in our basis as follows: The physics ofĤ is contained entirely within the real-valued hopping parameters t (δ ) αβ ≡ ⟨α 0 |Ĥ|β δ ⟩, which denote the energy cost of an electron transitioning from the orbital |β δ ⟩ to the orbital |α 0 ⟩ in the principal unit cell. The periodicity of the crystal ensures that all matrix elements can be written in terms of the hopping parameters: Note that hopping parameters between orbitals centered on atoms far apart (ie. large δ ) will tend to vanish. Thus, one typically adopts a "nearest-neighbor approximation", in which t (δ ) αβ is non-zero only for those unit-cells δ which hold the nearest images of |β ⟩ to the atomic center of |α 0 ⟩. Realness and hermiticity ofĤ guarantee t (δ ) β α . The number of independent hopping parameters can often be reduced further by exploiting additional crystal symmetries, but we take all hopping parameters as given for the purposes of this work.
The orbital register projection |α⟩⟨β | can be written as the annihilation of an electron in orbital |β ⟩ followed by the creation of an electron in orbital |α⟩, encouraging us to use the fermionic creation and annihilation operators a The fermionic creation and annihilation operators are typically mapped onto Pauli words using well-known transformations such as the Jordan-Wigner or Bravyi-Kitaev transformations, which enforce the fermionic anticommutation relations necessary for representing indistinguishable electrons. 23 However, as pointed out in our previous work, 13 this full machinery is unnecessary for band structure calculations constrained to single-electron states, and the following, simpler mappings suffice: where we have used the notationX α (Ŷ α ) to represent the Pauli word with anX (Ŷ ) on the qubit q α and the identityÎ on every other qubit. The site register projection |ν⟩⟨ν ′ | must bring the basis state |ν ′ ⟩ into the state |ν⟩, and it must annihilate every other state. This action is factorizable into each qubit: where ν q represents the q-th bit in the binary representation of the integer coordinate ν. The reader may verify the following mappings for ν q ν ′ q hold: More concisely, (Note the use of the ⊕ operator to denote bitwise addition, also known as the "XOR" operator.) Given all hopping parameters t (δ ) αβ , Eqs. 10-16 are sufficient to generate the mappingĤ → ∑ j c jPj . Generally,Ĥ will consist of O(M 2 N 2 ) non-zero Pauli words. Adopting a nearest-neighbor approximation reduces this to O(M 2 N). In the next section, we will show how these Pauli words can be partitioned into O(M log N) commuting groups to generate a cost function E = ⟨Ĥ⟩ which can be efficiently evaluated in a quantum computer.

Measurement Strategy
We now write out the full single-electron tight-binding Hamiltonian under a nearest-neighbor approximation. The methods below generalize to multiple dimensions and are easily implemented in code, but the equations become very unwieldy, so we will restrict ourselves in this section to the one-dimensional case. Under the nearest-neighbor approximation, t (δ ) αβ will be non-zero only if δ ∈ {0, ±1}, and we can expand the HamiltonianĤ from above as: where eachĤ (δ ) is defined as follows: Note that we impose periodic boundary conditions so that the projection |ν⟩⟨ν + 1| for ν = N − 1 is identified with |N − 1⟩⟨0|.

7/14
Consider first the orbital register factors ∑ α,β t (δ ) αβ a † α a β . Substituting Eq. 13, In the real partÂ Consider now the site register factors in Eq. 18. The sum ∑ ν |ν⟩⟨ν| is immediately recognized as the identity operatorÎ N acting on the site register. The sum ∑ ν |ν⟩⟨ν + 1| is less trivial, but it too can be decomposed into real and imaginary partsÂ N , In the Supplementary Information, we show thatÂ N andB N each consist of Θ(log N) commuting groups. We now rewrite the HamiltonianĤ of Eq. 17 in terms of theÂ,B operators.
We have used the symmetry relation t Each expectation value in E is measured by the quantum computer. Let Q n be the set of all n-qubit Pauli words consisting of onlyÎ andẐ operators. The procedure for simultaneously obtaining the expectation values ⟨Q j ⟩ of all Pauli wordsQ j ∈ Q n is well understood (see eg. our previous work 13 for a brief tutorial). The same procedure can be applied to any commuting group P n if a basis rotation circuit is applied prior to measurement which transforms each Pauli wordP j ∈ P n to an element of Q n . 24 For example, all Pauli words of the formX αXβ inÂ M can be measured simultaneously by applying the Hadamard gate to each qubit prior to measurement. Fig. 6 presents the measurement circuits required for each register to measure all of the commuting groups inĤ.
So far, the measurement strategy discussed in this section is agnostic of the quantum circuit presented above. Our quantum circuit treats the orbital and site registers as independent, introducing no entanglement between the two. As a consequence, any expectation value ⟨Ô MÔN ⟩ = ⟨Ô M ⟩ · ⟨Ô N ⟩, and we can estimate the expectation values of eachÂ andB operator independently, reducing the number of circuit measurements to Θ(M + log N). Going a step further, we see from Eq. 8 that ⟨Â N ⟩ and ⟨B N ⟩ depend only on the momentum quantum number p, and we can evaluate them a priori as the real and imaginary parts of ∑ ν ⟨Φ N |ν⟩ ⟨ν + 1|Φ N ⟩: Now we may write our cost-function as: In this perspective, k enters in (via p) as a classical parameter of the cost function rather than as an input into the quantum circuit, establishing the equivalence of this method and those of previous results. 9,12,13 Choosing between Eq. 22 and Eq. 24 is a matter of preference and logistical convenience. Figure 7 shows band structures calculated for several model systems, using Eq. 22 and the circuits presented in Figs. 2-6. Fig. 7a corresponds to a one-dimensional lattice consisting of alternating atoms. The unit cell consists of two distinct atoms, each contributing a single orbital with rest energies separated by 1eV. The hopping parameter between each adjacent atom is also 1eV. The lower band corresponds to the "bonding" orbital for individual electrons at each momentum, and the upper band corresponds to the "anti-bonding" orbital. The solid curve is calculated analytically with the standard classical algorithm. Our quantum solution requires two qubits for the orbital register, and we choose to use three qubits (N = 8) for the site register, for a total of five qubits. (Contrast this with a direct simulation of all sixteen atoms in the supercell, which requires sixteen qubits.) Square markers are obtained by simulating our quantum circuits with Google's cirq software package, using the COBYLA optimization protocol (as implemented in the scipy python package) and 8096 circuit evaluations per expectation value. When these results differ noticeably from the analytical solution, we use an X to denote the result of our algorithm in ideal conditions, with perfect optimization and no sampling noise. Fig. 7b corresponds to a two-dimensional hexagonal lattice, such as in graphene. The unit cell consists of two identical atoms, each contributing a single orbital. Each atom has three neighbors, each with a hopping parameter of 1eV. Our results use a supercell size of N = 8 for each dimension, requiring six qubits in the site register and two in the orbital register, for a total of eight qubits. One feature of special interest is the "Dirac cone" surrounding the degeneracy at the high-symmetry point labelled K. Characterizing the Dirac cone and identifying band gaps in materials which lift the degeneracy at K is one common objective of 2D materials science. However, a very high resolution in reciprocal space is required to obtain an accurate characterization. This highlights one significant limitation of band theory using direct-space orbitals: increasing resolution requires larger supercells, which require more orbitals. Our hybrid mapping ensures the number of qubits required scales only logarithmically with the size of the supercell. Nevertheless, obtaining useful resolutions on multi-dimensional lattices still requires many more qubits than are available on present-day quantum hardware, ensuring that the classical approach to band theory will remain dominant for the time being. Fig. 7c corresponds to a three-dimensional simple cubic lattice. The unit cell consists of a single atom with an s orbital and three p orbitals. Rest energies and hopping parameters are selected to qualitatively match the band structure of polonium, which forms a simple cubic lattice in standard conditions (on-site Coulomb interaction and relativistic corrections are required for a more accurate representation). 25,26 Our results use a supercell size of N = 8 in all three dimensions, requiring nine qubits in the site register and four in the orbital register, for a total of thirteen qubits. Note that the path we have shown through reciprocal space includes two redundant branches, ΓR and XM. In particular, note that the quantum results for the RΓ branch happen to be of much lower quality than those of the ΓR branch, despite considering the exact same momentum vectors. This highlights the probabilistic nature of quantum devices, and the urgent need for robust operator estimation and optimation protocols. This need is further exacerbated by the thermal noise and low-fidelity gate operations which plague present-day quantum hardware, and the prevalence of barren plateaus in most cost functions of interest. 27 The purpose of this paper is to showcase the algorithm, rather than obtain high-precision results, and so we have favored simplicity and ease of implementation over rigor. However, the interested reader should be aware that developing efficient and robust quantum protocols is a highly active and relevant area of research, especially in operator estimation, 24, 28-31 optimization, 32-36 and error mitigation. [37][38][39][40]

Conclusion
We have presented a novel approach to band structure calculations in a quantum computer by adopting a basis of local atomic orbitals. This has enabled us to prepare a single cost function which can be used for all k and for all bands. This is an improvement over previous approaches 9,12,13 , which require recalculating new matrix elements for each point in k space and add additional overlap terms in the cost function to explore higher bands. In exchange, we require an additional log N qubits, where N is the resolution in k explored by the band structure.
While band structures are obtainable through classical approaches, the surge of quantum algorithms developed in recent months heralds a new era of computational materials science. Quantum computing enables efficient electronic structure calculations of highly correlated systems for which the single-electron approximation fails. The most promising approaches 7,8,10 10/14 (a) Alternating chain lattice (one dimension). The orbital register requires two qubits, with one orbital per atom and two atoms per unit cell.
(b) Hexagonal lattice (two dimensions). The orbital register requires two qubits, with one orbital per atom and two atoms per unit cell.  assign a unique set of qubits for each periodic basis function, requiring O(N) qubits, although this can be reduced somewhat by tapering methods. 10,41 The hybrid first/second quantized qubit mapping for periodic systems presented in this paper may be adapted to accommodate multiple electrons, exponentially reducing the number of qubits required to express the system. The quantum circuit can be adjusted to express multi-electron states, long-range interactions can be accommodated by introducing entanglement between the orbital and site registers, and additional multi-body terms can be added to the cost function.
One particular difficulty in adapting our approach to multiple electrons is the need to enforce fermionic anticommutation relations, which we have omitted to fully exploit the single-electron approximation. Mapping multi-body fermionic integrals onto the hybrid first/second quantized registers is a subject of further research. Alternatively, one could develop a novel quantum circuit which enforces the appropriate antisymmetry relations on the ansatz directly. Such a strategy is consistent with the theme of this paper: simplifying the cost function and measurement complexity of the VQE algorithm by creatively introducing symmetries and constraints in the quantum circuit.