Hybrid Oscillator-Qubit Quantum Processors: Simulating Fermions, Bosons, and Gauge Fields
Eleanor Crane Kevin C. Smith Teague Tomesh Alec Eickbusch John M. Martyn Stefan Kühn Lena Funcke Michael Austin DeMarco Isaac L. Chuang Nathan Wiebe Alexander Schuckert Steven M. Girvin








Department of Physics, Co-Design Center for Quantum Advantage, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Joint Quantum Institute and Joint Center for Quantum Information and Computer Science,
University of Maryland and NIST, College Park, Maryland 20742, USA

Brookhaven National Laboratory, Upton, New York 11973, USA
Yale Quantum Institute, PO Box 208 334, 17 Hillhouse Ave, New Haven, CT 06520-8263, USA
Departments of Applied Physics and Physics, Yale University, New Haven, CT 06511, USA

Department of Computer Science, Princeton University, Princeton, NJ 08544, USA
Infleqtion, Chicago, IL 60604, USA

Present address: Google Quantum AI, Santa Barbara, CA
Yale Quantum Institute, PO Box 208 334, 17 Hillhouse Ave, New Haven, CT 06520-8263, USA
Departments of Applied Physics and Physics, Yale University, New Haven, CT 06511, USA

Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
The NSF AI Institute for Artificial Intelligence and Fundamental Interactions

CQTA, Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany

Transdisciplinary Research Area “Building Blocks of Matter and Fundamental Interactions” (TRA Matter) and Helmholtz Institute for Radiation and Nuclear Physics (HISKP),
University of Bonn, Nussallee 14-16, 53115 Bonn, Germany
Department of Physics, Co-Design Center for Quantum Advantage, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
The NSF AI Institute for Artificial Intelligence and Fundamental Interactions

Brookhaven National Laboratory, Upton, New York 11973, USA
Department of Physics, Co-Design Center for Quantum Advantage, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

Department of Physics, Co-Design Center for Quantum Advantage, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

Department of Computer Science, University of Toronto, Canada
Pacific Northwest National Laboratory, Richland WA, USA
Canadian Institute for Advanced Studies, Toronto, Canada


Joint Quantum Institute and Joint Center for Quantum Information and Computer Science,
University of Maryland and NIST, College Park, Maryland 20742, USA


Yale Quantum Institute, PO Box 208 334, 17 Hillhouse Ave, New Haven, CT 06520-8263, USA
Departments of Applied Physics and Physics, Yale University, New Haven, CT 06511, USA

†Corresponding authors: emc2@mit.edu, aschu@umd.edu, steven.girvin@yale.edu



§ ABSTRACT  
We develop a hybrid oscillator-qubit processor framework for quantum simulation of strongly correlated fermions and bosons that avoids the boson-to-qubit mapping overhead encountered in qubit hardware. This framework gives exact decompositions of particle interactions such as density-density terms and gauge-invariant hopping, as well as approximate methods based on the Baker-Campbell Hausdorff formulas including the magnetic field term for the U(1) quantum link model in (2+1)D. We use this framework to show how to simulate dynamics using Trotterisation, perform ancilla-free partial error detection using Gauss's law, measure non-local observables, estimate ground state energies using a oscillator-qubit variational quantum eigensolver as well as quantum signal processing, and we numerically study the influence of hardware errors in circuit QED experiments. To show the advantages over all-qubit hardware, we perform an end-to-end comparison of the gate complexity for the gauge-invariant hopping term and find an improvement of the asymptotic scaling with the boson number cutoff S from 𝒪(log(S)^2) to 𝒪(1) in our framework as well as, for bosonic matter, a constant factor improvement of better than 10^4. We also find an improvement from 𝒪(log(S)) to 𝒪(1) for the U(1) magnetic field term.  While our work focusses on an implementation in superconducting hardware, our framework can also be used in trapped ion, and neutral atom hardware. This work establishes digital quantum simulation with hybrid oscillator-qubit hardware as a viable and advantageous method for the study of qubit-boson models in materials science, chemistry, and high-energy physics.




Hybrid Oscillator-Qubit Quantum Processors:Simulating Fermions, Bosons, and Gauge Fields
    Steven M. Girvin0000-0002-6470-5494^†
    September 9, 2024
========================================================================================






  
  





§ INTRODUCTION

Underlying many mechanisms in the natural world is the interaction between fermions and bosons, such as phonons mediating electron attraction, leading to superconductivity or the down-conversion of light during photosynthesis. Examples of fermionic particles include electrons, holes, fermionic atoms, baryons, protons, neutrons, quarks, leptons, and examples of bosonic particles include photons, phonons, excitons, plasmons, Cooper pairs, bosonic atoms, mesons, W-bosons, Z-bosons, Higgs bosons, and gluons. When the interactions among these particles is strong, analytical and numerical methods fail, leading to a large range of open problems in high energy physics, chemistry and condensed matter physics. In particular, quantum Monte-Carlo simulations are effective only for static equilibrium properties (or imaginary-time dynamics) and then only in special cases where there is no fermion sign problem. Tensor network methods work well only for systems with limited entanglement.  Some of the hardest examples to simulate can be found in the lattice gauge theory (LGT) formulation of the standard model of particle physics <cit.>. LGTs constitute an especially hard problem because, further to containing fermionic or bosonic matter, they also contain gauge fields which require the implementation of N-body interactions with N>2 such as the gauge-invariant kinetic energy (N=3) and the plaquette magnetic field term (N=4).

Quantum simulators <cit.> offer an attractive alternative to study these problems as they can naturally host entangled states which can be strongly correlated <cit.> and they can be used to study real-time equilibrium <cit.> and non-equilibrium dynamics <cit.>. Analog simulations of LGTs have been proposed <cit.> and carried out using, for example, cold atom <cit.>, ion trap <cit.> and donors in silicon <cit.> platforms, but analog quantum simulators are limited in terms of the Hamiltonians that can be implemented and observables that can be measured. In particular, the multi-body interactions appearing in LGTs makes their analog simulation challenging. Separately, while proof-of-principle digital simulations of LGTs have been performed, the overheads encountered in mapping bosons and fermions to qubits make these problems challenging to implement on qubit-based discrete-variable (DV) hardware <cit.>. While digital qudit-based quantum processors offer a promising approach toward efficiently encoding the large Hilbert space of bosons and gauge fields <cit.>, the qudit dimensions thus far explored offer only a slight increase in Hilbert space dimension over qubits. More importantly, neither qubit nor qudit systems natively implement the bosonic operations necessary for creating bosonic Hamiltonian terms – in particular, synthesis of the occupation dependent square-root factors needed for bosonic raising and lowering operators in qubit-based encodings leads to a large overhead for simulating bosonic models <cit.>. Simulation of purely bosonic (1+1)D theories using continuous-variable hardware has been proposed, for example using microwave cavities <cit.> or photonic devices <cit.>, but such an approach does not easily facilitate the study of models involving interactions between bosons and fermionic matter. 

In this work, we investigate the promise of hybrid oscillator-qubit quantum processors <cit.> for solving the above challenges. In particular, we leverage this hardware to realize a straightforward mapping between model and computational degrees of freedom, as illustrated in Fig. <ref>. By encoding both the bosonic matter and gauge fields in native oscillator modes, we avoid the costly overheads incurred by boson-to-qubit encodings. More specifically, the primary goal of our approach is to explore the advantage given by the native availability of bosonic gate sets, a crucial component for simulation of Hamiltonians containing bosons, yet costly to implement in qubit-based hardware. While the universal control of oscillator modes is challenging due to the equal level spacing of linear oscillators, here this issue is resolved by leveraging a hybrid oscillator-qubit architecture with non-linearity provided by the qubits. To that end, we build on the universal oscillator-qubit gate sets recently developed in a companion paper <cit.>. In particular, we focus on a circuit quantum electrodynamics (cQED) architecture with high-Q 3D cavities, although many of the hybrid oscillator-qubit operations and techniques we discuss can be implemented in other platforms such as ion traps <cit.>, where a similar oscillator-qubit approach for gauge theories has been proposed <cit.>, or neutral atoms in tweezer arrays <cit.>.


We focus on the task of simulating dynamics, i.e., evolving an initial state under an in-general time-dependent Hamiltonian Ĥ(t), using a digital approach in which we compile evolution under that Hamiltonian into elementary, native oscillator-qubit operations. Many quantum algorithms can be formulated with (controlled) time evolution under the Hamiltonian <cit.> as a fundamental subroutine. The goal of our work is to find these compilations for the most important models in quantum many-body physics involving fermions, bosons and gauge fields and to show that this decomposition leads to advantages compared to compilations using qubit-only hardware.    
 

This work is structured as follows. We start by reviewing elementary concepts (Sec. <ref>). We then present our main results on compilation strategies. This is divided into fundamental primitives (Sec. <ref>), methods for matter fields (Sec. <ref>), and finally gauge-fields and their interactions with matter (Sec. <ref>). We then perform an end-to-end resource estimation of our methods and compare to all-qubit hardware (Sec. <ref>).  Finally, we introduce methods to measure observables, algorithms to perform dynamics simulation and ground-state preparation, and numerical benchmarks (Sec. <ref>), and conclude (Sec. <ref>).

As a guide to the reader, here is some more detail about each of the seven sections of this paper, in sequence.  In Sec. <ref>, i.e., the remainder of this introductory section, we give context by broadly introducing LGTs (Sec. <ref>) and the challenges associated to their simulation. We choose paradigmatic models <cit.> as a case study, the and the , which we later use in Sec. <ref> to benchmark oscillator-qubit algorithms. We then introduce the superconducting hybrid oscillator-qubit hardware comprising transmon qubits coupled to high-Q microwave cavity resonators on which we base our study (Sec. <ref>). Finally, we introduce Trotter algorithms (Sec. <ref>), which are the basic subroutines used in our oscillator-qubit algorithms.



In Sec. <ref>, we introduce a toolset of hybrid oscillator-qubit compilation strategies used throughout this work. We review previously introduced methods, but also develop exact methods for synthesizing bosonic parity-dependent and density-density operations, as well as oscillator-mediated qubit-qubit entangling gates that we extend to multi-qubit gates. 

In Sec. <ref> we leverage this toolset to develop compilation strategies for implementing the time evolution operator for common bosonic (Sec. <ref>) and fermionic  (Sec. <ref>) Hamiltonian terms, including single-site potentials, nearest-neighbor and onsite interactions, and hopping terms in hybrid superconducting hardware. We consider both (1+1)D and (2+1)D. For bosonic matter, we utilize the direct mapping of a bosonic matter site to a single mode of a microwave cavity. For fermionic matter, we leverage a Jordan-Wigner mapping using either transmon qubits or dual-rail qubits in microwave cavities.

In Sec. <ref>, we present compilation strategies for implementing pure gauge and gauge-matter Hamiltonian terms, specializing to the case of ℤ_2 and U(1) gauge fields in (2+1)D. 

In Sec. <ref> we introduce explicit Trotter circuits for the and and perform an analysis of the Trotter errors. From this, we derive the end-to-end asymptotic scaling of the number of gates required for a dynamics simulation, including the errors encountered in the compilation. We then perform an explicit gate count comparison between our oscillator-qubit and an all-qubit simulation of the bosonic hopping Hamiltonian. To do so, we develop a qubit algorithm in the Fock-binary encoding to simulate bosonic hopping. 

In Sec. <ref>, we adapt algorithms for ground-state preparation and dynamics to our oscillator-qubit approach, utilising the compilation methods from the previous sections. We use the (1+1)D and as example models for numerical benchmarks using Bosonic Qiskit <cit.>. In particular, we develop an oscillator-qubit Variational Quantum Eigensolver (VQE) (Sec. <ref>, <ref>) ansatz. Within the VQE approach and in the context of the model, we investigate the effects of hardware and shot noise on VQE ground-state preparation (Sec. <ref>) and demonstrate the utility of post-selection upon Gauss'ss law in improving ground state preparation. Following this, we show how to measure observables in these ground states, in particular non-local observables such as string order correlators and the superfluid stiffness (Sec. <ref>). Finally, we discuss ground state preparation with quantum signal processing (QSP) and compare the error in ground state preparation as a function of the time to solution between QSP and VQE (Sec. <ref>).

In Sec. <ref> we conclude with a summary of our main results, a discussion of the prospects for quantum advantage with our approach, and an outlook on some of the future research directions that are opened up by this work.

Before continuing, we briefly comment on vocabulary and notation used throughout this work. We will use the terms bosonic mode, qumode, oscillator, resonator, and cavity mode interchangeably. In general, we will refer to certain hybrid oscillator-qubit gates as being `controlled' meaning that the generator of the gate contains either a qubit or a mode number-operator, i.e. the phase is proportional to the number of excitations in the state (the gate does nothing to |0⟩). We generally call gates `conditioned' when they obtain a state-dependent phase, i.e. the generator contains a qubit Ẑ or mode Fock-state dependent phases, which could for example be implemented using a projector.



 §.§ Lattice gauge theories

In the following, we briefly introduce the LGTs we will use as paradigmatic examples, a and a , which we will show how to compile onto oscillator-qubit hardware in the rest of the study. Here, we introduce the Hamiltonian terms of these models in order to prepare our discussion on how to implement them. In Sec. <ref>, we will introduce their physics in the context of our numerical experiments, especially also emphasizing past work. Note that these models only serve as examples; due to the digital nature of our methods, our study enables the simulation of a wide range of models involving bosons, fermions and gauge fields.

In lattice gauge theories, matter fields reside on sites and gauge fields reside between sites, usually (for the case of 2+1-D) in a square lattice geometry.

. In this model, the gauge fields are spin 1/2 degrees of freedom. The Hamiltonian is written as follows <cit.>:

    Ĥ_ℤ_2   = - g ∑_⟨ i,j⟩X̂_i,j -B ∑_i,j,k,l∈□Ẑ_i,jẐ_j,kẐ_k,lẐ_l,i
       - J ∑_⟨ i,j⟩( m̂^†_iẐ_i,jm̂_j + h.c.) + U ∑_in̂_i^2,

where i,j indicates the link between lattice site i and j on a square lattice, ∑_⟨ i,j ⟩ indicates a sum over all combinations of i,j such that i,j are nearest-neighbours  and  is a plaquette of four sites (see Fig. <ref>a). The annihilation operator at matter site i is given by m̂_i and the density or occupation by n̂_i =m̂^†_i m̂_i. The matter can be either bosonic (m̂=â) or fermionic (m̂=ĉ), with corresponding commutation/anticommutation relations [â_i,â_j^†]=δ_ij and {ĉ_i,ĉ_j^†}=δ_ij. X̂_i,j, Ẑ_i,j are Pauli operators that act on the gauge field at the link between sites i and j. The first, second, and third terms in Eq. (<ref>) represent the electric field (coupling g), magnetic field (coupling B), and gauge-invariant hopping (coupling J), respectively. The last term is an on-site interaction of strength U which is the simplest possible interaction term for bosonic matter, but does not contribute as an interaction term for (spinless) fermions since n̂^2=n̂, due to the Pauli exclusion principle.

The defining feature of gauge theories is invariance under local gauge transformations. This leads to the presence of Gauss's law dividing up the Hilbert space of the theory into sectors which are not coupled by the Hamiltonian. This constraint is in analogy to the relation between the integral of the electric field flux through a surface and the charge enclosed by the surface in classical electrodynamics, with a product of field operators replacing the integral. (-1)^n_i takes the role of the charge. In quantum mechanics, Gauss's law is formulated as a constraint on the eigenstates |Ψ⟩; for all i, Ĝ_i|ψ⟩ = q_i|ψ⟩, where q_i∈Z and

    Ĝ_i = ∏_⟨ i,j⟩X̂_i,j(-1)^n̂_i,

in (2+1)D, where the product is over the four links ⟨ i,j⟩ that connect to site i, see Fig. <ref>a and c. We choose q_i=1 throughout this work. For open boundary conditions, terminated by matter sites, the Gauss's law constraints on the boundary are fixed by imagining virtual qubit sites around the system and fixing them in a X̂ product state. We set X̂=1 for these qubits when necessary.

. The U(1) LGT, which is equivalent to quantum electrodynamics in the limit of vanishing lattice spacing, hosts continuous gauge fields. The quantum link formalism <cit.> replaces the gauge fields with discrete spin-S degrees of freedom with commutation relations [Ŝ^z_i, Ŝ^+_j]=δ_ijŜ^+_i and (Ŝ^z_i)^2+1/2(Ŝ^+_iŜ_i^-+Ŝ^-_iŜ_i^+)=S(S+1) <cit.>, where S is the spin length leading to a Hilbert space cutoff of 2S+1. The physics of continuous U(1) gauge fields are recovered for S→∞ <cit.>. Choosing (large) finite S, this encoding approximately preserves the bosonic algebra. The Hamiltonian in this formalism reads <cit.>

    Ĥ_U(1)   = g^2/2∑_⟨i,j|⟩(Ŝ^z_i,j +τ/2π)^2  
       -1/4g^2(S(S+1))^2∑_i,j,k,l∈□(Ŝ^+_i,jŜ^+_j,kŜ^-_k,lŜ^-_l,i +h.c.)
         
        + J/2√(S(S+1))∑_⟨i,j|⟩(m̂^†_iŜ^+_i,jm̂_j +  h.c.) 
       + M∑_i(-1)^i n̂_i,

where the first term is the electric field term, the second term is the magnetic field term, and the third is the gauge-invariant hopping. g is the gauge-matter coupling strength and J is the hopping strength. In one spatial dimension, the constant background electric field τ corresponds to a topological term <cit.>. The last term is the staggered mass term with mass M, which comes from the mapping of the continuous Dirac fermion fields onto the lattice and is therefore only necessary for fermions <cit.>. We set the lattice spacing to unity here. Also note that in order to take the above Hamtiltonian to the continuum limit of QED, the phases in the hopping and staggered mass terms need to be slightly altered, c.f. <cit.>. This does not introduce any additional challenges in our implementation, so we choose the above convention for simplicity.

In addition, the physical states of the Hamiltonian, i.e. the gauge-invariant ones, have to fulfill Gauss'ss law: for all 𝐢, Ĝ_𝐢|ψ⟩  = q_𝐢|ψ⟩, where q_𝐢∈Z and 

    Ĝ_𝐢 = Ŝ^z_𝐢-𝐞_x+Ŝ^z_𝐢-𝐞_y - Ŝ^z_𝐢+𝐞_x-Ŝ^z_𝐢+𝐞_y - Q̂_𝐢.

To clarify the spatial structure of Gauss's law, we specify the lattice site as a two-dimensional vector 𝐢 in this paragraph only. 𝐞_x/y is a unit vector in the x/y direction of the lattice, respectively. In the above expression Q̂_𝐢 =  m̂^†_𝐢m̂_𝐢 - (1-(-1)^r_𝐢)/2 is the staggered charge operator, where r_𝐢 is the Manhattan distance of the site 𝐢 from the arbitrarily chosen origin. The fixed values q_𝐢 correspond to a choice of static background charge configuration, and for the rest of the paper we restrict ourselves to the sector of vanishing static charges, i.e., q_𝐢 = 0 for all 𝐢. 

In order to encode the spin S degrees of freedom representing U(1) gauge fields in the into bosonic modes, we use the Schwinger boson mapping. It consists of two bosonic modes, labelled a and b, which fulfill the constraint n̂^a_i,j + n̂^b_i,j = 2S.
The electric field operator in this encoding is represented by

    Ŝ^z_i,j = n̂^a_i,j - n̂^b_i,j/2,

and the electric field raising operator S^+_i,j becomes

    Ŝ^+_i,j   = â^†_i,jb̂_i,j
    Ŝ^-_i,j   = â_i,jb̂^†_i,j

 Compared to an encoding in just a single oscillator <cit.>, our encoding allows for a smaller occupation of the modes to represent electric field zero, reducing the impact of mode decay. 

 We show a possible encoding of this model in our architecture in Fig. <ref>d.



 §.§ Circuit QED platform


We now outline the main features of the platform we consider in the rest of the study – a lattice of microwave cavities coupled to superconducting transmon qubits (see Fig. <ref>b). In particular, we discuss the native Hamiltonian of this platform and the gates that can be implemented via microwave pulses.

Within the circuit QED framework, long-lived bosonic modes can be realized in three-dimensional superconducting microwave cavities with typical lifetimes on the order of one to ten milliseconds (ms) in aluminum <cit.>, Fig. <ref>a. Recently, niobium cavities have been used to improve coherence times, with observed single photon lifetimes on the order of 30 ms <cit.> to over 1 s <cit.>. Separately, recent work realising planar microwave resonators, which may be more scalable due to their smaller size and fabrication advantages, have reached 1 ms single photon lifetimes <cit.>. In this platform, it is common to utilize one mode per cavity, although recent work has demonstrated that multi-mode cavities are also feasible <cit.>.

Superconducting transmon qubits are routinely coupled to microwave cavities and used for control and readout via on-chip readout resonators <cit.>. The relaxation T_1 and dephasing times T_2 of transmons when coupled to microwave cavities are typically on the order of 250 μs <cit.>, with some recent coherence times in tantalum-based transmons reaching 500 μs <cit.>. When simulating operation in hardware in section <ref>, we take into account the dominant source of hardware error, qubit decay and dephasing, adopting a conservative value of T_1=T_2=200 μs if not stated otherwise.



The SNAIL (Superconducting Nonlinear Asymmetric Inductive eLement) <cit.> is an attractive choice  for realizing a programmable coupling between cavities. The SNAIL consists of three large junctions shunted by one small junction. When biased with a non-zero DC magnetic flux, the SNAIL (or SNAIL array) generates a three-wave mixing Hamiltonian that enables a parametrized beamsplitter interaction of strength g(t) between the cavities while minimizing unwanted nonlinear static interactions, such as the cross-Kerr <cit.>.

Our proposed architecture is an array of `oscillator-qubit pairs,' each comprising a high-Q mode of a cavity dispersively coupled to a transmon qubit (see Fig. <ref>b, top panel). These units are arranged into a 2D square lattice, where adjacent cavities are connected via a SNAIL, which can be used to implement a beamsplitter (see Fig. <ref>b, bottom panel). In the remainder of this work, we represent cavities (qubits) by circles (crosses). This symbolic representation will later prove useful for denoting the various possibilities for mapping systems of fermions, bosons and gauge fields onto this 2D hardware layout. 

While this architecture is physically large in size (each cavity is about 50 mm long <cit.>) and might therefore exhibit scaling limitations, planar resonators, while at present exhibiting shorter coherence times compared to 3D cavities, are smaller and therefore more scalable. Thus, the proposed architecture in Fig. <ref>b can also be implemented with planar resonators, and all results presented herein extend to this case (with the caveat that Hamiltonian parameters, loss rates and gate speeds are informed by recent experiments in 3D cavities, and would need to be adjusted to appropriately reflect the planar case). 


In this architecture, each transmon-cavity pair in the array is in the strong-dispersive coupling regime <cit.>. In the rotating-wave-approximation, the effective Hamiltonian of the 2D array can be written as:

    Ĥ = Ĥ_0 + Ĥ_1
       Ĥ_0 = ∑_i=1^Nχ_i/2Ẑ_in̂_i
        + ∑_i=1^N(ε_i(t)â_i^† + Ω_i(t)σ̂_i^+ + h.c.)
       + ∑_⟨i,j|⟩(g_i,j(t) â_i^†â_j + h.c.)
       Ĥ_1 = ∑_i=1^N  K_i n̂_i^2 + ∑_⟨i,j|⟩χ_i,jn̂_i n̂_j

Here, â_i (â_i^†) is the annihilation (creation) operator for the ith cavity mode, n̂_i = â_i^†â_i is the corresponding number operator, and σ̂_i^- (σ̂_i^+) is the lowering (raising) operator for the transmon coupled to the ith mode. Ĥ_0 includes the dispersive shift (χ_i), cavity drives (ε_i(t)), transmon Rabi drives (Ω_i(t)), and tunable beamsplitter couplings (g_i,j(t)) implemented via SNAILs. The time-dependence in all parameters can be used to engineer different couplings. Ĥ_1 includes additional unwanted static couplings such as self-Kerr (K_i) and cross-Kerr (χ_i,j). These nonlinear couplings incur contributions from the hybridization of each cavity mode with the transmon and SNAIL couplers. With careful engineering, it is possible to design the components such that the Kerr contributions from transmon and SNAIL cancel or nearly cancel, and χ_i,i+1≈ 0, K_i≈0 could be achieved at a single flux point, with small residual couplings not likely to contribute to the principal results of this work <cit.>. When the beamsplitter interaction is not activated (g_i,j = 0), each oscillator-qubit pair can be independently controlled. For a full discussion of the gates available in the proposed superconducting hybrid oscillator-qubit platform, see Ref. <cit.>. The native gates arising directly from the Hamiltonian of the system are collected in Tab. <ref>. Note that in some parameter regimes, certain transmon-conditioned bosonic operations, such as the conditional displacement (sometimes referred to as `controlled displacement') <cit.>,  can be implemented natively.

All phase-space rotations R_i(θ) can be realized `in software' by absorbing them into the phases of the subsequent microwave pulses. This is because all operations are carried out in the rotating frame of the free evolution of the cavity oscillator mode ωâ^†â (each transmon is only ever connected to a single cavity). Therefore, a phase-space rotation can be applied by changing the phase of the subsequent microwave drive carrying out a gate (such as a beamsplitter or displacement) on the corresponding mode.  Arbitrary qubit Z rotation gates can similarly be carried out in software.

Arbitrary initial states for each cavity can be prepared using optimal control <cit.>: both qubit and cavity are initialised to the ground state, the qubit prepares the cavity in the required Fock state and is then measured to validate this preparation <cit.>. In particular, Fock states |n⟩ with low photon number can be prepared with typical state-of-the-art Fock preparation infidelities of 10^-2 <cit.>.

As mentioned above, it is possible to effectively turn off the `always on' dispersive interaction using dynamical decoupling. Specifically, the transmons can be driven with Rabi rate Ω_i ≫χ_i ⟨â_i^†â_i|$⟩ during a timeT = 2πk/Ω_i, wherekis an integer, to decouple the cavity mode from the transmon and leave the transmons in their initial state. With this, the beamsplitter gateBS_i,j(φ, θ)is performed by simultaneously driving transmonsiandj, and driving the SNAIL coupler linking modesiandjsuch thatg_i,j > 0. Alternative schemes based on transmon echoing can also be used to decouple and perform transmon-state-independent beamsplitting <cit.>.

In addition to unitary control, the transmons can be leveraged for measurements of the cavity modes. In particular, it is possible to measure the boson occupation in the mode in a total duration that is logarithmic in the maximum Fock state. This is done by reading out the photon number bit-by-bit in its binary representation. Each bit is mapped onto the transmon qubit dispersively coupled to the cavity, which is measured and then reset  <cit.>.



 §.§ Trotter simulations


In this work, we focus on a Trotter decomposition of dynamics as our basic simulation algorithm. In such simulations, the time evolution operatorexp(-iĤ T)until timeTwith HamiltonianĤ=∑_γ=1^ΓĤ_γis divided intortime steps of durationt=T/r, i.e.,exp(-iĤ T)=(exp(-iĤ t) )^r. Iftis chosen sufficiently small, a product formula approximation toexp(-iĤt)can then be employed to decompose the Hamiltonian into unitaries  which only act on a few qubits and oscillator modes. For example, a(2p)^th-order Trotter formula in general can be written as a sequence ofN_expexponentials of Hamiltonian termsĤ_γ_idrawn from{Ĥ_1,…,Ĥ_Γ}and a sequence of timest_iwith|t_i|≤tthat satisfies

    exp(-iĤ t)=e^-iĤ_γ_1 t_1⋯ e^-iĤ_γ_N_exp t_N_exp +𝒪(t^2p+1).

There is great freedom in the choice of this decomposition; however, certain orderings can be cheaper, more accurate or more natural.
We choose our decomposition to be given by the local terms inside the sums over lattice sites in many-body Hamiltonians, c.f. Eq. (<ref>) and Eq. (<ref>). For example, the first-order Trotter decomposition of the  forU=0in one spatial dimension, whereB=0, is given by

    e^-i Ĥ_ℤ_2t=∏_⟨ i,j⟩ e^igt  X̂_i,j∏_⟨ i,j⟩ e^iJt ( m̂^†_iẐ_i,jm̂_j + h.c.)+ 𝒪(t^2).

One of the main goals of this paper is to show how to compile unitaries acting on two or more sites, such as the gauge-invariant hoppingexp [ iJt ( m̂^†_i Ẑ_i,j m̂_j + h.c. ) ], into the native operations available in circuit QED hardware or other hardware with equivalent instruction set architectures.

In section <ref>, we discuss the errors incurred from the Trotter approximation and their interplay with compilation errors. 



§ IMPLEMENTATION PRIMITIVES


In this section, we introduce the key compilation strategies we use to synthesize oscillator-oscillator, qubit-qubit, and oscillator-qubit gates from native operations. Throughout, we use the word “qubit” to refer either to a transmon qubit or a dual-rail qubit encoded in two cavities
<cit.>, and only distinguish between the two possibilities when clarification is necessary. A summary of the gates we obtain through these compilation strategies is provided in Tab. <ref>. These primitives will then be used in Sections <ref> and <ref> to compile Hamiltonian terms comprising fermions, bosons, and gauge-fields to native operations. 





Specifically, in Sec. <ref>, we discuss the bosonic SWAP gate, a particularly useful primitive that enables the synthesis of entangling gates between non-neighbouring modes. In Sec. <ref>, we demonstrate how one can synthesize qubit-conditioned bosonic operations using conditional parity gates and, in Sections <ref> and <ref>, develop strategies to realize bosonic operations conditioned or controlled on the Fock-space occupancy of another mode. Building upon these “Fock-projector-conditioned” and “Fock-projector-controlled” operations, we then discuss how one can apply them in sequence to implement more complex density-controlled gates in Sec. <ref>. We then discuss a technique for implementing oscillator-mediated, multi-qubit gates in Sec. <ref>, and follow this with a discussion on approximate methods for synthesizing multi-mode gates in Sec. <ref>. Finally, we discuss the dual-rail encoding in Sec. <ref> as an alternative scheme for realizing qubits in the proposed hybrid oscillator-qubit architecture.




 §.§ Bosonic SWAP

Our architecture in Fig. <ref> is not all-to-all connecting, and hence it is useful to be able to implement SWAP operations between the bosonic sites. In particular, bosonic SWAP gates enable quantum communication and entangling operations between remote bosonic modes, and additionally allow for entangling operations between qubits using the bosonic modes as a quantum bus <cit.>. 

To that end, a bosonic SWAP gate can be realized using a beamsplitter and a pair of phase-space rotations,

    SWAP_i,j = R_i(-π/2)R_j(-π/2)BS_i,j(0,π/2).

withBS_i,j(φ,θ)andR_j(θ)defined in Tab. <ref>. Here, the role of the phase-space rotations is to cancel the spurious phase obtained from the action of the beamsplitter,BS_i,j(0,π/2)|Ψ_i, Ψ_j⟩=e^-i π/2[n̂_i+n̂_j]|Ψ_j, Ψ_i⟩. As mentioned previously, all such phase-space rotationsR_i(θ)can be realized `in software' by absorbing them into the phases of the subsequent microwave drives. In the proposed architecture with SNAILs linking adjacent cavities, SWAP operations have been demonstrated with a duration of around 100 ns and a fidelity of∼0.999<cit.>.  

As discussed in Sec. <ref>, our proposed architecture assumes native entangling gates between each mode-qubit pair and beamsplitters are only available between adjacent sites. However, as discussed in Ref. <cit.>, it is possible to use a beamsplitter network (i.e., a sequence of bosonic SWAPs) to synthesize non-native oscillator-oscillator and oscillator-qubit gates, e.g.,BS_i,j(φ,θ)for any modesi≠j, andCR_i,j(θ)for any transmoniand modej. Thus, we will henceforth assume access to such non-native gates, abstracting away the underlying bosonic SWAPs.



 §.§ Qubit-conditioned bosonic operations


Certain bosonic operations, such as beamsplitters or displacements, can be conditioned on qubits by conjugating them with conditional parity gates,CΠ_i,j<cit.>. The intuition for this is provided by the relation

    CΠ_i,jâ_jCΠ^†_i,j = iẐ_i â_j,

where we used the Baker-Campbell-Hausdorff formula

    e^ÂB̂ e^-Â=B̂+[Â, B̂]+1/2 ![Â,[Â, B̂]]+…,

andCΠ_i,jis the conditional parity gate defined in Tab. <ref>, experimentally realized by evolving under the dispersive interaction for a timeT = π/χ(see Eq. (<ref>)). 

It is possible to implement the Hermitian adjointCΠ_i,j^†by simply appending a bosonic phase-space rotation:

    CΠ_i,j^†   = R_j(-π) CΠ_i,j = CΠ_i,jR_j(-π).

This is becauseCΠ_i,j^†= e^iπ/2Ẑ_i n̂_j = e^iπẐ_i n̂_je^-iπ/2Ẑ_i n̂_j, ande^iπẐ_i n̂_j = cos(πn̂)1 + i sin(πn̂)Ẑ= cos(πn̂) = e^iπn̂_j. In effect, applying aCΠ_i,j^†gate requires the application ofCΠ_i,jand a phase-space rotationR_j(-π)=R_j(π)which, as explained in Sec. <ref>, can be implemented `in software' via modification of the phases of the subsequent microwave drives.

The relation in Eq. (<ref>) allows us to “convert” the generator of certain bosonic operations into one that depends upon the state of a qubit. For example, the conditional displacement gate can be synthesized from an unconditional displacement and a pair of conditional parity gates via

    CD_i,j(α)    =e^Ẑ_i(αâ_j^† - α^*â_j)
       = e^-iπ/2Ẑ_i n̂_je^i(αâ_j^† + α^*â_j)e^iπ/2Ẑ_i n̂_j
       =CΠ_i,jD_i(iα)CΠ_i,j^†,

whereα∈C. Note the modified phase in the unconditional displacement. This is due to the factor ofiinherited from the transformation in Eq. (<ref>). Furthermore, we emphasize that while this sequence produces exactly the conditional displacement gate as defined in Table <ref>, this gate can more efficiently be realized natively using the techniques developed in Ref. <cit.>.




As another example, we can exactly compile the conditional beamsplitter gate as follows:

    CBS_i,j,k(φ,θ)    = e^-iθẐ_i(e^iφâ^†_j â_k + e^-iφâ_j â^†_k), 
       = e^-iπ/2Ẑ_i n̂_je^θ(e^iφâ^†_j â_k - e^-iφâ_j â^†_k)e^iπ/2Ẑ_i n̂_j
       = CΠ_i,j BS_j,k(φ+π/2,θ)CΠ_i,j^†,

where{φ,θ} ∈R. This gate will play an instrumental role in Sec. <ref> for the compilation of interaction terms between matter and gauge fields. Its duration is experimentally limited by the speed of each conditional parity gate, requiring roughly1μs each. As such, the conditional beamsplitter gate requires∼2.1μs, including the worst-case∼100ns required for a full SWAP beamsplitter <cit.>.

More generally, as shown in Fig. <ref>a, we can employ this strategy to realize a conditional unitaryCU^Z_i,jfrom its unconditional counterpart:

    CU^Z_i,j   = e^Ẑ_i(Θ̂â^†_j - Θ̂^†â_j)
       = CΠ_i,je^i(Θ̂â_j^† + Θ̂^†â_j)CΠ_i,j^†.

This decomposition holds for any choice of operatorΘ̂that satisfies[Θ̂,â_j] = [Θ̂,Ẑ_i]=0. For instance,Θ̂→α, corresponds to the conditional displacement in Eq. (<ref>) whileΘ̂→θe^iφ â_kyields the condition beamsplitter in Eq. (<ref>).

This strategy can be iterated to realize bosonic gates conditioned on multiple qubits. For example, one can realize a doubly-conditioned gate of the form

    CU^ZZ_i,j,k=e^Ẑ_i Ẑ_j (Θ̂â_k^† - Θ̂^†â_k )

for any operatorΘ̂that commutes withẐ_i,Ẑ_j, andâ_k. In this case, conditional parity gates are implemented on modekfor qubitsjandkin sequence with the aid of SWAP gates. Alternatively, one can apply conditional parity gates acting on different qubit-mode pairs in parallel as shown in Fig. <ref>b.

We note that, in the above, we use `qubit-conditional' and `qubit-controlled' interchangeably to refer to a gate that enacts a bosonic operation that depends upon theZbasis state of a qubit. In particular, distinct non-trivial operations are enacted for both qubit state|0⟩and|1⟩. However, in certain contexts, it is useful to synthesize a `1-controlled' bosonic operation, i.e., one that acts only if the qubit is in|1⟩. In this work, we will refer to this as ann-controlled gate, referring to the operator that appears as a multiplicative factor in the generator. For example, for ann-controlled mode rotation, we have

    e^-iθn̂_qubitn̂_mode = R_i(-θ/2)CR_i,i(θ),

wheren̂_qubit=(1-Ẑ)/2. More generally, we have
    
    e^iθn̂_qubitÔ_i = e^iθ/2Ô_ie^-iθ/2Ẑ_i Ô_i,

for anyÔ_iwhich commutes with operations on the qubit.





 §.§ Fock-projector-conditioned bosonic operations


An important primitive in the hybrid oscillator-qubit architecture of Fig. <ref> is the synthesis of nontrivial two-mode entangling gates. In this section, we develop a technique to realize bosonic operators that are conditioned on the Fock-space information of a second mode:

    CU_i,j^P̅ = e^-iẐ_ancP̅_iÔ_j,

whereP̅_iis an operator with eigenvalues±1defined below (hence the name `conditioned'), `anc' refers to an ancillary transmon qubit andÔ_jis an arbitrary (Hermitian) operation acting on modej. For an arbitrary subspace𝒫of the oscillator,P̅_igives states inside𝒫a phase of-1, and states outside𝒫a phase of1. Mathematically, 

    P̅_i   =∑_n∉𝒫|n⟩⟨n|_i-∑_n∈𝒫|n⟩⟨n|_i,

where each projector acts on modei. Note that we define the operatorP̅_iwithout a hat for simplicity of notation.

BecauseP̅_ihas eigenvalues±1, the core idea is to map this information to an ancillary qubit, which is then used to mediate this information to a second mode via a hybrid qubit-mode gate. This effectively realises the Fock-projector-conditioned bosonic gate in Eq. (<ref>). This requires the ancillary qubit to be initialized to an eigenstate ofẐ_anc. This ancillary qubit can correspond to the transmon qubit dispersively coupled to modeiorj, or any other qubit via appropriate use of bosonic SWAP gates. 

To explain further, we begin with the particular example of a parity-conditioned gate, defined as

    CU^Π̅_i,j   =e^-iẐ_ancΠ̅_iÔ_j,

where

    Π̅_i    = e^iπn̂_i

is the parity operator acting on modeiwith eigenvalues±1. Therefore, for an ancilla initialized to|0⟩_anc,CU^Π̅_i,japplies the unitarye^-iÔ_j(e^iÔ_j) to modejif modeihas even (odd) parity. To realize this gate, it is helpful to first note that 

    SQR_anc,i(π⃗_0,0⃗) = e^-iπ/2X̂_anc∑_n  odd|n⟩⟨n|_i,

whereπ⃗_0 = {0,π,0,π,…}and0⃗={0,0,0,…}, andΠ̂_i(as opposed toΠ̅_i) is the projector over odd Fock states. In other words, an SQR gate is simply a Fock-projector-controlled qubit rotation. By leveraging the Baker-Cambell-Hausdorff formula in Eq. (<ref>), it can be shown that

    SQR_anc,i(π⃗_0,0⃗)Ẑ_ancSQR^†_anc,i(π⃗_0,0⃗) = Ẑ_ancΠ̅_i.


Thus, in direct analogy to the strategy for synthesizing qubit-conditional bosonic operations in Sec. <ref>, we can leverage a pair of SQR gates to condition a bosonic operation on the parity of a second mode:

    CU^Π̅_i,j = SQR_anc,i(π⃗_0,0⃗)e^-iẐ_ancÔ_jSQR_anc,i(-π⃗_0,0⃗),

where we have used the fact thatSQR^†_anc,i(θ⃗,φ⃗)=SQR_anc,i(-θ⃗,φ⃗). This requires the realization of the intermediary qubit-oscillator gateexp(-iẐ_ancÔ_j), either natively or, e.g., using the technique in Sec. <ref>.

Returning to the general form in Eq. (<ref>), it is possible to generalize this technique by substitutingΠ̅_iwithP̅_iwhich allows for any arbitrary set of Fock states𝒫by choosing appropriate angles for the SQR gate in Eq. (<ref>). In particular, Eq. (<ref>) generalizes as

    SQR_anc,i(θ⃗_P,0⃗)Ẑ_ancSQR^†_anc,i(θ⃗_P,0⃗) = Ẑ_ancP̅_i,

whereθ⃗_P = {θ_0,θ_1,…, θ_n …}and

    θ_n =
        π   if |n⟩∈𝒫
    
          0    otherwise.

Akin to Eq. (<ref>), this leverages the fact thatSQR_anc,i(θ⃗_P,0⃗)realizes a projector-controlled qubit rotation:

    SQR_anc,i(θ⃗_P,0⃗) = e^-iπ/2X̂_anc∑_n∈𝒫|n⟩⟨n|_i.
 Thus, with this, one can realize the generalized Fock-projector-conditional gate,

    CU_i,j^P̅   = e^-iẐ_ancP̅_iÔ_j
       =SQR_anc,i(θ⃗_P,0⃗)e^-iẐ_ancÔ_jSQR_anc,i(-θ⃗_P,0⃗).

This sequence is illustrated in Fig. <ref>a.




 §.§ Fock-projector-controlled bosonic operations


In this section, we build on the technique developed in the previous section in order to realize bosonic operators that are controlled on the Fock-space information of another mode,

    CU_i,j^P̂   = e^-iẐ_ancP̂_iÔ_j,

whereP̂_iis a Fock-space projector onto the subspace𝒫, 

    P̂_i = ∑_n∈𝒫|n⟩⟨n|_i
        
    
        = (1 - P̅_i)/2.

The second equality reminds us that the projector and its reflection operator about the orthogonal complementP̅_iare reminiscent ofn̂_qubitandẐfor qubits, wheren̂_qubit=(1-Ẑ)/2. It follows that whileCU^P̅_i,japplies a distinct, nontrivial operation for each value of the modei∈𝒫, this gateCU^P̂_i,jacts withÔ_jonly if modei∈𝒫(hence the `controlled' nomenclature).

In fact, using (<ref>), we can easily construct the target gate (<ref>), similar to the construction in (<ref>):

    CU_i,j^P̂   = e^-iẐ_anc((1_i - P̅_i)/2)Ô_j
       = e^-iẐ_ancÔ_j/2e^iẐ_ancP̅_iÔ_j/2.


 §.§ Boson density-controlled operations


A powerful use-case of the generalized controlled gate in Eq. (<ref>) is to apply many of them in sequence to synthesize complex entangling gates between two modes, mediated by an ancillary qubit in a known eigenstate ofẐ_anc. In this section, we leverage this technique to exactly compile a bosonic density-controlled operation:

    CU^n̂_i,j=e^-iẐ_ancn̂_i Ô_j


To obtain this gate, we follow the iterative procedure illustrated in Fig. <ref>c. The core idea is to sequentially condition the oscillator-qubit gateU=exp(-iẐ_ancÔ_j)on each bit of the binary representation ofn̂_i. With each iteration, an appropriate rotation angle is chosen such that the sequence altogether produces the desired operationCU_i,j^n̂. This strategy is reminiscent of binary photon number readout <cit.>, and requires a circuit depth that is logarithmic in the boson number cutoffn_max.



Specifically, as a particular instance of Eq. (<ref>), we define the operator generated byÔ_jand conditioned on thekth superparity of modei,

    CU_i,j^Π̅_k(θ)    = e^-iθẐ_ancΠ̅_kiÔ_j
       = SQR_anc,i(π⃗_k,0⃗)e^-iθẐ_ancÔ_jSQR_anc,i(-π⃗_k,0⃗),

where we have made the variable angleθexplicit for reasons that will become clear shortly.
Similarly, we defineCU_i,j^Π̂_k(θ)in analogy to Eq. (<ref>), see Fig. <ref>b. Here,Π̂_kis projection operator that returns thekth bit ofn̂_i(andΠ̅_kits corresponding reflection operator, returning±1), andπ⃗_k = {θ_0, θ_1, θ_2, …}with

    θ_n =
        π   if ⌊ n/2^k⌋ mod 2 = 1
    
          0    otherwise.

For exampleπ⃗_0 = {0,π,0,π,…}is equivalent to the parity vector defined below Eq. (<ref>), converting the SQR gate into a qubit rotation controlled on the least-significant (zeroth) bit ofn̂_i. We note thatΠ̅_0≡Π̅, where the latter is defined in Eq. (<ref>).  Likewise, the choiceπ⃗_1 = {0,0,π,π,…}corresponds to qubit rotation controlled on the first bit ofn̂_i. In general,

    SQR_anc,i(π⃗_k,0⃗) = e^-iπ/2X̂_ancΠ̂_ki,

and the choiceπ⃗_kenables a qubit rotation controlled on thekth bit ofn̂_i.

Applying successive operations controlled on each bit ofn̂_iin sequence then yields the general form

    CU_i,j^f(Π̂_0, Π̂_1, Π̂_2,…)   = ∏_k=0^K-1CU_i,j^Π̂_k(θ_k) 
       = exp(-iZ_anc∑_k=0^K-1θ_kΠ̂_kiÔ_j),

requiringK = ⌈log_2(n_max+1)⌉iterations of (super)-parity controlled gates, and wheref(Π̂_0,Π̂_1,Π̂_2,…)is any linear function in the bit operatorsΠ̂_k. Returning to the initial goal of synthesizingCU_i,j^n̂in Eq. (<ref>), we note that

    n̂_i = ∑_k=0^K-1 2^kΠ̂_ki.

Consequently, choosingθ_k = 2^krealizes the desired gate. This sequence is shown in Fig. <ref>c.

As a particular example of this gate that will prove useful in Sec. <ref> for implementing density-denstiy interactions, one can chooseÔ_j→θn̂_jto realize the phase-space rotation of one mode controlled on the density of second,

    RR_i,j(θ) =  e^-iθẐ_ancn̂_in̂_j.


Finally, we emphasize that while our focus has been on synthesizing the density-controlled operatorCU_i,j^n̂, the general formCU_i,j^f(Π̂_0, Π̂_1, Π̂_2,…)enables further possibilities beyond this particular choice. Specifically, it allows for operations conditioned on any functionf(Π̂_0, Π̂_1, Π̂_2,…)that is linear in the bit operatorsΠ̂_ki. Morevover, we note that an arbitrary functionf(n̂_i)can be realized by iterating the primitiveCU^P̂_ija number of times linear inn_max. This possibility is discussed in our companion work, Ref. <cit.>.



 §.§ Oscillator-mediated multi-qubit gates


As the architecture illustrated in Fig. <ref> does not include direct couplings between transmons, their utility for encoding fermionic and gauge field degrees of freedom relies on the ability to synthesize transmon-transmon entangling gates using native gates. Here, we adopt an oscillator-mediated approach  proposed in our companion work Ref. <cit.> to realize such entangling gates that bear close similarity to Mølmer-Sørensen gates in trapped-ion platforms <cit.>, and also some parallels to mediated gates in silicon <cit.>. Crucially, this approach allows for the exact analytical compilation of multi-qubit gates, independent of (and without altering) the initial state of the mediating oscillators in the absence of noise. The core idea is to use a sequence of conditional oscillator displacements that form a closed phase-space path that depends upon the state of the qubits (see, for example, Fig. <ref>a). Upon completion, this closed path leaves the oscillator unchanged yet imparts a geometric phase on the system that depends upon the state of the qubits, thus enacting the target gate. 

To realize these oscillator-mediated multi-qubit gates, we require the ability to enact displacements of thejth oscillator conditioned on theith transmon,

    CD_i,j(α) = e^Ẑ_i(αâ_j^† - α^* â_j).

Fori=j, such conditional displacements can be realized either natively <cit.> or compiled using conditional-parity gates and unconditional displacements, as shown in Eq. (<ref>); furthermore, as discussed in Ref. <cit.> and summarized in Sec. <ref>, it is possible to synthesize conditional displacements for arbitraryiandjby additionally leveraging a beamsplitter network that links the two sites. 

As a particular example, the two-qubit entangling gateZZ_i,j(θ) = exp(-iθ/2Ẑ_i Ẑ_j)can be realized via a sequence of four conditional displacements,

    ZZ_i,j(θ)    = CD_i,j(-iα)CD_i,i(α)
       ×CD_i,j(+iα)CD_i,i(-α).

The corresponding phase-space trajectories and circuit diagram are shown in Fig. <ref>a and Fig. <ref>b, respectively. As explained in Ref. <cit.>, this sequence can be further optimized by leveraging conditional displacements on both modesiandjsuch that each accumulates a geometric phase in parallel, leading to a reduction in gate duration by a factor of√(α)compared to the above single-mode approach, while also reducing the potential impact of non-idealities such as self-Kerr interactions (anharmonicities). See Ref. <cit.> for more details.

This strategy can be extended to realize oscillator-mediatedN-qubit gates that are useful for studying lattice gauge theories. For example, as will be discussed in Sec. <ref>, the four qubit gate

    e^-i θ/2Ẑ_i Ẑ_j Ẑ_k Ẑ_l

is useful for implementing a magnetic field term in theZ_2LGT introduced in Eq. (<ref>). To realize this gate, we modify the sequence in Eq. (<ref>) such that each displacement is conditioned on the joint state of pairs of qubits, e.g.,Ẑ_i Ẑ_jorẐ_k Ẑ_l, as shown in Fig. <ref>c. To that end, we define the doubly-conditional displacement as

    CD^ZZ_i,j,r(α) = e^Ẑ_iẐ_j(αâ^†_r - α^* â_r).

This gate be implemented via the technique described in Sec. <ref>. As illustrated in Fig. <ref>d, the four-qubit term can be exactly realized by sequencing four such doubly-conditional displacements:

    e^-i θ/2Ẑ_i Ẑ_j Ẑ_k Ẑ_l   = CCD_i,j,r(-iα)CCD_k,l,r(α)
       ×CCD_i,j,r(+iα)CCD_k,l,r(-α),

withα= √(θ)/2.

Similar to the implementation ofZZ_i,j(θ), it is beneficial to use multiple modes in parallel to accumulate the necessary geometric phase and reduce the displacement parameterα. Moreover, this strategy can be extended to realize arbitrary operations of the forme^-iθŴ, where hereŴis an arbitrary-weight Pauli operator. See Ref. <cit.> for a more complete discussion.

Finally, we note that a similar strategy for oscillator-mediatedN-body entangling gates has been proposed <cit.> and demonstrated <cit.> in trapped-ion platforms. There, spin-dependent squeezing operations are used in place of conditional displacement operators, but the overall framework bears similarity to the strategy described here for the superconducting architecture described in Fig. <ref>.



 §.§ Approximate synthesis of multi-mode gates



To approximately implement arbitrary multi-mode operations that fall outside of the form in Sec. <ref> and Sec. <ref>, we use the method introduced in Ref. <cit.>. It relies on a variant of the Baker-Campbell-Hausdorff (BCH) formula, 

    e^Â t e^B̂ te^ -Â t e^-B̂ t=e^[Â, B̂] t^2+O(t^3).

The key insight is to choose the operatorsÂ, B̂such that they enact commuting bosonic operationsÔ_1,Ô_2(possibly acting on distinct sets of different modes), but conditioned on the same qubit  with respect to non-commuting Pauli operators, e.g., using the method described in Sec. <ref>. For example, we can leverage this idea as follows:

    Γ_1,2   = e^i X̂Ô_1 t e^i ŶÔ_2 t e^-i X̂Ô_1 t e^-iŶÔ_2 t
       =e^[X̂Ô_1, ŶÔ_2](-i t)^2+O(t^3),
       =e^-2iẐÔ_1Ô_2t^2+O(t^3) .

See also Fig. <ref>a for the circuit diagram. In this way, one can implement operator multiplication betweenO_1andO_2. Combining this with the Trotter formula to implement addition, following Ref. <cit.> and Fig. <ref>b, we can implement exponentials of non-linear functions in the creation and annihilation operators.

When the commutator obtained through the above BCH formula contains both desired and undesired high-order mode terms, we can cancel the latter by applying Eq. (<ref>) a second time with a different choice of operators. Via the Trotter formula we then obtain: 

    Δ_1,2   = e^-2iẐÔ_1Ô_2t^2e^-2iẐÔ'̂_̂1̂Ô'_2t^2+O(t^3)
       = e^-2iẐ(Ô_1Ô_2+Ô'_1Ô'_2)t^2+O(t^3),

where the error term will be different between the two lines and we rewrote the error as an additive error by Taylor expanding the exponential. Here, each term is conditioned on the same qubit, either natively or via SWAP operations between the qubits or the modes.

In Fig. <ref>c, we use this technique to approximately synthesize the four mode operator,

    e^-iẐ(â^†b̂ĉ^†d̂ + h.c.)(2θ)^2.

which will later be used in Sec. <ref> to realize the plaquette term in theU(1)quantum link model.



 §.§ Dual-rail qubits


While so far we have discussed primitives for mapping the model to the hardware with transmons being used as data or ancillary qubits and the modes being used to encode bosonic degrees of freedom, another possibility is to use the modes to encode qubits. In this section, we discuss such an encoding: the dual-rail qubit.
The core idea is to encode a single qubit using a pair of cavity modes that share a single photon. The location of the photon then represents the qubit state <cit.>
    |0⟩_DR   =|0,1⟩
    
        |1⟩_DR   =|1,0⟩.

Denoting the cavity photon numbers to be|n^a,n^b⟩, we see that the PauliZoperator has various equivalent representations,

    Ẑ^DR   =n̂^b - n̂^a
       =1-2n̂^a=2n̂^b-1
       =e^iπn̂^a=-e^iπn̂^b,

where we have used the fact thatn̂^a + n̂^b = 1. The remaining Pauli operators are given by,

    X    = a^† b + a b^†
    
        Y    = i[a^† b - a b^†].

Note that this is simply the Schwinger boson representation of a spin one-half.

The primary advantage of the dual-rail encoding is that photon loss in either the left or right cavity results in an erasure error that is detectable via joint-parity measurements <cit.>. Furthermore, single-qubit operations are straightforward: rotations about any axis in the azimuthal plane can be carried out with a beamsplitter, and rotations about theZaxis correspond to a phase-space rotation of one of the modes (up to a global phase). Furthermore, it has been shown that one can construct a universal gate set that enables detection of ancillary transmon dephasing and relaxation errors, the latter made possible by leveraging three levels of the transmon qubit <cit.>. Consequently, for models that require many qubit-qubit gates (such 2D fermionic problems that necessitate fermionic SWAP networks), the dual-rail encoding is potentially advantageous for mitigating errors with post-selection in near-term experiments.

ZZ_i,j^DR(θ). As an example of an entangling gate between dual-rail qubits, we consider the compilation ofZZ^DR_i,j(θ)=e^-i(θ/2)Ẑ^DR_i Ẑ^DR_j. Using the fact thatẐ^DR_i=1-2n̂_i^a, this gate can be expressed as

    ZZ_i,j^DR(θ)   =e^-iθ/2e^-i2θn̂_i^an̂_j^ae^iθn̂_i^ae^iθn̂_j^a
    
            =   e^-iθ/2(⟨0|_ anc⊗ )RR_i,j(2θ)(|0⟩_ anc⊗ )
       ×R_i(-θ)R_j(-θ),

where, for simplicity of notation, we have dropped the superscripta. Therefore, this gate can be realized by combining a pair of phase-space rotations with the technique introduced in Sec. <ref> to realizeRR(2θ)– see Eq. (<ref>). 




However, note that in the dual rail encoding, the boson number of modeiis restricted to either 0 or 1, which means that the SQR gate used to implement Eq. (<ref>) of Sec. <ref> can be reduced to a conditional phase-space rotation up to single qubit operations. For the choice ofθ⃗ = π⃗_0, this conditional rotation in the{ |0⟩,|1⟩}subspace corresponds to a conditional-parity gate:

    SQR_i,j(π⃗_0, 0⃗) ⇒H_i CΠ_i,jH_i =CΠ^X_i,j

whereH_iis a Hadamard gate on qubiti, and we have defined the shorthandCΠ^X_i,jto denote theX-conditional parity gate.
Therefore, we can re-write the decomposition ofZZ^DR_i,j(θ)as:

    ZZ_i,j^DR(θ)=   e^-iθ/2(⟨0|_ anc⊗ )CΠ^X_anc,iCR_anc,j(-2θ)CΠ^X_anc,i
       ×(|0⟩_ anc⊗ )R_i(π)R_j(θ)R_i(-θ)R_j(-θ),

where the top line, up to the global phase, corresponds toRR(2θ)in the truncated dual rail subspace. The complete circuit (with single qubit rotations simplified) is presented in Fig. <ref>a. We note that sequential ancilla-conditional gates using the same transmon (but distinct modes) requires implicit SWAP operations not shown. 

While exact, the circuit in Fig. <ref>a does not enable the detection of transmon dephasing errors. To that end, a scheme to synthesize error-detectable entangling gates was presented in Ref. <cit.>. It leverages an “exponentiation circuit” that interleaves ancilla-controlled unitaries (controlled-P̂) with ancilla rotations to construct any unitaryP(θ)=exp(-i θ/2 P̂), providedP^2=1̂. Returning to the case ofZZ_i,j^DR(θ), this gate can be constructed following this prescription in the previous paragraph, choosingP=Ẑ ⊗Ẑand implementing controlled-P̂using a pair of conditional-parity gates. We illustrate the full implementation in terms of our discrete gate set in Fig. <ref>b. Notably, this latter circuit can be understood as a further decomposition of its counterpart in Fig. <ref>a, where the middleCR(-2θ)gate is broken into conditional-parity gates and single qubit operations using the identities in Sec. <ref>. 

While both circuits are mathematically equivalent, the structure of Fig. <ref>b enables partial error detection of some of the dominant errors, including transmon dephasing. However, it comes at the cost of requiring four conditional-parity gates, each with a durationT_CΠ = π/χ. Thus, ignoring the cost of SWAPs and single-qubit and mode operations, the minimum duration ofZZ_i,j^DR(θ)using the method in Fig. <ref>b isT_b,min = 4π/χ, independent ofθ. In contrast, the method in Fig. <ref>a requiresT_a,min = 2(π+  θ)/χ, a reduction for even the maximal entangling caseθ=π/2. Thus, there is a tradeoff in error detectability and gate duration, and the approach in Fig. <ref>a may therefore be beneficial in certain contexts, particularly for repeated applications ofZZ_i,j^DR(θ)for smallθ(e.g., for Trotterized circuits).



§ IMPLEMENTATION OF MATTER FIELDS


In this section, we leverage the compilation strategies of the previous section to realize the unitary time evolution operatore^-iĤ tfor a range of common Hamiltonian terms involving bosonic and fermionic matter. For the purely bosonic matter terms, in Sec. <ref>, we show that some of the common dynamics terms correspond to native gates, and others to terms that can be realized using the compilation techniques in the previous section. For the purely fermionic matter terms, in Sec. <ref>, we demonstrate how to encode fermionic SWAP networks into various qubit encodings possible in hybrid oscillator-qubit hardware. For the fermion-boson interaction terms, in Sec. <ref>, we present novel compilation strategies specific to hybrid oscillator-qubit hardware. We summarize these results in Tab. <ref> and Tab. <ref>. All terms involving gauge fields will be treated in Sec. <ref>.



 §.§ Bosonic matter


In this section, we discuss the most common Hamiltonian terms for bosonic matter fields (which typically involve small powers of the creation and annihilation operator). The results are summarized in Table <ref>. We first discuss hopping as this is the most elementary term coupling bosonic modes. We then discuss terms related to powers and products of number operatorsn̂: the onsite potential (linear inn̂), the onsite interaction (quadratic inn̂), and the intersite interaction (bilinear,n̂_in̂_j). We map the bosonic matter fields to the native bosonic modes supported in the hybrid oscillator-qubit hardware in Fig. <ref>.



  §.§.§ Hopping

The time evolution operatore^-iĤ tfor the hopping interactionĤ = ∑_⟨i,j|⟩J_i,j ( â^†_i â_j + h.c. )between bosonic matter sites can be implemented via Trotterization using a sequence of beamsplitter gatesBS_i,j(φ,θ)(defined in Tab. <ref>) by choosingθ= Jt.
Furthermore, a static gauge field (vector potential) appears as a complex phase in the hoppingJ→J e^iφ<cit.>,  and is important to simulate models of the fractional quantum Hall effect <cit.>. To implement this term, one can simply adjust the phaseφof the beamsplitter on each link such that a net phase is accumulated in moving around each plaquette.



  §.§.§ Onsite potential


The onsite potential isĤ = ∑_i μ_i n̂_i, withn̂_i = â_i^†â_iμ_i ∈ℝ. The time evolution of this term can be implemented on each site separately using a phase space rotation gate (defined in Tab. <ref>) as is shown in Tab. <ref>.

As mentioned previously in Sec. <ref>, it is possible to absorb phase space rotations into operators that include creation/annihilation operators such as beamsplitters.  In particular, a constant site energy shift can be created with a fixed frequency detuning of the microwave tones that activate the beamsplitters.



  §.§.§ Onsite interaction


The bosonic onsite interaction is ∑_i U_i n̂_i^2withU_i ∈ℝ. To implement the time evolution operator for this term, one can use the Selective Number-dependent Arbitrary Phase (SNAP) gate <cit.>.SNAP_i(θ⃗), as defined in Tab. <ref>, operates in the strong-dispersive coupling regime where the qubit frequencies depend on the mode occupation such that the SNAP gate effectively imparts an independently chosen phase to the system comprising qubit and mode, for each vaue of the mode occupation. Choosingθ_k=U_i k^2t, wherekis photon number corresponding to a particular Fock state implements the onsite interactionn̂^2_ion theith mode. This procedure additionally requires an ancillary qubit initialized to the state|0⟩.



  §.§.§ Intersite interaction


The time-evolution of a nearest-neighbor density-density interactionĤ =  ∑_⟨i,j|⟩ J_i,j n̂_i n̂_jbetween two modes in separate cavities, which is useful for example in simulation of the extended Hubbard model or dipolar interacting systems, can be implemented with the method proposed in Sec. <ref> via theRR(θ)gate defined in Eq. (<ref>):

    e^-i J_i,jn̂_in̂_j t =    (⟨0|_ anc⊗ )RR_i,j(J_i,jt)(|0⟩_ anc⊗ ).

Here, we have used an ancillary qubit that begins and (deterministically) ends in the state|0⟩_anc. As discussed in Sec. <ref> and shown in Table <ref>, this gate requires a native gate count that is logarithmic in the bosonic cutoffn_max. Long-range interactions between sitesiandjcan be implemented at the cost of bosonic SWAP gate count that is linear in the Manhattan-metric distance between the sites.




 §.§ Fermionic matter


In this section, we show that it is possible to treat purely fermionic terms within our proposed hybrid qubit-oscillator architecture. While this is not an area in which to expect an explicit advantage over qubit-only hardware, we show that no significant cost is added compared to existing qubit platforms. 

We treat fermionic matter by mapping the fermions to qubits and rely on FSWAP networks <cit.> of Jordan-Wigner (JW) strings for a low-depth implementation. For an overview of the set of operations required to implemente^-i Ĥtfor a2D fermionic HamiltonianĤin qubit-based hardware,  see App. <ref>. In summary, the two-dimensional lattice is virtually mapped to a one-dimensional chain. Nearest-neighbour hoppings in the2D lattice become either nearest-neighbour hoppings in the1D JW chain (which we call “JW-adjacent”) or long-range hoppings (which we call “non-JW-adjacent”). The former hoppings only require iSWAP operations, whereas the latter also require FSWAP operations. Gates that involve the fermionic number density are mapped to gates involving theẐoperator and do not add any overhead associated with the Jordan-Wigner encoding.

We first summarize the possibilities for encoding fermions (i.e., implementing iSWAP and FSWAP operations) either in transmons (Sec. <ref>) or in cavity dual-rail qubits (Sec. <ref>). While the former is more hardware efficient, the latter enables error-detection and, as a result, is a particularly promising avenue for near-term simulations on noisy hardware.  The high-level results for both encodings are summarized in Table <ref>. In Sec. <ref>, we discuss the fidelity of gates that underly the simulation of fermions in the dual-rail encoding. 



  §.§.§ Transmon qubits
 
In the proposed 3D post-cavity circuit QED setup, encoding the fermion into the transmon enables the use of the cavities for a separate purpose, such as for representing the phonon modes in the Hubbard-Holstein model, or for encoding bosonic gauge fields as will be demonstrated in the context of aℤ_2lattice gauge theory later in this work.

As the architecture illustrated in Fig. <ref> does not include direct couplings between transmons, we rely on the oscillator-mediated compilation strategy discussed in Sec. <ref> to realize entangling gates,ZZ_i,j(θ)in particular. We reemphasize that this approach does not require that the mediating oscillators are in a known, particular state. Thus, they can simultaneously be used to encode matter or gauge degrees of freedom.



  §.§.§ Dual-rail qubits 


We now discuss the separate possibility to map fermions to dual rail qubits. As discussed in Sec. <ref>, dual rail qubits present a number of advantages, including the possibility for detection of photon loss and ancilla relaxation and dephasing errors <cit.>.

FSWAP. For the non-JW-adjacent hopping terms, a fullSWAP_i,joperation is required between the dual-rail qubits as part of the FSWAP operation. This can be implemented by separately applying aSWAP_i,joperation between the cavities representing the|1⟩and|0⟩states, as represented in Fig. <ref>.  As described in App. <ref>, the FSWAP also requires aCZ_i,j(θ)gate which is equivalent to theZZ_i,j(θ)gate up to single qubit phase gates. As discussed in Sec. <ref> (c.f. App. <ref>), these single qubit phase gates correspond to phase-space rotations on the mode which hosts the boson when the dual rail qubit is in the state|1⟩. The corresponding circuit is shown in Fig. <ref>a.

iSWAP(θ). For the JW-adjacent hopping terms, following Eq. (<ref>), we only need to implement a variable-angle iSWAP gate defined in Fig. <ref>b. This requires onlyZZ_i,j(θ)gates andHandSsingle-qubit rotations. An error-detectable implementation ofZZ_i,j(θ)between dual-rail qubits was presented in Ref. <cit.> and is summarized here in Sec. <ref>. We show the circuit for implementinge^­iθ( X̂_i X̂_i+1+Ŷ_i Ŷ_i+1)at the level of the modes supporting the dual-rail qubits in Fig. <ref>b.





  §.§.§ Comparing advantages of different platforms for encoding fermions


We note that different platforms have different advantages for implementing non-JW-adjacent hopping. 

In a planar transmon qubit array with tunable couplers, FSWAP operations can be implemented using the FSIM gate using FSIM(θ=π/2,ϕ=0)with high fidelity (≈0.9930) <cit.>.e^iθ(X̂ X̂ + Ŷ Ŷ )operations can be implemented using FSIM(θ,ϕ=0), also with high fidelity (currently≈0.9962) <cit.>.

In the ion trap QCCD architecture <cit.> FSWAP operations can be implemented using a SWAP operation and aZZ(π/2). Because the ions can be physically moved around the circuit to achieve all-to-all connectivity, the fidelity of the FSWAP operation reduces to the fidelity ofZZ(π/2)which is currently 0.9980 <cit.>, yielding a fidelity advantage compared to the planar transmon implementation. Furthermore,e^iθ(X̂ X̂ + Ŷ Ŷ )can be implemented with twoZZ(θ)which currently has a fidelity of≈0.9960<cit.>.

In the circuit QED setup using dual rail qubits, as discussed, we can use theZZ(θ)gates to implement fermionic SWAP networks similar to ion trap QCCD architectures. Beamsplitter operations in the high-Q cavity setup required for a SWAP operation between dual-rail qubits have very high fidelity of0.9992<cit.> (comparable to single-qubit gates on transmons). Importantly, we can use the error-detection capability to increase the fidelity of the operations with post-selection, as discussed in <cit.>. This results inZZ(π/2)gates with a theoretically calculated fidelity in the presence of noise and post-selection of≈0.9999<cit.>, for a pure dephasing channel with timeT_ϕ=200μs. The error detection comes at the price of a post-selection requirement, leading to a shot overhead. For the above value ofT_ϕ, it was estimated that an error would be detected in2%<cit.> of shots. While this percentage increases exponentially with circuit depth, even for100ZZ_i,i+1(θ)gates, only1-(1-0.02)^100≈87%of shots would have to be discarded due to a detected error. This leads to a shot overhead of a factor of1/(1-87%)≈7.7. 



 §.§ Fermion-boson interactions


Here we compile the unitary time evolution operator for two paradigmatic fermion-boson models: phonon-electron interactions in solids and the coupling of strong light fields with electrons. Employing the techniques of the previous sections, we discuss the use of either transmons or dual rail qubits to encode the fermionic matter, and employ modes to represent the bosonic matter.




  §.§.§ Phonon-matter interactions


The Holstein model <cit.> describes (dispersionless) optical phonons in the solid state coupled to the density of electrons. It consists of a kinetic energy term for the electrons (c.f. first row in Tab. <ref>), an onsite potential term for for the phonons (c.f. Sec. <ref>) as well as an electron-phonon interaction term 
    Ĥ= ∑_i g_i ĉ^†_i ĉ_i (â_i^† + â_i).

To implement the latter, note that in the Jordan-Wigner encoding, independent of the dimension, this Hamiltonian term becomes equivalent to the Spin-Holstein coupling <cit.>,

    Ĥ= ∑_i g_i/2 (Ẑ_i+1) (â_i^† + â_i).

   Time evolution under this Hamiltonian can be easily implemented since all terms in the sum commute and each has a simple implementation; the term proportional toẐ_igenerates a conditional displacement gate and the remaining term can be eliminated by a simple coordinate frame shift for each oscillator (resulting in a constant chemical potential shift for the fermions).
    
Transmon qubits. If the qubits are mapped to transmon qubits, conditional displacements can be compiled using conditional paritiesCΠ_i,j(Tab. <ref>) and a displacementD_i(g_i t/2). This control was discussed in Sec. <ref> and is illustrated in Fig. <ref>a.

Dual-rail qubits. In the case of dual-rail qubits,Ẑ_i^DR = e^iπn̂_i,1, where the subscriptsi,1indicates the mode which hosts the boson when the dual rail qubit is in state|1⟩. This means that the required dual rail qubit-conditional displacement will take the forme^-iθe^iπn̂_i,1(â_j^†+ â_j), whereθ= g_i t/2. This is a displacement conditioned on the parity of another mode{i,1}and can be realized byCU^Π̅_i,j, whose compilation is presented in Sec. <ref>. The full operation illustrated in Fig. <ref> is as follows:

    e^-iθẐ_i^DR (â_j^† + â_j)= 
       ⟨0|_ ancCΠ^X_anc,iCD_anc,j(iθ)CΠ^X_anc,iR_i(π)|0⟩_ anc,

whereCD(α)corresponds to a conditional displacement gate as defined in Table <ref>, andCΠ^Xdenotes the Hadamard-conjugated conditional parity gate introduced in Eq. (<ref>). The subscript `anc' refers to the ancilla transmon qubit which is initialized to|0⟩_ancand guaranteed to end in|0⟩_ancat the end of the sequence unless an error occurred – a fact which could be used for error-detection purposes as previously mentioned.



  §.§.§ Light-matter interactions
 When light interacts with charged matter on a lattice, for instance when a material is placed in an optical cavity, the Hamiltonian is given by <cit.>
    Ĥ= ωâ^†â-J∑_⟨ i,j⟩ĉ^†_i ĉ_j e^ig/√(L)(â^† + â) +h.c.

    for a single light modeâ, whereωcorresponds to the free evolution of the photons,Jis the strength of the hopping,gis a coupling constant which depends on the specifics of the system, such as the geometry and material composition of the cavity, andLis the number of lattice sites. Within the Peierls substitution approximation, the exponential term,Ã=g/√(L)(â^†+ â), is given by the particle charge multiplied the line integral of the electromagnetic vector potential along the link from sitejto sitei.
    
    For smallg/√(L)the exponent in the second term can be expanded:
    
    Ĥ   = ωâ^†â-J/2∑_⟨ i,j⟩(ĉ^†_i ĉ_j + h.c.) 
       +g J/√(L)(â^† + â)∑_⟨ i,j⟩ĉ^†_i ĉ_j + O(Jg^2/L).

    The first term can be implemented using a bosonic phase-space rotation as discussed in Sec. <ref>. The second term corresponds to fermionic hopping, and can be implemented using the strategy of App. <ref> for fermions encoded in either transmons or dual rail qubits. The second term can be reexpressed asexp(-itg J/√(L)(â^†+ â)X̂_i X̂_i+1)exp(-itg J/√(L)(â^†+ â)Ŷ_i Ŷ_i+1)for Jordan-Wigner adjacent sites. This corresponds to two doubly-conditional displacements, i.e., conditioned on two qubits, with single qubit rotations.
    
    Transmon qubits. If the qubits are mapped to transmon qubits, doubly-conditional displacements can be compiled using the technique described in Sec. <ref> and illustrated in Fig. <ref>b. It requires a pair of conditional parity gatesCΠ(Tab. <ref>) and a conditional displacementCD(α)(or alternatively, two pairs of conditional parity gates and an unconditional displacement).
    
   Dual-rail qubits. If the qubits are mapped to dual-rail qubits, doubly-conditional displacements can be synthesized by iterating the compliation strategy for the singly-conditional displacement in Eq. (<ref>), there developed in the context of phonon-matter interactions. The full sequence is as follows:

    e^-iθẐ^DR_i Ẑ^DR_j (â_k^† + â_k)=
       (⟨0|_ anc⊗ )CΠ^X_ anc,iCΠ^X_ anc,jCD_anc,k(iθ)
       ×CΠ^X_ anc,jCΠ^X_ anc,iR_i(π)R_j(π)(|0⟩_ anc⊗),

where, similar to previous gates, we have projected out the ancillary qubit, which is initialized to|0⟩_ancand guaranteed (in the absence of noise) to be disentangled at the end of the sequence. With the above gate realized, single dual rail qubit gates (i.e., beam splitters) can be used to rotate the operatorẐ^DR Ẑ^DRtoX̂^DR X̂^DRorŶ^DR Ŷ^DR, enabling the implementation of all desired gates.






    
 

§ IMPLEMENTATION OF GAUGE FIELDS AND GAUGE-MATTER COUPLING

In this section, we present implementations ofe^-iĤ tfor pure gauge Hamiltonian terms – i.e., for the electric and magnetic fields – and for the gauge-invariant hopping which couples the gauge fields to matter. In particular, we discuss the implementation of these terms for the two illustrative examples in(2+1)D: aℤ_2gauge field and aU(1)gauge field, each playing a vital role in the paradigmatic LGTs discussed in Sec. <ref>. However, we emphasize that due to the digital nature of our implementations, more complex models can implemented beyond those in Sec. <ref> using techniques in this section, for instance models that contain bothℤ_2andU(1)fields coupled to the same matter fields. 

We present the implementation of pureℤ_2andU(1)gauge fields in Sec. <ref> and Sec. <ref>, respectively. 
 Following this, we briefly discuss the possibility to realize static gauge fields (which do not posses any of their own dynamics) in Sec. <ref>. We then return to our two examples and describe the implementation of gauge-invariant hopping terms in Sec. <ref> and Sec. <ref> forℤ_2andU(1)fields, respectively. Finally, we present a summary of the implemented terms  in Tab. <ref> and Tab. <ref>.



 §.§ Pure gauge dynamical Z2 fields


In this section, we present two strategies for representingℤ_2gauge fields in our architecture – one that uses the transmon qubits, and another the relies on the dual-rail encoding  discussed in Sec. <ref>. Also see Sec. <ref> for the introduction of theℤ_2Hamiltonian and the notation, in particular how we label the links on which the gauge fields reside.



  §.§.§ Electric Field

The time evolution operator of the electric field defined in Tab. <ref> is e^-i g X̂_i,j t , where the coupling strength isgas defined in Eq. (<ref>). This corresponds to a single-qubit rotation.

Transmon qubits.
This term is directly realized through on-resonant driving of each transmon using theR_i^ϕ=0(θ)gate defined in Tab. <ref>.

Dual-rail qubits.
This term is a single-qubit rotation. As discussed in Sec. <ref>, single-qubit rotations (defined in Tab. <ref>) around theX̂^DRandŶ^DRaxis  are carried out using beamsplittersBS_i,j(φ,θ)between the two cavities composing the dual-rail qubit.



  §.§.§ Magnetic field


The time evolution operator associated with the magnetic field term defined in Tab. <ref> ise^(- i t (-B) Ẑ_i,j Ẑ_j,k Ẑ_k,l Ẑ_l,i), corresponding to a four-qubit entangling gate.

Transmon qubits. The exact gate sequence that implements this term is discussed in Sec. <ref> and illustrated in Fig. <ref>. Multi-body gates in the transmon qubits in our proposed circuit QED architecture have a compact form: the required displacements can all be carried out using a single ancilla mode. This removes the requirement for a costly Pauli-gadget [100], which implements this four-qubit gate using 6 two-qubit gates, though instead requires bosonic SWAP gates to mediate the multi-qubit interaction. However, as previously mentioned, the bosonic SWAP operation has a very high fidelity in the proposed architecture (see Sec. <ref>).

Dual-rail qubits.
In Sec. <ref> we discussed how to obtainZZ_i,j(θ)gates using the exponentiation gadget presented in Ref. <cit.> which implementse^iθP̂for any operatorP̂ifP̂^2 = 1. Therefore, using the same method we can extend this to include higher-weight Pauli strings inP̂. This circuit is illustrated in Fig. <ref>.
As mentioned previously, ancilla relaxation and dephasing errrors, in addition to photon loss, can be detected using this general gate scheme <cit.>.


 §.§ Pure gauge dynamical U(1) fields


Next, we turn to the  described by the Hamiltonian in Eq. (<ref>). As described in Sec. <ref>, we utilize the Schwinger-boson representation to encode the link degrees of freedom, requiring two bosonic modes per link. We note that the Schwinger-boson representation may be understood as an extension of the dual-rail qubit encoding to spinS > 1/2<cit.>.



  §.§.§ Electric field


  In the Schwinger-Boson representation, the electric field energy can be written as 

    Ĥ_E = g^2/2∑_⟨i,j|⟩(n̂^a_i,j - n̂^b_i,j/2 - τ/2π)^2,

wheregis the strength of the electric field andτhere is the topological parameter. When multiplying out the square, the cross-Kerr terms-g^2/2 n̂^a_i,j n̂^b_i,jappear. Although we can in principle implement these terms, c.f. Sec. <ref>, they are harder to implement than self-Kerr terms on the same site. Luckily, we can remove the cross-Kerr terms by using the Schwinger-Boson constraintn̂^a_i,j + n̂^b_i,j = 2S . Due to this identity, we may add a termg^2/2∑_i,j (n̂^a_i,j + n̂^b_i,j/2)^2, which due to the constraint is just a constant, to the Hamiltonian in Eq. (<ref>) to obtain

    Ĥ_E   = g^2/2∑_⟨i,j|⟩( 1/2((n̂^a_i,j)^2 + (n̂^b_i,j)^2) -τ/2π (n̂^a_i,j - n̂^b_i,j))
       ≡∑_⟨i,j|⟩ĥ_E,i,j

up to a constant, thereby removing the cross-Kerr term.

The time evolution underĤ_Efulfillsexp(-iĤ_E t)=∏_⟨i,j|⟩ exp(-iĥ_E,i,jt). Because theâandb̂modes commute,exp(-iĥ_E,i,jt)=e^-ig^2/2(1/2n̂_i,j^a^2-τ/2πn̂^a_i,j)te^-ig^2/2(1/2n̂_i,j^b^2+τ/2πn̂^b_i,j)t. This can be implemented exactly using twoSNAP(θ⃗)gates defined in Tab. <ref> on modesaandb:

    exp(-iĥ_E,i,j t)=SNAP_a(θ⃗^a)SNAP_b(θ⃗^b)

whereθ^a_n = g^2/2(1/2n^2-τ/2πn)tandθ^b_n=g^2/2(1/2n^2+τ/2πn)t.




  §.§.§ Magnetic field

In(2+1)D, the magnetic field Hamiltonian, illustrated in Fig. <ref>a, acts on the four gauge field raising operators linking a square (also called plaquette) of fermionic sites labelledi,j,k,l. We discuss the compilation of a single plaquette exponential termexp(-iĤ_□t)here, where

    Ĥ_□=    -1/4g^2(S(S+1))^2(Ŝ^+_i,jŜ^+_j,kŜ^-_k,lŜ^-_l,i +h.c.).

Using our Schwinger-boson mapping of the gauge fields in Eq. (<ref>), we obtain a product of creation and annihilation operators which we relabela,bfor the modes of link⟨i,j⟩,c,dfor link⟨j,k⟩,e,ffor link⟨k,l⟩andg,hfor link⟨l,i⟩for simplicity of notation, such that:  

    Ĥ_□=    -1/4g^2(S(S+1))^2(â^†b̂ĉ^†d̂f̂ê^†ĥĝ^† + h.c.).

We show here how to realize the plaquette term using the method introduced in Ref. <cit.> and summarized in Sec. <ref>, which makes use of the group commutator relations given by the Baker-Campbell-Hausdorff (BCH) formula, specifically the relation 

    e^iÂθ e^iB̂θe^-iÂθ e^-iB̂θ=e^-[Â, B̂] θ^2+O(θ^3).

Through appropriate choice of hybrid mode-qubit operatorsÂandB̂, this gives us the means to implement multiplication between commuting bosonic operators. We also use Trotter formulas to implement addition.

We rely on the above method and present a possible implementation that leverages the conditional beamsplitter gate introduced in Eq. (<ref>) as the fundamental building block. This approach has the advantage of relying on simple hardware primitives, at the cost of Trotter error. Following this, we briefly comment on the possibility to instead leverage a three-wave mixing term realized on the hardware level, c.f. App. <ref>. This latter approach would enable the compilation of the desired term without Trotter error, but uses a gate which is more complicated to engineer on the hardware level. 

The key idea for this synthesis is to apply the BCH formula to conditional beamsplitters acting on the four pairs of modes around the plaquette, conditioned on the same ancillary qubit. Intuitively, this leads to multiplication of creation and annihilation operators for these four pairs of modes, yielding the eight-mode term in Eq. (<ref>).

We first show how to synthesise a product of four mode operators. Choosinge^-i Â θ = CBS^Y_a,b(-π/2,θ)= e^-iθ(iâ^†b̂ - iâb̂^†)Ŷande^-i Ĉ θ=CBS^Z_c,d(0,θ)= e^-iθ(ĉ^†d̂ + ĉd̂^†)Ẑ, whereĤis a Hadamard gate on the qubit andCBS(φ,θ)is defined in Eq. (<ref>), we define the primitive 

    V(Aθ,Cθ) = e^i Âθ e^i Ĉθe^-i Âθe^-i Ĉθ.

This set of operations is illustrated in Fig. <ref>b. In order to reduce the error sufficiently, we use this primitive in a higher-order product formula <cit.>
    Γ_A,C(θ)
       = V(Aγθ,Cγθ) V(-Aγθ,-Cγθ) V(Cγ^2 θ, Aγ^2 θ) 
       × V(-Cγ^2 θ, -Aγ^2 θ)V(Aγθ,Cγθ) V(-Aγθ,-Cγθ)
       = e^-2θ^2( â^†b̂ĉ^† d +â^†b̂ĉd̂^†- âb̂^†ĉ^†d̂- âb̂^†ĉd̂^†) ⊗X̂+𝒪(θ^5),

whereγ=1/√(4-2√(2)).





This expression contains two terms which are not part of the plaquette operator. In order to remove the unwanted terms, we apply the same four-mode synthesis again with different phases. This yields a new termΓ_A',C'obtained using the primitiveV(A'θ,C'θ)withe^-i Â' θ = CBS^Y_a,b(0,θ)ande^-i Ĉ' θ =CBS^Z_c,d(-π/2,θ). CombiningΓ_A,CandΓ_A',C'in a second-order Trotter formula, we remove the unwanted terms, 

    Δ_A,C   = Γ_A,C(θ/2)Γ_A',C'(θ)Γ_A,C(θ/2)
       = e^((-2i â^†b̂ĉ^†d̂+2i âb̂^†ĉd̂^†) ⊗ 2i X)(-i θ)^2+𝒪(θ^5),

This set of operations is illustrated in Fig. <ref>c and we dropped the𝒪(θ^6)Trotter error as the𝒪(θ^5)BCH error dominates. 
Using a single-qubit gate on the qubit, one can also obtain(Δ_A,C)^†. The same expression can be obtained on the other four modes which are also involved in the plaquette term: 
    Δ̃_E,G= e^((-2i ê^†f̂ĝ^†ĥ+2i êf̂^†ĝĥ^†) ⊗ 2i Ŷ)(-i θ)^2+𝒪(θ^5).
 Note thatΔ̃_E,Gis chosen such thatŶis carried out on the qubit rather thanX̂– this is to enable the further use of the BCH formula which we explain next.

In order to obtain the eight-wave mixing required to realise the plaquette term, we concatenate four four-mode terms according to 
    Ξ_{A,C},{E,G} = Δ_A,CΔ̃_E,GΔ_A,C^†Δ̃_E,G^†.
 This set of operations is illustrated in Fig. <ref>d.
Similarly to the four-mode term, equation (<ref>) yields terms which are not part of the plaquette term. In order to remove them, we multiplyΞ_{A,C},{E,G}with a term which is obtained in the same way but with a different choice of phases in the beamsplitters:

    Ξ_2,{A_2,C_2},{ E_2, G_2}= Δ_A_2,C_2Δ̃_E_2,G_2Δ^†_A_2,C_2Δ̃^†_E_2,G_2.

The constituent operations part of this term are the following, where we use the same notation with primes and tildes as in the discussion of the synthesis ofΞ_{A,C},{E,G}. On modesa,b,c,dwe use 

    A_2   = Ĥ CBS_a,b(π,1)  Ĥ, 
    C_2   =CBS_c,d(0,1),
    A'_2   = Ĥ CBS_a,b(π/2,1)  Ĥ, 
    C'_2   =CBS_c,d(π/2,1).

On modese,f,g,hwe use 

    Ẽ_2   = Ĥ CBS_a,b(0,1)  Ĥ,
    G̃_2   =Ĥ Ŝ CBS_c,d(0,1)Ĥ Ŝ^†, 
    Ẽ'̃_2   = Ĥ CBS_a,b(π/2,1)  Ĥ, 
    G̃'̃_2   =Ĥ Ŝ CBS_c,d(-π/2,1)Ĥ Ŝ^†.

In total, we obtain

    Ξ_{A,C},{E,G}Ξ_2,{A_2,C_2},{ E_2, G_2}=
       e^-i(16( â^†b̂ĉ^†d̂f̂ê^†ĥĝ^† + âb̂^†ĉd̂^†f̂^†êĥ^†ĝ) ⊗Ẑ) θ^4 +O(θ^5),

illustrated in Fig. <ref>e. Choosingθ= (1/64g^2(S(S+1))^2t)^1/4we get time evolution under a single plaquette Hamiltonian for a  timestept.
Each CBS gate can be implemented using2conditional parity operations and1beamsplitter, c.f. Eq. (<ref>). We count conditional parities in the following. EachΓcan be implemented using6Voperations (see Eq. (<ref>)), which require8conditional parities to implement and hence eachΓrequires48conditional parities. EachΔoperation requires3Γoperations using the symmetric Trotter formula implying that it requires144conditional parities. EachΞis built from4Δterms and hence each requires576conditional parities. Then the final approximation,V_□requires1152conditional parities as it is composed of twoΞoperations.

We close with a discussion of the error term. Because each pair of mode operators has norma^†b=√(S(S+1))in the Schwinger-boson encoding of a large spin and the error term contains5such pairs, our result implies that the distance between the thus-implemented unitaryV_□and the exact plaquette exponential is

    V_□ - e^-i(16( â^†b̂ĉ^†d̂f̂ê^†ĥĝ^† + âb̂^†ĉd̂^†f̂^†êĥ^†ĝ) ⊗Ẑ) t
         = 𝒪(t^5/4/g^5/2).

In particular, this error is independent ofS.

Another way of synthesising this term would be to use the three-wave mixing term introduced in App. <ref>. Because the additional terms appearing in the beamsplitter approach mentioned above would not appear using three-wave mixing, we could save on Trotter error and reduce the number of native operations needed by a prefactor. However, because the Trotter error is sub-leading, this will not change the asymptotic error scaling shown above.



 §.§ Static gauge fields coupled to matter


Static gauge fields appear when a magnetic field <cit.> is coupled to matter and are crucial when, for example, considering the fractional quantum Hall effect <cit.>. In this case, the hopping term of the Hamiltonian has a fixed complex amplitude on each link, i.e.,

    -J∑_⟨i,j|⟩(e^iφâ^†_iâ_j + h.c.),

where we have specialized to the case of bosonic matter. This term can be directly implemented by a beamsplitterBS_i,j(φ,θ)defined in Tab. <ref> via the appropriate choice for the phaseφfor each link.

For fermionic matter, the hopping term generalizes as

    -J∑_⟨i,j|⟩(e^iφĉ^†_iĉ_j + h.c.).

The phase can be implemented as a part of the gates necessary for the Jordan-Wigner encoding described in App. <ref>, for either  choice of qubit type. TheFSWAPoperation would remain unchanged, but the phase can be incorporated by conjugating the fermionic hopping bye^ iφẐon one of the qubits.




 §.§ Dynamical Z2 fields coupled to matter


Here we discuss  gauge fields coupled to bosonic and fermionic matter. We first discuss the bosonic gauge-invariant hopping in Sec. <ref>, mapping the gauge fields to transmon and then dual rail qubits. We then discuss the fermionic gauge-invariant hopping in Sec. <ref>, mapping the fermions and the gauge fields to transmon and dual rail qubits.



  §.§.§ Bosonic gauge-invariant hopping



The time evolution operator of the  gauge-invariant bosonic hopping mediated by aℤ_2gauge field takes the form 
    e^-i t Ẑ_i,j(â^†_i â_j + â_i â^†_j).

In the following, we discuss the implementation of this operator using either transmons or dual rail (DR) qubits. For both cases, we use modes to directly encode the bosonic matter.

Field: transmons | Matter: modes. For the situation where we map the gauge fields onto transmon qubits, Eq. (<ref>) corresponds to a conditional beamsplitter, discussed in Sec. <ref>.


Field: dual rail | Matter: modes. For the case where we map the gauge fields onto dual-rail qubits, we leverage the fact  that the parity of the photon numbern̂in one of the rails is equivalent to the Pauli operatorẐ^DR(see Eq. (<ref>)). Therefore the time evolution of the hopping between sitesiandjtakes the form

    e^-i t e^iπn̂_i,j(a_i^† a_j+ â_i â^†_j).

and can be compiled using the methods described in Sec. <ref> and shown in Fig. <ref>.



  §.§.§ Fermionic gauge-invariant hopping


The fermionic SWAP networks used to implement fermionic gauge-invariant hopping require the exact same FSWAP operations as are discussed in App. <ref>, where  the gauge field also needs to be swapped. To make the hopping gauge invariant, we replace the iSWAP operation with 

    e^-i t (-J) Ẑ_i,j(X̂_i X̂_j + Ŷ_i Ŷ_j)
       =e^-i t (-J) Ẑ_i,j(X̂_i X̂_j)e^-i t (-J) Ẑ_i,j(Ŷ_i Ŷ_j)
    
        =   R^π/2_i(π/2)R^π/2_j(π/2)e^-i t (-J) Ẑ_i,j(Ẑ_i Ẑ_j)R^π/2_i(-π/2)R^π/2_j(-π/2)
    ×   R^0_i(π/2)R^0_j(π/2)e^-i t (-J) Ẑ_i,j(Ẑ_i Ẑ_j)R^0_i(-π/2)R^0_j(-π/2),

using the fact thatX̂ = R^π/2(π/2) Ẑ R^π/2(-π/2)andŶ = R^π(π/2) Ẑ R^π(-π/2) , whereR^φ(θ)is defined in Tab. <ref>. We note that this gate sequence can be further compressed by combining successive single-qubit gates into a single rotation.

We do not discuss the all-transmon qubit implementation of this model as this case does not leverage the advantages of our platform.

Field: transmon | Matter: dual rail.
Following Eq. (<ref>), single-qubit rotations can be performed on the dual rail qubits using beamsplitter gates, whilee^-i t (-J) Ẑ_i,j (Ẑ_i^DR Ẑ_j^DR)is a transmonZ-conditional dual railZZ^DRoperation. This latter gate is compiled in Eq. (<ref>) using an ancillary transmon initialized to|0⟩_ancthat is ultimately traced out. However, we can adapt this same sequence to realize the desired gate by replacing the ancilla with the transmon encoding the gauge field in an unknown state (and removing the matrix element projection in Eq. (<ref>)).

Field: dual rail | Matter: transmon. This term requires the implementation ofe^-i t (-J) e^iπn̂_i,j (X̂_i X̂_j + Ŷ_i Ŷ_j)=e^-i t (-J) e^iπn̂_i,j X̂_i X̂_je^i t (-J) e^-iπn̂_i,j Ŷ_i Ŷ_jwhich can be implemented using additional single-qubit rotations frome^-i θe^iπn̂_i,j Ẑ_i Ẑ_j. The latter is an instance of the parity-controlled gateCU^Π̅_i,jintroduced in Sec. <ref>. Specifically, the intermediate gate in Eq, (<ref>) would be replaced by the two transmon qubit gateZZ(θ)implemented in our architecture via the technique in Sec. <ref>.

Field: dual rail | Matter: dual rail. Following Eq. (<ref>), aside from single-qubit rotations, we require implementation ofe^-i t (-J) Ẑ_i,j^DR (Ẑ_i^DR Ẑ_j^DR). This type of multi dual-rail gate can be realized using the Pauli-exponentiation gadget as discussed in Sec. <ref> and shown in Fig. <ref>, there for the slightly more complex case of four dual-rail qubits. Alternatively, one can iterate the parity-control synthesis technique of Sec. <ref>, similar to the doubly-conditional displacement realized in Eq. (<ref>).






 §.§ Dynamical U(1) fields coupled to matter


Next, we discuss the implementation of dynamicalU(1)fields coupled to bosonic (Sec. <ref>) and fermionic matter (Sec. <ref>). We show how to implement the dynamics of the hopping described in Eq. (<ref>) which has the form: 

    exp(-i Jt/2√(S(S+1))∑_⟨i,j|⟩(m̂^†_iŜ^+_i,jm̂_j +  h.c.)).
 

In the remainder of the section, we will assume largeSand write√(S(S+1))≈S. As in Sec. <ref>, we leverage the Schwinger-boson encoding the represent theU(1)gauge fields. Furthermore, for the case of fermionic matter, we separately consider transmons and dual-rail qubits.


  §.§.§ Bosonic gauge-invariant hopping
 

For a bosonic matter siteidenoted by operatord̂, the gauge-invariant hopping between siteiand sitejis given by

    e^-i t J/2S(d̂^†_iâ^†_i,jb̂^†_i,jd̂^†_j + h.c.),

whereâ^†_i,j b̂^†_i,jis the Schwinger boson representation of the link operatorŜ^+_i,j. This unitary is generated by a four-mode mixing term, and can be synthesized using the BCH formula, see Eq. (<ref>). 






  §.§.§ Fermionic gauge-invariant hopping


For fermions in(2+1)D, we first use FSWAP networks to move the sites to neighbouring sites within the JW encoding. To then perform the gauge-invariant hopping, we also SWAP the corresponding gauge fields as depicted in Fig. <ref>. The remaining operation is then given by the nearest-neighbour hopping

    e^-i t J/2S(σ̂^+_iâ^†_i,jb̂^†_i,jσ̂^-_j + h.c.),

wherejis a neighbouring site toi. In(1+1)D, SWAP and FSWAP operations are not necessary. We now discuss how this term is synthesized for the cases where the fermions are encoded in transmons or dual-rail qubits.


Fermionic matter: transmons.
We present two schemes, a naive implementation which contains Trotter error, and a scheme which avoids Trotter error (for this term alone), but requires a particular experimental adaptation of the hardware.

Scheme 1– We rewrite the gauge-invariant hopping term as

    e^-i Jt/2S(Ĥ_1+Ĥ_2)≈ e^-iJt/2SĤ_1e^-iJt/2SĤ_2+𝒪((Jt/2S)^2),

where we used Trotterization and defined

    Ĥ_1   =(X̂_iX̂_j +Ŷ_iŶ_j )(â^†_i,jb̂_i,j+h.c.)
    Ĥ_2   =(Ŷ_iX̂_j -X̂_iŶ_j )( e^iπ/2â^†_i,jb̂_i,j+h.c.).

Noting thatX̂_i X̂_jandŶ_i Ŷ_jcommute, as doŶ_i X̂_jandX̂_i Ŷ_j, we are then left with implementing doubly conditional beamsplitters.  This is done by sandwiching a beamsplitter with  conditional parity operations conditioned on the qubitsiandj, as shown in Eq. (<ref>). However,[Ĥ_1, Ĥ_2]≠0and this term will therefore acquire Trotter error. 

Scheme 2– This scheme does not introduce Trotter error, but requires additional engineering at the hardware level. We define the unitary ase^-iθĤwithĤ=σ_i^+ a^† b σ_j^-+σ_i^- a b^† σ_j^+where, for simplicity we have dropped the indicesi,jon the modesaandband we definedθ= Jt/2S.
We start by breaking upĤintoĤ=Ĉ+D̂, where

    Ĉ   =X̂_i(X̂_j(a^† b+a b^†)+i Ŷ_j(a^† b-a b^†)) / 4 
       ≡X̂_iÊ_j 
    D̂   =Ŷ_i(-i X̂_j(a^† b-a b^†)+Ŷ_j(a^† b+a b^†)) / 4 
       ≡Ŷ_iF̂_j.

Indeed,[Ĉ, D̂]=0because{Ê_j, F̂_j}=0.

It is helpful to rewrite:

    Ê_j    =(a b^†σ_j^++a^† b σ_j^-) / 2 
    F̂_j    =i(a b^†σ_j^+-a^† b σ_j^-) / 2.

With this, we can now express the desired time-evolution operator as 

    e^-i θĤ   =e^-i θD̂ e^-i θĈ
       =e^-i θŶ_i  F̂_j e^-i θX̂_i  Ê_j
       = R^0_i(-π/2) C Π_i,j e^-i θF̂CΠ_i,j R^0_i(π/2) 
       × R^π/2_i(-π/2) C Π_i,j e^-i θÊCΠ_i,j R^π/2_i(π/2),

whereR^φ(θ)(a single transmon qubit rotation gate) andCΠgates are defined in Tab. <ref>.
This reduces the problem of realizing four-body interactions to the more manageable problem of realizing the three-body interaction Hamiltonians – the latter being more natural to realize in the circuit QED architecture given the primary source of non-linearity is the four-wave mixing property of Josephson junctions. We derive one possible implementation of these interactions in App. <ref>.

Fermionic matter: dual-rail.
In this case, we must implement a many-body gate that involves the two Schwinger-boson modes representing the gauge fieldsa_i,jandb_i,j, and the four modes composing the two dual-rail qubitsc_i, d_iandc_j, d_j:

    e^-iθ(ĉ_i^†d̂_iâ^†_i,jb̂^†_i,jĉ_jd̂_j^† + h.c.).

To implement such a gate, we use the BCH formula synthesis methods presented in Sec. <ref> for pureU(1)gauge fields. In particular, we apply Eq. (<ref>), where we sete^i X̂ Ô_1 θ = e^((-2i c_i^† d_i a_i,j^† b_i,j+2i c_i d_i^† a_i,j b_i,j^†) ⊗2i X)(-i θ)^2+O(θ^3)for which the compilation was shown in Eq. (<ref>), and sete^i Ŷ Ô_2 θ = e^-i Ŷ_k ( c_j^†d_j + h.c.) θ. This yields unwanted terms which we cancel using the Trotter formula in Eq. (<ref>) through use of BCH with appropriate phases in the beamsplitters.








§ GATE COMPLEXITY AND COMPARISON TO ALL-QUBIT ALGORITHMS 


In this section, we analyse the gate complexity of our qubit-boson Trotter approach. In particular, we discuss the interplay of Trotter errors with errors introduced by some of our compilation schemes. We first review the general principles of Trotter error analysis. We then discuss a particular Trotter implementation of the 
 and  and its error. While all compilations are exact for the , the compilation of the magnetic field term in the  (c.f. section <ref>) introduces additional error. We therefore discuss the interplay of Trotter error with this compilation error. While these sections focus on the asymptotic scaling of the whole Trotter algorithm, we close the section with a comparison of the explicit gate counts of a single Trotter step in our qubit-boson approach and compare it to an all-qubit approach. To do so, we introduce a qubit algorithm for simulating beamsplitters between bosonic modes.



 §.§ General Trotter error analysis




We write the Hamiltonian as a sum of termsĤ=∑_γ=1^ΓĤ_γ. We first decompose the time evolutionexp(-iĤ T)until timeTintortime steps of durationt=T/r, i.e.exp(-iĤ T)=(exp(-iĤ t) )^r. 
An estimate of the error 

    ϵ_Trotter≡(𝒮(t))^r-e^-iĤT

of a2p-th order Trotter formula𝒮(t)as introduced in section <ref> is <cit.>
    ϵ_Trotter=𝒪((∑_γĤ_γ)^2p+1T^2p+1/r^2p).


The number of Trotter steps for a target Trotter error is hence

    r   =𝒪((∑_γĤ_γ)^2p+1/2pT^2p+1/2p/ϵ_Trotter^1/2p)
       =𝒪((∑_γĤ_γT)^1+1/2pϵ_Trotter^-1/2p).

From this expression, we get the total number of gates as 
    N_gates,2p   =r5^p-1 N_gates,2p=2
       =r5^p-1 2 ∑_γ N_gates,γ,
 
whereN_gates,γis the number of gates required for the individual exponentialsexp(-iĤ_γt).  We used that a Trotter formula of order2prequires5times the number gates of a Trotter formula of order2p-2when defined recursively <cit.>. The second-order Trotter formula requires twice as many gates as the first-order Trotter formula. Finally, the number of gates required for the first-order Trotter formula is the sum of the number of gatesN_γrequired for each exponential. In cases where the individual exponentials can only be synthesized with error (dependent ont),N_gates,γcontributes to the overall complexity. Note that in general, the complexity of Trotter formulas is evaluated at fixedp. Hence, factors such as5^p-1above are dropped for estimating the asymptotic complexity in𝒪notation.


Approximating the Trotter error in terms of the norms of the Hamiltonian terms overestimates the error because most terms in many-body Hamiltonians commute with each other. We therefore expect the Trotter error to be much smaller in practice, see e.g. <cit.> for a tighter estimate for the(1+1)D fermionic .

Bosonic Hamiltonians have an in-principle unbounded norm and therefore naively exhibit an unbounded product-formula error. However, in practice, infinitely high Fock states are not occupied during the dynamics such that the error stays bounded. In particular, in the , the spectral norm of bosonic operators is bounded by√(S). Similarly, in the , using an initial state with a total number ofNbosons in the system, bosonic operators are bounded by√(N).

We detail specific Trotter decompositions of the  and  and their errors next. We specifiy to open boundary conditions on anL ×Lsquare lattice and use the triangle inequalityĤ_0+Ĥ_1≤Ĥ_0 +Ĥ_1 as well asaĤ_0=|a|  Ĥ_0. For bosonic matter, we assume evolution within a sector of total boson numberNto estimate the Trotter error.





 §.§ Trotter decomposition and gate complexity of the Z2-Higgs model
 
The electric field and on-site interaction term commute with each other such that we can choose 
    Ĥ_0=-g∑_⟨i,j|⟩X̂_i,j + U∑_i=1^L^2n̂_i^2,

where⟨i,j|$⟩ indicates summation over the 2L(L-1) bonds on a square lattice with L^2 sites. It follows from the triangle inequality and X̂_i =1, n̂_i=N in a fixed particle sector with N particles that 
    Ĥ_0 ≤ 2gL(L-1)+|U|L^2N^2.

Next, we split the gauge-invariant hopping into four non-commuting terms ξ=1,2,3,4, where 1(2) denote the even(odd) hoppings in the horizontal direction and 3(4) the even(odd) hoppings in the vertical direction, c.f. Fig. <ref>a).  We write this as 
    Ĥ_1^ξ=-J∑_⟨i,j|_⟩ξ(m̂^†_iẐ_i,jm̂_j +  h.c.).
 Next we find m̂^†_iẐ_i,jm̂_j +  h.c.≤2m̂^†_iẐ_i,jm̂_j=2N, where we used that m̂^†_i= √(N)=m̂_i. From this, we get 
    ∑_ξ=1^4Ĥ_1^ξ≤ 4JNL(L-1),

where N=1 for fermions.
Finally, we divide up the plaquette terms into the two non-overlapping checkerboard patterns of the square lattice, c.f. Fig. <ref>b), which we label with ξ=1,2, such that 
    Ĥ_2^ξ=-B∑_i,j,k,l∈□_ξ(Ẑ_i,jẐ_j,kẐ_k,lẐ_l,i),


and from the triangle inequality,

    ∑_ξĤ_2^ξ≤ B(L-1)^2,

where we used that there are (L-1)^2 plaquettes.

Inserting Eqs. (<ref>), (<ref>), (<ref>) into Eq. (<ref>),  we find a number of Trotter steps given by

    r=
       𝒪(((2g+|U|N^2+4JN+B)L^2T)^1+1/2p1/(ϵ_Trotter)^1/2p).

where N=1 and U=0 for fermionic matter.

In particular, we find that for bosonic matter the interaction term dominates the error. If U=0, instead the hopping term is the leading source of error. For fermionic matter, both hopping and field term contribute equally to the Trotter error.

All our compilations of the terms in the are exact and therefore N_gates,γ=𝒪(1). There are (L-1)^2 plaquette terms on a square lattice. Therefore, we find that the  total number of gates is given by

    N_gates,2p=r(L-1)^2 5^p-1 2 ∑_γ N_gates,γ=
       𝒪(((2g+|U|N^2+4JN+B)L^2T)^1+1/2pL^2/(ϵ_Trotter)^1/2p)

for fixed p. This means that for a large value of p and N=𝒪(1) in L, the number of gates approximately scales as L^4, i.e. the square of the total number of sites. This makes sense as there is a number of terms in the Hamiltonian which is linear in the number of sites and in our crude estimate, the number of Trotter steps r for a certain target error is linear in the number of terms in the Hamiltonian. We note that in most interesting situations, N∝ L^2, such that the scaling in L will be higher and in particular, the on-site term dominates the scaling with L in that case, i.e. approximately L^8 for large p. However, we expect this scaling to be much improved when considering tighter bounds on the Trotter error <cit.>.



 §.§ Trotter decomposition of the U(1) quantum link model


The electric field and mass terms commute with each other and between lattice sites. Hence, we choose 

    Ĥ_0=g^2/2∑_⟨i,j|⟩(Ŝ^z_i,j+τ/2π)^2 + M∑_i(-1)^i n̂_i.


Using that spin S operators fulfill Ŝ^z =S and n̂_i=N, we get from the triangle inequality that

    Ĥ_0 ≤g^2/2L(L-1) (S+τ/2π)^2+MNL^2.

Further, we split the gauge-invariant hopping in the same way as for the , i.e. into the two checkerboard patterns of the square lattice, writing 

    Ĥ_1^ξ=J/2√(S(S+1))∑_⟨i,j|_⟩ξ(m̂^†_iŜ^+_i,jm̂_j +  h.c.).
 Because Ŝ^+ =Ŝ^- =√(S(S+1)), the S-dependence of the norms is cancelled by the prefactors and we find that 

    ∑_ξĤ_1^ξ≤ NJL(L-1).

Finally, akin to the , we use the checkerboard decomposition of the plaquette term, which we write as 

    Ĥ_2^ξ=-1/4g^2(S(S+1))^2∑_i,j,k,l∈□_ξ(Ŝ^+_i,jŜ^+_j,kŜ^-_k,lŜ^-_l,i +h.c.).
 Again, there is no S-dependence of the norm. Because there are are (L-1)^2 plaquettes, the norm is

    ∑_ξĤ_2^ξ≤1/2g^2 (L-1)^2.
 Inserting Eqs. (<ref>), (<ref>),(<ref>) into Eq. (<ref>), we find

    r=𝒪(((1/2g^2S^2+MN+NJ+1/2g^2)L^2T)^1+1/2p
                   ×1/(ϵ_Trotter)^1/2p),

where we dropped some subleading contributions in S and L. Again, N=1 for fermionic matter. For large S, we find

    r=𝒪(((gSL)^2T)^1+1/2p1/(ϵ_Trotter)^1/2p).


Note that a further decomposition of the Ĥ_n^(ξ) into terms that only act locally yields no further Trotter error as the constituent Hamiltonians commute with each other. For example, exp(-iĤ_2^ξ t)=∏_i,j,k,l∈□_ξexp(-it/4g^2(S(S+1))^2(Ŝ^+_i,jŜ^+_j,kŜ^-_k,lŜ^-_l,i +h.c.)) with no error.

Contrary to the , not all our compilations are exact for the . Specifically, the magnetic field term is compiled using the BCH formula, hence leading to error dependent on the timestep T/r. Because N_gates,γ for the magnetic field term therefore contributes to the overall gate count and we need to include this error in the analysis. This is the topic of the following subsection.






 §.§ Gate complexity of the
U(1) quantum link model





We focus on the simulation of the pure gauge theory, i.e. only including electric and magnetic field terms to keep the equations simpler to read. The gauge-invariant hopping can be synthesized without error and the Trotter error is in any case dominated by the electric field term for large S; therefore, the analysis below generalizes trivially to the 
with matter. We first discuss the BCH synthesis error encountered in the plaquette term and then combine this with the Trotter error. 

First, consider the error of the synthesis of a single plaquette term, which we showed in section <ref> to be given by 𝒪(g^-5/2t^5/4). This error can be further reduced by splitting the evolution in r_□ steps of time t/r_□ (c.f. Fig. <ref>c), yielding 

    ϵ_□   =r_□𝒪(g^-5/2t^5/4r_□^-5/4)
       =𝒪(g^-5/2t^5/4r_□^-1/4)

for the error ϵ_□ to evolve under the Hamiltonian of a single plaquette until time t. Solving for r_□, we get

    r_□   =𝒪(g^-5/2t^5/4r_□^-5/4)
       =𝒪(g^-10t^5ϵ_□^-4).

Because the BCH error depends on the timestep t chosen in the Trotterization scheme, r_□ depends on the relation chosen between the target Trotter error and the BCH error. We make the choice that the total BCH error incurred to synthesize all plaquette terms for all Trotter steps scales the same as the Trotter error ϵ_Trotter. Because there are (L-1)^2 plaquette terms in the magnetic field Hamiltonian, and a 2p Trotter formula calls the magnetic field term 2× 5^p-1 times in each of the r Trotter steps, this total error is given by

    ϵ_BCH=2 5^p-1 r(L-1)^2 ϵ_□.

Solving this equation for ϵ_□ and demanding ϵ_BCH=ϵ_Trotter, we find 

    ϵ_□ =𝒪(ϵ_Trotter/rL^2).

Inserting this equation into Eq. (<ref>), we find

    r_□   =𝒪(t^5r^4L^8/ϵ_Trotter^4g^10)
       =𝒪(T^5L^8/rϵ_Trotter^4g^10),

where in the second step we inserted t=T/r.

There are 𝒪(L^2) plaquette and  electric field terms in each Trotter step. The electric field term can be synthesized with 𝒪(1) gates, while the plaquette term requires 𝒪(r_□) gates. Hence,  the total gate complexity for implementing the U(1) gauge theory is given by

    N_gates,2p=𝒪(rL^2(1+r_□)).
    
Inserting Eq. (<ref>) for r and Eq. (<ref>) for r_□, we find

    N_gates,2p=
    𝒪(L^2(gSLT)^1+1/2p/(ϵ_Trotter)^1/2p+T^5L^10/ϵ_Trotter^4g^10),

where the first(second) term inside the bracket is the electric field (magnetic field) contribution. 

This expression shows that the contribution to the total number of gates coming from the magnetic field term (the second term) is independent of the cutoff S. This is in contrast to qubit methods in the Fock-binary encoding, where this contribution scales as 𝒪(log(S)) <cit.>. This improvement comes at the cost of a high polyonmial scaling with ϵ_Trotter, T and L. However, this scaling can be systematically improved by using higher-order BCH formulas in the synthesis of the magnetic field term <cit.>.



 §.§ Gate-complexity advantage over all-qubit hardware


 While in the previous sections, we were considering end-to-end asymptotic gate counts for the whole algorithm, here we compare our hybrid qubit-oscillator compilation to an all-qubit compilation by calculating explicit gate counts for implementing single Hamiltonian terms used in a Trotter decomposition. As most Hamiltonian terms only require a handful of qubit-boson gates, this is straight-forward in our qubit-boson approach. By contrast, compiling qubit-boson operations onto qubits is highly non-trivial and hence, most of this section is about constructing a suitable qubit algorithm and counting its gate requirement. To do so, we must first choose a mapping of bosons to qubits and calculate the gate complexity within this qubit approach. A number of different boson-to-qubit mappings exist: the unary encoding <cit.> which is wasteful in number of qubits but straightforward for gate application, the Jordan-Lee-Preskill encoding for phase-space-formulated bosonic Hamiltonians <cit.>, the Gray encoding which requires minimal numbers of bit flips when transferring between Fock states <cit.>, and the binary encoding <cit.> of the Fock-state basis. See also Ref. <cit.> for a discussion of the complexity of simulating bosons in gauge theories.

The Fock-binary encoding – for which the action of a displacement operator is explained in App. F of <cit.> – is efficient in terms of qubit number and operations because native Fock space operations usually contain diagonal elements which can be transformed into diagonal qubit operations such as qubit rotations. The encoding of the state in binary requires a register of qubits as large as n=⌈log_2(N_max+1)⌉ where N_max is the maximum Fock state which can be accessed during the simulation. The protocol entangles this register with an ancilla qubit. The ancilla qubit stores the relative amplitude of Fock states 2j and 2j+1 and the register state 2j+1 hold the information needed to convert the relative amplitudes to absolute amplitudes.

One severe limitation of simulating bosonic operations with qubits is the prefactor coming from the action of the creation and annihilation operators: â^†|k⟩=√(k+1)|k+1⟩ and â|k⟩=√(k)|k-1⟩, where |k⟩ is the Fock state. When performing quantum simulation in regimes intractable with classical methods, for models with bosonic matter (such as the we study here) it is essential to calculate the square-root factors as they are crucial for understanding the physics of the model. One way of calculating a square-root with qubit-only hardware on-the-fly is to first calculate the inverse square root using Newton iterations, and then multiplying the result with the initial value to obtain the square-root. This results in <cit.>

    N_CNOT / SQRT=
       (270m+126)n^2 + (228m+96)n-12m

 operations, where n is the number of qubits required to represent the maximum Fock state in binary and m is the number of Newton iterations required for a set precision.  As Newton's method converges quadratically, m=2 or m=3 is usually sufficient to obtain errors less than 10^-4 <cit.>. Therefore, at least when using a naive arithmetic approach, the square-root requirement leads to an extremely large overhead of thousands of CNOT gates when only using qubits as we show below.

We calculate here the number of CNOT gates required to perform a single Trotter step of the in (1+1)D on all-qubit hardware using the Fock-binary encoding, keeping the gauge fields. We do not consider here the magnetic field term which appears in higher dimensions. However, for ℤ_2 gauge fields which can be mapped to qubits, the qubit-only and the qubit-oscillator circuits to implement this term will be the same. Therefore no additional Fock-binary gate count calculation and comparison would need to be made in higher dimensions for this model. The CNOT gate count originates from expressing the gauge-invariant hopping term of the bosonic matter in the qubit representation, because the electric field term for ℤ_2 gauge fields corresponds only to single qubit gates (see Sec. <ref>). We do not treat the bosonic Hubbard on-site interaction, however it could be implemented by imprinting the phase corresponding to n^2 in binary using an oracle (e.g., QRAM). The gauge-invariant hopping corresponds to a qubit-conditional beamsplitter term. We convert this to an all-qubit circuit using the Fock-binary encoding for which the gate counts are explained in App. <ref>. We obtain the following scaling in terms of CNOT gates for the qubit-conditional beamsplitter term (with m=2 Newton iterations), i.e. the bosonic in (1+1)D, to be:

    N_CNOT / ℤ_2   = 12 (2673 n^2 + 1160 n - 34).

 
 

In Fig. <ref>, we compare gate counts for the and in both types of hardware. The jumps in the qubit-only lines are due to the logarithmic scaling of the qubit register required to represent increasing Hilbert space cutoffs. The in (1+1)D corresponds to the Schwinger model, for which we report the number of CNOT gates required for one Trotter step to be <cit.>:

    N_CNOT / U(1) = 9n^2-7n+34.

In comparison, the number of entangling gates required to implement these models in hybrid qubit-oscillator hardware is 𝒪(1). The number of conditional parities CΠ required for the and is:

    N_CΠ / ℤ_2   = 2
    
        N_CΠ / U(1)   = 8.




Finally, we consider the corresponding circuit fidelities for the execution of one Trotter step for a system size of 2 sites. This is of utmost importance in the NISQ era where entangling gate fidelities rarely reach above 99.99% <cit.>. In the fidelity calculation, we assume single qubit gates and beamsplitters to have ideal fidelity. We plot the circuit fidelity with regards to the entangling gate fidelity which corresponds in the all-qubit case to the CNOT gate fidelity and in the hybrid oscillator-qubit case to the CΠ fidelity. Both types of hardware are based on the same oscillator-qubit circuits and will therefore encounter the same Trotter error (the bosonic mapping is exact). Therefore, Trotter error is not taken into account in the circuit fidelity.
 
 In the fault-tolerant era, ideal fidelities can only be approached with large distance codes which directly increase the space-time volume of the computation. Therefore, even in the fault-tolerant era for all-qubit hardware, entangling gate fidelities will be of concern, and non-error corrected hybrid oscillator-qubit hardware may provide advantages. This is especially true considering the extremely large circuit depths required when mapping gauge fields or bosonic modes to qubits. In the case of the which contains bosonic matter, the square root factors lead to the requirement for extremely high fidelity operation of a qubit-only device, unless we see substantial algorithmic improvements in arithmetic. In this case, the experimental improvements in hybrid oscillator-qubit hardware may enable progress to be made on bosonic models even in relatively large system sizes such as the ones presented here. Conversely, however, we must note that the harmonic approximations used for native bosonic hardware may break down if large occupation cutoffs are required. For these reasons, this work suggests the existence of a regime of advantage but experimental and algorithmic realities need to be considered in greater detail to verify this advantage.
 
These results suggest that for Trotter-based simulations an advantage may be seen for hybrid hardware; however, a full comparison of qubit to oscillator-qubit devices is challenging to make. For one, a full comparison between the two requires that we compare the best implementations of the best algorithms for both hardware in an agnostic fashion. This is especially significant because alternative methods of simulation, such as qubitization or LCU <cit.> do not require square roots at the price of an increased number of ancillary qubits.  A full comparison would need to consider optimized versions of all available algorithms for both hardware platforms. 
Second, a more detailed discussion of the approximation errors incurred in the approaches (such as the large S limit or the use of BCH approximations and the different errors in the Trotter formula based on the form of the Hamiltonian) is needed to properly compare different approaches. For these reasons, we cannot claim that this work suggests a universal advantage for simulating ℤ_2 or U(1) LGTs, but it does strongly suggest advantages are likely to exist in certain regimes and subsequent work is needed to understand the shape and size of these regions of advantage.



§ MEASUREMENTS, ALGORITHMS AND NUMERICAL BENCHMARKS

In the previous sections, we have introduced all of the necessary compilation methods to simulate both the and (with either bosonic or fermionic matter) in (2+1)D. As an explicit demonstration of these compilation methods in practice, in this Section we apply them to the (1+1)D variants of these two models. In particular, we present schemes to measure short- and long-range observables, and near-term (Trotter and variational approaches) and longer-term (quantum signal processing) algorithms for dynamics and ground-state preparation, which we benchmark numerically.  In particular, to enable approximate ground state preparation using the variational quantum eigensolver (VQE) approach, we show how to measure the expectation value of the Hamiltonian. We also develop an ancilla-free method for using post-selection on the Gauss's law gauge constraint for partial error detection.

We use Bosonic Qiskit for our simulations <cit.>. Each mode is represented with ⌈log_2(N_b+1)+1⌉ qubits, with N_b is the total number of bosons in the system (which is conserved for our Hamiltonians) in order to make sure that the commutation relation for the truncated mode operators obeys [a,a^†]=1 within the accessible Hilbert space. We also include single qubit amplitude and phase damping error channels using Qiskit Aer <cit.> and neglect photon loss and dephasing of the cavities. For transmon qubits, T_2∼ T_1 <cit.>. In tantalum transmons it was shown that T_1>500 μs is possible <cit.>; however, we choose a conservative value of T_1=T_2=200 μs <cit.>.  This implies a pure dephasing time of T_ϕ= (1/T_2-1/2T_1)^-1=400 μ s. Cavities have very little dephasing and photon loss times exceeding 1ms <cit.>, and most of our numerical examples do not access high Fock states, justifying us neglecting these processes.

This section is structured as follows. In Sec. <ref>, we discuss the (1+1)D variants of the ℤ_2-Higgs and U(1) quantum link models, their physics and challenges in previous approaches to solve them. In Sec. <ref> we investigate time dynamics, numerically benchmarking a simulation on the while including the dominant source of noise. In Sec. <ref>, we show that ground states of these models can be prepared using a variational quantum eigensolver (VQE), for which we introduce an optimization technique better adapted to noise. In Sec. <ref>, we analyze in detail how to mitigate the influence of shot noise and decoherence. In Sec. <ref>, we show how to measure observables relevant to identifying phase-transitions, such as the pair-gap, the superfluid density, and the string order correlator with respect to the VQE ground states. In Sec. <ref> we present a perspective on how to implement a ground-state preparation technique with provable success probabilities, known as quantum signal processing.



 §.§ Models and Physics

To begin, we present the two models that we use for explicit demonstration throughout this section: the containing bosonic matter sites, and the U(1) quantum link model containing fermionic matter sites, both in (1+1)D. These models are special cases of the more general LGT Hamiltonians presented in Section <ref>; we discuss the details and implementation of each below. In particular, we can use the architectures presented in Fig. <ref>c),d) for one spatial dimension by winding a “snake” pattern through the square lattice. This way, every cavity and qubit represents either a gauge field or matter site.



  §.§.§ Z2-Higgs model

The ℤ_2 lattice gauge theory was introduced by Wegner <cit.> as a simple model showing signatures of confinement and deconfinement at the critical point <cit.>, making it an ideal testbed for early testing of novel methods for LGT simulations. In two spatial dimensions in the absence of an electric field, its ground state encodes a spin-1/2 degree of freedom and can therefore be used for quantum error correction <cit.>. ℤ_2 lattice gauge theory with bosonic matter is reminiscent of the Higgs sector coupled to the non-Abelian SU(2) electroweak gauge sector (W and Z bosons) of the (3+1)D Standard Model of particle physics
with the two simplifications of a lower dimensionality and a ℤ_2 gauge field instead of the non-Abelian SU(2) electroweak gauge fields. The Hamiltonian in (1+1)D is written as follows:

    Ĥ_ℤ_2   = - g ∑_i=1^L-1X̂_i,i+1  + U ∑_i=1^Ln̂_i^2 
         - J ∑_i=1^L-1( â^†_iẐ_i,i+1â_i+1+ h.c.) ,

where the symbols are the same as described in the introduction (see Section <ref>), and L is the number of sites in the chain.

The ℤ_2 model in (1+1)D in the hardcore boson limit U→∞ has previously been investigated in depth <cit.>, revealing intriguing phase transitions between insulating and Luttinger liquid phases <cit.>. A digital simulation of this model for hardcore-bosons has been performed in Ref. <cit.>. A basic building block of two matter and one gauge field site has been implemented in an experiment with ultracold atoms in a Floquet approach <cit.>. However, this approach is challenging to generalize away from the single particle limit in the presence of on-site interactions U due to Floquet heating, which is important when considering parameter regimes away from the hardcore-boson limit. Our digital implementation enables a fully flexible study of this model, including regimes away from the hardcore limit. While we focus here on the implementation aspects of the model, enabling a study away from the hardcore boson limit, we study its physics and phase diagram in a separate publication <cit.>.








  §.§.§ 


The constrained to (1+1)D and fermionic matter  <cit.> reads


    Ĥ_U(1)   = g^2/2∑_i=1^L-1(Ŝ^z_i,i+1 +  
        τ/2π)^2 + M∑_i=1^L(-1)^i ĉ^†_i ĉ_i 
        + J̃/2S∑_i=1^L-1(ĉ^†_iŜ^+_i,i+1ĉ_i+1 + h.c.),

where compared to Eq. (<ref>), the definition of J differs by J̃=J/2 and we set √(S(S+1))→ S in the denominator of the hopping term, i.e. we consider the limit S→∞.

This model can be connected to a continuum field theory by introducing a lattice spacing a by J→ J/a and g→ g√(a) and taking the limit a→ 0.

This model shares many properties with quantum chromodynamics in (3+1)D: chiral symmetry breaking, confinement, a U(1)_A quantum anomaly, and a topological τ-term, which gives rise to a non-trivial topological vacuum structure <cit.>. Due to the sign problem, this τ-dependence of the model cannot be studied with quantum Monte Carlo methods. One possible way to circumvent the sign problem is using the tensor network (TN) approach, in particular matrix product states (MPS) in (1+1)D. Using MPS, the spectrum of the has been computed, and the model has been studied at non-zero temperature, non-zero chemical potential, and with a non-zero τ-term <cit.> (see Ref. <cit.> for a review). In particular, Refs. <cit.> studied the first-order phase transition at  τ=π, in which the CP symmetry is spontaneously broken. The first-order phase transition line terminates at a second-order critical point <cit.>. This model has been extensively simulated in superconducting qubits <cit.>.

Despite these successes of studying sign-problem afflicted regimes with TN, TN approaches have several limitations. First, highly entangled states such as those with volume law entanglement can occur in the evolution after a quench or in thermal states, and cannot be efficiently described using matrix product states <cit.>, the most widely used TN. This is particularly apparent in multi-stage thermalization dynamics such as after instabilities <cit.>. Second, the computational cost of TN algorithms increases with increasing dimensionality of the lattice field theory under consideration. For MPS, a particular kind of one-dimensional TN, the leading-order computational cost of the variational ground-state optimization in the case of open boundary conditions scales as 𝒪(D_b^3) <cit.>, where D_b, called the bond dimension, indicates the tensor size in the MPS. For the generalization to (2+1)D, called projected entangled pair states (PEPS) <cit.>, computational costs can scale up to 𝒪(D_b^10) <cit.>. These costs are still polynomial in the bond dimension, but substantially limit the values of D_b that can be reached in practice.  There have been TN studies of gauge theories in (2+1)D and (3+1)D <cit.>, but designing efficient TN algorithms in higher dimensions remains a fundamental challenge. 

For this Hamiltonian (Eq. (<ref>)), Refs. <cit.> explored whether the gauge degrees of freedom can be efficiently represented with finite-dimensional systems. For intermediate values of the coupling, g√(a)∼ 0.1, Ref. <cit.> demonstrated that the results converge rapidly, implying that a reasonable accuracy can be obtained even with small cutoffs N_max of the gauge-field truncation. To be more precise, the truncated yielded a fast convergence to the exact ground state for N_max ranging from 3 to 9, reaching sub-permille precision in the energy density for N_max=9 <cit.>. However, when approaching the continuum limit, g√(a)→ 0, the increasing electric field fluctuations require an increasingly large number of states and therefore an increasingly large (up to exponentially large) cutoff. This strongly motivates applying hybrid oscillator-qubit quantum devices to U(1) lattice gauge theory.





 §.§ Dynamics


To illustrate the power of our digital approach to qubit-boson quantum simulation of lattice gauge theories we show that some gauge-constrain-induced dynamics can be observed in near-term noisy hardware. Here, we focus on the ℤ_2-Higgs model for U=0.

We first need to prepare an initial state of interest. As discussed in the exposition of the cQED platform in Section <ref>, it is straightforward to prepare the modes in a product state of Fock states. We consider a L=5-site chain prepared in the states |ψ_0⟩=|00100⟩ and |ψ_1⟩=|00200⟩. The limitation to five sites is not restrictive for short timescales – because the interactions are short-range, excitations spread with a finite velocity and as long as the “lightcone” does not reach the boundary, the results obtained here are the same as would be in a larger system. The qubits representing the gauge fields must then be initialised according to Gauss's law, which requires consistency between the qubits and the parity in each mode: choosing the leftmost gauge field qubit to be in |+⟩ and the sites to be in Fock state |00100⟩, the gauge field qubits situated on the links between the bosonic sites need to be initialized in |++–⟩. Likewise, for bosonic sites in |00200⟩, the corresponding gauge field state is |++++⟩. The fact that the state of the gauge fields is fully determined by the state of the sites shows us that the gauge fields could be integrated out, which is usually the case for lattice gauge theories in (1+1)D. The reverse is not true however: the state of the gauge fields does not fully determine the state of the bosonic sites but only their parity. This is in contrast to hard-core bosons, where the parity fully determines the occupation number <cit.>. See also Ref. <cit.> for a more in-depth consideration of these dynamics in the context of hybrid simulation in trapped ions.

We evolve the initial states |ψ_0⟩ and |ψ_1⟩ under the with U=0 using a first-order Trotterization e^-iĤ_ℤ(2)t≈ (Û_1 Û_2)^r and timestep Δ t=t/r with

    Û_1    = ∏_i^L e^-i g X̂_i,i+1Δ t,
    Û_2    =  ∏_i=1^L-1 e^-i J Ẑ_i (â^†_i â_i+1 + â_i â^†_i+1)Δ t.

We have discussed in Section <ref> how the time evolution operator for the individual Hamiltonian terms in Eq. (<ref>) can be implemented using native gates in our proposed architecture. 

At the end of the simulation, we measure the state of the qubits and modes. The measurement of the mode occupations are carried out bit-by-bit in a binary representation, first mapping the boson number parity onto the qubit using an SQR gate, making a mid-circuit measurement, resetting the qubit, then measuring the super-parity and so on, as has been demonstrated experimentally <cit.>. This requires a circuit with a depth that is logarithmic in the maximum boson occupation cutoff. If the qubits representing the gauge fields are read out first, they can be reused to iteratively read out the Fock state of the cavity mode to which they are dispersively coupled.


To show that these dynamics can be probed in near-term hardware, we implement the (1+1)D gate sequences in Bosonic Qiskit <cit.>, assuming an implementation in which the ℤ_2 gauge fields are encoded in transmons.

 In order to simulate the dynamics in presence of the noise, we need to specify the gate durations we assume for the potential experimental implementation using the hardware in 
 Sec. <ref>. We assume a value of χ≈ 3 MHz for the dispersive interaction strength. The duration of each Trotter step is then dominated by the conditional beamsplitter operation, which is on the order of 2.1 μs, consisting of two conditional parity gates which take 1 μs each and one beamsplitter which has a duration on the order of 100 ns <cit.>. Single qubit gates take on the order of 10 ns. The Trotter step size in our numerical simulation, Δ t=0.1/J, was chosen to be small enough to recover the expected physical behavior in comparison to exact diagonalisation. We chose 1000 shots per circuit, which leads to both low shot-noise errors as well as a realistic experimental runtime: because each Trotter step has depth 2, a depth one circuit takes <2.2μs and we need 85 Trotter steps to evolve to Jt=8.5. Each circuit therefore takes at most 2× 85× 2.2μs=374μs. We measure 85 circuits and therefore, 1000 shots per circuit leads to a total projected time to take the experimental data of Fig. <ref>b) for g=5J of 1000×85×374μs≈ 32s.




On the one hand, in the regime of infinite hopping strength (i.e., g=0), starting from the Fock state |00100⟩ with open boundary conditions, we find a linear light-cone-like spreading of the boson, c.f. Fig. <ref>. On the other hand, large field strengths prohibit the gauge fields from flipping during hopping events which leads the boson to remain confined to the central site, as predicted and observed in Ref. <cit.> in cold atoms. Strikingly, for large field strengths, when starting from the Fock state with an even initial number of bosons such as |00200⟩, we again find a linear light-cone even in the large field strength regime. We find that the hopping takes place on a much slower time scale. This hopping is enabled by a perturbative energy-conserving pair-hopping process of strength J^2/2g in which a qubit which is flipped by a first boson hopping and then flipped back to its original value by the second boson hopping <cit.>. Despite the transmon decoherence and the much longer time scale (the circuit to evolve to the latest time point for g=5J takes around 2×85×2.2μs≈ 374μs experimental time), these coherent dynamics are observable and match exact diagonalisation. 





 §.§ Qubit-boson ground state preparation with VQE


In this section, we show how the variational quantum eigensolver algorithm (VQE) <cit.> (see Refs. <cit.> for important follow-up works) can be implemented not only on oscillators/modes <cit.>, but also on our proposed hybrid qubit-oscillator architecture to prepare ground states of the ℤ_2-Higgs and U(1) quantum link model. VQE is a variational ansatz for the ground state wavefunction in which layers of unitaries Û_k act on some initial state. The unitaries in our approach, sometimes called the quantum approximate optimisation algorithm, are related to the Hamiltonian 

    Ĥ=∑_k,iĥ_k(i),

where ĥ_k(i) are terms of the Hamiltonian with i the site index on which they are applied. We choose 

    U_k(θ,i)=e^iθĥ_k(i),

where θ is a variational parameter which is optimized. Each gate in the ansatz has an independent angle; we denote the vector of angles with θ.

In the VQE, the classical and quantum processors  work together to reach the ground state by executing an optimization algorithm which updates the vector of ansatz angles θ→θ ', and uses the quantum computer to measure the expected energy of the anstatz characterised by those angles ⟨ψ(θ ') | Ĥ | ψ(θ ')|$⟩. The intuition behind this VQE procedure comes from the fact that for our choice of ansatz, a sufficiently large number of layers and a particular choice ofθ, VQE realizes adiabatic state preparation (note the parallel between choosing the gate variable angles to be proportional to time for dynamics, and choosing them to be the optimisation parameter for ground state search). Therefore we know that in principle, VQE is able to prepare the ground state if the Hamiltonian has a large enough gap such that the necessary evolution time is within the coherence time of the quantum computer. In the absence of first-order phase transitions, the gap is always polynomial in the system size and hence an efficient state preparation path can be found. In the gauge theory setting, constructing a variational ansatz this way also has the advantage of fulfilling Gauss's law. This is exactly the case if every term in the Hamiltonian conserves Gauss's law individually (which will be the case for us below). If this is not true, Gauss's law can be fulfilled approximately by choosing the angles small enough to essentially realize the Trotter limit of adiabatic time evolution. 

The hope behind VQE is that by minimizing the energy with respect to theθangles, a circuit could be found which prepares the ground state in a shallower circuit than would be required for adiabatic state preparation, which in this digital setting would require Trotterisation. In the following, we discuss VQE in the cQED architecture for theℤ_2-Higgs model andU(1)quantum link model in(1+1)D.



  §.§.§ Z2-Higgs model


In the case of the , theÛ_krequired for the ansatz discussed in the first paragraph of (Sec. <ref>), Eq. (<ref>), are

    Û_1(θ,i)    = e^-θẐ_i,i+1(â^†_i â_i+1 - â_i â^†_i+1)
    Û_2(θ,i)    = e^iθẐ_i,i+1(â^†_i â_i+1 + â_i â^†_i+1)
    Û_3(θ,i)    = e^-iθX̂_i,i+1
    Û_4(θ,i)    = e^-iθn̂_i^2.


While technically not necessary (all other gates together are already able to prepare the ground state by approximating an adiabatic evolution), we find that includingÛ_1leads to a lower energy error for the same number of layers than when it is not included. All of the above unitaries can be realized with our native gate set using techniques previously discussed; in particularÛ_1andÛ_2are conditional beamsplitters (see Eq. (<ref>)),Û_3is a single-qubit rotation, andÛ_4can be implemented using a SNAP gate (see Section <ref>, or alternatively App. <ref> for an ancilla-free implementation). 

Using these resource unitaries, we construct anM-layer ansatz withNfield sites andN+1matter sites as

    |ψ_M(θ)⟩ = ∏_l=1^M    [
        ∏_i=1^NÛ_3(γ^i_l, i)∏_i=1^N+1Û_4(δ^i_l, i)
         ∏_i ∈even^N Û_2(α^i_l, i)Û_1(β^i_l, i)
         ∏_i ∈odd^N Û_2(α^i_l, i)Û_1(β^i_l, i)
        ]|ϕ⟩,

where|ϕ⟩is an easily-preparable state,θis a vector of the(4N+1)Mreal-numbered angles in this ansatz, and the products overi∈even (odd)include only the field qubits with even (odd) index. This means that every single gate in each layer and at each site has a different angle to be optimised. Fig. <ref>a) shows one layer of the VQE ansatz.



Since theℤ_2resource unitaries in Eq. (<ref> - <ref>) exactly preserve Gauss's law (in the absence of decoherence), the entire circuit evolution will take place within the Gauss's law subspace chosen by the initial state. Restricting the evolution of VQE circuits within a particular subspace is a powerful technique introduced in prior work <cit.>, with many applications to constrained optimization problems such as maximum independent set <cit.>. Here, we choose the initial state|ϕ⟩to be a simple-to-prepare state which obeys Gauss's law, specifically, a product state of Fock states with unit occupation.

 Because we aim to provide scalable methods adapted for large system sizes, global optimizers such as `DIRECT' which were used in previous studies of lattice gauge theories <cit.> are not well suited here. We compared optimizers (with Python Scipy names `COBYLA', `Nelder-Mead', `CG', `trust-constr', `BFGS', `L-BFGS-B', `SLSQP') and find that, among the optimizers which were tested, `SLSQP' (Sequential Least Squares Programming) has  the fastest convergence rates for theℤ_2model  without qubit decoherence or shot noise, measuring the energy by calculating the expectation value of the Hamiltonian matrix directly with the circuit state vector without shot noise.

We find that the VQE circuit (without noise and therefore no sampling error) can estimate ground state energy with relative error with respect to the true ground state energy of around10^-3for two, three, and four sites using one, four, and seven layers, respectively, as can be seen in Fig. <ref>b. The VQE requires on average around10^2,10^3, and10^4circuit evaluations for two, three, and four sites respectively, an increase which is due to the commensurate increase in the number of Fock states required to represent a larger Hilbert space for unit filling. Increasing the number of layers increases the total number of variational parameters, which can lead to a higher final fidelity; however, this presents a tradeoff – in the presence of noise, reliable results can be extracted only from shallow depth circuits. Thus, the overall fidelity will decrease when a large number of layers is used <cit.>.

Furthermore, we find that the higher the boson filling, the larger the number of layers required, with up to four layers required for four bosons in two sites to reach a relative error below10^-3. This may be related to the increased size of the Hilbert space, i.e., due to the increase in possible combinations for distributing the bosons among a set number of sites.

To measure the expectation value of the Hamiltonian as required for the VQE optimisation procedure, we discuss how to measure each term of the Hamiltonian individually <cit.>. The field term⟨- g ∑_i=1^L-1 X̂_i,i+1|$⟩ is a simple sum of X-basis qubit measurements. The onsite interaction term on a single site i can be measured by reading out the photon number of the cavity using the binary search readout <cit.>. This enables the single-shot readout of n̂_i, which can then be used to reconstruct ⟨ n_i^2 ⟩. For the gauge-invariant hopping term - J ∑_i=1^L-1⟨Ẑ_i,i+1( â^†_iâ_i+1 + h.c.)|$⟩, we are faced with the difficulty that we can not directly read out the beamsplitter term by measuring the cavities in the number basis. 
Instead, we note the following identity:

    ⟨ψ|â_i^†â_i+1  + h.c.|ψ⟩
       =⟨ψ̃|BS_i,i+1(π/2,π/4)(â_i^†â_i+1  + h.c.)BS_i,i+1^†(π/2,π/4)|ψ̃⟩
       =⟨ψ̃|( â_i^†â_i - â_i+1^†â_i+1)|ψ̃⟩
       = ⟨ψ̃|n̂_i|ψ̃⟩ - ⟨ψ̃|n̂_i+1|ψ̃⟩,

where we have inserted1 = BS_i,i+1(π/2,-π/4) BS_i,i+1(π/2,π/4)on either side of the hopping, and have cast the expression in terms of the transformed state|ψ̃⟩=BS_i,i+1(π/2,π/4)|ψ⟩. Thus, much like a Hadamard gate facilitates anX-basis measurement for qubits, we can measure the expectation value of a bosonic hopping term by `rotating' the modes into the Fock state basis via a 50:50 beamsplitter, followed by a Fock-basis measurement. By correlating the result with that of the gauge field linkZ-basis measurement result, one can reconstruct the expectation value of the gauge-invariant hopping.

Only commuting terms can be measured simultaneously; for theℤ_2model, we therefore perform three separate measurement circuits to estimate the energy. In the first circuit in Fig. <ref>c, the onsite interaction is measured on all sites simultaneously.  In the same circuit, it is possible to measure the field terms, though this step must precede measurement of the onsite interaction such that the gauge link qubits can then be reset and used for binary readout. Next, we measure half of the gauge-invariant hopping terms, preceeded by a measurement of the electric field (X̂_i) at every other gauge link. Finally, the other half of the gauge-invariant hopping terms are measured along with the other half of the field sites.

In summary, in this section we showed that a VQE ansatz for the  which is specifically tailored to the types of operations available in hybrid qubit-oscillator hardware indeed only requires very few layers to determine the ground state for the small system sizes we can access numerically for proof of principle demonstrations. The implementation of this ansatz directly follows from the methods presented in the first part of this paper (Secs. <ref> - <ref>), which implies that the hybrid qubit-oscillator hardware will require shallower circuits than qubit hardware as discussed in Sec. <ref>. However, more work needs to be done beyond these proof of principle results to determine the circuit depths and the number of shots which will be required as a function of system size.



  §.§.§ 


In the following, we describe how to carry out a hybrid oscillator-qubit VQE of the  in Eq. (<ref>) in(1+1)D. Furthermore, we demonstrate how to detect the first-order phase transition atτ=πin our proposed architecture. While the correspondence between the model in Eq. (<ref>) and quantum electrodynamics in(1+1)D is only applicable for an even number of lattice sites due to the staggering of the fermions, we also study an odd number of sites here for benchmarking purposes.




Following the same scheme as for the , we choose resource unitaries that are generated by the individual terms of the Hamiltonian in Eq. (<ref>). Furthermore, we leverage the Schwinger-boson mapping of the gauge fields as described in Sec. <ref>, and map the fermions to qubits using the Jordan-Wigner encoding (which is particularly simple in(1+1)D). With these choices, we construct our VQE ansatz using the following resource unitaries:

    Û_S1(θ,i)    = e^i θẐ_i
    Û_S2(θ,i)    = e^-i θ (n̂_i,i+1^a-n̂^b_i,i+1)
    Û_S3(θ,i)    = e^-i θ((n_i,i+1^a)^2+(n_i,i+1^b)^2)
    Û_S4(θ,i)    = e^-iθ(σ̂^+_i+1â^†_i,i+1b̂_i,i+1σ̂^-_i + h.c.)
    Û_S5(θ,i)    = e^-θ(σ̂^+_i+1â^†_i,i+1b̂_i,i+1σ̂^-_i - h.c.).

Each of these unitaries can be implemented using techniques described in prior sections. In particular,Û_S1andÛ_S2are native gates – see Table <ref>. As explained in Sec. <ref>,Û_S3can be synthesized using a SNAP gate. Finally, the gauge-invariant hopping unitariesÛ_S4andÛ_S5can be implemented with the strategies described in Section <ref>.

Our ansatz for the wavefunction is

    |ψ_M(θ)⟩ = ∏_l=1^M    [ ∏_i ∈even^N Û_S5(α_l, i)Û_S4(β_l, i) 
          ∏_i ∈odd^N Û_S5(α_l, i)Û_S4(β_l, i)
         ∏_i=1^NÛ_S3(γ_l, i)  ∏_i=1^NÛ_S2(γ_l, i)
       ∏_i=1^N+1Û_S1(δ_l, i)] |ϕ⟩,

forN+1matter sites. Note that contrary to the ansatz for the , each site has the same parameter such that there are only4Mtotal parameters.

The expectation value of the Hamiltonian, can be measured using correlated measurements. To evaluate the gauge-invariant hopping term, we first useσ̂^+=1/2(X̂+iŶ), to split it into

    σ̂^+_i+1â^†_i,i+1b̂_i,i+1σ̂^-_i   + h.c.
    =
        1/4(   X̂_i+1(â^†_i,i+1b̂_i,i+1 +h.c.) X̂_i
    +   Ŷ_i+1(â^†_i,i+1b̂_i,i+1+h.c.) Ŷ_i
    
        +   X̂_i+1(-iâ^†_i,i+1b̂_i,i+1 +h.c.) Ŷ_i
    +   Ŷ_i+1(iâ^†_i,i+1b̂_i,i+1+h.c.) X̂_i)

Using our method in Eq. (<ref>) we can map a measurement of the bosonic terms in above equation to a measurement of the number densities of the two modes. Correlating those number densities with the appropriate qubit measurements, we can therefore measure all four terms in the above expansion. Because we require two circuits to measure all possible lattice sitesi, we  require eight separate measurement circuits for the gauge-invariant hopping. All the other terms are straightforwardly accessible by measuring the values of the boson numbers in the modes <cit.> and the state of the qubits.

Using Bosonic Qiskit <cit.>, we implement this ansatz and test the convergence of the VQE with respect to the number of layers forS=1and show the results in Fig. <ref>. For a single gauge field site and two matter sites, we initialize the system in the “vacuum state”|↓0 ↑⟩, where the middle label indicates the eigenstate ofS^z, and assume perfect energy measurements, i.e. no shot noise and decoherence. We find that one layer is sufficient for obtaining a relative energy error of10^-5. This is related to the fact that the ground state for this system is given by1/√(2)(|↑1 ↓⟩ + |↓0 ↑⟩). Three sites, initialised in|↓0 ↑0 ↓⟩, we find from our numerical evaluations that  two layers are required to attain low relative errors of10^-8, with only a small increase in total number of iterations from the two-site case.

In the regime of largeg, the  features a first order phase transition <cit.> at half-integer values ofτ/2π. This can be seen trivially by takingg→∞, in which case the electric term dominates the Hamiltonian, pinningŜ^zat the integer closest toτ/2π. This phase persists for finiteguntil the mass reachesM/g≈1/3(see definitions ofMandgin Hamiltonian in Eq. (<ref>)). The energy gap between the ground and first excited state is generally exponentially small in system size at first order phase transitions. Therefore, we expect that this transition will be relatively sharp and visible even for small systems of two or three sites. Hence, this is a good proof-of-principle target for small-scale experiments. 

For various values ofτused to find ground states with the VQE ansatz discussed above (Eq. (<ref>)-<ref>), we measure the expectation value of the electric field (in the Schwinger-boson representation)⟨S^z_i,i+1|=⟩ 1/2⟨n̂^a_i,i+1 - n̂^b_i,i+1|$⟩. In Fig. <ref> we indeed find sharp drops corresponding to the first order phase transition at τ = ±π. The Bosonic Qiskit VQE results match our exact diagonalization results to relative errors of around 10^-7.

Our results show that our VQE approach using hybrid oscillator-qubit hardware can prepare ground states of the in (1+1)D for two and three site systems. Although in (1+1)D the can be probed with qubit-only hardware by integrating out the gauge fields, this would not generalize to higher dimensions. In that case, each gauge field link would be encoded in many qubits and Hamiltonian operations would be highly costly as discussed in Section  <ref>. Another possible advantage of our approach is the usage of the Schwinger-boson constraint n̂_a+n̂_b=2S as a (single) photon loss error detection scheme. This can be implemented by measuring the joint-parity of modes a and b; to what extent this improves results from noisy hardware is an interesting direction for future exploration.





 §.§ VQE with hardware and shot noise


Next, we test the resilience of our methods to decoherence and shot noise in the energy measurements. To do so, we simulate the full measurement schemes that would be employed in the experiment and also include qubit decay and dephasing, which are the leading decoherence mechanism <cit.> for the small system sizes and unit filling we consider here (i.e., we neglect photon loss). For systems where large photon numbers are possible, mode decay will become appreciable. The model studied in this section is the , for which we defined the ansatz in Sec. <ref>. We choose U/J=g/J=1, which we expect to be among the most difficult regimes to simulate due to the balanced competition between all energy terms. We consider three matter sites and two gauge links at unit filling, requiring three modes and two qubits. 

We start by explaining how to measure the expectation value of the Hamiltonian in experiment. Following this,  we first benchmark the simultaneous perturbation stochastic approximation (SPSA) optimiser <cit.> optimizer (as opposed to SLSQP used so far) without shot-noise, as this is the optimizer of choice for optimisations with shot noise. We investigate how to use averaging over previous iterations to improve convergence of SPSA (still without any noise). Third, we run our improved optimisation procedure with shot noise. Upon finding that there is a certain independence to the number of shots in the initial stages of optimisation where the gradients are large, we further adapt our procedure to incrementally increase the number of shots during the optimisation. Fourth, we include hardware noise (qubit T_1 decay and additional pure T_ϕ dephasing) into the numerical simulation. Finally, we error-mitigate our results by implementing a procedure to post-select on preservation of Gauss's law.



  §.§.§ Introducing shot noise in the energy measurement


We simulate a binary search routine <cit.> for measuring the cavity in the Fock basis. This naively requires an additional readout qubit per cavity for the measurement (i.e. separate from the transmon qubits that encode the gauge field links). However, for the we can remove this requirement by first measuring the gauge field links, and then reusing the qubits as ancillas for mode readout, c.f. Figs <ref>c. To do so, we measure and reset the field qubits prior to the cavity measurement. In our VQE simulations, we simulate the full binary search, which enables us to simulate the effective shot noise on the mode measurement that results from the qubit measurements. Due to shot noise, each experiment requires 𝒪(1/ϵ^2) circuit executions to estimate ⟨Ĥ_ℤ(2)|$⟩ to additive precisionϵ.





  §.§.§ Using averaging in the stochastic VQE optimisation


In the presence of shot noise error, VQE optimisation is far more challenging <cit.>. We use the SPSA, which is particularly well suited for VQE because it only evaluates the Hamiltonian twice per iteration regardless of the dimension of the optimization problem.

In Fig. <ref> we show the energy found in each iteration using the SPSA optimiser. The gray dots indicate an optimisation without any shot noise, and we see that it does not converge within7000iterations.

It has been found that averaging variable parameters can increase the performance of stochastic optimisation <cit.>. We found in our case that the convergence of the SPSA optimiser, which is stochastic, can be improved by periodically choosing parameters for the next iteration which are averages of previous iterations. Therefore, we adapt our implementation of the SPSA optimiser in order to perform a parameter averaging step throughout the optimization: everyniterations, we use the average of the variable parameters of the lastmiterations. It is important to not choosenso small that the averaging is too frequent, resulting in a small sample size, or choosemtoo large such thatm≳n, which may constrain the optimisation too much and lead to trapping in local minima. In Fig. <ref>, we show data corresponding ton=15andm=7in dark blue. We find notable improvements using this strategy over the one without averaging, i.e.n=1andm=1.



  §.§.§ Adaptive shot budget in the VQE optimisation


We now include shot noise in the estimation of the energy. The number of shots required for good convergence of the VQE is determined by the size of the energy gradient. Generally, gradients are larger in the first iterations and gradually decrease as the minimum is approached. This is evident in the large fluctuation between iterations for the early iterations<2000in Fig. <ref> (black dots). Therefore, it is natural to increase the number of shots used in the energy measurement of later iterations <cit.>. Specifically, we increase the shots per iteration by an order of magnitude every1500iterations (starting at10shots per iterations). We show the resulting energy trajectory as light blue points in Fig. <ref>. We find that (for this specific system size) this method yields similar performance to the scheme without shot noise (dark blue). In general, we expect the total number of shots necessary to prepare the ground-state with fixed fidelity to increase exponentially with the system size, see e.g. Ref. <cit.> for empirical evidence.



  §.§.§ VQE results including qubit decoherence, averaging and adaptive shot budget


Finally, we add the dominant source of experimental noise in the hardware – qubit decay and dephasing – to our simulations in addition to the shot-based measurement of the energy discussed above. This requires specifying the duration of each of the gates in the VQE ansatz. The ansatz is explained in Sec. <ref>– see Fig. <ref>a for a circuit diagram. Below Eqs. (<ref>-<ref>), we explain how every resource unitary can be realized using native gates. Altogether, the depth of one ansatz layer consists of four conditional beamsplitters and aSNAPgate. We assume a value ofχ=πMHz for the dispersive coupling. In this case the maximal duration of a conditional beamsplitter is≈2.25μs, consisting of two conditional parities which takeT=π/χ=1μs each, and one beamsplitter of250ns duration (100ns has been shown to be feasible  <cit.>). We do not include duration cancellations, a range of which are possible. For example, two successive conditional beamsplitters acting on the same oscillator modes and qubits can cancel each others' conditional parity operations, leading to a total time of2.25μs. We also note that in principle, the duration depends on the gate parameters, which in experiment could lead to a shorter duration than estimated here. TheSNAPgate takesT=2π/χ=1μs <cit.>. The total duration of one VQE layer is roughly(4×2.25 + 1)=10μs. Four layers (T=40μs) are therefore well within the decay time of the qubit (T_1=200μs) and the mode (T_1 ≈1ms). Therefore, the probability for the|1⟩state of a mode to decay is 4% and for a qubit to decay during the four-layer ansatz cicuit is 20%.  While in most regions of the phase diagram, the occupation of the higher-lying mode states is small, and therefore the effect of mode decay (and other non-idealities such as self-Kerr <cit.>) minimal, this is not the case in the superfluid and clump phase whenU/Jis small. In these cases, mode decay could contribute significantly. We leave the study of this effect for future work and neglect photon loss throughout.

We show the optimisation including averaging, adaptive shot budget, and qubit decay in black dots in Fig. <ref>. Comparing this result to the optimisation without qubit decoherence in light blue, we find that qubit decoherence leads to a significant energy error of about0.86J. However, we still find that the energy is well below the first excited state. Therefore, we still expect a reasonable overlap with the ground state. In order to test this expectation, we plot the fidelity with the ground state in the inset of Fig. <ref>, using matching colors to the main figure. We see that for the simulations excluding decay and dephasing errors (light and dark blue), the fidelity remains  high, following the same trend as the energy.

We also note that total runtime of the optimisation in the presence of noise is small: we find an energy error below the energy of the first excited state after around3000iterations with a total shot count of<2×10^5, yielding a total quantum runtime of<10s, if this experiment were to be run on current superconducting hardware that is limited only by gate times. 



  §.§.§ Post-selection on Gauss's law

The conservation laws obeyed by lattice gauge theories, such as Gauss's law and particle number conservation, can be used to discard unphysical shots in the presence of hardware noise, enabling partial error detection.

In order to use Gauss's law for error-mitigation, we must measure the operatorĜ_i(see Eq. (<ref>)) for each shot of the experiment. For theℤ_2-Higgs model in(1+1)D Gauss's law reads

    Ĝ_i = X̂_i-1,i e^iπn̂_iX̂_i,i+1 = 1,

for alli. Because Gauss's law commutes with all terms in the Hamiltonian, it is possible to measure Gauss's law for allias well as the energy of the state in the same shot. This requires mid-circuit measurements and is most often done using the Hadamard test <cit.>. However, the Hadamard test involves ancillary qubits which we do not have at our disposal. Therefore, we will present a method below that avoids ancillary qubits. This method requires additional gates and measurements, compared to the measurement circuits presented in Fig. <ref>c). We also discuss that even for the measurement circuits as described in Fig. <ref> (i.e. without additional gates and measurements), Gauss's law can be partially checked. 

The first scheme is inspired by the schemes presented in Ref. <cit.> for measuring Toric-Code stabilizers. The key observation is that the individual parts of the Hamiltonian commute with Gauss's law and therefore we can check Gauss's law for each measurement individually. For the measurement of the gauge-invariant hopping, replacing the beamsplitters in Eq. (<ref>) with conditional beamsplitters, a measurement of the number operator of the modes directly yields a measurement of the gauge-invariant hopping. More explicitly,

    CBS_i,i+1(π/2,-π/4)(n̂_i-n̂_i+1)CBS_i,i+1(π/2,π/4) 
       =Ẑ_i,i+1(â_i^†â_i+1  + h.c.).

Hence, conjugating a number measurement of modesiandi+1with conditional beamsplitters yields a measurement of the gauge-invariant hopping without having to measure the qubit, as shown in the green part of the circuit shown in Fig. <ref>a). Crucially, one can check that also the operators effectively measured byn̂_iandn̂_i+1individually commute with Gauss's law. BecauseX̂_i-1,icomutes with both the gauge-invariant hopping and Gauss's law, it can be measured simultaneously. Gauss's law is then measured by the circuit shown in black and white in Fig. <ref>a. It relies on the fact that conjugation of a qubitẐ_i,i+1operator with SQR gates yieldsẐ_i,i+1 e^iπn̂_i, c.f. Eq. (<ref>). Similarly, conjugation with Hadamard gates and CNOTs convertsẐ_i,i+1intoX̂_i-1,i X̂_i,i+1. Taken together, this converts a measurement ofẐ_i,i+1to a measurement of Gauss's law operatorĜ_i. Hence, these two measurements allow for a full check of Gauss's law while measuring the value of the gauge-invariant hopping in the same shot. This scheme requries additional gates. 

In the second scheme, we do not require additional gates and instead  propose a form of `partial' error-mitigation based on the information we already have access to in the context of the VQE: the measurements of the Hamiltonian terms shown in Fig. <ref>c. This does not enable a local check of Gauss's law for every matter site (which is why we call it `partial'), and hence only slightly improves the results. These three measurement circuits do not commute, and therefore only one can be used per shot of the experiment. The measurement circuit of the number density (shown in Fig. <ref>c) carries out a photon number measurement on each mode, and could thus be combined with the measurement ofX̂on each qubit to check Gauss's law at each site (see Fig. <ref>a). However, the other measurement circuits (see Fig. <ref>c) are not compatible with a measurement of Gauss's law at each site, as some qubits are individually measured alongẐ(which alone does not commute withĜ_ifor an overlapping region) and the mode occupations are not resolved for each site. However, in these circuits, we measureX̂for every other gauge link, and can furthermore resolve total photon number in the intermediary pairs of modes. Thus, we can conduct checks of Gauss's law for blocked pairs of sites:

    Ĝ_iĜ_i+1 = X̂_i-1,i e^iπ (n̂_i + n̂_i+1)X̂_i+1,i+2.

This possibility is illustrated in Fig. <ref>b. Thus, we can post-select on valid Gauss'-law preserving states for all three measurement circuits needed for estimating the energy. This scheme captures only single Gauss's law violations, but not errors on two consecutive sites.





We simulate full Gauss's law post-selection (the first scheme) and focus on its influence on the energy measurement. To do so, we take the final VQE parameters found in a noisy simulation (i.e. the final iteration of the black dots in Fig. <ref>) and run a simulation including hardware noise in the VQE circuit, but neglect noise in the measurement circuit. We then check whether an error happened in every quantum trajectory of the noisy simulation. The error-free trajectories exhibit high overlap≈0.96%with the ground state (c.f. blue line in Fig. <ref>c). For trajectories with an error, we project onto the correct Gauss's law sector. The norm of the part of the wavefunction we project out gives us an estimate for the number of shots that would need to be discarded. We find that13%of shots need to be discarded. We show in red the overlap with the eigenstates in the correct Gauss's law sector for some such states in  Fig. <ref>c. Those states have still large overlap with the ground state, but also some overlap with excited states. To calculate the average post-selected energy, we renormalize the projected states and average over all trajectories. We find a post-selected energy of0.36J, compared to a result of0.88Jwithout post-selection and an exact ground-state energy ofE/J=0.02. This corresponds to approx.20%of the difference between ground and first excited state energy.

In our simulations, we find that84%of quantum trajectories have no error. This is roughly in agreement with the dephasing survival probability of a single qubit, whereT_2=200μs, only one of the qubits starts in|+⟩, the other in|-⟩, and the total circuit duration isT=40μs, which isexp(-T/T_2)≈82%. These trajectories with no error have energy0.16J. The1-84%=16%of trajectories in which an error happened therefore have average energy4.66J. Gauss-law post-selection reduces this average energy of the noisy trajectories to1.41J, where the remaining energy error results from noise processes within the same Gauss's law sector. We therefore find that removing errors that bring us out of the correct Gauss's law sector drastically reduces the energy error introduced by the noise.

More intricate mitigation methods can be engineered from this Gauss's law check. For example, we note that one can cross-checkĜ_iat different sites to not only detect that an error occured, but also resolve its position. This opens up the possibility to not only post-select on the case of no errors, but additional develop novel error-mitigation techniques where the error can be partially corrected. In addition, if the error is spatially disconnected in such a way that one can ascertain it did spread to a region of interest, then it may be desirable to retain this shot for the computation of observables local to this region. In addition to post-selection using gauge symmetries <cit.>, VQE implementation for the  proposed here is amenable to other commonly used noise mitigation techniques such as zero-noise extrapolation <cit.> and probabilistic error cancellation <cit.>. An alternative to post-selection is to introduce an additional term in the Hamiltonian which enforces the symmetry constraint by prethermalization <cit.>, or pseudo-generators for gauge protection <cit.>. 



 §.§ Measuring observables


Having characterised the hybrid oscillator-qubit VQE state preparation in terms of the energy error and overlap with the true ground state, we now show how to extract relevant observables in the proposed cQED architecture. While we benchmark the measurement of the observables for the ground state, the measurement schemes can in principle be applied to any state. 



  §.§.§ Detecting phases of the Z2-Higgs model in VQE-prepared ground states

In Fig. <ref>, we calculate different gauge-invariant observables in the  for varying onsite interaction, with respect to VQE ground states of a system with four bosonic sites and three gauge field sites found using noiseless simulations employing the SLSQP optimizer. We compare to exact diagonalisation results. We find that all VQE observables agree very well with the exact results, as expected from the low relative energy error found in Fig. <ref> for this small system size.

Importantly, despite the small system size, we can see the hallmarks of three qualitatively different states asU/Jis scanned <cit.>. The three states are indicated by the blue (rightmost), red (middle) and yellow (leftmost) color of the lines, representing a Mott-insulating phase, pair superfluid phase and the “clump phase”– a phase in which the gauge field perturbatively generates a nearest-neighbour attraction between the bosons, leading to the bosons “clumping” together on the same site <cit.>. For four sites, we define the “clump order parameter” as the probability of all four bosons to occupy the same site (defined in caption of Fig. <ref>). This observable indeed becomes large for smallU/J<cit.>. A (Mott) insulating state is characterized by a lack of spreading of particles in the system. We observe signatures of this effect in the string order correlator and pair hopping, which are both small for largeU/J. Moreover, on-site fluctuations are also small in the insulating state as expected. Finally, the parity tends towards-1, indicating that on average there is an odd occupation at each site. By contrast, in a superfluid phase, particles spread far. In the case of the , single particles stay confined due to the presence of the gauge field while two particles can spread <cit.>, which we observe in the large pair hopping and fluctuations. Moreover, contributions from states with a single boson are suppressed, which we see from the parity tending towards+1.





Having shown that signatures of different phases can be found from VQE simulations, we now discuss how the different observables can be measured in our cQED architecture.



  §.§.§ Measuring the string order correlator
 

The gauge-invariant string order correlator is large in a deconfined phase and small in a confined phase. It is defined as⟨ â_i^† Ẑ_i,i+1...Ẑ_j-1,jâ_j + h.c.|$⟩. To measure the value of this correlator on the ground states prepared via VQE, we need to correlate a measurement of the gauge field qubits with the measurement of the beamsplitter operator of two bosonic sites. We note that one possibility is to adapt the scheme used in Sec. <ref> to measure a long-range hopping through the use of bosonic SWAPs. Ordinarily, such a measurement would require a time linear in the distance |j-i|. Alternatively, here we propose a separate strategy to measure this nonlocal observable using homodyne measurements: first, we expand â_i = (x̂_i + ip̂_i)/2, where x̂ and p̂ are the bosonic position and momentum operators respectively, such that 

    ⟨â_i^†Ẑ_i,i+1...Ẑ_j-1,jâ_j + h.c.|⟩
       =1/2(⟨x̂_i Ẑ_i,i+1...Ẑ_j-1,jx̂_j|+⟩⟨p̂_i Ẑ_i,i+1...Ẑ_j-1,jp̂_j|⟩)

We must therefore perform correlated measurements of x̂_i, Ẑ_i,i+1 and x̂_j (and, separately, p̂_i, Ẑ_i,i+1 and p̂_j). In order to measure x̂_i, we use a Hadamard test implemented with a displacement conditioned on an ancilla, defined as CD(α) = e^( αâ^† - α^*â)Ẑ^anc (see  Table <ref>). Choosing an arbitary real angle θ such that α=iθ, we get CD(iθ)=e^iθx̂Ẑ^anc. A measurement of the ancilla after applying the Hadamard test circuit with a small angle θ, shown inside Fig. <ref>, then yields 

    ⟨Ẑ^anc|_⟩α=iθ   =⟨( sin(2θx̂)|
    ⟩   ≈ 2θ⟨x̂|+⟩ O(θ^3).
 
To measure the first part of the string order correlator we perform two Hadamard tests simultaneously in order to correlate sites i and j, i.e., after applying the Hadamard circuits on both sites i and j we measure

    ⟨Ẑ^anc_i Ẑ_i,i+1…Ẑ_j-1,jẐ^anc_j|_⟩α_i=α_j=iθ
       ≈ 4 θ^2⟨x̂_i Ẑ_i,i+1…Ẑ_j-1,jx̂_j|.⟩
 Similarly, we measure the second part of the string order correlator, i.e., those involving p̂ correlations, by using conditional displacements with α=θ. 

In summary, we get 

    ⟨â_i^†Ẑ_i,i+1...Ẑ_j-1,jâ_j + h.c.|⟩
       ≈1/8θ^2(⟨Ẑ^anc_i Ẑ_i,i+1…Ẑ_j-1,jẐ^anc_j|_⟩α_i=α_j=iθ
         +⟨Ẑ^anc_i Ẑ_i,i+1…Ẑ_j-1,jẐ^anc_j|_⟩α_i=α_j=θ).


To test that we can be in the small θ regime while still being able to measure the string order correlator accurately with a reasonable number of shots, we simulate the measurement using Bosonic Qiskit, c.f. Fig. <ref>. We find that a choice of θ≈ 0.1 suffices, requiring approximately 1/θ^4=10^4 shots, as expected from the statistical uncertainty of shot noise. Note that for this measurement, we do not need additional ancillae as one can recycle gauge link qubits that are not within the support of the string order correlator for this purpose. Lastly, we remark that one could alternatively use the method of quantum signal processing (see App. <ref>) to extract ⟨x̂⟩ similarly to Eq. (<ref>), but with a larger value of θ. This could be achieved by enacting a target function 2/πarcsin(·) on the signal sin(2θx̂); if 2θx̂≤ 1/2, this outputs θ/π⟨x̂⟩ + O(ϵ) while requiring O(log(1/ϵ)) queries to CD(iθ) <cit.>. 





  §.§.§ Measuring the superfluid density




The superfluid density ρ_s in 1D is given by <cit.>  

    ρ_s = L/J∂^2 E_TPB(ϕ)/∂ϕ^2,

where E_TPB(ϕ) is the ground state energy in the presence of “twisted” periodic boundary conditions and J is the hopping matrix element between neighbouring sites. The phase ϕ is given by the phase difference between two edges, for example for the , we have

    Ĥ^TPB_ℤ(2)(ϕ) = Ĥ_ℤ(2) - J (e^iϕâ_L^†Ẑ_n â_1 + e^-iϕâ_L Ẑ_n â_1^†),

where the entirety of the accumulated phase ϕ is incorporated by the hopping term between the two edge sites, where Ĥ_ℤ(2) is defined in (<ref>) and L is number of sites. To measure ⟨Ĥ^TPB_ℤ(2)|$⟩, we can follow the general procedure outlined in section <ref>. The procedure to estimate the stiffness is as follows: we first prepare the ground state for a range of small values ofϕ. We do so by extending the variational ansatz in Section <ref> to incorporate edge-to-edge hopping. Once prepared, we estimate Eq. (<ref>) via finite-difference methods or, alternatively, through more sophisticated, noise-resilient techniques such as parameter shift rules <cit.>.

Experimentally, periodic boundary conditions are not required in the hardware to implement edge-to-edge hopping. It is sufficient to synthesize a conditional beamsplitter between the two end modes and leverage this for an edge-to-edge hopping term in the Hamiltonian. We can synthesize a beamsplitter between the edges of the chain using the bosonicSWAPgate defined in Tab. <ref> to move the mode information from both edges of the chain to adjacent sites, carry out the required operation between the modes, andSWAPthe information back.




 §.§ Qubit-boson ground state preparation with quantum signal processing


While above we employed the heuristic VQE approach to approximately prepare the ground state, new developments in an algorithm known as quantum signal processing (QSP) provide an alternative approach to ground state preparation with provable accuracy and guarantees. Here we provide a brief comparison of VQE and QSP for ground state preparation of theℤ_2-Higgs model.

Initially pioneered in Refs. <cit.>, QSP provides a systematic method to apply a nearly arbitrary polynomial transformation to a linear operator embedded in a unitary matrix <cit.>, thus furnishing a unifying framework for quantum algorithms <cit.>. In short, QSP works by interleaving the unitary matrix with a sequence of parameterizable SU(2) rotations; see Fig. <ref> for an illustration. The output of this sequence is an embedding of a polynomial transformation of the initial linear operator, parameterized by the chosen SU(2) rotations. For a more detailed introduction to QSP, see App. <ref>. 


In App. <ref>, we illustrate how QSP furnishes an algorithm for ground state preparation, as initially presented in Ref. <cit.> for implementation on early fault-tolerant quantum computers. In this scenario, one considers a HamiltonianHwith spectral gap bounded below byΔ, and seeks to prepare its ground state with fidelity≥1 - ϵ. As input to this algorithm, one assumes access to an easily preparable initial state that has overlap≥γwith the ground state. One also assumes access to the time evolution operatorse^±iH, implemented for instance with Trotterization. To avoid alisiaing, it is taken that the eigenvalues ofHlie in the range(0, π); othwerwise can implement a rescaled time evolution operatore^±i (HΔt + ϕ)to meet this constaraint. 

With this setup, the algorithm uses QSP to approximate a projector onto the ground state, which is then applied to the initial state to project out its component parallel to the ground state. This algorithm ultimately requires a circuit depthO(1/Δ log(1/ϵγ) ), consisting of this many coherent queries to the controlled time evolutione^±iH. This circuit is repeatedO(1/γ^2)times to prepare the ground state with high probability. The circuit of this algorithm is illustrated in Fig. <ref>.




Recent work has shown that QSP is naturally implementable on hybrid CV-DV quantum procesors <cit.>, rendering it applicable to the LGT simulations studied here. To exemplify the requirements of ground state preparation with QSP, let us consider theℤ_2model withn=3sites and couplingsJ=g=U=1, as investigated in Sec. <ref>. Let us also use the same initial state as above, which numerically has overlapγ≈0.36with the ground state. In this scenario, VQE produces the ground state with infidelityϵ≈2·10^-2(see Fig. <ref>). Aiming to achieve a similar performance with QSP, our calculations in App. <ref> indicate that the QSP approach requires∼50queries to the controlled time evolutione^±iH. Using even a first order Trotterization as outlined in Sec. <ref>, a single such Trotter step would require∼3controlled qubit-boson gates per site, for a rough total of≳10^4qubit-boson gates, which is prohibitively more expensive than VQE.

This prohibitive circuit depth suggests that QSP-based ground state preparation will only become practical and competitive with VQE once qubit-boson gate infidelities≪10^-4can be achieved, in contrast to the current qubit-boson gate infidelities∼10^-3<cit.>. And moreover, even if such a gate fidelity can be achieved, the number of repetitions of the QSP circuit will generally scale exponentially in the system size without prior knowledge of the ground state (see App. <ref>). However, this scaling is to be expected because generic ground state preparation is QMA-hard <cit.>, and of course the VQE ground state preparation is also expected to require resources exponential in the system size, without further prior knowledge.


In Fig. <ref>, we compare the performance of VQE to QSP, as applied to the above scenario of theℤ_2model. For both methods, we plot the dominant source of error vs. time to solution, assuming no gate errors in the underlying circuit. The dominant source of error for VQE is the infidelity with the ground state, and for QSP it is the probability of failure (see App. <ref> for details); the time to solution is taken to be the circuit depth times the number of shots. For VQE, increasing the depth (i.e., number of layers) can increase the fidelity with the ground state, but also leads to a high dimensional optimisation problem that suffers from barren plateaus and requires a high shot count to facilitate optimisation <cit.>. On the other hand, QSP requires a relatively deeper circuit to estimate a projector onto the ground state, and a modest number of shots to successfully project out the ground state. Incidentally, the empirical results of Fig. <ref> suggest that both dominant sources of error scale similarly as a function of the time to solution, implying that their total resource costs are comparable. This arises because our VQE circuit has a shallow depth yet requires many shots, whereas QSP has a deeper circuit yet requires fewer shots. We speculate that this seemingly universal behavior arises from the fact that the complexity of ground state preparation is fundamentally governed by the size of the spectral gap of the Hamiltonian, agnostic to the algorithm used to estimate the ground state. In any event, this analysis implies that VQE is better suited to settings with low gate fidelities and short coherence-times, whereas QSP is more favorable if gates are fast yet circuit repetitions are expensive. 






§ CONCLUSIONS AND OUTLOOK


To conclude, in this paper we have developed a resource-efficient framework for implementing fermion-boson quantum simulation of a wide variety of models including lattice gauge theories in(1+1)D and(2+1)D using oscillator-qubit quantum computation. To do so, we developed a range of compilation tools and subroutines to realize multi-body interaction terms using the natively available oscillator-qubit gate set in circuit QED hardware. In particular, we developed parity-controlled, and density-density bosonic gates and oscillator mediated multi-qubit gates as foundational subroutines, which are useful beyond our simulation motivation. In the context of LGTs, we showed how to implement gauge-invariant hopping terms for transmon and dual-rail qubit fermionic encodings. We also presented a method implementing the magnetic field term forU(1)gauge fields.

By performing an end-to-end gate complexity analysis, including gate synthesis errors, we showed that our approach yields an asymptotic improvement from𝒪(log(S)^2)to𝒪(1)scaling in terms of the boson cutoffSfor the gauge-invariant hopping term and an𝒪(log(S))to𝒪(1)improvement for the magnetic field term. Moreover, in order to also compare the prefactors of the compilation, we developed an efficient all-qubit algorithm in the Fock-binary encoding for simulating bosonic hopping terms. We show that even this efficient approach yields gate counts that are10^4times higher than our oscillator-qubit approach. 

We then applied our compilation strategies to some of the most pertinent quantum simulation tasks – ground state preparation and time evolution, for which we developed a  oscillator-qubit variational quantum eigensolver method. We benchmarked our hybrid oscillator-qubit approach through numerical simulations of a  and  in(1+1)D in Bosonic Qiskit, and demonstrated that low-depth oscillator-qubit circuits can prepare states with a high overlap with the ground state for small systems. In particular, in the  we found relative energy errors of<10^-3for just a two-layer ansatz for three matter sites. Separately, we simulated both the time dynamics and VQE for the  in the presence of the leading source of noise in circuit QED architectures 
 – transmon decay and dephasing. We then developed an approach for post-selecting on Gauss's law in VQE circuits, using the fact that each term in the Hamiltonian individually conserves Gauss's law, and showed that energy errors of20%of the difference between ground and first excited state can be reached with this method. We also showed how non-trivial observables can be measured within our approach, including the superfluid stiffness and string order correlators. Finally, we discussed how quantum signal processing – an alternative approach to eigenstate preparation with performance guarantees – can be implemented in oscillator-qubit systems, showing that its time to solution for a fixed error is comparable to VQE. 




 §.§ Prospect for quantum advantage


For simulations in the NISQ-era, it is important that errors do not completely destroy the signal to be observed. In cavities, the probability of a Fock state|N⟩to decay to|N-1⟩scales linearly withN, therefore posing a challenge to the achievable fidelity for largeN, in particular in the context of quantum computation (see e.g. section IIE2 in Ref. <cit.>). By contrast, when mapping bosons to qubits using the Fock-binary encoding, the probability to be in an encoded Fock state|N⟩has a decay-rate which scales asΩ(1)and𝒪(log_2(N)). However, this apparent advantage of the Fock-binary encoding does not directly materialize when considering observables. As we discuss in App. <ref>, for both cavities and qubits, even non-linear observables such as⟨n̂^2|$⟩ have a decay rate which is 𝒪(1) in N. This is in analogy to the fact that for many-qubit systems, local observables can have low error even when the global fidelity of the state is essentially zero <cit.>. Therefore, the linearly growing decay rate of number states in cavities need not lead to particular difficulties for quantum simulation when considering such local observables. Moreover, the effective errors induced by qubit decay in the Fock-binary encoding are highly non-linear in N and also lead to large jumps in occupations. For example, Fock states with N=2^n, n∈ℕ, decay immediately to N=0, leading to an effective loss of 2^n bosons at once. Hamiltonian evolution with non-linear terms such as an on-site interaction in the presence of Fock-state-dependent loss can lead to highly correlated errors. This is not the case for bosonic hardware, where the decay rate of ⟨n̂^2|$⟩ for an initial Fock state is a smooth function of the Fock state number and in a single decay event, only a single boson is lost. Finally, in cavity QED, boson loss is far slower (decay rateκ≈1kHz, with some cavities shown to exhibitκ<1Hz) than qubit decay (decay rateγ≈5-10kHz). The excellent coherence of cavities is illustrated by the recent experimental preparation of coherent (i.e., definite parity) photon cat states with as many as1000photons <cit.>. Therefore, the differing nature of errors in oscillators, which is closer to physical errors appearing in bosonic systems, can pose advantages for qubit-boson hardware compared to all-qubit simuluations of bosonic systems.
 
As to the prospect of quantum advantage, we believe that non-equilibrium dynamics in regimes with large boson number fluctuations, starting in Fock states, will be challenging to simulate classically: the widely-used time-evolving-block-decimation algorithm for tensor networks scales cubically with the boson number cutoff <cit.>. Monte-Carlo methods struggle with non-equilibrium and real-time dynamics, while semi-classical methods such as the truncated Wigner approximation <cit.> with initial states that are non-Gaussian. Finally, brute-force sparse diagonalization methods struggles with large systems – we expect dynamics of system sizes ofL>20bosonic sites with an average of two bosons per site to be challenging to simulate when considering initial states leading to large boson fluctuations such as|0⋯0LL0 ⋯0⟩, in which case the cutoff needs to be chosen to be essentially given by the number of bosons in the system,2L, leading to a Hilbert space size of3L-12L. The largest simulations to date have reached Hilbert space sizes of<2^40<cit.>. Hence, we expect to surpass classical simulability using sparse matrix methods forL≈16. This is in contrast to qubit-only simulations, where systems ofL≳50are needed to reach beyond-classical regimes <cit.>.



 §.§ Future directions

Our study opens a promising pathway towards using native oscillator-qubit hardware to simulate boson-fermion-gauge field systems. Compared to qubit-based hardware, this approach will be most advantageous when simulating systems with strong boson number fluctuations, for example in chaotic field theory evolution <cit.>, as such a case requires large cutoffs and, consequently, the synthesis of  square-root factors (when acting with creation and annihilation operators) that are expensive to synthesize in all-qubit hardware. A key advantage of our digital approach compared to classical or analog methods is its versatility – our approach can simulate topologies that do not match the hardware topology. For example, a3D model can in principle be mapped to the2D array at the expense of needing to implement long-range entangling gates with SWAPs.

Extending our approach to non-Abelian Gauge theories, most importantlySU(N)to simulate the weak and strong forces, are most pertinent. To this end, our Schwinger-Boson approach can be relatively straightforwardly generalized <cit.> and we expect similar scaling advantages to materialize as we found forU(1)fields.

Many applications of qubit-oscillator quantum simulation such as vibronic excitations of molecules or phonon-electron interactions in condensed-matter physics show similar decoherence processes to the one present in qubit-oscillator hardware. This opens up the direction of hybrid qubit-oscillator digital simulation of open-system dynamics, which would be an interesting topic to explore next. In noisy qubit hardware, this would not be immediately possible as the physical loss processes do not resemble boson loss in Fock-binary encodings as discussed above.

Looking into the future beyond noisy simulations, our work motivates the study of bosonic error-correcting codes, i.e., codes that directly encode a bosonic logical operator. Furthermore, it would be interesting to explore, whether our proposed hybrid oscillator-qubit approach also yields advantages for equilibrium state preparation away from the ground state. To that end, quantum signal processing <cit.>, methods based on filtering <cit.>, or thermalization <cit.> could be employed. Another interesting direction is to develop explicit benchmarks of dynamical quantum simulations of, e.g., string breaking and particle scattering problems <cit.> to better compare hybrid oscillator-qubit and qubit-based approaches.

While challenging in cavity QED hardware, a purely native approach for fermion-boson problems could be explored by combining our approach with digital fermionic quantum processing in neutral atom arrays <cit.>, where bosons can be encoded in (the mechanical) harmonic oscillator degrees of atoms in tweezers <cit.>. 




§ ACKNOWLEDGMENTS


This research was funded by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (C2QA) under contract number DE-SC0012704.  C2QA led this research. Support is also acknowledged from the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator. S.K. is supported with funds from the Ministry of Science, Research and Culture of the State of Brandenburg within the Centre for Quantum Technologies and Applications (CQTA). J.M.M. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. 2141064. L.F. is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as part of the CRC 1639 NuMeriQS – project no. 511713970.

We acknowledge helpful discussions with Marko Cetina, Zohreh Davoudi, Alexey Gorshkov, Fabian Grusdt, Mohammad Hafezi, Or Katz, and Torsten Zache. 

External interest disclosure: SMG is a consultant for, and equity holder in, Quantum Circuits, Inc.