Almost optimal measurement scheduling of molecular Hamiltonian via finite projective plane

We propose an efficient and almost optimal scheme for measuring molecular Hamiltonians in quantum chemistry on quantum computers, which requires $2N^2$ distinct measurements in the leading order with $N$ being the number of molecular orbitals. It achieves the state-of-the-art by improving a previous proposal by Bonet-Monroig et al. [Phys. Rev. X 10, 031064 (2020)] which exhibits $\frac{10}{3}N^2$ scaling in the leading order. We develop a novel method based on a finite projective plane to construct sets of simultaneously-measurable operators contained in molecular Hamiltonians. Each measurement only requires a depth-$O(N)$ circuit consisting of $O(N^2)$ one- and two-qubit gates under the Jordan-Wigner and parity mapping, assuming the linear connectivity of qubits on quantum hardwares. Because evaluating expectation values of molecular Hamiltonians is one of the major bottlenecks in the applications of quantum devices to quantum chemistry, our finding is expected to accelerate such applications.


I. INTRODUCTION
One of the most promising, practical, and industry-relevant applications of quantum computers is quantum chemistry calculation [1,2].Measuring expectation values of operators (relevant to a physical or chemical system of interest) is an important and often indispensable subroutine in such application.For example, the variational quantum eigensolver (VQE) [3,4], which has been intensively studied for utilizing noisy quantum computers called the noisy intermediate-scale quantum devices (NISQ devices) [5], is a method to prepare an approximate ground state of a given Hamiltonian by iteratively minimizing the expectation value of the Hamiltonian for a trial quantum state.In VQE, measuring the expectation value of the Hamiltonian is essential and consists of the most important part of the algorithm.Even for more advanced algorithms in the quantum chemistry applications of quantum computers, such as quantum phase estimation [6,7], measuring expectation values of operators is important.It is because the phase estimation algorithm allows us to sample eigenvalues and eigenstates of the Hamiltonian but it does not provide properties of the eigenstates.For example, expectation values of operators with respect to the ground state are often required to investigate properties of simulated molecules or materials, such as the force acted on molecule [8,9] or transition amplitude with respect to some perturbation [10,11].Considering that the operators arising in these applications share the same form with the Hamiltonian, we believe that a strategy for measuring the expectation value of the Hamiltonian will remain to be an essential tool even in the future when the quantum phase estimation is executable on fault-tolerant quantum computers.
The number of copies of quantum states needed to estimate the expectation value of the Hamiltonian with sufficient accuracy becomes very large even for relatively small systems, due to the need to estimate expectation values of O(N 4 ) distinct fermionic operators contained in the Hamiltonian for a system with 2N fermions (see Eq. ( 1)).Even for small systems corresponding to ∼ 20 qubits, the required number of copies of quantum states can be over 10 9 using naive strategies [12].Gonthier et al. [13] recently pointed out that several days may be needed to evaluate the expectation value of the Hamiltonian just one time to analyze the combustion energies of small organic molecules with sufficient accuracy, assuming the current speed of NISQ devices.Improvement for the fast evaluation of the expectation values is therefore highly demanded for practical applications of quantum computers to quantum chemistry.We note that the expectation value of the Hamiltonian can be evaluated by the fermionic one-particle and two-particle reduced density matrices (fermionic 1,2-RDMs, Eq. (37)) and that such fermionic RDMs are utilized in several algorithms of quantum computational chemistry such as orbital-optimization [14][15][16].
There have been various studies to tackle the problem of a large number of measurements for evaluating the expectation value of the Hamiltonian [9,[17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].Especially, efforts have been put to reduce the number of the sets of the simultaneously-measurable operators in the Hamiltonian.One can use commutativity or anti-commutativity of the Pauli operators in the qubit representation of the Hamiltonian to group those operators into simultaneously-measurable sets.Bonet-Monroig, Babbush, and O'Brien [29,33] showed that the number of such simultaneously-measurable sets to evaluate the fermionic 2-RDM, which is sufficient to determine the expectation value of the Hamiltonian, is lower bounded by Ω(N 2 ) for 2N fermion systems.They also presented an algorithm to construct the 10  3 N 2 + O(N ) sets of simultaneously-measurable operators for evaluating the fermionic 2-RDM (see Sec. V C), which is optimal up to a constant factor, by using a technique based on a binary partitioning of 2N integers.Hereafter, we call their algorithm BBO algorithm after the initials of the authors.
Contribution.In this work, we provide an algorithm based on a finite projective plane in mathematics to construct 2N 2 sets of simultaneously-measurable fermion operators which are sufficient to determine the expectation value of generic molecular Hamiltonian in quantum chemistry [Eq.( 1)].By using the symmetry of the coefficients of the quantum chemistry Hamiltonian [Eq.( 4)], we classify terms in the Hamiltonian into several types and associate each of them with sets of simultaneously-measurable operators and quantum circuits to measure the operators in the set.The crucial and novel contribution of our algorithm is the use of a finite projective plane to construct the sets of simultaneously-measurable operators for evaluating expectation values of the products of four distinct fermion operators (Sec.III B 4).We formulate the problem of finding the sets of simultaneously-measurable operators as a minimal edge clique cover problem of a certain graph, and find the almost optimal solution of it by using the finite projective plane.The total number of the sets of simultaneously-measurable operators for the whole Hamiltonian scales as 2N 2 + O(N ) for 2N fermion systems in our algorithm when N = Π K for a prime Π and an integer K > 0. It improves the BBO algorithm which scales as 10N 2 /3 + O(N ) when applied in the same setup.We also show that the quantum circuits for measuring those sets can be constructed by O(N 2 ) one-and two-qubit gates with O(N ) depth in the case of Jordan-Wigner mapping [34] and the parity mapping [35].We perform numerical simulation of our method for molecular Hamiltonians of hydrogen chains up to N = 30 orbitals.Our work gives a simple and efficient protocol to measure the expectation value of the molecular Hamiltonian in quantum chemistry so it can accelerate various applications of quantum computers to quantum chemistry calculation.Moreover, since the minimum edge clique cover problem for general graphs is known to be NP-hard, our algorithm using the finite projective plane to create an almost optimal solution for the specific graph may be of an independent interest in the graph theory (we stress that our algorithm does not solve the minimum edge clique cover problem for general graphs; it does solve for the specific graph corresponding to the quantum chemistry Hamiltonians).
This article is organized as follows.We define our problem of measuring the expectation value of the Hamiltonian in Sec.II.Section III describes our main results: an algorithm to measure the expectation value of the Hamiltonian with 2N 2 distinct quantum circuits.We discuss our result and compare it with the previous studies in Sec.V. Summary and outlook are presented in Sec.VI.

II. PROBLEM DEFINITION
In this section, we describe a concrete setup of the problem which we study.

A. Electronic structure Hamiltonian in quantum chemistry
A central task in applying quantum computers to quantum chemistry calculations is to estimate an expectation value of the electronic structure Hamiltonian with 2N spin orbitals (fermions) in the form of where E n is a nuclear repulsion energy (scalar), a † pσ and a qσ are fermion creation and annihilation operators, respectively, satisfying the fermion anti-commutation relation {a † pσ , a † qτ } := a † pσ a † qτ +a † qτ a † pσ = 0, {a pσ , a qτ } = 0, {a † pσ , a qτ } = δ pq δ στ .Two-electron integral g pqrs,στ is defined as and h pq,σ is defined as where ϕ pσ (r) is a one-particle wavefunction for the orbital p with spin σ, N a is the number of nuclei in the molecule, and R A (Z A ) is the coordinate (charge) of the nucleus A, respectively.When the orbital ϕ pσ (r) is real-valued, which is the case for most calculations in quantum chemistry1 , there are symmetries in the coefficients of the Hamiltonian: Using these symmetries, we can rewrite the Hamiltonian as where A pq,σ is defined as We can therefore estimate the expectation value of Hamiltonian for a given state |ψ⟩ by measuring ⟨A pq,σ ⟩ := ⟨ψ|A pq,σ |ψ⟩ and ⟨A pq,σ A rs,τ ⟩ := ⟨ψ|A pq,σ A rs,τ |ψ⟩ for all possible combinations of p, q, r, s, σ, τ .In the rest of this article, we describe how to measure such expectation values.

B. Measurement cliques
Mutually-commuting operators are simultaneously measurable.We call a set of mutually-commuting operators a measurement clique.For a given measurement clique, we can associate it with one quantum circuit to perform projective measurement on simultaneous-eigenstates of all the operators in the measurement clique.Repetitive execution of that circuit yields an estimate of an expectation value of any product of the operators in the measurement clique.For example, when a measurement clique consists of two commuting operators We want to construct a set of measurement cliques that covers all operators contained in the Hamiltonian [Eq.( 5)], namely, where C i is i-th measurement clique, M is the total number of the measurement cliques, O j is j-th operator in C i , and L i is the number operators in C i .We require U to satisfy the condition that all A pq,σ and A pq,σ A rs,τ in the Hamiltonian [Eq.(5)] can be written as a single operator or a product of operators contained in some measurement clique C i .If that condition is satisfied, we can measure all the expectation values of ⟨A pq,σ ⟩ and ⟨A pq,σ A rs,τ ⟩ by using M distinct quantum circuits.Therefore, the problem of estimating the expectation value of the Hamiltonian with a small number of distinct quantum circuits is reduced to find a set of measurement cliques whose number is as small as possible.We construct such a set with M ∼ 2N 2 + O(N ) in the next section.

III. MAIN RESULT: AN EFFICIENT MEASUREMENT SCHEME FOR MOLECULAR HAMILTONIANS
We explain our main results in this section.First we classify the terms in the molecular Hamiltonian [Eq.( 5)] into several types.We then construct measurement cliques to cover all the types of the terms.There are O(N 4 ) terms in the Hamiltonian, but the number of measurement cliques can be reduced to O(N 2 ) by the strategy that the terms in the Hamiltonian are reconstructed by a product of two operators in the clique.Finally, we explain the quantum circuits associated with those measurement cliques and discuss their number of gates.
For the later use, we define the particle number operator of fermion as n pσ := a † pσ a pσ .We also define σ representing the opposite spin of σ, i.e., ↑ =↓, ↓ =↑.

A. Classification of terms
The Hamiltonian (5) consists of A pq,σ and A pq,σ A rs,τ .Independent terms coming from A pq,σ can be classified into two types by considering two cases p = q and p ̸ = q: (1-1) n pσ for p = 0, • • • , N − 1 and σ =↑, ↓.
We note that we have only to consider p > q when p ̸ = q because A pq,σ = A qp,σ .

B. Construction of measurement cliques
Here we construct a set of measurement cliques to cover all the types in the Hamiltonian discussed above.The resulting measurement cliques are summarized as Tab.I.

Measurement clique for particle number operators
The first measurement clique is for the particle number operator: This is a measurement clique because all the particle number operators commute: [n pσ , n qτ ] = 0. C part can determine the expectation values of the terms of the types (1-1), (2-1), and (2-4).
TABLE I. Summary of measurement cliques, the number of cliques for each type (when the number of molecular orbitals N satisfies N = Π K + 1 for a prime Π and an integer K > 0), and corresponding terms in the Hamiltonian.

Measurement clique for Apq,σ
Next, we consider a measurement clique that can evaluate ⟨A pq,σ ⟩ with 0 ≤ p < q ≤ N − 1.By observing that A pq,σ and A rs,σ commute if p, q, r, s are mutually different, a set of pairs of integers, Therefore, it suffices to find the sets satisfying Eq. ( 9) and M [1]   i=1 to estimate all expectation values of ⟨A pq,σ ⟩ with p < q.For any given p < q, there is some I i such that (p, q) ∈ I i so that we can evaluate ⟨A pq,σ ⟩ by the measurement clique C [1] i,σ .One of the ways to find such sets was proposed in the paper of BBO algorithm [29] but we can also use various scheduling algorithm for round-robin tournaments.The number of measurement cliques is M [1] = N − 1 (M [1] = N ) when N is even (odd).
U [1] U [1]  σ Note that we add the particle number operators of the opposite spin to C [1] i,σ to cover the type (2-2).The total number of cliques contained in U [1] is 2M [1] .
4. Measurement clique for Apq,σArs,τ with the same spin Expectation values ⟨A pq,σ A rs,σ ⟩ and ⟨n pσ A rs,σ ⟩ (the same spin case, the types (2-6) and (2-7)) are most non-trivial to construct measurement cliques.Here, we propose a novel algorithm to create the measurement cliques for them by using a projective finite plane.This part is one of the main contributions of this study.
First, observe that the following equations hold if p, q, r, s are mutually different.We want to find sets of pairs of integers , satisfying three conditions, • r • For any mutually distinct integers 0 ≤ p, q, r, s ≤ N − 1 with p < q and r < s, there exists some J i that contains (p, q) and (r, s).
• For any mutually distinct integers 0 ≤ p, r, s ≤ N − 1 with r < s, there exists some J i that contains (p, p) and (r, s).
Note that the set J i allows a pair of the same integer like (p, p) in contrast with I i in Eq. ( 9).Once such sets are found, the measurement cliques can evaluate all the expectation values ⟨A pq,σ A rs,σ ⟩ and ⟨n pσ A rs,σ ⟩ for any mutually distinct integers 0 ≤ p, q, r, s ≤ N − 1 with p < q and r < s.We note that A pp,σ = 2n pσ .Formulation as a edge clique cover problem.We then focus on how to obtain the sets J i for a given N .We formulate the problem of finding the sets J i as a clique cover of a graph.Let us define a graph G with a set of vertices A vertex v = (p, q) corresponds to an operator A pq,σ .Edges of the graph G between two vertices v a = (v a,1 , v a,2 ) and ) are placed if either of the following three cases is met: Moreover, the edge between v a and v b can be associated with the operator O va O v b , and all the operators ⟨A pq,σ A rs,σ ⟩ and ⟨n pσ A rs,σ ⟩ with any mutually distinct integers 0 ≤ p, q, r, s ≤ N − 1 with p < q and r < s are associated to the edges of the graph with one-to-one correspondence.Therefore, finding the sets of J i for measurement cliques is reduced to the problem of finding the edge clique cover of the graph G. Finding an edge clique cover of a given graph with the minimal number of the cliques is NP-hard in general.However, in our case, we know the property of the graph well and we can explicitly construct the edge clique cover having an almost optimal number of cliques with the classical computational time of O(N 3 ).We note that our graph G is different from those in the previous studies [19][20][21][22][23][24][25][26] where the Pauli terms of the Hamiltonian are vertices; rather, edges correspond to the terms in the Hamiltonian in our graph similarly to the BBO's approach [29].Construction of cliques using the finite projective plane.Our algorithm works when N = Π K + 1 for a prime Π and an integer K > 0. When this is not the case, it suffices to choose the smallest Π K larger than N − 1.We consider a finite projective plane of the order Π K as the Desarguesian plane using the finite field F Π K [36].Below, for simplicity, we restrict ourselves to the case of K = 1 where the finite field is Z Π and all the arithmetics (addition, multiplication, and their inverse operations) on the elements of the field are those of integers modulo Π.All the discussion described here applies to the case of arbitrary K > 1 by considering the corresponding arithmetics on the finite field.
The projective plane is defined by P, a set called points, and L, a family of the subset of P which is called lines.The set P consists of three types of points: a single point P α , Π points {P β (y)} Π−1 y=0 , and Π 2 points {P γ (x, y)} Π−1 x,y=0 .
(see also Fig. 2).The lines also consists of three types of the subset of points: , defined as The total number of lines is Π 2 + Π + 1. Intuitively, the point {P γ (x, y)} is interpreted as a two-dimensional plane with coordinates x, y = 0, • • • , Π − 1, and the points {P β (x)} and P α are points at infinity.P and L are known to constitute a finite projective plane.This means that for any two points in P there exists only one line that contains (passes) both.Also, for any two lines in L, there exists only one point that is contained by both simultaneously.
Our idea to find a edge clique cover of G begins by placing Π + 1 = N points on the plane.Namely, we define S = {S(0), • • • , S(Π)}, a subset of P, as (see Fig. 2) We then associate a line L v of the plane to each vertex v of the graph G as follows.For v = (l, l) ∈ V 1 (Eq.( 18)), the line L v is defined as the line that passes S(l) and does not pass other points in S.There always exists only one such line on the plane (see Appendix A 2 for proof).For v = (l, l ′ ) ∈ V 2 (Eq.( 18)), the line L v is defined as the line that passes both S(l) and S(l ′ ).The property of the finite projective plane assures the existence and uniqueness of such line.Moreover, distinct vertices are always associated to distinct lines, i.e., v because there is no line that passes three distinct points in S [37,38] (for completeness, we prove this fact in Appendix A 1).
The complementary set of S in P, K := P \ S, consists of (Π 2 + Π + 1) we associate a set of vertices C k defined as (see Fig. 3), We claim that C 1 , • • • , C Π 2 are cliques of the graph G and constitute a edge clique cover of G.
k=1 is a edge clique cover of G.That is, for any edge between two distinct vertices v and v ′ on the graph G, there exists some associated to these vertices intersects at p k ∈ K by definition.Suppose that there does not exist an edge between v and v S( 2) S(3) S(4) P γ (4, 0) FIG. 2. The projective plane used in our algorithm for N = 6 (Π = 5).Each grid represents a point of the projective plane.
As an example of a line, Lγ(1, 0) is represented by shaded grids.We label N points {S(k)} Π k=0 such that a line passes either only one S(k) or two S(k)'s.For each point that are not labeled S(k), our algorithm lists all lines that pass it.The set of all lines that pass a specific point specifies a measurement clique, that is, the sets of operators that are simultaneously measurable.Examples are given in Fig. 3.
intersects at the point S(v 1 ) ∈ S other than p k .This contradicts with the uniqueness of the intersection between two lines.Similarly, if v, v ′ ∈ V 2 , the non-existence of the edge implies that L v and L v ′ intersects at some point in S other than p k .Again this contradicts with the uniqueness of the intersection between two lines.Therefore, there must exist an edge between v and v ′ on the graph G, which proves (a).
To prove (b), let us take an edge between v and v ′ on the graph G (note that v ̸ = v ′ ).There are two distinct lines L v and L v ′ , and it has exactly one intersection P ∈ P on the plane.The point P must be in K due to the same argument to prove (a).Therefore, there is some k satisfying p k = P , and C k contains v and v ′ .
Definition of measurement clique.Now that we have a set of pairs of integers {C k } Π 2 k=1 by using the finite projective plane, which satisfies Eq. ( 16) and the conditions described there, the measurement clique to evaluate all the terms ⟨A pq,σ A rs,σ ⟩ and ⟨n pσ A rs,σ ⟩ for both σ =↑, ↓ can be explicitly constructed by taking J k = C k in Eq. ( 17).The number of the measurement cliques is Π 2 = (N − 1) 2 .

Summary for measurement cliques
To summarize, the measurement cliques we need to evaluate all the terms in the Hamiltonian [Eq.( 5)] are U = {C part } ∪ U [1] ∪ U [2,diff] ∪ U [2,same] . ( The number of the total measurement cliques is 1 + 2 The result is summarized in Tab.I.Note again that, when N = Π K + 1 for a prime Π does not hold, we can take the smallest Π K to construct U [2,same] .This strategy, in fact, upper-bounds the number of cliques contained in U [2,same] by N 2 + O(N (log N ) 2 ) for any practical N .It is because (c.f.Cramér's conjecture) Π n+1 − Π n < (ln Π n ) 2 holds for 7 < Π n ≤ 5 × 10 16 , where Π n denotes nth prime [39].This guarantees that, for a given (practical) N , we can choose a prime Π > N such that Π = N + O((log N ) 2 ) and run the algorithm to obtain Π 2 = N 2 + O(N (log(N )) 2 ) cliques.

C. Quantum circuits for measuring operators in cliques
After constructing the measurement cliques, we consider how to measure all the operators A pq,σ contained in each measurement clique simultaneously.Strategies for measuring a single A pq,σ on a quantum computer depends strongly on the choice of the fermion-to-qubit mapping.However, typical mappings (e.g., Jordan-Winger, parity, or Bravyi-Kitaev; see also [40]) transform fermionic states with the fixed particle number, pσ ) fpσ |vac⟩ with f pσ = 0, 1 and |vac⟩ being the fermionic vaccum state, into a computational basis state of qubits and thus a particle number operator n pσ of fermions into a projection operator on computational basis states.For such cases, measurement of A pq,σ reduces to a computational basis measurement after applying a unitary U pq,σ such that U † pq,σ A pq,σ U pq,σ can be written by the number operators.In the language of fermions, U pq,σ can, for example, be which leads to The challenge is how to efficiently implement such U pq,σ for all operators in a clique on a quantum computer with possibly limited connectivity.One choice is to use the product of Eq. ( 23) for all A pq,σ in the clique.The overall unitary in this case becomes a spin-conserving fermionic orbital rotation.It is known that such a rotation can be performed using a depth-(2N −3) quantum circuit on linearly connected qubits under the Jordan-Wigner encoding [41].This indeed is an efficient strategy, but below we will present even more efficient strategy which uses parallel application of nearest-neighbor fermionic swap gates for at most N times.Our construction is essentially parallel to that of BBO [29].Suppose that we want to measure m operators in a clique, The measurement cliques other than C part can be written in this form.The operators in C can always be reordered in up-then-down order as, where m σ is the number of σ-spin operators in the clique.For simplicity, we explain the algorithm for the cases where C contains only A kk ′ ,σ with k ̸ = k ′ .The extension of the algorithm to the cases with A kk,σ is straightforward.
Our strategy for measuring the operators in the clique C consists of two main steps.
1. Apply nearest-neighbor fermionic swap gates [Eq.( 29)] to sort the indices of the operators in Eq. ( 26) into The total number of the fermionic swap gates for this sorting is at most N 2 and the depth (counting a parallel fermionic swaps as one layer) is at most N .
2. Transform the operators in Eq. ( 27) to number operators.For example, this can be done by the parallel application of orbital rotations given in Eq. ( 23).However, as we will see later, it is not the only choice.By inspecting the concrete forms of A j,j+1,σ under specific fermion-qubit mappings, we can design simpler quantum circuits for diagonalizing them.We find that the Bell measurement circuit suffices our purpose for Jordan-Wigner mapping and a single layer of Hadamard gates does for parity and Bravyi-Kitaev mappings.
We explain these two steps in order.

Sorting operators in the clique by fermionic swaps
We want to find a unitary U swap such that which is equivalent to the transformation from Eq. ( 26) to Eq. ( 27).To this end, we utilize the fermionic swap gate, which is a Hermitian and unitary operator, i.e., (f swap ij,σ ) 2 = I.Its action can be described as We construct U swap as a product of M ↑ +M ↓ fermionic swap gates which indicates the position of the integer l in the sequence (k 0 , k 1 , ..., k 2m ↑ −1 ).p l = −1 corresponds to the case where l is not found in the sequence.Note that p l is well-defined and there is one-to-one correspondence between (k 0 , k 1 , ..., k 2m ↑ −1 ) and (p 0 , p 1 , ..., p N −1 ) because we assume k 0 , k 1 , • • • , k 2m ↑ −1 are mutually different.Applying fermionic swap gate f swap ij,↑ to the operators in the clique C invokes a change in p l as an exchange between p i and p j and corresponding change in k.Therefore, the problem reduces to finding a sequence of swaps for the elements of p l which results in because this p ′ l corresponds to This problem can be solved by applying the odd-even sort algorithm [42] to p l .This algorithm first compares the odd adjacent pairs (p 1 , p 2 ), • • • , (p N −3 , p N −2 ) and exchanges two components in each pair if necessary.Second, it does the same for the even adjacent pairs (p 0 , p 1 ), • • • , (p N −2 , p N −1 ).Each of these steps requires N/2 operations at most but they can be executed in parallel.Iterating this process, it successfully sorts the integer with N steps (i.e., with N/2 even and N/2 odd steps) at most.Hence we can construct U swap ↑ by interpreting exchanges in the sort algorithm as the fermionic swap gates.U swap ↓ can be constructed in the same manner and can be executed in parallel with U swap ↑ .Therefore, U swap in this construction consists of at most N layers of at most N parallel fermionic swap gates.
The actual representation of the femionic swap gate f swap ij,σ on qubits depends on the fermion-qubit mapping, but it becomes typically local when the swap is occurred at adjacent sites, i = l, j = l + 1 for some integer l.To see this, let us first define fermimonic operators ãj with j = 0, 1, • • • , 2N − 1 labeling both the orbital and the spin simultaneously by ãj := a j,↑ (0 This is the so-called up-then-down convention for ordering spin-orbitals.The Jordan-Wigner [34], parity [35], and Bravyi-Kitaev [35,43,44] transformation respectively maps ãk to where U (j) and P (j) are set of indices defined in Appendix B, and the notation like X U (j) means l∈U (j) X l .From the definition of f swap i,i+1,σ (Eq.( 29)), f swap i,i+1,σ becomes a two-and three-qubit operator in the Jordan-Wigner and parity mappings, respectively, because the chain of Pauli Z operators and X operators cancels out.For the Bravyi-Kitaev mapping, it becomes a O(log N )-qubit operator since the number of elements in the sets U (j) and P (j) is O(log N ).
2. Diagonalizing Aj,j+1,σ in specific fermion-to-qubit mappings Note that the sorted clique (Eq.( 27)) only contains Ã2j,2j+1 and therefore we only need to consider how to diagonalize them.Under three different fermion-to-qubit mappings, Ã2j,2j+1 becomes It is easy to see that, for the Jordan-Wigner mapping, we can diagonalize it by applying the Bell measurement circuit to qubits 2l and 2l + 1.Also, it can be seen that the set P (2l + 2) \ {2l + 1} and P (2l) contains odd integers i < 2l by inspecting Algorithm 1 in Appendix B. This means that, for the parity and Bravyi-Kitaev mappings, the Hadamard gate to all even-site qubits suffices our purpose.The total number of local gates are O(N ) and the depth is O(1) in all of three cases.

IV. NUMERICAL SIMULATION FOR HYDROGEN CHAINS
In this section, we apply our method to the molecular Hamiltonians for hydrogen chains.We compute the number of measurement cliques for these Hamiltonians and estimate the number of measurement shots to realize the standard deviation of the energy typically required by precise quantum chemistry calculations.

A. Setup
We consider hydrogen chains H 4 , H 6 , H 8 , • • • , H 30 where the atomic distance between two adjacent hydrogens is 1 Å.We employ the STO-3G basis set to perform the Hartree-Fock calculation with the numerical library PySCF [45,46], and the Hartree-Fock orbitals are used to construct the molecular Hamiltonian [Eq.(1)].Note that the number of orbitals N is the same as the number of hydrogens with this setup (e.g., N = 4 for H 4 ).The Jordan-Wigner transformation [34] is used to map the fermionic Hamiltonian into the qubit one, implemented by the library Open-Fermion [47].
We compare our method with two existing grouping techniques: qubit-wise commuting (QWC) [18,21] and general commuting (GC) grouping [22,26].Both methods divide the terms in the Hamiltonian (1) into simultaneouslymeasurable Pauli operators.Note that the methods based on the factorization of fermionic Hamiltonians is more efficient than these [27].However, since such methods require us to apply fermionic orbital rotation operations which are non-Clifford and need more gates depending on hardware, they cannot be directly compared to our methods.We thus omit comparison with them in this work.For QWC grouping, the quantum circuits for measurement consist of O(N ) one-qubit Clifford gates with depth one.For GC grouping, the number of two-qubit gates (CZ or CNOT, both are Clifford) for the measurement circuits is O(N 2 / log N ) [26,48].We implemented the "sorted insertion" algorithm for both groupings [26].The details of both grouping methods and their implementations can be found in Appendix C 1.
After the grouping, we estimate the number of shots required to suppress the standard deviation of the energy expectation value down to 10 −3 Hartree.This value (10 −3 Hartree) is comparable to the so-called chemical accuracy that is typically targeted at in precise quantum chemistry calculations.The concrete procedures to estimate the number of shots are in Appendix C 2. 4. Computational time required for performing each grouping methods.We show two plots for our method.The blue one shows the time required for generating grouped qubit operators, that is, generating {Gi} with Gi = {P i j } where P i j is a 2N -qubit Pauli operators such that [P i j , P i j ′ ] = 0 for all i, j, j ′ .This is for a fair comparison with GC and QWC methods which generate such lists of groups.The light blue plot illustrates the time required to construct all cliques, represented by U, defined as lists of integer pairs.

B. Numerical results
We performed grouping of the hydrogen chain Hamiltonians by our method up to N = 30, while we did it by QWC and GC grouping up to N = 20.This is because the classical computational time is much longer for the two (QWC and GC) methods as illustrated in Fig. 4. Our implementations of grouping methods are pure-Python strongly depending on OpenFermion [47], works as a single-thread, and is not optimized for speed.The CPU used in this experiment is AMD EPYC 7252. Figure 4 shows that our method is highly efficient compared to the existing methods.
Numerical results are shown in Fig. 5.In Fig. 5(a), we compare the number of groups (measurement cliques) obtained by each grouping method.To see the scaling of the number of groups, we also fit the results by the function aN b , where a, b are the fitting parameters, by the linear regression of the data points with N ≥ 10 on the log-log plot.This fitting yields b = 2.03 (7) for our method, b = 2.34(3) for GC grouping, and b = 4.017(3) for QWC grouping.We observed that GC grouping exhibits the slightly smaller number of groups than that of our method for N ≤ 20.However, the predicted scaling exponent of GC grouping is larger than that of our method.We expect that our method beats GC grouping for larger N in terms of the number of groups (see the fitting lines in the figure).We also note that the number of groups of our method seems to follow the theoretical scaling 2N 2 − 2N + 1 (the total number of all cliques [Eq.( 22)]).
In Fig. 5(b), we compare the number of shots required to make the standard deviation of the energy 10 −3 Hartree.We again performed the fitting of the data by the function aN b in the same way as we did for the number of groups, which results in b = 3.51(5) for our method, b = 3.12(7) for GC grouping, and b = 4.542 (7) for QWC grouping.We obeserve that GC grouping shows the smaller number of shots than that of our method.In contrast to the number of groups, the scaling exponent is smaller in GC grouping, indicating that it would performed well than our method if we could execute GC grouping for larger N > 20.This is because "sorted insertion" algorithm performs grouping considering the magnitudes of coefficients of Pauli operators while our method are not aware of such information.We stress that, however, the huge classical computational cost will prevent GC grouping to be applied for larger N .Our method can be seen as a practical and effective way to group the Pauli operators in the Hamiltonian and reduce the number of shots for the desired accuracy of the energy expectation value.It might be possible to further reduce the number of shots of our method by taking the information about coefficients into account.

V. DISCUSSION
We here discuss several points of our study and compare our result with the previous studies.

A. Optimality of the obtained solution
Our construction of cliques for the graph G in Sec.III B 4 is almost optimal from the viewpoint of the minimal edge clique cover problem.To show this, we give a lower bound to the number of cliques to cover edges of G. Consider a subgraph G ′ where we remove all of vertices in V 1 and corresponding edges in G.The number of cliques to cover all edges in G ′ is clearly less than or equal to that of G, and it follows that a lower bound for the number of cliques for the edge cover of G ′ is also a lower bound for G. First, we claim that a clique in G ′ can have at most N/2 vertices and thus N 2 ( N 2 − 1)/2 edges.This is because there is an edge between two vertices (p, q) and (r, s) in V 2 if and only if p, q, r, s are mutually different, and thus a clique corresponds to a pairing of N integers.Second, the number of edges in G ′ is N (N − 1)(N − 2)(N − 3)/8.To see this, for every tuple of four integers (p, q, r, s) such that p < q < r < s, there are three edges each of which connects the vertices {(p, q), (r, s)}, {(p, r), (q, s)}, and {(p, s), (q, r)} in a graph G ′ .The number of combinations of such integers is N  4 , and therefore the statement follows.Combining these two facts, we conclude that the number of cliques to cover all edges in G ′ must be at least 3).Therefore, our algorithm, where the number of cliques is (N − 1) 2 , provides an optimal solution to the minimal edge clique cover problem of G up to the subleading order terms.

B. Classical computational cost for constructing measurement cliques
The classical computational cost for constructing the measurement cliques in our method is O(N 3 ), which comes from creating a set C k in Eq. ( 21).For each p k ∈ K (the number of k is Π 2 = O(N 2 )), we add a vertex v on G to C k by inspecting the line passing p k and S(i) (the total number of i is Π + 1 = N ): when the line passes S(i) and S(j), we add (i, j) or (j, i) to C k , and when the line passes only S(i), we add (i, i) to C k .The exponent observed in Fig. 4 is slightly worse than this theoretical scaling, which we attribute to additional logarithmic costs (the overhead happening when N ̸ = Π + 1, memory access, integer arithmetics, etc.).Note that this computational cost is optimal when we seek to find O(N 2 ) measurement cliques, because the output of the algorithm must be of size O(N 3 ) (O(N 2 ) cliques with O(N ) vertices) by the discussion in Sec.V A. BBO's paper [29] does not explicitly give the classical computational cost, but the size of the output of their algorithm is also O(N 3 ) and thus its classical computational cost cannot be less than O(N 3 ).
This O(N 3 ) scaling of the classical computational cost is much smaller than the methods to create the measurement clique consisting of the terms of the Hamiltonian as is [17][18][19][20][21][22][23][24][25].In those methods, one should calculate the commutativity or anti-commutativity of pairs of the O(N 4 ) terms in the Hamiltonian, which results in as huge as O(N 8 ) classical computational cost.The target problem of the application of quantum computers in quantum chemistry lies for more than N ∼ 50, and such scaling can be problematic in practice.
In this study, we have proposed a measurement scheme for general molecular Hamiltonians in quantum chemistry [Eq.(1)].We have used the symmetries of the coefficients in the Hamiltonian and reduced the problem of evaluating the expectation value of the Hamiltonian into the evaluation of the terms ⟨A pq,σ ⟩ and ⟨A pq,σ A rs,τ ⟩.We have classified all the possible terms appearing in the Hamiltonian and proposed a measurement clique for each type of the terms.Especially, the measurement cliques for the terms ⟨A pq,σ A rs,σ ⟩ with mutually-different p, q, r, s have been constructed by finding the edge clique cover of the specific graph using the novel method based on the finite projective plane.The total number of the distinct measurement cliques (or measurement circuits) is 2N2 + O(N ) with 2N being the number of spin orbitals (fermions), which exhibits better scaling compared with the previous method (BBO algorithm).We have shown explicit quantum circuits to measure operators in the measurement cliques, and those circuits consist of O(N 2 ) one-or two-qubit gates with depth O(N ) for Jordan-Wigner and parity mappings.We have also performed numerical simulation for molecular Hamiltonians of hydrogen chains and evaluated the number of groups of simultaneously-measurable operators generated by our method as well as the number of measurement shots required to estimate the energy expectation values with sufficient accuracy.Evaluation of the expectation value of the Hamiltonian is one of the most fundamental subroutines among various algorithms for the applications of quantum computers to quantum chemistry, so our method can be utilized in broad studies on quantum algorithms.
As future work, it is intriguing to apply the proposed method to actual Hamiltonians in quantum chemistry and numerically investigate a standard deviation of the estimated energy expectation value with using our measurement cliques.Another interesting direction is to generalize our method using the finite projective plane to higher order fermionic RDMs, such as ã † p ã † q ã † r a s a t a u .
x 0 + x 1 + • 2 are mutually different.When an operator corresponding to a vertex v = (p, q) is denoted O v , i.e., O v := A pq,σ , this graph G has an edge between v a and v b if and only if [O va , O v b ] = 0.As an example, Fig. 1 shows G for N = 6.It means that a clique C (a subset of vertices with edges between all pairs of vertices contained in the set) of G can define a measurement clique.

↑
and U swap ↓ can be constructed exactly in the same manner, so we consider how to construct U swap ↑ below.Define the following integers p l for l = 0, 1, • • • , N − 1 as FIG.4.Computational time required for performing each grouping methods.We show two plots for our method.The blue one shows the time required for generating grouped qubit operators, that is, generating {Gi} with Gi = {P i j } where P i j is a 2N -qubit Pauli operators such that [P i j , P i j ′ ] = 0 for all i, j, j ′ .This is for a fair comparison with GC and QWC methods which generate such lists of groups.The light blue plot illustrates the time required to construct all cliques, represented by U, defined as lists of integer pairs.

FIG. 5 .
FIG. 5. Numerical results for the molecular Hamiltonians for hydrogen chains HN .(a) The number of groups for each grouping method.Fitting lines for the data points with N ≥ 10 are also shown.The Black dotted line denotes the theoretical scaling for our method, 2N 2 − 2N + 1.(b) The number of shots required to achieve the standard deviation of the energy lower than 10 −3 Hartree, estimated by Haar random averaging (see main text).

we can estimate the expectation values such as ⟨O 1 ⟩ , ⟨O 2 ⟩, and ⟨O 1 O 2 ⟩ by the projective measurement on the simultaneous eigenstates of O 1 and O 2 .
• • + x 15 x 0 + x 1 + • • • + x 7 x 0 + • • • + x 3 x 8 + • • • + x 11 x 0 + x 1 x 4 +x 5 x 8 + x 9 x 12 + x 13 0 y 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y 10 y 11 y 12 y 13 y 14 y 15 FIG. 6. Illustration of Fenwick tree of size n = 16 which is a basis of Bravyi-Kitaev mapping.An n-dimensional vector x can be transformed to the Fenwick tree y by taking appropriate range sums defined in Eq. (B1).Arrows represent what is stored as yi.
y Algorithm 1 Generation of P (l) i ← i − LSB(i) P ← P ∪ {i − 1} return P