Scalable Circuits for Preparing Ground States on Digital Quantum Computers: The Schwinger Model Vacuum on 100 Qubits

The vacuum of the lattice Schwinger model is prepared on up to 100 qubits of IBM's Eagle-processor quantum computers. A new algorithm to prepare the ground state of a gapped translationally-invariant system on a quantum computer is presented, which we call Scalable Circuits ADAPT-VQE (SC-ADAPT-VQE). This algorithm uses the exponential decay of correlations between distant regions of the ground state, together with ADAPT-VQE, to construct quantum circuits for state preparation that can be scaled to arbitrarily large systems. These scalable circuits can be determined using classical computers, avoiding the challenging task of optimizing parameterized circuits on a quantum computer. SC-ADAPT-VQE is applied to the Schwinger model, and shown to be systematically improvable, with an accuracy that converges exponentially with circuit depth. Both the structure of the circuits and the deviations of prepared wavefunctions are found to become independent of the number of spatial sites, $L$. This allows for a controlled extrapolation of the circuits, determined using small or modest-sized systems, to arbitrarily large $L$. The circuits for the Schwinger model are determined on lattices up to $L=14$ (28 qubits) with the qiskit classical simulator, and subsequently scaled up to prepare the $L=50$ (100 qubits) vacuum on IBM's 127 superconducting-qubit quantum computers ibm_brisbane and ibm_cusco. After introducing an improved error-mitigation technique, which we call Operator Decoherence Renormalization, the chiral condensate and charge-charge correlators obtained from the quantum computers are found to be in good agreement with classical Matrix Product State simulations.


I. INTRODUCTION
Quantum simulations of physical systems described by the Standard Model [1][2][3][4][5][6], and descendant effective field theories (EFT), are anticipated to provide qualitatively new predictions about matter under extreme conditions; from the dynamics of matter in the early universe, to properties of the exotic phases of quantum chromodynamics (QCD) produced at the LHC and RHIC (for overviews and reviews, see Refs.[7][8][9][10][11][12][13][14][15]).One of the major challenges facing quantum simulations of physical systems is the preparation of initial states on quantum computers that can be used to determine important quantities that are inaccessible to classical high-performance computing (HPC), i.e., the problem of state preparation.While simulating the dynamics of any given initial state is known to be efficient for an ideal quantum computer [16], residing in the BQP complexity class, preparing an arbitrary state generally requires quantum resources that asymptotically scale super-polynomially with increasing system size [17], residing in the QMA complexity class. 1 However, states of physical systems are not the general case, and are often constrained by both local and global symmetries. 2 In some instances, these symmetries allow observables to be computed by perturbing around states that can be efficiently initialized [9].In the foreseeable future, quantum simulations will be far from asymptotic in both system size and evolution time, and the resources required for both time evolution and state preparation will be estimated by direct construction and extrapolations thereof.Furthermore, successful quantum simulations will require specialized quantum circuits and workflows that are optimized for specific quantum hardware.
Many systems of physical interest, including QCD, have translational symmetry and possess an energy (mass) gap Λ between the unique ground state and first excited state.The gap defines a characteristic length scale of the system ξ = 1/Λ, and parameterizes the decay of the longest distance correlations in the ground state wavefunction, falling as ∼ e −r/ξ /r α for regions separated by r ≫ ξ, for some α.A natural way to encode a lattice QFT onto a register of a digital quantum computer is by identifying subsets of qubits (or qudits) with spatial points of the lattice that align with the connectivity of the quantum computer.A realization of the ground state on the register of a quantum computer should reflect the localized correlations between these subsets of qubits separated by r ≫ ξ [105,106].In the absence of topological order, one way to establish the ground state is to initialize the quantum register in a state without correlations between qubits, e.g., a tensor product state, and then systematically introduce correlations through the action of quantum circuits.A crucial point is that the localized correlations imply that the state preparation circuits need to have structure only for qubits spatially separated by r ≲ ξ [105,106].This is sufficient to obtain exponentially converged accuracy in the prepared state.Due to translational invariance, the ground state for an arbitrarily large lattice can be prepared by repeating these circuits across the entire register.
To study the dynamics of physically relevant systems in a quantitative way, with a complete quantification of uncertainties, simulations of large volumes of spacetime are typically required.Motivated by the discussion in the previous paragraph, we introduce Scalable Circuits ADAPT-VQE (SC-ADAPT-VQE), a new method for quantum state preparation that uses the hierarchies of length scales present in physical systems; see Fig. 1 for an illustration.In SC-ADAPT-VQE, quantum circuits that (efficiently) prepare a given state to a specified level of precision are determined on modest-sized lattices that are large enough to contain the longest correlation lengths.As long as ξ is not too large, these circuits can be determined using classical computers.This avoids the challenging task of optimizing circuits on a quantum computer with both statistical uncertainty and device noise [107,108].Once determined, (discrete) translation invariance is used to scale these circuits up to the full lattice.Since the quality of the prepared state becomes independent of the spatial lattice length L, with O(e −ξ/L ) corrections, this is a potential path toward quantum simulations of lattice QFTs that are beyond the capabilities of HPC.
In this work, SC-ADAPT-VQE is applied to the Schwinger model and is used to prepare the vacuum on up to 100 qubits on IBM's Eagle quantum processors.Underlying the development is the ADAPT-VQE algorithm [109] for quantum state preparation, which is modified to generate scalable circuits.After the necessary Trotterized circuits have been built, SC-ADAPT-VQE is performed using the qiskit classical simulator on system sizes up to L = 14 (28 qubits).It is found that both the energy density and the chiral condensate converge exponentially with circuit depth to the exact results.Importantly, both the quality of the prepared state and the structure of the associated circuits are found to converge with system size.This allows the state preparation circuits, determined on small lattices using classical computing, to be extrapolated to much larger lattices, with a quality that becomes independent of L. The scaled circuits are used to prepare the L ≤ 500 vacua using qiskit's Matrix Product State (MPS) circuit simulator, and to prepare the L ≤ 50 (100 qubits) vacua on the registers of IBM's superconducting-qubit quantum computers ibm brisbane and ibm cusco.An improved and unbiased error mitigation technique, Operator Decoherence Renormalization (ODR), is developed and applied to the quantum simulations to estimate error-free observables.The results obtained from both the MPS circuit simulator and IBM's quantum computers are found to be in excellent agreement with Density Matrix Renormalization Group (DMRG) calculations.
FIG. 1. Pictorial description of the SC-ADAPT-VQE algorithm.Once a pool of scalable operators { Ôi} has been identified, ADAPT-VQE is performed using classical computers to determine a quantum circuit (parameterized by {θi}) that prepares the vacuum up to a desired tolerance.ADAPT-VQE is repeated for multiple lattice sizes, {L1, L2, ...}, and the circuit parameters are extrapolated to the desired L, which can be arbitrarily large.The extrapolated circuits are executed on a quantum computer to prepare the vacuum on a system of size L.

II. THE LATTICE SCHWINGER MODEL
The Schwinger model [21] has a long history of study in the continuum and using numerical lattice techniques.In the continuum, it is described by the Lagrange density Electrically-charged fermions are described by the field operator ψ with mass m ψ , the electromagnetic gauge field by A µ with field tensor F µν , and the covariant derivative is defined as D µ = ∂ µ − ieA µ .It is the Hamiltonian lattice formulation, first developed and studied by Banks, Kogut and Susskind [110,111], that is relevant for quantum simulations.One feature of gauge theories in 1 + 1D, which distinguishes them from theories in higher dimensions, is that the gauge field is completely constrained by the distribution of fermion charges through Gauss's law.In axial gauge, the spatial gauge field is absent, and the effects of the time-component of the gauge field are included by non-local (Coulomb) interactions [81,95].With open boundary conditions (OBCs), using the staggered fermion discretization [110] of the electron field, and applying the Jordan Wigner (JW) transformation to map fermion field operators to spins, the Schwinger model Hamiltonian is (for a derivation, see, e.g., Ref. [23]) Here, L is the number of spatial lattice sites, corresponding to 2L staggered (fermion) sites, m and g are the (bare) electron mass and charge, respectively, and the staggered lattice spacing a has been set to one. 3"Physical" quantities are derived from the corresponding dimensionless quantities by restoring factors of the spatial lattice spacing.Even (odd) sites correspond to electrons (positrons), as reflected in the staggered mass term and charge operator. 4A background electric field can be included straightforwardly, equivalent to a θ-term, but will be set to zero in this work.
In the sector with vanishing total electric charge, Ĥel in Eq. ( 2) can be rewritten in a way that reduces the number of gates required in quantum circuits for time evolution, and is less demanding on device connectivity, see App. A. Due to the confinement, the low energy excitations are hadrons and the mass gap is given by Λ = m hadron .For our purposes, m hadron is defined to be the energy difference in the Q = 0 sector between the first excited state (single hadron at rest) and the vacuum.

A. Infinite Volume Extrapolations of Local Observables
Central to the development of state preparation circuits is the scaling of expectation values of local observables in the ground state, with both the correlation length ξ = 1/m hadron , and the volume L. Due to the exponential suppression of correlations in the ground state between regions separated by r > ξ, it is expected that, locally, the wavefunction has converged to its infinite volume form, with corrections of O(e −ξ/L ).As a result, expectation values of local observables will be exponentially converged to their infinite volume values.However, near the boundaries of the lattice, the wavefunction is perturbed over a depth proportional to ξ, causing local observables to deviate from their infinite volume values.Equivalently, boundary effects cause deviations in volume averages of local observables that are O(ξ/L).This scaling of observables is responsible for the SC-ADAPT-VQE prepared vacuum converging exponentially in circuit depth, and enables the circuits to be systematically extrapolated to larger system sizes.I and DMRG calculations (orange squares) given in Table V are extrapolated to L → ∞, as shown by the darker points, respectively.The solid lines correspond to linear extrapolations and the dashed lines correspond to quadratic extrapolations, and are found to overlap (see inset).The difference between the L → ∞ values of these two extrapolations defines the (fitting) uncertainties associated with the darker points.
Two quantities associated with the ground-state wavefunction (vacuum) that we focus on are the chiral condensate χ, and the energy density ε.The chiral condensate 5 is an order parameter of chiral symmetry breaking, and in the JW mapping is The energy density is defined as ε = ⟨ Ĥ⟩/L, and in axial gauge is not a local observable because the contribution from the electric-field term in the Hamiltonian, Ĥel , involves all-to-all couplings.However, this is an artifact of using axial gauge and enforcing Gauss's law.In Weyl gauge, with explicit (local) gauge degrees of freedom, the Hamiltonian is manifestly local, and therefore the energy density is a local observable.These quantities are computed for m = 0.5, g = 0.3 using exact diagonalization for L ≤ 14 (Table I) and DMRG for L ≫ 14 (Table V).As anticipated, a linear extrapolation in 1/L is found to be consistent with these results, as seen in Fig. 2. Additional details, along with results for m = 0.1 with g = 0.3 and g = 0.8, can be found in App.B.

III. SC-ADAPT-VQE FOR THE LATTICE SCHWINGER MODEL
Underlying SC-ADAPT-VQE is ADAPT-VQE [109], a quantum algorithm for state preparation that has been applied to spin models [113], systems in quantum chemistry [109,[114][115][116][117][118][119] and nuclear structure [120,121].It builds upon the Variational Quantum Eigensolver (VQE) [122], in which parameterized quantum circuits are optimized to minimize the expectation value of a Hamiltonian.The parameterized circuits are constructed step-wise (or equivalently in layers), where the incrementally-improved ansatz states converge to the ground state with successive iterations.At each step, the unitary transformation that maximally decreases the energy of the ansatz state is identified from a pre-defined set ("pool") of unitaries.The quantum circuit corresponding to this unitary is then appended to the state preparation circuit.The (initial) state from which the algorithm starts will often be chosen to be a tensor product or an entangled state that can be efficiently prepared on a quantum computer, such as a GHZ-state.If the target state is the ground state of a confining gauge theory, e.g., the Schwinger model, the strong-coupling (trivial) vacuum, can be a good choice for such an initial state as it has the correct long-distance structure in the gauge fields.The ADAPT-VQE algorithm can be summarized as follows: 1. Define a pool of operators { Ô} that are constrained to respect some or all of the symmetries of the system.
2. Initialize the register of the quantum computer to a strategically selected state, |ψ ansatz ⟩, with the desired quantum numbers and symmetries of the target wavefunction.
3. Measure the expectation value of the commutator of the Hamiltonian with each operator in the pool, ⟨ψ ansatz |[ Ĥ, Ôi ]|ψ ansatz ⟩.These are estimators of the associated decrease in energy from transforming the ansatz wavefunction by |ψ ansatz ⟩ → e iθi Ôi |ψ ansatz ⟩, for an arbitrary parameter θ i .
4. Identify the operator, Ôn , with the largest magnitude commutator with the Hamiltonian.If the absolute value of this commutator is below some pre-determined threshold, terminate the algorithm.If it is above the threshold, update the ansatz with the parameterized evolution of the operator |ψ ansatz ⟩ → e iθn Ôn |ψ ansatz ⟩.
6.Return to step 3. 5 In the continuum, the chiral condensate is defined as χcont = ⟨ψψ⟩, which on the lattice becomes χ lat = 1 L j ⟨ψ j ψ j ⟩, where j labels the spatial site.To have a positive quantity, we have added a constant to the definition of χ, χ ≡ χ lat + 1.This counts the average number of electrons and positrons on a spatial site.
For a given pool of operators, it is a priori unknown if this algorithm will furnish a wavefunction that satisfies the pre-determined threshold for the observable(s) of interest, but it is expected that the pool can be expanded on the fly to achieve the desired threshold.The systems that have been explored with this algorithm show, for a fixed pool, exponential convergence with increasing numbers of iterations [109,114,117,119,120].
Generally, different terms contributing to operators in the pool do not commute with each other.Constructing quantum circuits that exactly implement the exponential of a sum of non-commuting terms is challenging, and in practice approximations such as first-order Trotterization are used.This introduces (higher-order) systematic deviations from the target unitary operator in each case, and defines the pool of unitary operators, These Trotterized unitary operators correspond to the quantum circuits that are implemented in state preparation.In optimization of the quality of the state prepared on a given quantum computer, particularly a NISQ-era device, there are tradeoffs between the gate-depth of a particular circuit implementation, the coherence time, the errors associated with gate operations, and the associated Trotter errors.This is explored in App. C. Typically, ADAPT-VQE is a hybrid classical-quantum algorithm that evaluates matrix elements of the Hamiltonian in trial wavefunctions on a quantum computer, with parameters that are optimized classically.One disadvantage of this is that the evaluation of expectation values of the Hamiltonian requires a large number of measurements (shots) on quantum computers.A novel part of SC-ADAPT-VQE is the use of a classical simulator to determine the ADAPT-VQE state preparation circuits.As shown in Sec.V, these circuits can be scaled and used to prepare the vacuum on arbitrarily large lattices.

A. A Scalable Operator Pool for the Lattice Schwinger Model
A successful application of SC-ADAPT-VQE to the preparation of the lattice Schwinger model vacuum requires choosing an efficient and scalable pool of operators (first step in Fig. 1).These operators are used to systematically improve the ansatz vacuum wavefunction, and are (only) constrained to be charge neutral, symmetric under chargeconjugation and parity (CP) and, as a consequence of the CPT theorem [123][124][125], invariant under time reversal. 6 Ideally one wants to find the smallest pool of operators that is expressive enough to converge rapidly toward the vacuum.For a lattice with OBCs, the system has translational symmetry in the volume that is broken by the boundaries (surface).In the vacuum, the effects of the boundaries are expected to be localized, with penetration depths set by the mass gap.Therefore, the pool of operators should contain translationally invariant "volume" operators, and "surface" operators that have support only near the boundaries.In addition, a hierarchy is anticipated in which one-body operators are more important than two-body operators, two-body more important than three-body, and so on. 7Note that because wavefunctions are evolved with exp(iθ i Ôi ), arbitrarily high-body correlations are built from n-body operators (analogous to connected vs disconnected Feynman diagrams).For the Schwinger model, we observe that one-body operators are sufficient.
With the above discussion as guidance, it is convenient to define two classes of one-body operators, one containing volume operators, and the other containing surface operators: 6 In the total charge Q = 0 sector, there is a CP symmetry corresponding to the composition of a reflection through the mid-point of the lattice, exchanging spatial sites n ↔ L − 1 − n, and an interchange of an electron and a positron on each spatial site.In terms of spins on staggered sites this is realized as σi n ↔ σi Unlabelled Ẑs act on the qubits between the leftmost and rightmost operators (e.g., X0 Ẑ2 X3 = X0 Ẑ1 Ẑ2 X3 ).The first two operators in Eq. ( 7) are translationally invariant, ΘV m is the mass term in the Hamiltonian, and ΘV h (d) is a generalized hopping term that spans an odd-number of fermion sites, d, connecting electrons and positrons at spatial sites separated by ∆L = (d − 1)/2.Only d-odd operators are retained, as the d-even operators break CP.The second two operators in Eq. ( 7) correspond to surface terms, of the form of a mass-density and of a hopping-density at and near the boundaries.For ΘV h (d), d ∈ {1, 3, . . .2L − 3}, and for ΘS h (d), d ∈ {1, 3, . . .2L − 5}, preventing hopping between boundaries (which is found to improve convergence).
Time reversal symmetry implies that the vacuum wavefunction can be made real up to an overall phase.The SC-ADAPT-VQE ansatz is built from unitaries of the form e iθi Ôi , and furnishing a real wavefunction requires that the operators in the pool are imaginary and anti-symmetric.The operators in Eq. ( 7) are real and are therefore disqualified from being members of the pool.Instead, consider a pool comprised of their commutators,8 While the contributions to extensive quantities from the volume operators, ÔV , typically scale as O(L), the surface operators, ÔS , make O(1) contributions as they are constrained to regions near the boundaries. 9When acting on the strong-coupling vacuum, the exponential of an operator in the pool creates and annihilates e + e − pairs separated by distance d.As the operators that are being considered are one-body, the variational algorithm is essentially building a coupled cluster singles (CCS) state (see, e.g., Refs.[126,127]).

IV. SCALABLE QUANTUM CIRCUITS FROM CLASSICAL COMPUTING
Integral to the application of SC-ADAPT-VQE is performing ADAPT-VQE on a series of systems that are large enough to enable a robust scaling of the parameterized circuits.These scalable circuits can either be determined with classical computing, or by use of a smaller partition of a larger quantum computer.In this section, SC-ADAPT-VQE is implemented using the qiskit noiseless classical simulator [128,129].

A. Trotterized Quantum Circuits for the Scalable Operator Pool
As discussed above, implementing the unitary operators in the pool, i.e., Eq. ( 5), on classical simulators or quantum computers requires mapping them to sequences of quantum gates.For the individual terms in Eq. ( 8), we have chosen to do this using Trotterization.The optimal gate decomposition is less important for implementation using a classical simulator, but is crucial for successful simulations on a quantum computer.With the goal of using IBM's superconducting-qubit quantum computers [128,129], our circuit designs aim to minimize two qubit gate count and circuit depth and require only nearest-neighbor connectivity.
As can be seen in Eq. ( 8), each term in a given operator in the pool is of the form ( X Ẑd−1 Ŷ − Ŷ Ẑd−1 X) for some odd value of d.The construction of circuits implementing the corresponding unitary operators follows the strategy outlined in Ref. [130].First, consider the Trotterization of terms with d = 1, i.e., constructing a circuit corresponding to e iθ/2( X Ŷ ± Ŷ X) ≡ R ± (θ).There is a known 2-CNOT realization of this unitary operator [130], shown in Fig. 3a.For terms with d > 1, this circuit can be extended in an "X" pattern as shown in Fig. 3b and 3c for d = 3 and d = 5, respectively. 10Terms with larger d are constructed by extension of the legs of the "X".Compared with the traditional CNOT staircase-based circuits, there is a reduction by two CNOTs, and a reduction by a factor ×2 in CNOT-depth. 11owever, the primary advantage of these circuits is that they allow for an efficient arrangement of terms leading to cancellations among neighboring R + (± π 2 ) gates.As depicted in Fig. 4, this is made possible by arranging the circuit elements so that sequential terms are offset by d − 1 qubits, i.e., start on qubit {0, d − 1, 2(d − 1), . ..}.This allows the outermost gates to cancel (using the identity in the upper left of Fig. 4).Also, for d ≥ 5, the next layer should start (d − 1)/2 qubits below the previous one, as the circuit depth can be reduced by interleaving the legs of the "X".Further optimizations are possible by noting that distinct orderings of terms, while equivalent up to higher order Trotter errors, can have different convergence properties; see App. C.

B. Building Scalable State Preparation Quantum Circuits using SC-ADAPT-VQE with Classical Computing
In this section, SC-ADAPT-VQE is used to prepare approximations to the vacuum of the lattice Schwinger model on up to L = 14 spatial sites (28 qubits) using classical simulations of the quantum circuits developed in the previous section (second step in Fig. 1).I, and the sequencing of the corresponding Trotterized operators and the variational parameters are given in Table II.The corresponding results for m = 0.1, g = 0.3 (7 th step) and m = 0.1, g = 0.8 (6 th step) can be found in App.F.
In addition to the energy density and chiral condensate introduced in Sec.II A, the infidelity density, is also studied, where |ψ exact ⟩ is the exact vacuum wavefunction on a lattice with L spatial sites.An infidelity density that is constant in L corresponds to constant deviations in local observables evaluated in the prepared state.
To investigate the interplay between L and ξ = 1/m hadron , three sets of parameters are considered: m = 0.1, g = 0.3 (ξ L=14 = 2.6), m = 0.1, g = 0.8 (ξ L=14 = 1.3) and m = 0.5, g = 0.3 (ξ L=14 = 0.9).The ξ are determined with exact diagonalization, and are found to weakly depend on L. Note that increasing either m or g decreases the correlation length.To make systematically improvable predictions of observables from the QFT that emerges from a given lattice model, extrapolations to the continuum (lattice spacing to zero) and infinite-volume (L → ∞) limits must be performed.This requires that the relevant correlation length(s) are all much greater than the lattice spacing, ξ ≫ 1 in lattice units, but are well contained in the lattice volume, L ≫ ξ.We primarily focus on extrapolation to large lattices, and therefore only require L ≫ ξ.As a result, the parameter set m = 0.5, g = 0.3 is used as the primary example throughout this work.
The values of ε, χ and I L obtained at the 7 th step of SC-ADAPT-VQE with m = 0.5, g = 0.3 are given in Table I, while their deviations from the exact values are shown in Fig. 5, as a function of increasing number of SC-ADAPT-VQE steps.The corresponding numerical values obtained from the other parameter sets are presented in App.F. 12  As seen by their approximately linear behavior in the log-plots in Fig. 5, the error in each of these quantities decreases exponentially with algorithm step, indicating convergence to the target wavefunction.This exponential trend is demonstrated out to 10 steps, reaching a convergence comparable to the systematic errors introduced in the Lextrapolations below.This provides evidence that this choice of initial state and operator pool does not suffer from "barren plateaus" or local minima.For a given step in the algorithm, the error is seen to become independent of system size.This indicates that extrapolations of the circuits to arbitrarily large systems will have errors that are independent of L. As discussed above, it is expected that SC-ADAPT-VQE will converge more rapidly for systems with smaller correlation lengths.This is indeed seen in Fig. 5, where the correlation length decreases from left to right, while the convergence improves.Also included in Table I is the number of CNOTs per qubit in the SC-ADAPT-VQE circuit.It is seen to scale as a constant plus a subleading O(1/L) term, leading to an asymptotic value of 48 CNOTs per qubit.This scaling is due to there being (2L − d) terms in each volume operator.
The structure of the SC-ADAPT-VQE state preparation circuit and the corresponding variational parameters for m = 0.5 and g = 0.3 are given in Table II.Notice that initially localized operators are added to the wavefunction (small d), followed by increasingly longer-range ones, as well as surface operators.Systems with longer correlation lengths require larger d operators (e.g., compare Table II and Table VI), in line with previous discussions on the exponential decay of correlations for d > ξ.It is also seen that the surface operators become less important (appear later in the ansatz structure) for larger lattices.For example, as shown in Table II, the 5 th step of SC-ADAPT-VQE transitions from being a surface to a volume operator at L = 10 (causing the jump in convergence at the fifth step in the right column of Fig. 5).This is expected as they contribute O(1/L) to the energy density, whereas volume operators contribute O(1).
Importantly, Table II shows that the order of operators, and the corresponding variational parameters are converging with increasing system size (third step in Fig. 1).This is due to exponentially decaying correlations for d ≫ ξ, and it is expected that the variational parameters will also converge exponentially, once L is sufficiently large to contain ξ, and we assume the following form: Table II shows that this convergence sets in for L > 7, 13 and the variational parameters extrapolated to L = ∞ are given in the last row of Table II.These are used in the next section to initialize the vacuum on lattices up to L = 500.An example of extrapolating the variational parameters is shown in Fig. 6 for the parameter θ 1 , associated with ÔV mh (1).The exact results obtained for L ≤ 14 are well reproduced and extrapolated with the exponential functional form in Eq. (10). 14A more complete discussion of the parameter extrapolations, along with examples for m = 0.1 and g = 0.3 and for m = 0.1 and g = 0.8, can be found in App.D.  8).An entry of " -" means that the operator does not contribute.The bottom row corresponds to an extrapolation to L = ∞ as detailed in Eq. (10).

V. PREPARING THE VACUUM OF THE SCHWINGER MODEL ON LARGE LATTICES
The vacuum preparation circuits, determined for L ≤ 14 with SC-ADAPT-VQE using an exact (statevector) classical simulator, are scaled to prepare the vacuum on much larger lattices.These scaled circuits are used to prepare the vacuum on lattices of up to L = 500 (1000 qubits) using a classical MPS circuit simulator and up to L = 50 (100 qubits) using IBM's Eagle-processor quantum computers (fourth step in Fig. 1).We emphasize that this scaling requires no further optimization of the circuits.The chiral condensate and energy density determined from the classical simulator are found to be consistent with DMRG calculations.On the quantum computers, the chiral condensate and chargecharge correlators are measured to probe the quality of one-and two-qubit observables.The results are in agreement with those from the classical MPS simulator, within statistical uncertainties.II, are extrapolated to L = ∞ by i) use of a 3-parameter fit given in Eq. ( 10), as shown by the blue line, with an asymptotic value shown by the blue region, and by ii) the forming of effective-θ (orange diamonds) defined in Eq. (D3), with the maximum and minimum values shown as the orange shaded region.Results for large lattices with m = 0.5, g = 0.3 through 7 steps of SC-ADAPT-VQE using circuits scaled from L ≤ 14.The superscript "SC-MPS" denotes the results obtained from the scaled circuits using the qiskit MPS classical simulator, and the superscript "DMRG" denotes results obtained from DMRG calculations.

A. Classical Simulation
Very large quantum circuits that do not generate long-range entanglement can be efficiently simulated using the qiskit matrix product state classical simulator.Here it is used to simulate the preparation of the Schwinger model vacuum on L ≫ 14 lattices, applying the scalable circuits determined in the previous section from 7 steps of SC-ADAPT-VQE on L ≤ 14 lattices.The values obtained for the chiral condensate and energy density up to L = 500 are compared with DMRG results, and are presented in Table III.The deviations in the energy density (∼ 1 × 10 −4 ) and chiral condensate (∼ 1 × 10 −3 ) are in good agreement with what was found for smaller L; see Table I.This demonstrates that the systematic errors in the vacuum wavefunctions prepared with the scaled quantum circuits are (approximately) independent of L over this range of lattice volumes. 15The scaled circuits corresponding to m = 0.1, g = 0.3 and m = 0.1, g = 0.8 have also been used to successfully prepare the vacuum.However, due to the larger correlation lengths, MPS calculations with L ≳ 100 required excessive classical resources, and were not performed.See App.F for more details.
It is worth summarizing what has been accomplished in this work with classical simulations: • In Sec.II, the vacuum energy density and chiral condensate were determined exactly for L ≤ 14 (28 staggered lattice sites) using exact diagonalization, and for L ≤ 10 3 using DMRG.The results for L ≥ 9 were (consistently) extrapolated to L → ∞, with 1/L scaling.
• In Sec.IV, SC-ADAPT-VQE, based on the scalable operator pool determined in Sec.III, was performed on L ≤ 14 lattices.Intensive quantities were found to converge exponentially with circuit depth, and the errors in these quantities, as well as the structure of the state preparation circuits, were found to become independent of L. This enabled the variational parameters defining the state preparation circuits to be consistently extrapolated to arbitrarily large L.
• In this section, the quantum circuits corresponding to 7 steps of SC-ADAPT-VQE were scaled and applied to large lattices using the qiskit MPS circuit simulator.The deviations of the energy density and chiral condensate computed from these wavefunctions were found to be independent of L, i.e., consistent with L ≤ 14.
These main points indicate that the quantum circuits determined classically with SC-ADAPT-VQE can be used to prepare the vacuum of the Schwinger model on quantum computers at scale with a precision that is independent of system size. .This is compared with the expected results obtained from the qiskit MPS circuit simulator (black dashes).Averaging χj over all of the qubits (including at the boundaries) gives the chiral condensate presented in Table IV.The layout of the qubits used on the processor is shown on the right.These results were obtained by performing 150 Pauli twirls, each involving 8 × 10 3 shots for the physics circuits and the corresponding mitigation circuits.The blue icon in the upper right indicates that this calculation was done on a quantum device [138].
are used.Fewer steps equates to shallower circuits, and a preliminary study of the performance of the computer with more steps can be found in App.G.The variational parameters extrapolated to the chosen range of L for 2 steps of SC-ADAPT-VQE are given in Table XII in App.G.The large number of qubits and two-qubit gates involved in these simulations make error mitigation essential to obtain reliable estimates of observables.Specifically, this work uses readout-error mitigation (REM), dynamical decoupling (DD), Pauli twirling (PT), and decoherence renormalization.The qiskit Runtime Sampler primitive is used to obtain readout-corrected quasi-distributions via the matrix-free measurement mitigation (M3) from Ref. [132].Also included in the primitive is DD, which is used to suppress crosstalk and idling errors [133][134][135].Crucial to the error mitigation is decoherence renormalization [95,96,104,136], modified in this work for simulations on a large number of qubits, which we call Operator Decoherence Renormalization (ODR).Underpinning decoherence renormalization is PT [137], which turns coherent two-qubit gate errors into incoherent errors, which can be inverted to recover error-free expectation values.Unlike previous applications of decoherence renormalization, which assume a constant decoherence across the device, ODR estimates the decoherence separately for each operator.This is done by running a mitigation circuit, which has the same operator structure as the one used to extract the observables, but with the noise-free result being known a priori.We choose the state preparation circuits with the variational parameters set to zero for mitigation, and in the absence of noise this prepares the strong coupling vacuum, |Ω 0 ⟩ in Eq. ( 4).Naively, it could be expected that post-selecting results on states with total charge Q = 0 would eliminate the leading bit-flip errors [98].However, when post-selection is combined with ODR, which accommodates single-qubit decoherence, undesirable correlations between qubits are introduced.We find that performing both mitigation techniques (post-selection and ODR) degrades the quality of two-qubit observables, and post selection is not used in this work as it is found to be less effective.More details about ODR and post-selection can be found in App.E.
The local chiral condensate, χ j in Eq. ( 3), obtained from ibm cusco for L = 50 is shown in Fig. 7, where the subscript "j" denotes the qubit index. 16Three different sets of results (in different stages of error mitigation) are shown: with only DD applied (squares), with DD and PT applied (diamonds), and after ODR (circles).Looking at the results with only DD (squares), it is seen that the noise is not uniform across the device, signaling a significant contribution of coherent noise.After PT (diamonds), this coherent noise is averaged out, and is transformed into incoherent (depolarizing) noise, seen by the almost-constant shift of the results compared with the MPS simulation.Finally, ODR removes this shift by mitigating the effects of depolarizing noise.More details on the interplay between these methods can be found in App.G.
With the statistics and twirlings gathered, the 1σ uncertainties in each point are ∼ 15% of their mean, and each χ j is within 3σ of the MPS simulator result (the individual values of χ j can be CP averaged to reduce the uncertainty, as shown in Fig. 14 in App.G).It is expected that these uncertainties will reduce with increased statistics and twirlings.
L Qubits CNOTs χ (SC−IBM) before ODR χ (SC−IBM) after ODR χ (SC−MPS) χ (DMRG) for large lattices with m = 0.5, g = 0.3 using the scaled circuits from two steps of SC-ADAPT-VQE.The values before and after application of ODR are given in columns four and five.Column six gives results obtained from running the two step SC-ADAPT-VQE circuits on an MPS classical simulator (the noiseless result), while column seven gives the results from DMRG calculations.
Notice that the expected values of χ j deviate from the volume average for only a few qubits near the boundaries.This is because the boundaries perturb the wavefunction only over a few correlation lengths, leaving the rest of the volume unaffected.The chiral condensates for L = 14, 20, 30, 40 and 50 are given in Table IV.This is an average over the whole lattice, Eq. ( 3), and therefore the uncertainty decreases with increasing L due to increased sampling.Despite having smaller uncertainties, the results remain within 3σ of the MPS simulator result.Also given in Table IV is the number of two-qubit CNOT gates.The number of CNOTs is seen to grow linearly with L, without affecting the quality of the result, and 788 CNOTs over 100 qubits is well within the capabilities of the quantum computer.This is in line with other quantum simulations that have been performed with large numbers of qubits and CNOTs using IBM's quantum computers [139][140][141].This highlights the fact that it is not the total number of CNOT gates in the quantum circuit that is limiting the scale of simulations, but rather it is the number of CNOT gates per qubit.This, of course, assumes that the CNOT gates in a single layer of the circuit can be enacted in parallel.Due to this, increasing L actually improves volume-averaged quantities by ∼ 1/ √ L due to statistical averaging.In a similar vein, since scalable circuits repeat structures of size ξ many times over the whole lattice, the number of Pauli-twirls being sampled is effectively multiplied by L/ξ.
To further probe the quality of the prepared wavefunctions, correlations between electric charges on the spatial sites are considered.The charge on a spatial site is defined to be the sum of charges on the two associated staggered sites, Qk = Q2k + Q2k+1 , where k is an integer corresponding to the spatial site.Of particular interest are connected correlation functions between spatial charges,17 defined as These correlations decay exponentially for |j − k| ≳ ξ due to confinement and charge screening.Unlike the chiral condensate, which is a sum of single qubit observables, ⟨ Qj Qk ⟩ c is sensitive to correlations between qubits, i.e., requires measurement of ⟨ Ẑj Ẑk ⟩.The results from ibm cusco for L = 50 are shown in Fig. 8.The correlations are symmetric under j ↔ k, and only the lower-triangle of the correlation matrix is shown.Each measured value is within 3σ of the MPS simulator result, consistent with statistical fluctuations.Also shown in Fig. 8 are the spatial charge-charge correlations as a function of distance, averaged over the lattice volume, To reduce the effects of the boundaries, this sum omits the first and last two spatial lattice sites.As anticipated, this correlation function decays exponentially, with a characteristic length scale proportional to ξ = 1/m hadron . 18For d > 2, the correlations are consistent with zero within 2σ (note that the log scale distorts the error bars), and increased numbers of shots and twirlings are needed to distinguish additional points from zero.The local chiral condensate and charge-charge correlations corresponding to the other values of L are given in App.G.

VI. SUMMARY AND OUTLOOK
In this work, the vacuum of the lattice Schwinger model was prepared on up to 100 qubits of IBM's 127-qubit Eagleprocessor quantum computers, ibm brisbane and ibm cusco.This was accomplished with SC-ADAPT-VQE, an algorithm for identifying systematically improvable state preparation quantum circuits that can be robustly scaled to operate on any number of qubits The utility of scalable circuits relies on physically relevant systems often having a (discrete) translational symmetry, and a finite correlation length set by the mass gap.Together, these imply that the state preparation circuits have unique structure over approximately a correlation length [105,106], which is replicated across the lattice.The lattice Schwinger model with OBCs was chosen to explore these ideas as its vacuum has (approximate) translational invariance and, due to confinement, has a mass gap.By performing SC-ADAPT-VQE on a classical simulator, state preparation circuits for lattices of L ≤ 14 (28 qubits) were built from an operator pool containing both translationally invariant terms and ones localized to the boundaries.Exponential convergence in the quality of the prepared state with both system size and circuit depth enabled the extrapolation of circuits that can be scaled to arbitrarily large lattices.This methodology was successfully demonstrated by preparing the Schwinger model vacuum on up to 100 superconducting qubits of IBM's quantum computers.Both the charge-charge correlators and the chiral condensate were measured, and were found to agree with results from an MPS simulator, within statistical uncertainty.Vital to the success of these quantum simulations involving a large number of qubits was the development of an improved error mitigation technique, which we have called Operator Decoherence Renormalization (ODR).
Due to its generality, we expect that the scalable circuit framework embodied by SC-ADAPT-VQE can be applied to other gapped theories with translationally-invariant ground states.Of particular importance is QCD, for which the initialization of ground states for quantum simulations continues to be a daunting prospect.It is likely that many of the ideas used to construct efficient state preparation circuits for the Schwinger model can be applied to the initialization of the ground state of QCD.Of course, the operator pool that informs the state preparation circuits will be more diverse since the gauge field is no longer completely constrained by Gauss's law.Local quark-field operators, extended quark operators with associated gauge links, and closed loops of gauge links will need to be included in the pool.It is also expected that the variational parameters defining the ground-state preparation circuits will converge exponentially, once the simulation volume can completely contain the pion(s).
The utility of SC-ADAPT-VQE is that it provides a straightforward prescription for determining low-depth quantum circuits that prepare the ground state on systems of arbitrary size with only classical computing overhead.This not only allows for the quantum simulation of ground state properties, but will be important for future simulations of dynamics, where preparing the initial state is a crucial first step.Scalable circuits can likely be used to prepare single-and multi-hadron states, for example, a vector meson in the Schwinger model or a baryon in QCD.Once these states are initialized, they can be used to simulate scattering, electroweak processes and probe properties of dense matter.As an example, localized e + e − pairs on top of the Schwinger model vacuum could be prepared and evolved forward in time.Such localized distributions can have high energy components, whose propagation through the lattice leaves behind showers of particles.These processes probe the dynamics of fragmentation, confinement, and hadron production, and lead to long-range correlations that entangle distant regions of the lattice (see, e.g.,    V for L ≥ 9.The solid lines correspond to linear extrapolations and the dashed lines correspond to quadratic extrapolations, and are found to overlap (see the insets).The darker points show the L = ∞ extrapolated value, with an uncertainty determined by the difference between the linear and quadratic extrapolations.To initialize large quantum registers, the variational parameters defining the state-preparation quantum circuits need to be extrapolated with high precision.In volumes large enough to contain the longest correlation length, the variational parameters are expected to be exponentially close to their infinite-volume values.Therefore, we assume that the form of the volume dependence for practical purposes is that given in Eq. ( 10), and check the self-consistency of this form. 25While there could be a polynomial coefficient of the exponential, we find that this is not required.Fitting exponential functions can be challenging; however, with results over a sufficient range of L, algebraic techniques, such as effective masses, have proven useful in lattice QCD calculations to eliminate "uninteresting" parameters, while at the same time mitigating correlated fluctuations in measurements [151][152][153][154][155].
With the goal of initializing large lattices, it is the θ ∞ i that are of particular interest.Assuming the volume dependence given in Eq. (D1), it is useful to form four relations These relations can be combined to isolate θ ∞ i , providing an L-dependent "effective-θ ∞ i ", denoted as θ ∞ i,eff : For a sufficiently large set of results, θ ∞ i,eff (L) will plateau for large L if the functional form in Eq. (D1) correctly describes the results.This plateau can be fit by a constant, over some range of large L, to provide an estimate of θ ∞ i .This method is similar to using varpro (variable projection) in a multi-parameter χ 2 -minimization.
As an example, the results for θ ∞ 1 from a 3-parameter fit of θ 1 to Eq. (D1) are compared with a determination using θ ∞ 1,eff (L) from Eq. (D3).Results obtained with these two methods for m = 0.1, g = 0.3 and for m = 0.1, g = 0.8 are shown in Fig. 12.The result obtained from fitting a constant to θ ∞ 1,eff (L) is consistent with the asymptotic result  VIII are extrapolated to L = ∞ by i) use of a 3-parameter fit given in Eq. ( 10), as shown by the blue line and shaded region, with an asymptotic value shown by the gray region, and by ii) forming of effective-θ (orange diamonds) with the maximum and minimum values shown as the orange shaded region (where possible).

Appendix E: Operator Decoherence Renormalization (ODR)
To mitigate the effects of noise, the decoherence renormalization technique [95,96,104,136] is modified for use with larger systems.In its original form, decoherence renormalization assumes that each qubit decoheres at the same rate under a depolarizing noise channel.When working with a small number of qubits, this is a reasonable approximation, but for larger systems, it is necessary to consider the rate of decoherence of each qubit individually.After Pauli twirling, the qubit errors are well described by a Pauli error channel [156], which maps the N qubit density matrix to where Pi is a tensor product of Pauli operators ( Î, X, Ŷ or Ẑ) acting on N qubits, and the set of η i characterizes the error channel.It is important to understand the effect of this error channel on observables.Generic observables can be written as a sum over tensor products of Pauli operators, so it suffices to consider an observable, Ô, that is a tensor product of Pauli operators.Under a Pauli error channel, the measured (noisy) expectation value, ⟨ Ô⟩ meas , is given by Note that Pi Ô Pi = ± Ô, depending on whether or not Ô and Pi commute or anti-commute.Using this fact, the measured (noisy) expectation value ⟨ Ô⟩ meas , can be seen to be directly proportional to the predicted (noiseless) expectation value, ⟨ Ô⟩ pred = Tr( Ôρ), i.e., The ODR factor η O is, in general, distinct for each operator, and can be estimated by running a mitigation circuit that has the same structure as the physics circuit, but where ⟨ Ô⟩ pred is already known.In this work, the mitigation circuit was taken to be the state preparation circuit with variational parameters set to zero, which is the identity in the absence of noise.This mitigation circuit will have the same noise channel as the physics circuit provided that the noise is dominated by errors in the two-qubit gates and is independent of the single qubit rotation angles in the circuit.Without noise, the mitigation circuit prepares the strong coupling vacuum, where ⟨ Ô⟩ pred is known, and therefore η O can be computed.Once η O is determined, Eq. (E3) is used to estimate the value of the noiseless observable from the results of the physics circuits.
An added benefit of ODR is that it reduces the need for other error mitigation techniques.For example, readout errors are partially mitigated since the measured observables are affected by both gate and measurement errors.This is convenient as current measurement mitigation techniques require a large classical computing overhead.It also reduces the need for post-selection, which in our work could have been performed on states with total charge Q = 0.This post-selection removes single-qubit errors, but introduces further correlations between qubits.These correlations effectively increase the size of the single-qubit errors (making observables sensitive to errors anywhere on the register).This reduces the efficacy of the Pauli error model, making post-selection incompatible with ODR. 26 Another desirable feature of ODR is that it allows simulations to retain the results of a much larger fraction of the ensemble.This is because the probability of a single-qubit error increases with system size, and therefore much of the ensemble is lost with naive post selection.Further, such errors have little effect on local observables that are summed across the entire qubit register.

Appendix F: Additional Results From Classical Simulations
The results corresponding to Fig. 5 for m = 0.1, g = 0.3 are given in Table VI IX.The 6 th step of the algorithm is chosen for m = 0.1, g = 0.8 because the operator structure through L = 14 has converged, allowing a consistent extrapolation of the circuits to large L.This can be seen by comparing the operator structure in Table VIII (6 steps) and Table X (7 steps).An interesting observation is that the sum of parameters for a particular operator in the ansatz remains approximately unchanged when an additional insertion of the operator is added.For example, compare the sum of parameters for ÔV mh (1) between L = 8 and 9 in Table VI.Using the same method as for m = 0.5, g = 0.3 in Sec.IV, scalable circuits for m = 0.1, g = 0.3 and m = 0.1 and g = 0.8 were also determined.The results of running these circuits on qiskit's MPS simulator for m = 0.1, g = 0.3 and m = 0.1 and g = 0.8 are given in Table XI.Due to the longer correlation lengths for these parameters, it was not possible to go to L = 500 with the available computing resources.In these MPS simulations, qiskit's default settings were used, where the bond dimension increases until machine precision is achieved.The details of the qiskit MPS simulator can be found on the qiskit website [157].Again, the energy density and chiral condensate are found to have precision comparable to that found on smaller systems.This shows that, despite the longer correlation lengths for m = 0.1, g = 0.3 and m = 0.1, g = 0.8, it is still possible to accurately extrapolate the state preparation circuits to large lattices.Note that stabilization of operator ordering for the different m and g (see Tables II, VI and X) does not follow the hierarchy in correlation lengths.This is because larger ξ increases both the contribution of the volume ∼ e −d/ξ and surface ∼ ξ/L terms to the energy density.
To emphasize the advantage of performing SC-ADAPT-VQE using a classical simulator, we give an estimate of the number of shots required to perform SC-ADAPT-VQE on a quantum computer.For m = 0.5, g = 0.3, L = 14 performing 10 steps of SC-ADAPT-VQE required ∼ 6000 calls to the optimizer, in addition to about 500 evaluations of ⟨[ Ĥ, Ôi ]⟩ for pool operators Ôi .Each one of these calls required roughly 10 −3 precision in the measured observable, corresponding to about 10 6 shots on a noiseless device.Therefore, SC-ADAPT-VQE for L = 14 would require ∼ 10 10 shots on a noiseless device.Factoring in the effects of device noise would increase this estimate by at least a factor of 10, and probably close to a trillion shots would be required to perform SC-ADAPT-VQE on a quantum computer.This is infeasible on current hardware.In this appendix, we provide additional details about how our results are obtained from IBM's quantum computers, together with the additional figures not shown in Sec.V B. All measurements are performed on ibm brisbane (L ≤ 40) and ibm cusco (L = 50) by sending the state preparation circuits, with measurements in the computational (z) basis, via the qiskit Runtime Sampler primitive.The values of the variational parameters obtained from fitting to the exponential form in Eq. ( 10) for 2 steps of SC-ADAPT-VQE are given in Table XII.The different qubits used for each lattice size can be seen in the insets in Figs.7 and 13. χ j , obtained from ibm brisbane for L = 14, 20, 30 and 40, is shown in Fig. 13, and the charge-charge correlation functions are shown in Fig. 15.In Fig. 14, the CP symmetry relating χ j = χ 2L−1−j is used to effectively double the number of shots, resulting in statistical error bars that are smaller by a factor of √ 2.
In an effort to explore the limitations of the quantum computer, the 3-step SC-ADAPT-VQE state preparation circuits for L = 30 and L = 50 were implemented on ibm brisbane and ibm cusco, respectively.The structure of the ansatz wave function and corresponding variational parameters can be found in Table XII.The local chiral condensate and charge-charge correlators obtained from 80 (L = 30) and 40 (L = 50) twirled instances, each with 8 × 10 3 shots, are shown in Figs.16 and 17.Despite the factor of three increase in the number of CNOTs relative to 2 layers (1254 versus 468 for L = 30, and 2134 versus 788 for L = 50), the results are consistent with those obtained from the qiskit MPS circuit simulator.Note that qubit 0 and 2 have decohered for both volumes, and in principle could be removed from volume averaged quantities, such as the chiral condensate.By sending the circuits with the Sampler primitive, several error mitigation techniques are applied during runtime, as mentioned in Sec.V B. Specifically, the readout mitigation technique used (for L ≤ 40) is M3 [132].This method is based on correcting only the subspace of bit-strings observed in the noisy raw counts from the machine (which usually include the ideal ones plus those with short Hamming distance, introduced by the noise in the measurement), and using Krylov subspace methods to avoid having to compute (and store) the full assignment matrix.
Unlike the other works that have utilized ≥ 100 superconducting qubits [139][140][141], which used zero-noise extrapolation (ZNE) [158][159][160] in conjunction with probabilistic error correction (PEC) [159,161] to remove incoherent errors, we use Operator Decoherence Renormalization (ODR), as explained in App.E. Both methods require first transforming coherent errors into incoherent errors, which is done via Pauli twirling.However, the overhead in sampling using ZNE and PCE, compared with ODR, is substantial.For ZNE, one has to add two-qubit gates to increase the noise level, and then perform an extrapolation to estimate the noiseless result.In the minimal case, this leads to running only another circuit, like in ODR, but with a circuit depth that is three times as large as the original circuit (e.g., replacing each CNOT with 3 CNOTs).However, this leads to a large uncertainty in the functional form of the extrapolation, and ideally the circuit is run with multiple noise levels to have multiple points from which to extrapolate.For PEC, the overhead is even larger, as it involves learning the noise model of the chip, by running multiple random circuits with different depths (see Ref. [161]).For ODR, as explained in App.E, only the same "physics" circuits are run, but with all rotations set to zero, meaning the sampling overhead is only doubled.
To generate the different twirled circuits, the set of two-qubit Pauli gates G 2 and G ′ 2 that leave the (noisy) two-qubit gate invariant (up to a global phase) must be identified.For the quantum processors used in this work, the native two-qubit gate is the echoed cross-resonance (ECR) gate, which is equivalent to the CNOT gate via single qubit rotations.Explicitly, .

(G1)
Using the functions from the package qiskit research [162], together with the two-Pauli gate set shown in Table XIII, a total of 40 (150) twirled circuits for both mitigation and physics were generated for L ≤ 40 (L = 50), each with 8 × 10 3 shots.From Fig. 7, the effects of each error mitigation method can be seen.The first set of results shown are semi-raw, obtained directly from the quantum computer.They are not raw since DD is integrated into the circuits that are run on the machine (REM is also included for L ≤ 40).To check the effect that DD has, several runs were performed without it, and a degradation of the signal was visible when qubits were idle for long periods (the effects of not using DD were more evident when the deeper 3-step circuit was run).Regarding REM, while the final fully-mitigated results for L = 50 (no REM applied) and L ≤ 40 (REM applied) systems are similar in quality, a larger statistical sample for L = 50 was required to achieve an equivalent level of precision (2.4 × 10 6 vs 6.4 × 10 5 shots).The second set shows the effects of applying PT (the results for no Pauli twirling corresponded to one twirled instance).It is seen that all the coherent noise on the different qubits has been transformed into uniform incoherent noise.The last set shown is after ODR has been used to remove the incoherent noise.

m
FIG.2.L-extrapolations of the vacuum energy density ε (left panel) and chiral condensate χ (right panel) for m = 0.5 and g = 0.3.The results of exact diagonalization calculations for L ≥ 9 (blue circles) given in TableIand DMRG calculations (orange squares) given in Table V are extrapolated to L → ∞, as shown by the darker points, respectively.The solid lines correspond to linear extrapolations and the dashed lines correspond to quadratic extrapolations, and are found to overlap (see inset).The difference between the L → ∞ values of these two extrapolations defines the (fitting) uncertainties associated with the darker points.

10 − 2 3 10 − 2 m 3 10 − 2 m 14 FIG. 5 .
FIG.5.Deviations from the exact values of the energy density δε, chiral condensate δχ, and wavefunction infidelity density IL obtained with SC-ADAPT-VQE.The deviation in quantity "x" is defined as δx = |(x (aVQE) − x (exact) )/x (exact) |, where x (exact) denotes the exactly calculated value at the same L. Results are shown for L = 6 to L = 14 as a function of step number for m = 0.1, g = 0.3 (left column), m = 0.1, g = 0.8 (center column) and m = 0.5, g = 0.3 (right column).The numerical values for m = 0.5, g = 0.3 for the 7 th step (highlighted with the dashed box) are given in TableI, and the sequencing of the corresponding Trotterized operators and the variational parameters are given in TableII.The corresponding results for m = 0.1, g = 0.3 (7 th step) and m = 0.1, g = 0.8 (6 th step) can be found in App.F.

FIG. 6 .
FIG.6.Example of fitting the asymptotic L-dependence of a parameter defining the SC-ADAPT-VQE state-preparation circuit.The results for θ1, corresponding to evolving by ÔV mh (1) (blue circles), determined from classical simulations, for m = 0.5 and g = 0.3 given in TableII, are extrapolated to L = ∞ by i) use of a 3-parameter fit given in Eq. (10), as shown by the blue line, with an asymptotic value shown by the blue region, and by ii) the forming of effective-θ (orange diamonds) defined in Eq. (D3), with the maximum and minimum values shown as the orange shaded region.

B. Quantum Simulations on 100 FIG. 7 .
FIG. 7.Local chiral condensate χj for L = 50, as obtained from IBM's Eagle-processor quantum computer ibm cusco after different steps of error mitigation: DD (squares), PT (diamonds), and ODR (circles).This is compared with the expected results obtained from the qiskit MPS circuit simulator (black dashes).Averaging χj over all of the qubits (including at the boundaries) gives the chiral condensate presented in TableIV.The layout of the qubits used on the processor is shown on the right.These results were obtained by performing 150 Pauli twirls, each involving 8 × 10 3 shots for the physics circuits and the corresponding mitigation circuits.The blue icon in the upper right indicates that this calculation was done on a quantum device[138].

FIG. 8 .
FIG. 8. Left: Connected contributions to the spatial charge-charge correlation functions, ⟨ Qj Qk ⟩c, for L = 50 (the inset shows the number of standard deviations by which the results obtained from ibm cusco deviate from the MPS simulator results).Right: Volume averaged correlation functions as a function of distance d, ⟨ Q Q⟩c(d), with the points following the same color map as in the left main panel (error bars show 1σ standard deviations).

FIG. 10 .
FIG. 10.L-extrapolations of the vacuum energy density ε (top) and chiral condensate χ (bottom) for m = 0.1 and g = 0.3 (left) and m = 0.1 and g = 0.8 (right).Each panel shows extrapolations of the exact values given in Tables VII and IX (blue circles) and of the results from DMRG calculations (orange squares) given in TableVfor L ≥ 9.The solid lines correspond to linear extrapolations and the dashed lines correspond to quadratic extrapolations, and are found to overlap (see the insets).The darker points show the L = ∞ extrapolated value, with an uncertainty determined by the difference between the linear and quadratic extrapolations.

FIG. 12 .
FIG. 12.Further examples of fitting the asymptotic L-dependence of the variational parameters defining the state-preparation quantum circuit, determined from classical simulations using SC-ADAPT-VQE.The results for θ1 = ÔV mh (1) (blue points) for m = 0.1 and g = 0.3 (left panel) given in Table VI and m = 0.1 and g = 0.8 (right panel) given in Table VIII are extrapolated to L = ∞ by i) use of a 3-parameter fit given in Eq. (10), as shown by the blue line and shaded region, with an asymptotic value shown by the gray region, and by ii) forming of effective-θ (orange diamonds) with the maximum and minimum values shown as the orange shaded region (where possible).

FIG. 13 .FIG. 16 .
FIG.13.Expectation values of χj for L = 14, 20, 30 and 40 (top to bottom), obtained from simulations using ibm brisbane.They are compared with the expected results obtained by using qiskit's MPS circuit simulator (black dashes).Averaging χj over all of the qubits provides the chiral condensates presented in TableIV.The layouts of the qubits used on the chip are shown in the insets.

FIG. 17 .
FIG. 17. Results for the L = 50 system obtained with use of three steps of SC-ADAPT-VQE, obtained from simulations using ibm cusco with 40 twirled instances.The top panel shows χj, the middle panel shows the CP averaged χj, and the bottom panels show the connected contribution to the spatial charge-charge correlation functions, ⟨ Qj Qk ⟩c (the first two spatial sites are not shown due to the errors on qubits 0 and 2), and the averaged correlation functions as a function of distance d, ⟨ Q Q⟩c(d), with the points following the same color map as in the left main panel (error bars show 1σ standard deviations).

TABLE I .
Energy density, chiral condensate and wavefunction infidelity for the vacuum of the Schwinger model with m = 0.5, g = 0.3.Both the results obtained from 7 steps of the SC-ADAPT-VQE (aVQE) algorithm using qiskit's classical simulator and the exact values are given.The last column shows the number of CNOTs/qubit in the state preparation circuit.

TABLE II .
Structure of the ansatz wavefunction with m = 0.5 and g = 0.3 through 7 steps of the SC-ADAPT-VQE algorithm obtained from a classical simulation using qiskit.For a given L, the order that the operators are added to the ansatz is displayed from left to right, with the associated parameter, θi, given as the entry in the table.The operators, ÔV mh (d h ) and ÔS mh (dm, d h ), are defined in Eq. (

TABLE IV .
Chiral condensate in the Schwinger model vacuum obtained from ibm brisbane (L ≤ 40) and ibm cusco (L = 50)

TABLE V .
Additional results for the energy density ε and chiral condensate χ obtained from DMRG calculations, and used in the extrapolations in Figs.2 and 10.Appendix D: Volume Extrapolations of the SC-ADAPT-VQE Variational Parameters: An "effective-θ ∞ i " and Table VII, and for m = 0.1, g = 0.8 are given in Table VIII and Table

TABLE VIII .
Same as TableIIexcept with m = 0.1 and g = 0.8 and through 6 steps of the SC-ADAPT-VQE algorithm.

TABLE IX .
Same as TableIexcept with m = 0.1 and g = 0.8 and through 6 steps of the SC-ADAPT-VQE algorithm.

TABLE X .
Same as TableIexcept with m = 0.1 and g = 0.8 and through 7 steps of the SC-ADAPT-VQE algorithm.

TABLE XI .
Same as TableIIIexcept with m = 0.1, g = 0.3 and m = 0.1, g = 0.8.Appendix G: Additional Details and Results From Simulations using IBM's Quantum Computers