Cooper-Pair Condensates with Non-Classical Long-Range Order on Quantum Devices

An important problem in quantum information is the practical demonstration of non-classical long-range order on quantum computers. One of the best known examples of a quantum system with non-classical long-range order is a superconductor. Here we achieve Cooper pairing of qubits on a quantum computer to represent superconducting or superfluid states. We rigorously confirm the quantum long-range order by measuring the large $O(N)$ eigenvalue of the two-electron reduced density matrix. The demonstration of maximal quantum long-range order is an important step towards more complex modeling of superconductivity and superfluidity as well as other phenomena with significant quantum long-range order on quantum computers.


I. INTRODUCTION
Phenomena like superconductivity and superfluidity arise from a Bose-Einstein-like condensation of fermion pairs into a quantum state with large non-classical longrange order . Recently, quantum computers have emerged as potentially powerful calculators of correlated quantum systems [22][23][24][25][26][27][28][29][30][31][32], which foreshadows the potential emergence of a significant advantage of quantum computers over classical computers for certain classes of problems-a phenomenon known as quantum advantage [33][34][35]. Here, we prepare and measure Cooper-like pairs of qubit particles on a transmon-qubit quantum computer. As the distinction between bosonic and fermionic statistics is lost as a result of a Jordan-Wigner mapping, condensations of Cooper-like qubit pairs can be interpreted as fermion-pair condensations, which can represent superconducting (or superfluid) states.
In this study, qubit particles-which are hard-core bosons-are entangled into Cooper-like bosonic pairs (see Fig. 1) to form superconducting-like states-extreme antisymmetrized geminal power (AGP) wave functions [12,[36][37][38][39][40][41][42][43][44][45][46]. As originally shown by Yang [37] and Coleman [39,41,47], such states are extreme in the set of two-electron reduced density matrices (2-RDM), exhibiting the largest possible eigenvalue of the 2-RDM on the order of the number N of electrons O(N ) that represents the maximum possible number of Cooper pairs in a common two-electron (geminal) eigenfunction of the 2-RDM. We use tomography on the quantum computer to measure a subblock of the 2-RDM [41,48,49] (see Eq. (9)) containing the large eigenvalue. Diagonalization of this subblock on a classical computer produces the large eigenvalue and confirms the preparation of the extreme states with maximal non-classical (off-diagonal) long-range order (ODLRO) [37,47]. Even though the extreme AGP functions are expressible as projections of product wave functions [39,41], they have contributions from an exponentially scaling number of orbital-product * damazz@uchicago.edu configurations (see Fig. 2). Moreover, the measurement of the large eigenvalue of the 2-RDM is applicable to confirming non-classical long-range order in a much richer set of quantum states. Because a necessary criterion for the modeling of superconductors (superfluids) on the quantum computer is the ability to capture the ODLRO, its demonstration provides a first step towards modeling more-complex superconducting (superfluid) materials. The superconducting-like state on the quantum computer is prepared by entangling pairs of qubits into Cooper-like bosonic states. Consider the creation of a state with a Cooper pair of electrons, an extreme geminal [2,[39][40][41], from the vacuum state FIG. 2: A schematic demonstrating the possible configurations (i.e., each row) for a given number r of qubits where a filled circle indicates the |1 state which corresponds to an occupied orbital and an unfilled circle represents the |0 state which corresponds to an unoccupied orbital.
where j andj are the indices of the paired orbitals φ j and φj, the sum over j is taken with respect to all pairs, and θ is an arbitrary global phase. If we represent each orbital by a qubit with the |0 state representing an unfilled orbital and the |1 state representing a filled orbital, we can use a specific Klein transformation [50] known as the Jordan-Wigner mapping [51] a † j = e iπ j−1 to map the fermionic operators in Eq. (1) to qubit operators to obtain |g = j e iθ e iπ j −1 If the paired orbital indices j andj are selected to be consecutive integers in the range [1, r] where r is the total number of orbitals, then the Jordan-Wigner mappings in Eq. (3) simplify to a negative global phase which we can cancel by selecting θ = π to obtain Hence, the extreme geminal of the Cooper pair |g jj of electrons can be represented as two-qubit excitations without approximation. The difference between the fermion and qubit statistics, typically included through an explicit many-qubit Jordan-Wigner mapping , disappears from the pairing of the orbitals to generate an extreme geminal. Moreover, the explicit details of the pairing of the particles is contained within the unspecified orbitals φ j and φj. Consequently, the extreme geminal can physically represent Cooper pairing of electrons in a superconductor or a superfluid in addition to representing even the Cooper-like pairing of qubit particles (hard-core bosons) [52] which are paraparticles [53,54].
The N -electron extreme AGP wave function |Ψ N AGP for even N can be generated from the wedge product of the extreme geminal with itself N/2 times [39][40][41]55] |Ψ N AGP = |g(12) ∧ |g(34) ∧ ... ∧ |g((N − 1)N ) (5) where the wedge ∧ denotes the sum of all products resulting from the antisymmetric permutation of the particles. We can also consider a wave function |Ψ AGP , also known as a Bardeen-Cooper-Schrieffer (BCS) wave function [1], that is a linear combination of the |Ψ N AGP for all N which is expressible as a product state Using the Jordan-Wigner transformation and simplifying as above, we can generate the AGP state in Eq. (6) with the qubit excitation operators which can also be cast as the tensor multiplication of r/2 distinct extreme geminals (or the |Φ + Bell states [56]) where j specifies the pair index and adjacent qubits with qubit indices 2j − 1 and 2j paired by definition. This state, which is composed of substates with all possible, paired, even-numbered excitations, can be prepared on a quantum device according to the general gate sequence given in Eq. (12). Figure 3 shows the specific r = 4 preparation.
On the quantum computer, the extreme AGP state is physically composed of Cooper-like pairs of qubits. Because the phase changes from the fermionic statistics are lost in the pairing-as seen in the above fermionic encoding of the qubits, the extreme state rigorously represents not only entangled pairs of qubits but also Cooper pairs of electrons that are entangled in superconducting or superfluid states. Moreover, the state can represent any physical model for pairing as the precise nature of the pairing (i.e. the pairing of electrons in physical space or in momentum space) is contained in the paired orbitals φ j and φj, which are left unspecified. All pairing states have fundamental entanglement and order properties that are independent of the physical definition of the orbitals. The non-classical long-range order of the extreme AGP state can be assessed from the number of Cooper pairs in the same extreme geminal, which is determinable from the largest eigenvalue of the 2-RDM [37][38][39]. In order to measure whether the experimentallyprepared quantum state and/or the number-conserving substates-which are all possible, even eigenstates of the number operator (Eq. (27)) that can be projected out from the overall ensemble quantum state according to the methodology presented in Appendix C-demonstrate off-diagonal long-range order, we conduct quantum tomography (see Appendix B and Appendix C) to probe directly the presence and extent of ODLRO . To determine the presence and degree of this long-range order for a specified quantum state, it is useful to establish a calculable, characteristic property [12,[36][37][38][39][41][42][43]. Such a signature of ODLRO is a large eigenvalue in the 2-RDM, which we denote as λ D [37,38]. As this large eigenvalue corresponds to the number of Cooper-like qubit pairs occupying the same two-qubit geminal (which is directly analogous to the one-qubit orbital), any λ D value exceeding the Pauli-like limit of one is indicative of ODLRO [37][38][39].
While the entire 2-RDM can be measured by quantum tomography, only the following subblock of the 2-RDM [41,48,49] is required due to the block diagonal structure of the 2-RDM of the AGP wave function where j/j and k/k represent paired fermions, which are given by j = 2m/j = 2m − 1 and k = 2n/k = 2n − 1 for integers m, n in the framework of the AGP wavefunction. After Jordan-Wigner transformation to the qubit representation, we can equivalently represent this block of the 2-RDM in terms of the qubit excitation operators as For fixed number N of electrons, if the 2-RDM is normalized to N (N − 1) as in second quantization, the maximum eigenvalue for even N is bounded from above by N as shown by Yang [37] and Sasaki [38]. Moreover, for a finite rank of r orbitals, this bound can be further tightened [39,41] to While the thermodynamic limit is not reached until r → ∞, even for finite r, as long as N ≥ 4, the 2-RDM exhibits a large eigenvalue that represents the non-classical long-range order associated with Cooper pairing. The 2-RDM from the non-number conserving extreme AGP state |Ψ AGP also exhibits a large eigenvalue, representing an average of the Cooper pairs in each of the fixed-N extreme AGP states (i.e., the number-conserving substates). The number-conserving blocks of the 2-RDM with even particle numbers-i.e., 2-RDMs of zero, two, four, ..., r − 2, and r particles for an r-qubit system-can be determined from the non-number conserving state via post-measurement analysis (see Appendix C). Analysis on the presence and extent of ODLRO(measured via λ D ) of both the overall entangled state (|Ψ ) and the numberconserving substates is conducted for various numbers r of total qubits in the following sections.

III. RESULTS
The extreme non-number conserving AGP state is prepared for both simulation utilizing IBM's QASM simulator (ibmq qasm simulator) and an experimental quantum device for all even-numbered qubit systems from r = 0 to r = 14. Post-measurement computation of the quantum signature of off-diagonal long-range order (λ D ) is then employed to probe the presence and extent of ODLRO for these overall states. As can be seen in Figure 4, the signature of ODLRO increases as the number r of qubits comprising the system is increased, and-for QASM simulation-ODLRO is observed (i.e., λ D > 1) for all prepared states with r ≥ 8. While the experimental results deviate from QASM simulation due to the noisy nature of near-term quantum devices [57] (see Appendix D and Appendix F), experimental systems with r = 12 and r = 14 qubits did demonstrate ODLRO. Further, the trend of the extent of ODLRO increasing as the number of qubits comprising the system increases holds for the experimental results, which is promising for future benchmarking of quantum computers through the preparation of extreme AGP states with larger number of qubits as well as efforts to probe more macroscopically-scaled materials demonstrating ODLRO on quantum devices.
As the non-number conserving extreme AGP state is an ensemble composed of number-conserving substates, the large eigenvalue associated with the ODLRO of the ensemble state is the ensemble average of the substates. By definition, then, the long-range order of the ensemble is less than that of the substate with the largest degree of ODLRO, which is expected to occur around the center of the number distribution N ≈ r/2. Additionally, realworld materials demonstrating ODLRO such as superconductors should conserve particle number. It is hence beneficial to probe the number-conserving substates that comprise the overall entangled state in order to both isolate the ODLRO behavior of the number-conserving substates and to more-closely model real-world materials. By projecting out a specific number of particles from the results obtained for overall entangled state (see Appendix C), we can probe the behavior and properties of the number-conserving substates. Specifically, as is shown in Fig. 5, the extent of ODLRO (λ D ) for each number-conserving state can be isolated from the overall r-qubit preparation described in Eq. (8). As can be seen from the QASM simulation results, all numberconserving substates with 2 < N < r demonstrate ODLRO (λ D > 1) where N = 2 fails to demonstrate condensation behavior as the maximum signature of condensation is N/2 for even N -particle systems [37,38] and where N = r fails to demonstrate condensation behavior as this substate describes the state in which all orbitals are fully occupied with no entanglement-i.e., the single Slater determinant |1 ⊗r -and is hence expected to have a maximum eigenvalue of λ D = 1. Further, the signature of condensation seems to follow a bell curve centered around (r+2)/2 such that maximum ODLRO is observed at half filling for N = (r +2)/2 if (r +2)/2 is even and for both N = (r+2)/2−1 and (r+2)/2+1 if (r+2)/2 is odd. Again, the extent of ODLRO is lesser for the experimental results for all particle-conserving states due to experimental error; however, the qualitative trends described for QASM simulation hold in general although the bell curve does demonstrate a slight negative (right-modal) skew, implying that the quantum computer does not exactly treat the particle and hole statistics symmetrically. Importantly, ODLRO is clearly observed for r = 14 experimental results for particle numbers N ≥ 6. Note that although only results for the largest-qubit preparation r = 14 are shown, all data is included in Table II; the trends in the r = 14 data hold for the lower-qubit results, and additionally, the r = 14 qubit data demonstrates the largest signature of off-diagonal long-range order as the largest λ D value for a fixed N increases as the number r of qubits is increased.

FIG. 5:
The λ D values for the number-conserving substates of rank r = 14 for QASM simulation and experimental melbourne results.

IV. CONCLUSIONS
Here we prepare superconducting-like states from the Cooper pairing of qubits on a transmon quantum computer-where each qubit is composed of a microwave phonon in an anharmonic well potential. Using the Jordan-Wigner mapping between fermions and qubits, we rigorously show that the prepared states are equally valid representations of condensations of Cooper pairs of fermions, bosons, or qubits. Hence, such Cooper-pairbased condensations and their associated non-classical long-range order are independent of the particle statistics. Moreover, the prepared states are also independent of the physical details of the paired orbitals, and consequently, are representative of superconducting, superfluid, or other pairing states. The studied states are known as extreme AGPs because they exhibit the maximum degree of non-classical (off-diagonal) long-range order as determined by the number of Cooper pairs in the same geminal state, which is equal to the largest eigenvalue of the 2-RDM [2, 37-39, 41, 44, 47]. We measure a subblock of the 2-RDM [41,48,49] on the quantum computer and compute its largest eigenvalue on a classical computer. We observe large eigenvalues both for the nonparticle conserving extreme AGP state and the particleconserving extreme AGP substates. The large eigenvalues confirm the preparation of these extreme AGP states, which are the only states to exhibit the largest possible eigenvalues [41], as well as the generation of maximum non-classical long-range order.
The upper bound on the largest eigenvalue of the 2-RDM, λ D = N , is technically only reached in the thermodynamic limit of r → ∞. However, as seen in Eq. (11), the large eigenvalue is rapidly approached with increasing r as the fraction of Cooper pairs that are removed from the condensate due to finite size effects scales as 1/r. Consequently, the quantum long-range order as well as its associated entanglement begin to appear for the range of r (r ≤ 14) explored in the present study. On both IBM's QASM simulator and an IBM quantum computer, we observe that the large eigenvalue follows the expected bell curve with respect to r. The known value for the maximum eigenvalue of the extreme AGP state provides a clear metric for not only confirming the presence of the extreme state and its long-range order but also benchmarking the fidelity with respect to noise of both current and future quantum computers.
An aspirational goal of quantum computing is to achieve a quantum advantage over traditional classical computing for the solution of a significant problem. One such area of chemistry and physics, which traces back to the original proposal of Feynman [58], is the simulation of molecules on quantum computers. The construction of the wave function on a classical computer scales exponentially in the number of orbital-based configurations. In principle, the quantum computer offers the possibility of preparing and measuring quantum states with nonexponential scaling. In the present case, the extreme AGP wave function is a product state composed of a product of extreme geminals (or Bell states). Consequently, the maximum degree of non-classical long-range order, at least as measured by the largest eigenvalue of the 2-RDM, can be achieved with polynomial cost on both classical and quantum computers. The extreme AGP states, nonetheless, provide an intriguing reference for the exploration of more complicated states demonstrating large ODLRO (i.e., demonstrating large eigenvalues) that cannot be easily expressed as product-state wave functions. From this perspective, the present work of preparing and measuring superconducting-like states from Cooper pairs of qubits on a quantum computer provides an initial step towards investigating more complicated condensates of Cooper pairs-with potential future applications to the study of both superconducting materials and simulation.
Certain electronic structure methodologies-such as Hartree Fock (HF), density functional theory (DFT) and coupled cluster (CCSD) theory-are incapable of demonstrating off-diagonal long range order (ODLRO) indicative of particle-particle condensation. As such, in order to accurately explore complex, real-world superconducting states on quantum devices, it is necessary to establish that noisy intermediate quantum devices are both capable of exhibiting ODLRO as well as the features of state preparations that yield such ODLRO. In the strongly-correlated regime, HF-based and multi-reference approaches on quantum devices are not suitable starting points for state preparations for exploring ODLRO. Instead, we show that numberprojected BCS wavefunctions-alternatively termed AGP wavefunctions-are ideal starting points for build-ing more-complicated wavefunctions that may model real-world superconducting states as we demonstrate that these AGP wavefunctions created by the pairing of adjacent qubits are capable of demonstrating ODLRO on the noisy near-term quantum devices. Such insight may additionally prove useful in the development of an appropriate ansatz for hybrid quantum-classical algorithmssuch as the variational quantum eigensolver (VQE)-in order to accurately describe strongly-correlated systems demonstrating ODLRO.

APPENDIX
We include details on the quantum algorithms used to prepare the qubit states presented in the article; the quantum tomography of the particle-particle reduced density matrix for both the overall non-number conserving state as well as number-conserving substates; a description of noise on near-term quantum devices; an analysis of errors by comparing simulated and experimental joint probabilities of occupation; relevant details on the experimental quantum backends employed; full results obtained from ibmq qasm simulation, ibmq 16 melbourne, ibmq 5 yorktown, ibmq santiago, and ibmq rochester; and information regarding the state preparation fidelity.

Appendix A. State Preparations
The overall quantum state is composed of all pairwise even excitations of the r/2 individually-paired qubits. This preparation is accomplished by where |0 ⊗r is the initial quantum state in which all qubits are in their ground state (i.e., all orbitals are unoccupied), p represents the index of each of the possible r/2 adjacent, paired qubits, H i is the Hadamard gate acting on Qi, and C j i is a CNOT gate with Qi and Qj acting as the control and target qubits, respectively. Application of the gate sequence given in Eq. (12)-represented pictorially for r = 4 in Fig. 3-produces the AGP wavefunction described by Eq. (8).
Appendix B. Quantum Tomography for the Particle-Particle RDM While the full particle-particle RDM has elements given by for all possible combinations of one-boson spin orbitals indexed by i, j, k, and l in a finite basis state with rank r, the large eigenvalue, λ D , is contained within the subblock of the 2-RDM given by [39,41] and, hence, only this portion of the matrix is constructed. Note thatâ † i andâ i are creation and annihilation operators for orbital i (and thereby qubit Qi), which can be represented in matrix form aŝ such that each creation operator creates a particle in an empty orbital i-takes a qubit i from |0 to |1 -and each annihilation operator kills a particle in a filled orbital itakes a qubit i from |1 to |0 -where and After construction of the subblock of the 2-RDM corresponding to the large eigenvalue ( 2 D s.b. ), the signature of off-diagonal long-range order-λ D -is then obtained by solving the eigenvalue equation with the signature corresponding to the largest D i value.

a. Tomography via Pauli Expectation Values
The subblock of the 2-RDM containing the large eigenvalue (i.e., Eq. (14)) can be obtained via translation of each of its elements into the bases of Pauli matrices, the expectation values of which can be directly probed on a quantum device. Specifically, the creation and annihilation operators can be rewritten aŝ with diagonal elements being given bŷ and off-diagonal elements being given by Therefore, each 2-RDM matrix element can be obtained by directly probing the expectation values of four-qubit tensor products of Pauli matrices.
As all wavefunctions prepared in this study are real, the 2-RDM should consist of only real-valued elements; hence, all imaginary components of 2-RDM elements should approach zero within a small range dictated by randomness inherent to quantum systems as well as error on the device employed. Therefore, only eight of the sixteen four-qubit expectation values corresponding to real contributions to a given 2-RDM element are nonzero and hence necessary for the determination of the subblock of the 2-RDM; thus, to lower computational expense, only these real values are used to construct the 2-RDM subblock where tomography via Pauli expectation values is conducted (as is consistent with previous analysis conducted in a manner similar to that described in Ref. 31).

b. Tomography via Direct Computation of the Wavefunction
As can be observed from Eq. (8), the phase angle of all qubits in the AGP wavefunction are known to be uniformly zero. As such, knowledge of the probabilities with which each of the possible 2 r basis states (|η i ) for an rqubit calculation are sampled out of the 81,920 times (8192 per trial) a given state is prepared and probed is sufficient information to completely construct the wavefunction (|Ψ ). Specifically, the wavefunction takes the form of a vector (using traditional qubit vector notation [59]) with each element of the wavefunction |Ψ i corresponding to the basis state |η i being given by which is the positive square root of the probability (p(η i )) with which the corresponding qubit basis state is measured. Each individual element of the matrix shown in Eq. (14) is then computed from the appropriate expectation value Ψ|â † jâ † j+1â k+1âk |Ψ for the wavefunction in vector form obtained for a given state preparation. The operatorâ † jâ † j+1â k+1âk can be represented as a 2 r × 2 r matrix according tô which is the tensor product of the creation and annihilation operators in matrix form (Eqs. (15) and (16)) corresponding to the appropriate qubits (j, j + 1, k, k + 1) and identity matrices for each spectator qubit.
On the error-prone melbourne experimental device, the Euclidean distances for the r = 2 and the r = 8 calculations are [r = 2 : 0.015, r = 8 : 0.338], which is likely due to a large degree of error on the quantum device, which increases with the number of qubits (and hence the number of two-qubit gates) involved. In fact, for the same number of qubits, the Euclidean distances between 2 D P s.b. [i, j] and the ideal, expected 2 D s.b. matrix is [r = 2 : 0.063, r = 8 : 0.541] and that between 2 D Ψ s.b. [i, j] and 2 D s.b. is [r = 2 : 0.048, r = 8 : 0.372], showing that the Euclidean distance between the two is on the same order or significantly smaller than the distance between each and the ideal. If anything, the direct wavefunction methodology seems to provide a slight mitigation of errors, likely due to the smaller number of circuits run as only one circuit per trial for the direct wavefunction is necessary while multiple circuits per trial must be run in order to obtain the expectation values of the Pauli matrices.
Note that to decrease computational expenseespecially in terms of the isolation of the numberconserving substates as explained in the following section-in this paper, we implement the tomography via the direct computation of the wavefunction.

Appendix C. Isolation of Number-Conserving Components
The number-conserving substates are eigenfunctions of the number operator, where Z q is the Pauli Z gate corresponding to qubit q. Therefore, each substate is composed of a definite number of particles (i.e., qubits in the |1 state) such that, for example, an overall r = 6 AGP wavefunction which has contributions from basis states with N = 0, N = 2, N = 4, and N = 6 particles can be decomposed into the following substates that are eigenfunctions of the number operator and hence have a definite number of particles: (31) and In general, the overall AGP wavefunction is constructed from its number-conserving substates according to The density matrices associated with the numberconserving substates ( 2 D r,N ) are able to be isolated as shown in the following sections such that the signature of ODLRO, λ D , can be attained for subblocks of the 2-RDM specific to each substate in addition to the signature of the overall non-particle-conserving AGP state with contributions from each of the number-conserving substates.

a. Tomography via Pauli Expectation Values
As can be extrapolated from Eq. (33), the full density matrix ( 2 D r,full ) obtained as shown using tomography via Pauli expectation values in Appendix B a can be represented as a sum of the density matrices for each of the individual substates ( 2 D r,N ) with elements given according to This relationship alone is insufficient to solve for the number-conserving density matrices in terms of the full density matrix. Instead, we introduce the number operator from Eq. (27) to define a new matrix with elements given by The elements of this new matrix can be obtained by again directly computing the expectation values of Pauli matrices in a manner directly analogous to that described in Appendix B a with the creation and annihilation operators in the basis of Pauli matrices being written as shown in Eqs. (15) and (16) and the number operator being written as shown in Eq. (27). Thus, each element can be expressed as so that the expectation values of four-and fivequbit Pauli expressions-such as X j Y j+1 X k+1 Y k and X j Y j+1 X k+1 Y k Z q -can be directly probed on the quantum device in order to compute the each element of the P 1 matrix. Each element of the P 1 matrix can additionally be represented as the following linear combination of numberconserving density matrices ( 2 D r,N ) Similarly, other matrices can be defined with elements which can be additionally be represented according to Obtaining 2 D r,full , P 1 , P 2 , etc. directly from probing expectation values of Pauli matrices and solving the system of linear equations described by Eqs. (34), (37), and (39) for the individual number-conserving density matrices ( 2 D r,N ) allows for the computation of the signature of ODLRO (λ D ) corresponding to each of the numberconserving substates.
As can be readily observed, this methodology for obtaining the number-conserving substates quickly becomes prohibitively expensive as system size is increased. As such, this study relies on the methodology presented in the following section.

b. Tomography via Direct Computation of the Wavefunction
Wavefunctions corresponding to each individual number-conserving substate can be determined from the probability information obtained from a quantum device in a manner analogous to how the full AGP wavefunction is prepared as described in Appendix B b. Specifically, for a specific possible, even value of N ∈ {0, 2, 4, . . . , r}, each element of the number-conserving wavefunction (|Ψ r,N i ) in the basis of the 2 r possible stated (|η i ) can be obtained according to whereby if the basis state |η i contains N particles (i.e., has N qubits in the |1 state), the element of the numberconserving wavefunction (|Ψ r,N i ) corresponding to that basis state is identical to that from the full wavefunction (|Ψ i ); however, if the basis state |η i doesn't contain N particles, |Ψ r,N i is set to zero. The resulting numberconserving wavefunction is then normalized to one.
The signature of ODLRO (λ D ) can then be obtained directly from the number-conserving wavefunction in the manner described in Appendix B b. Note that this projection can act as a form of error mitigation as contributions from bases corresponding to odd-numbered basis states-and indeed all basis states not corresponding to the number of particles of interest-are omitted. Additionally note that for QASM simulation, this methodology for isolating the number-conserving substates from the data obtained from the quantum device yields identical results within sampling error to the methodology presented in Appendix C a while being significantly less computationally expensive. As such, tomography via direct computation of the wavefunction is employed in this study.

Appendix D. Description of Noise on Near-Term Quantum Devices
Three main classes of errors lead to the deviation of physical qubits from the idealized logical qubits: namely, gate noise, readout noise, and decoherence. Quantum gate noise/error refers to a situation where the application of a unitary gateÛ to a quantum state |Ψ yields a result that deviates fromÛ |Ψ . This class of error is caused by either imprecisely calibrated control of the qubits and/or imperfect isolation of qubits from their environment, and the overall gate error increases with the number of gates applied. Readout noise/error refers to transmission line noise that makes the |0 state appear to be |1 or vice versa; it can be caused by the probability distributions of the measured physical quantities that correspond to the |0 and |1 states overlapping and/or the qubit decaying during readout. Decoherence involves interactions with external systems (vibrations, temperature fluctuations, electromagnetic waves, etc.) leading to the degradation of the quantum state prepared on quantum devices. As both the number of gates applied to a system-and hence gate noise-and decoherence tend to increase with system size (r), larger-scale quantum computations often involve more and more error [57]. See the Supplemental Information of Ref. 31 for a more thorough exploration of error on near-term quantum devices.

Appendix E. Analysis of Errors Via Joint Probabilities of Occupation
Comparing the results obtained from the simulated and experimental Melbourne data illustrates the error associated with the Melbourne device. The probability of a given orbital (i.e. qubit) being occupied if the orbital with index 0 (i.e., Q0) is occupied was computed for all combinations of total orbitals (qubits, r) and particles (|1 qubit states, N ), and the results for N = 2 are shown in Tab. I in order to comment on relative error based on system size (r). Note that all other joint probabilities for Melbourne and QASM simulation are given in Tabs. II and III, respectively.
As can be seen from the joint probability data in Tab. I, simulated results exactly demonstrate the orbital (qubit) pairing that we program into the system as the only possible two-particle orientation with Q0 being occupied should be to have Q1 simultaneously occupied. However, due to errors on Melbourne, this ideal behavior is not exactly recreated on the experimental quantum device. Specifically, the joint probability of occupying Q0 and Q1 is not unity and seems to decrease with increasing system size. Additionally, the joint probability of Q0 and Qi where i = 0, 1 should be zero as is seen in QASM simulation; however, the experimental data demonstrates that other double-excitations are contributing to the overall two-particle substate, indicating error in either state preparation or measurement. Overall, the error associated with noisy near-term quantum computers decreases the signature of ODLRO, indicating less ODLRO character for the experimentally-prepared states than predicted by QASM simulation. In order to best construct and probe entangled states on quantum computers, then, errors on real-world devices need to be minimized.

Appendix F. Experimental QASM Simulator and Quantum Device Specifications
Throughout this work, we have employed the ibmq qasm simulator [60] and the ibmq 16 melbourne [61] IBM Quantum Experience device, which are available online. The QASM simulator is a general-purpose simulator that emulates execution of quantum circuits in either an ideal manner (i.e., with only sampling error) or subject to highly-configurable noise modeling; in this study, all reported QASM results are ideal. The ibmq 16 melbourne device is composed of fixedfrequency transmon qubits with co-planer waveguide resonators [62,63]. Experimental calibration data and connectivity for this device-as well as that for other devices employed in obtaining supplemental data-is given in Tablels IV-VII.

Appendix G. State Preparation Fidelity
To provide a metric on which to judge the degree to which the expected state preparation was prepared on the quantum devices employed, we include the state preparation fidelity given by [64] where |Ψ ideal is the wavefunction corresponding to the result of applying the unitary obtained from the matrix form of the state preparation given in Eq. (12) to the all-zero initial state and is hence the ideal expected outcome for a given state preparation on a device with no error and where |Ψ exp. represents the wavefunction obtained from the quantum device. The state preparation fidelities for ibmq 16 melbourne and QASM simulation are reported in Tabs. II and III.

Appendix H. Additional Device Data
While only data from QASM simulation and the ibmq 16 melbourne quantum device are presented, additional data for ibmq 5 yorktown, ibmq santiago, and ibmq rochester are provided in Tabs. VIII-X.   II: All eigenvalue (λ D ) information for the non-number-conserving overall state (all) and the numberconserving substates are given with state preparation fidelities (F) and joint probability of occupation numbers of other orbitals (qubits) if the first orbital (Q0) is filled for ibmq 16 melbourne.   Version:

53.3
Qubit:    TABLE IX: All eigenvalue information for the non-number-conserving overall state (all) and the numberconserving substates are given with joint probability of occupation numbers of other orbitals (qubits) if the first orbital (Q0) is filled for ibmq santiago.
TABLE X: All eigenvalue information for the non-number-conserving overall state (all) and the number-conserving substates are given with joint probability of occupation numbers of other orbitals (qubits) if the first orbital (Q0) is filled for ibmq rochester.