Quantum algorithm for electronic 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 first 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 first/second-quantized qubit mapping which allows us to construct a single Hamiltonian, and a single cost-function, suitable for solving the entire electronic 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 modifications 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.


Scientific Reports
| (2022) 12:9867 | https://doi.org/10.1038/s41598-022-13627-x www.nature.com/scientificreports/ 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 word P 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 electronic 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 electronic 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: 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 (2) |α k � = r e ik·r |α r �.
(3) |α ν � = |α� ⊗ |ν�.  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. Equation (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.

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. www.nature.com/scientificreports/ 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.
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

Cost function
Our objective in this section is to present a cost-function suitable for electronic 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 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 † α , 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 singleelectron states, and the following, simpler mappings suffice: where we have used the notation X α ( Ŷ α ) to represent the Pauli word with an X ( Ŷ ) 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: 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. www.nature.com/scientificreports/ 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 onedimensional 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|.
Consider first the orbital register factors α,β t (δ) αβ a † α a β . Substituting Eq. (13), αβ (X αXβ +Ŷ αŶβ ); M , all terms of the form Ŷ αXβ>α and X αŶβ>α are each commutative for fixed α , so that B (δ) M contains �(M) Pauli groups. 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 , B N : In the Supplementary Information, we show that Â N and B 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 words Q 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 word P j ∈ P n to an element of Q n 24 . For example, all Pauli words of the form www.nature.com/scientificreports/ can be measured simultaneously by applying the Hadamard gate to each qubit prior to measurement. Figure 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 Â and B 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 Eqs. (22) and (24) is a matter of preference and logistical convenience. Figure 7 shows electronic band structures calculated for several model systems, using Eq. (22) and the circuits presented in Figs. 2, 3, 4, 5, 6. Figure 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 1 eV. The hopping parameter between each adjacent atom is also 1 eV. 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 16 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. Figure 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 1 eV. 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. Figure 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 electronic 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][33][34][35][36] , and error mitigation [37][38][39][40] .

Conclusion
We have presented a novel approach to electronic 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. www.nature.com/scientificreports/ 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 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 . With further research, we believe 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.