Simulating gauge theories with variational quantum eigensolvers in superconducting microwave cavities

Quantum-enhanced computing methods are promising candidates to solve currently intractable problems. We consider here a variational quantum eigensolver (VQE), that delegates costly state preparations and measurements to quantum hardware, while classical optimization techniques guide the quantum hardware to create a desired target state. In this work, we propose a bosonic VQE using superconducting microwave cavities, overcoming the typical restriction of a small Hilbert space when the VQE is qubit based. The considered platform allows for strong nonlinearities between photon modes, which are highly customisable and can be tuned in situ, i.e. during running experiments. Our proposal hence allows for the realization of a wide range of bosonic ansatz states, and is therefore especially useful when simulating models involving degrees of freedom that cannot be simply mapped to qubits, such as gauge theories, that include components which require infinite-dimensional Hilbert spaces. We thus propose to experimentally apply this bosonic VQE to the U(1) Higgs model including a topological term, which in general introduces a sign problem in the model, making it intractable with conventional Monte Carlo methods.

1 Introduction Quantum computers have the potential to solve currently intractable problems in fundamental and applied sciences.In particular, hybrid quantumclassical approaches bring interesting applications closer in reach by lowering the requirements for quantum hardware.In this paper, we pursue this direction by considering a variational quantum eigensolver (VQE) [1][2][3][4].This algorithm takes advantage of a quantum processor to efficiently generate a parametrized ansatz state.An appropriate cost function is measured on the state and fed into a classical processor.The classical processor uses this information to minimize the cost function by optimizing the parameters that control the ansatz state (see App.A for more details).VQEs exploit the quantum processor to efficiently evaluate a cost function that is hard to calculate classically, while the variational nature of the optimization algorithm ensures resilience to certain types of errors and lowers requirements for the quantum hardware.This approach has been shown to be extremely successful [5][6][7][8].
While VQEs have been widely implemented with qubit-based quantum hardware [9][10][11][12][13][14], we propose in this work a bosonic VQE using superconducting microwave cavities.Our approach is especially useful when simulating models involving degrees of freedom which are not well-described as qubits.Important examples are gauge theories that include quantum fields to be simulated, which are defined in infinitedimensional Hilbert spaces.Photonic systems are apt candidates for simulating such models, but the advantage of the large photonic Hilbert space usually comes with the disadvantage that strong nonlinear interactions are not easily available.By contrast, our proposal builds on recent advances in superconducting microwave cavities that allow for strong nonlinearities between photon modes [15][16][17][18][19].In su-perconducting microwave cavities, Josephson junctions are an important source of strong nonlinear dynamics: through parametric processes it is possible to easily select and enhance a variety of interactions between the cavity modes.This flexibility in the available interactions is an important resource in a VQE, since they determine the form of the ansatz state and are fundamental in the success of the algorithm.In microwave cavities, nonlinearities have been demonstrated to be highly customisable, and can be programmed in situ, i.e. during running experiments, by changing classical microwave fields generated at room temperature.Using these programmable microwave fields for realizing VQE schemes allows, therefore, for the realization of a wide range of bosonic ansatz states.
As a concrete example, we propose in this paper to simulate topological terms in gauge theories.Topological terms are particularly challenging for classical numerical studies of lattice gauge theories since they give rise to a sign problem [20][21][22].However, they play an important and fascinating role in quantum field theory.It would be very interesting and important to determine the amount of CP violation emerging from QCD itself to shed light on the the strong CP problem [23] and, in turn, on physics beyond the standard model [24].In analogy to QCD, also in the electroweak sector of the standard model, a CP violating topological term can be introduced which can lead to enhanced baryon number violating processes [25].This mechanism can thus provide one of the basic ingredients for generating the matterantimatter asymmetry of the universe.A topological term can already be realized in one spatial dimension using a compact formulation of a U(1) gauge theory.This allows us to study the effect of such a term, the strength of which is given by the parameter θ, in a much simpler setup than the full electroweak sector and constitutes a very first step towards addressing the complete gauge Higgs sector of the standard model.In addition to this perspective within the high-energy physics context, a gauge-Higgs model with a topological term is also closely related to interesting models in condensed matter physics [22,26,27].We note that compact 1+1 dimensional U(1) gauge theories in the presence of a nonzero topological term are also of interest by themselves.In particular, the topological term can lead to new phenomena [28][29][30][31][32][33], and an enriched phase structure of the model [20,27,34].The rich physical content found in this model can in general not be explored with conventional Markov Chain Monte Carlo (MCMC) lattice methods because the topological term induces a sign problem in Euclidean time (it has been realized though that in the special case of 1+1 dimensions the model can be written in a dual formulation which is sign-problem free [20,35]).Importantly, the topological term studied here can be extended to three spatial dimensions [36], showing a promising path to explore new physics with quantum simulators.
More specifically, we present an experimental proposal for simulating the U(1) Abelian Higgs lattice gauge theory in one spatial dimension with a topological term and open boundary conditions [37][38][39].We show how to use a VQE in order to study the phase diagram of the model.We also explore the physics of the model using matrix product states (MPS) [40][41][42].This allows us to diagonalize the Hamiltonian on lattices of various sizes and various truncations of the Higgs field's Hilbert space.We find that, despite the finite lattice size and truncated operators, we are able to capture the physics of the main features of the phase diagram from analytic results, and conclude that it can be simulated on near-term quantum devices.
The rest of the paper is organized as follows.In Sec. 2 we introduce the gauge theory under consideration, and explain how it can be mapped to a model suitable for simulations with parametric cavities.We then describe the experimental system, and explain how the experimental resources can be adjusted to perform a bosonic VQE for the U(1) Higgs model with topological term.In Sec. 3 we discuss in detail the VQE protocol and results.To put theses results into context, we summarize the main features of the phase structure of the U(1) Higgs model in Sec. 4. Subsequently we numerically demonstrate, using MPS, that we can capture the relevant physics of the model, even for the limited system sizes that can be simulated on current quantum hardware, and we systematically explore finite-size effects in Sec. 5. We summarize our findings and offer some possible perspectives in Sec. 6.

VQE with parametric cavities
In this section, we discuss how the physics of a gauge theory can be studied using a superconducting microwave cavity, using the one dimensional Higgs model with topological term as concrete example.As outlined in the introduction, we propose a VQE approach, which forgoes the need to implement the complicated interactions that appear in the Higgs model on the quantum simulator.Instead, the VQE protocol exploits a set of resource Hamiltonians that can be realized on a given platform, and allows for the preparation of an output state that approximates the targeted ground state using the limited set of resource interactions available (see App.A for a more detailed discussion about VQEs).We target quantum simulations with parametric microwave cavities [15,17,18], which have not been used for VQE so far, but are a promising candidate system as we explain below.In particular, we first introduce in Sec.2.1 the Hamiltonian of the model, which involves both scalar Higgs and gauge fields.We then describe the microwave platform and the features that make it suitable for our simulation in Sec.2.2.In Sec.2.3, we give the mapping from the fields of the original model to the photonic modes experimentally available.We present the specific resource Hamiltonian that can be engineered in the system, and the measurement scheme necessary to run the VQE in Sec.2.4 and Sec.2.5 respectively.We conclude by discussing the experimental imperfections and their role for the proposed simulation in Sec.2.6.

The Hamiltonian
Contrary to conventional lattice methods, for the purpose of quantum simulations, it is advantageous to work with a Hamiltonian lattice formulation.Here, we describe the U(1) Higgs lattice Hamiltonian with a topological term in one spatial dimension with open boundary conditions.
The Higgs fields are defined on the lattice sites, and the gauge fields are defined on the links between the lattice sites.In particular, the Higgs field on site n is φn and it has a canonically conjugate operator Qn , called the charge operator.Considering the case of fixed length of the Higgs field, these operators satisfy the commutation relation [37][38][39] [ Qn , φ † and as a result the Higgs field φn acts as a lowering operator for the eigenstates of Qn : The operators Ûn and Ên (called electric field operator) are associated with the gauge field on link n, which joins lattice sites n and n+1.These operators satisfy the commutation relation [37][38][39] [ Ên , and Ûn is a descending operator for the eigenstates of the electric field Ên : The U(1) Higgs Hamiltonian for one spatial dimension is given by [39] where N is the number of lattice sites, β = 1/g 2 , g is the coupling strength, and R 2 is inversely proportional to the mass of the Higgs field.Additionally, ε 0 is the background electric field.In one spatial dimension, the topological term is proportional to the background electric field, with the θ given by 2πε 0 .[43].In Eq. ( 5) we have fixed the lattice spacing to be a = 1 and throughout the paper we use natural units ℏ = c = 1.The first and third term in the Hamiltonian describe the Higgs field energy and electric field energy respectively.The second term is referred to as the kinetic term; it allows charge to be transferred between adjacent lattice sites, at the expense of changing the electric field between the sites.The presence of the gauge field operator in this term ensures the local gauge symmetry of the model is conserved.
In the Hamiltonian formulation of the model, physical states (i.e.gauge-invariant states) have to obey Gauss's law where we have defined for each lattice site Note that the Ĝn operators commute with the Hamiltonian and are the generators of timeindependent gauge transformations.The eigenvalues G n take integer values, and can be interpreted as static charges that can be introduced at every lattice site.For the rest of this paper we focus on the sector of vanishing static charges, G n = 0 ∀n.The Hamiltonian in Eq. (5) conserves the total charge As a result, the Hilbert space can be divided into subsectors, each containing states with a definite total charge.For the remainder of the paper, we focus on the ground-state properties of the theory in the Q total = 0 subsector.In order to ensure that the ground state is in the correct subsector during numerical calculations, we can add a penalty term to the Hamiltonian When the weight ℓ of this term is large enough, the states outside of the Q total = 0 sector are penalized and removed from the low-lying spectrum.Note that while this term is necessary for matrix product state results (discussed in Sec.5), it is not needed for the cost function of the VQE, since we choose an initial state in the correct sector, and a variational circuit that preserves the total charge.In one spatial dimension and with open boundary conditions, it is possible to integrate out the gauge field's degrees of freedom and work with an effective Hamiltonian described by the Higgs degrees of freedom only [5,[44][45][46].This procedure ensures that the eigenstates of the effective Hamiltonian are gauge invariant, and therefore allows for a contained quantum simulation, that takes place only in the gauge invariant (physical) subspace.
As described in more detail in App.B, the effective Hamiltonian is given by Notably, the elimination of the gauge fields introduces long-range interactions between the Higgs fields.It also allows one to simulate the model using fewer modes, as the redundant degrees of freedom have been removed, making it more accessible to quantum hardware.
In order to study the phase space of the model we consider as order parameter the electric field density (EFD) of the model, which is defined as where the expectation value is taken w.r.t. the ground state.

Microwave-photon cavity
Here we detail the microwave-photonic hardware that we consider for our VQE of the U(1) Higgs Hamiltonian with a topological term.As shown schematically in Fig. 1, the quantum processor of our VQE scheme consists of a multimode coplanar waveguide resonator terminated by a SQUID at one end [15,17].The SQUID is coupled both to a microwave pump mode for classical control and to the total flux of the cavity.The SQUID -consisting of a superconducting loop with two Josephson junctions -provides a high degree of dissipationless, nonlinearity in the system.This basic element is the key to the success of superconducting computing architectures [19,[47][48][49] and is used here to control microwave fields.θ determine how the SQUID is pumped and consequently which interactions are generated between the cavity modes to generate the ansatz state ψ( ⃗ θ) .A quantum-limited amplifier allows for performing heterodyne measurements on the cavity modes to obtain the cost function C( ⃗ θ).In our case the cost function is the expectation value of the Hamiltonian for which we want to find the ground state, as defined in Eq. (18).C( ⃗ θ) is used by the classical processing unit (CPU) to find the optimal variational parameters that minimizes it.
The several resonant modes available in the cavity are each mapped to one of the Higgs fields as described in Sec.2.3, so that the original spatial lattice of the model is simulated on a lattice in synthetic dimensions.Using microwave photons as quantum degrees of freedom in a VQE has the advantage that a larger Hilbert space can be accessed compared to traditional, qubit-based VQE protocols.This property is a naturally good fit for the simulation of gauge or Higgs fields, which are defined on infinite-dimensional Hilbert spaces.As a result, using bosonic degrees of freedom requires fewer experimental modes, and avoids the translation of twobody terms into three-and four-body interactions, which usually occurs when when encoding multilevel gauge fields using qubits [50].
As explained in Sec.2.4 and shown in Refs.[15,17,18,51], the SQUID in the experimental setup is an important source for nonlinear interactions, with its cosine potential providing access to a range of high-order nonlinearities.The pump mode provides a high degree of classical control over which of the available interactions are activated at a given time.Furthermore, the flexibility of circuit fabrication, allows for the design of other types of Josephson circuits, related to the SQUID, that can enhance or suppress certain nonlinearities [52].This means one can create highly customizable gates by controlling the interacting modes, the type of interaction, its strength and relative phases.In App.C, we explain in detail how to use the nonlinear operations to create interactions involving the product of four photonic operators, which are employed in our simulation of the Higgs model, but in principle these ideas can be extended much further.This flexibility is particularly suited for special purpose VQEs, including the one considered in this paper.In addition, it is natural for this platform to generate effective couplings that are nontrivial (see Sec. 2.4 and App.C for details), which can be beneficial to other VQEs of lattice gauge theories.Overall, our approach is complementary to existing VQE schemes in the sense that it has strikingly different features than qubitbased protocols, but also uses resource Hamiltonians that are different from other bosonic platforms, such as optical photons [53,54] or ultracold atoms [55].
For a microwave cavity the density of available modes can be increased by increasing the physical length of the cavity.Within the 8GHz of measurement bandwidth, we can reach approximately 10 modes with a frequency spacing of 100MHz within 1GHz, while keeping the same order for the strength of the necessary interactions.In order to increase further the number of available modes, more microwave cavities can be coupled parametrically to each other, as has been already demonstrated with hundreds of cavities in Ref. [56].We consider this a viable path to scale the proposed system to simulate larger lattices.

HOBM mapping
In order to simulate the Higgs Hamiltonian in Eq. (10) using a microwave cavity, we map its operators Qn and φn to photonic degrees of freedom using the Highly Occupied Boson Model (HOBM) [57][58][59], which is given by In the expression above, N 0 ∈ N is a constant such that the zero charge state in the Higgs model maps to the Fock state with N 0 photons (see Fig. 2).N and â are the photonic number and lowering operators, respectively.With this mapping, the U(1) Higgs Hamiltonian becomes The Hamiltonian given in Eq. ( 13) is now expressed in terms of the microwave modes available in the superconducting cavity.Note that the conservation of the total charge given in Eq. (8) corresponds to the conservation of the total number of photons N n=1 Nn .
FIG. 2: Schematic of HOBM mapping.The charge basis state (left) |Q⟩ with charge Q is mapped to the Fock state (right) with Q + N 0 photons |Q + N 0 ⟩ under the HOBM mapping.In the HOBM, â is proportional to the Higgs field operator φ according to Eq. (12b).

Resource Hamiltonians
In this section we discuss the resource Hamiltonians that are used in the VQE and how they can be generated by the microwave-photon platform's interaction Hamiltonian.The resources considered in the present work conserve the total number of photons, since it is a symmetry of the considered model, as discussed more in detail in Sec.2.3.However, we note that in general it is possible to realize also non photon number conserving interactions between photonic modes [17].By driving the pump at a frequency equal to the difference in the cavity mode frequencies of two modes c and d [17], we can create a c − d beam splitter interaction (see Ref. [51] for details) where g is the strength of the interaction, e iλ is a phase in the beam splitter interaction, and ⃗ Ω = (Ω 1 . . ., Ω N ) are the rotation frequencies of the modes obtained as the difference between the natural frequency and the frequency of a suitable rotating frame.This interaction induces the unitary transformation Û (t) = e −it Ĥ(c,d) BS , and we can use the time t, the phase e −iλ and the rotation frequencies ⃗ Ω as classical variational parameters.All three parameters can be optimized by the classical algorithm during the VQE, and while the first two can be controlled experimentally, the choice of rotation frequency can be made through a suitable change of frame of reference.Note that this interaction closely resembles the kinetic term in Eq. (13).
The experimental platform considered here has proven to be able to generate higher order non linear evolutions, such as the cubic interactions demonstrated in Ref. [17] and the underlying principle can be extended.In App.C we show that by going into a suitable rotating frame, and by driving the SQUID at a frequency significantly lower than the other characteristic frequencies of the system, we can obtain the following effective interaction The basic working principle to create such an interaction has been experimentally demonstrated; based on those, it is expected that this interaction can be realized with a coupling strength up to g ′ = 10 7 Hz [17,60].When applying the time evolution induced by the resource Hamiltonian in Eq. ( 16), we use the interaction time t and the frequencies ⃗ Ω as variational parameters.Note that the resource Hamiltonians Ĥ(c,d)

BS
and ĤNN are used to prepare the variational state for the VQE as indicated in Eq. (21).Importantly, these resources contain terms of the form Nn and N 2 n , both of which appear in the target HOBM Higgs Hamiltonian given in Eq. ( 13).Since ĤHOBM lies in the Lie algebra spanned by the resource Hamiltonians, the VQE is able to reach the ground state of the theory [61,62], without the need to exactly realize it.

Measurements of microwave photons
In contrast to the optical domain, a natural measurement for microwave photons is linear amplification.Here we consider, in particular, a well-established measurement available for microwave photons: the linear phase-insensitive amplifier, which simultaneously amplifies both quadratures of the input field [51,63].All the cavity modes are coupled out to the measurement line through a single coupling element that can be modulated [64] to adjust the coupling strength of all modes simultaneously.The resonator fields are then absorbed by the amplifier, as part of the measurement process.Quantum mechanics dictates that amplifying the field in a phaseinsensitive manner adds at least one unit of vacuum noise to the signal, once the signal is amplified.It is therefore natural to simultaneously measure two conjugate quadratures, since this does not affect the signal-to-noise ratio.The signal obtained is equivalent to a heterodyne measurement of the cavity modes.The recent development of ultra broadband quantum-limited amplifiers (JPAs and TW-PAs) improved significantly the measurement speed and signal-to-noise ratio [65] Additionally, higherorder coherence functions can also be calculated [66] to obtain statistics of nonlinear measurements that cannot be performed directly (for example, photonnumber measurements).
The quantum state emitted from the microwave cavity is amplified in the large-gain limit (see Ref. [17] for details).Measuring the field quadratures at the two output ports corresponds to simultaneously measuring the self-adjoint part XS = ( Ŝ + Ŝ † )/2 and the anti self-adjoint part PS = ( Ŝ − Ŝ † )/(2i) of the signal operator Here â refers to the cavity output state ψ( ⃗ θ) to be measured (see Fig. 1) and ĥ is the amplifier noise mode.Note that Ŝ is a normal operator, meaning [ Ŝ, Ŝ † ] = [ XS , PS ] = 0. Furthermore, the eigenvalues of XS and PS correspond to the real and imaginary parts of the eigenvalues of Ŝ, respectively.Repeated preparation and detection of the output state of the cavity yields a measurement distribution D(S) (shown in Fig. 1), from which any statistical moment of Ŝ and Ŝ † can be calculated [66].
During the VQE of the HOBM Higgs Hamiltonian ĤHOBM in Eq. ( 13), the cost function is calculated from the data comprising the measurement histogram D(S) as follows.First, we express the expectation value of the HOBM Higgs Hamiltonian in terms of the statistical moments of the measured signal operator Ŝ, The above mapping of the cost function from the target Hamiltonian ĤHOBM can be obtained using which assumes that the noise modes are in the vacuum state |0⟩ (i.e. this is a quantum-limited measurement) and the noise is uncorrelated with the cavity output modes.Finally, one can experimentally extract the value of C( ⃗ θ) from the last part of Eq. (18).In fact, expectation values of moments of the Ŝ operators can be obtained from the measured data distributions D(S) using the following expressions [67] Note that the indirect measurement of C( ⃗ θ) through Eq. (18) will cause the outcome's statistical variance to be larger compared to a direct measurement of the target Hamiltonian ĤHOBM .The effect of a finite number of repeated measurement will be discussed in the next section.

Measurement budget and experimental errors
In a VQE, the cost function C( ⃗ θ) must be computed many times by the quantum device over the course of the experiment (see App. A).Each cost function evaluation involves the preparation of the initial state, applying a sequence of gates, and the measurement of the final VQE state (together, we refer to this as a shot).As statistical and amplification noise are present in the cost function evaluation, the number of repeated measurements of the state M determines how accurately we know the cost function.In particular, for a large number of measurements the variance of the cost function is given by σ 2 H /M , where σ 2 H is the intrinsic variance of the Hamiltonian for the indirect measurement obtained through Eq. (18).Therefore, increasing the number of measurements reduces the statistical noise of the cost function, which means the number of measurements that can be performed and the repetition rate are crucial quantities when designing an experiment.
In the following, we discuss the measurement budget for the considered microwave photon setup in more detail, as it is very different from that of common VQE platforms.The time required to prepare the initial state can be very short, for example the fastest time to prepare the Fock state |2222⟩ is estimated to be on the order of 10 −8 s [68].The upper limit on the time required to apply all the gates in the circuit, including the beam splitter in Eq. (15) and the ĤNN interactions in Eq. (16), is given by the cavity decay time, which ranges between 10 −7 s and 10 −4 s [64].However, some gates can be applied considerably faster than that.For a cavity with fixed coupling, the time required to measure the state is also limited by the cavity decay time, as we must allow the state to escape the cavity before measuring it.There is however the possibility to employ cavities with tunable coupling to the environment [64]; the cavity can be closed to have a long cavity lifetime time and opened for a fast measurement process.Assuming a tunable coupling, the readout time is negligible compared to the gate application time.
As a result, the time budget for a single shot is dominated by the coupling strength, which is on the order of 10 7 Hz for the interactions considered in Sec.2.4 [17,60].Assuming the VQE circuit consists of about ten gates as considered in Sec. 3, and the time of application t of the Hamiltonians given in Eqs.(15)(16) is of order ≈ 1 g ≈ 1 g ′ , as is the case of our VQE (see Sec. 3), the time for a single shot is therefore on the order of 10 −6 s.One advantage of the microwave-photon platform over other potential VQE platforms is its ability to run continuously without the need for human intervention.As a result, ignoring the time required for classical computing in the VQE feedback loop, up to 10 11 measurements can be collected in a single day.
This discussion has so far ignored the effects of experimental imperfections present on this platform.The sources of error on the microwave-photon platform are similar to other gate-based superconducting processors and as a result we can expect a similar performance [69].We can also expect similar improvements to the platform as large, commercial entities continue to improve and develop gate-based quantum computers.In the three steps comprising a shot -state preparation, gate application, and measurements -we therefore expect the measurement to pose the biggest concern.The main issue is that measuring higher moments of the field operators, as required in Eq. (20) for instance, requires averaging times that increase polynomially in the experimental signal-to-noise ratio [66].That is, measurements of the higher moments are increasingly sensitive to measurement noise.To understand the impact of this, we present in Sec. 3 a realistic simulation of a VQE, which takes the measurement errors fully into account.

Microwave VQE for the U(1) Higgs model
In this section we propose and classically simulate a concrete and realistic microwave VQE protocol, that allows for the study of the U(1) Higgs theory with topological term as described in Sec.2.1.As a specific application example, we show how to experimentally detect the first-order quantum phase transition that appears in the model for values of R 2 below the critical value R 2 c .We further discuss how future microwave-based quantum simulators can study the smoothing of the phase transition if the parameter R 2 approaches the critical point.As discussed in the following, below the critical point, a realistic experiment would be able to find the ground state by repeating the measurement shot up to M = 10 7 times for each cost function evaluation, leading to VQE computation time of up to few hours.On the other hand, closer to the critical point the noisy measurement is more of an issue and the proposed VQE would require repetitions of at least M = 10 9 , therefore will benefit an increased experimental repetition rate.For a detailed description of the full phase diagram of the model will be given in Sec. 4.

Simulation and noise modelling
The detection of the considered quantum phase transition is performed by using a VQE for an approximate preparation of the ground state of Eq. ( 13), followed by the subsequent measurement of a suitable order parameter, namely the EFD introduced in Eq. (11).
To analyse the experimental feasibility of the observation of the targeted phase transition, we carry out classical simulations of the proposed VQE protocol.Our simulations include the main relevant source of imperfections, namely the noise added during amplification of the microwave signal that leads to imperfect measurement of the cost function, and a finite precision with which it is known.As detailed in Sec.2.5, for a given variational state ψ( ⃗ θ) prepared in the superconducting cavity, we have access to a high number of repeated measurements M .This allows us in our classical simulation to approximate the probability distribution of the measurement outcome as a Gaussian distribution with mean ψ( ⃗ θ) ĤHOBM ψ( ⃗ θ) and variance H can be obtained without considering the redundant noise modes in our classical simulation by suitably inverting Eqs.(19).
In the following, we consider a lattice consisting of N = 4 lattice sites and an HOBM mapping with N 0 = 2.For an analysis of the effects of the finite lattice size N and truncation effects due to the finite value of N 0 , see Sec. 5. To simulate the infinitedimensional photonic Hilbert space realized in the experiment in our simulation, we truncate the local Hilbert space of each microwave photon mode to a five-dimensional subspace.For more details on the truncation effects of our numerical description of the proposed experiment, see App.D (note that the values of N and N 0 affect the capability of the experiment to observe the phase transition, while the truncated numerical description of the photon modes by a five-dimensional subspace only affects our ability to predict the outcome of the experiment).

VQE detection of the phase transition
With N 0 = 2, the vacuum state in the model corresponds to the product state with two microwave photons in each mode |2222⟩, which we choose as the initial state for the VQE.As explained in Sec.4.2 and in App.G, the main contributions to the ground states in the two limiting cases corresponding to small and larger values of the background field ε 0 can be determined to be |2222⟩ and |1223⟩ respectively.These extremal cases motivate the choice of the variational quantum circuit: where N l layers are applied that each involve the beam splitter interaction and the N N -interaction, introduced in Sec.2.4.Optimizing over ⃗ θ and ⃗ Ω, the VQE protocol involves 2N l + 4 variational parameters (the relative beam splitter phase λ in Eq. (15) provides a possible additional parameter, which is not used and instead is zero in this example).
Figure 3 shows that the first-order phase transition that the model undergoes for the value R 2 = 0.3 < R 2 c can be studied in a microwave VQE experiment by observing a discontinuity in the order parameter F as a function of the background field ε 0 [F is here defined as the ground state expectation value of the EFD of the model given by Eq. (11)].The optimization required for each point is on the order of 10 3 cost function evaluations.We found that before (after) the discontinuity at ε 0 ≈ 1.6, M = 10 3 (M = 10 5 ) measurements per cost function evaluation were required for the used approach.
Points after the discontinuity are more challenging.Since we adopt a stochastic optimization algorithm, for each value of ε 0 the optimization is repeated fifteen times, and we post-select the successful optimization as the one that has reached the lowest energy.Experimentally, this corresponds to an order of necessary shots of 10 10 , which is well within the feasibility limit identified in Sec.2.6.Fig. 4 shows an example of a successful post-selected op-timization for ε 0 = 1.9.As optimizer we employ the Bayesian Adaptive Direct Search [70], which combines a mesh-based search strategy [71] with local Bayesian optimization steps [72].
FIG. 3: EFD of the ground state for different values of ε 0 , R 2 = 0.3 and β = 1.The ground state is calculated with a classical simulation of a VQE, including the statistical noise of the cost function, assuming M = 10 3 (M = 10 5 ) repeated measurements for points before (after) the discontinuity at around ε 0 ≈ 1.6.The discontinuity is a sign of the first-order transition in the ground state of the model discussed more in detail in Sec.4.2.Points and error bars are respectively the mean and standard deviation obtained with a sample of ten successful optimizations.Points after the discontinuity are post-selected as discussed in the main text.The continuous line is obtained from exact diagonalization results.

VQE close to the critical point
In the following, we explore the possibility to probe the quantum phase transition of the model in the vicinity of the critical point R 2 ≈ R 2 c , where the quantum simulation becomes most challenging.In this parameter regime, the ground state becomes highly entangled and hard to reach.For R 2 > R 2 c , the first-order quantum phase transition disappears, as explained in Sec. 4. This behavior can in principle be probed by a VQE protocol.To illustrate this, we provide in App.H a quantum circuit that is experimentally more demanding than the example shown in Sec.3.2, but allows for the VQE preparation of the ground state of the model for the difficult parameter area around R 2 = 1.In Fig. 5, we show the classically simulated VQE result for this case including measurement noise only, i.e. in the absence of statistical errors (corresponding to M → ∞).We also show how the VQE simulation with finite measurement budget approaches this limit for increasing values of the number of measurements per cost function evaluation M = 10 6 and M = 10 9 .The plot shows that while the proposed strategy can in principle obtain the ground state in the most challenging areas of the phase diagram, noise of the cost function has a more detrimental effect than in the case of the clear firstorder phase transition.For an experimental study of the ground state around R 2 = 1, it would therefore be necessary to devise a suitably adjusted measurement scheme or to increase the experimental repetition rate.
We note that in the parameter regime around R 2 = 1 , the sharp discontinuity that is visible in Fig. 3 is expected to become smooth.The analysis in Sec. 5 confirms that the disappearance in the discontinuity is indeed a genuine feature of the phase diagram of the model and not an artifact of finitesize or truncation effects (see Fig. 8).

Phase structure and ground-state properties
In the following, we provide more information on the phase transition that is considered as concrete application example for our VQE approach in Sec. 3. Despite its simplicity, the U(1) Higgs model shows a rich phase diagram, especially in the presence of a topological term.For the reader's convenience in this section we briefly review the phase structure of the model, before systematically exploring how limitations of current small-scale quantum hardware affect this picture in Sec. 5. We first focus on the case without a topological θ-term in Sec.4.1 before moving on to the phase structure in the presence of a θ-term in Sec.4.2.In particular, we discuss the signature of detecting a quantum phase transition when ε 0 in Eq. ( 5) is varied and which can be observed in a quantum simulation even for small system sizes.This quantum phase transition is of particular interest for the microwave-based VQE experiments proposed in the previous section, as it is inaccessible to conventional MCMC methods due to a sign problem.

Phase structure in the absence of a topological term
In the absence of a topological θ-term, the Higgs model with fixed length of the field can be studied numerically with conventional MCMC lattice methods [37,38,73].An intuition for the phase structure can also be obtained by examining the Lagrangian using a simple semiclassical approach (see App. E for details).In this picture, the scalar field φ is assumed to fluctuate only slightly around the minimum of the potential term V (|φ|) in the Lagrangian for the ground state of the model.Depending on the value of the mass and the coupling, the shape of V (|φ|), and consequently the nature of the ground state, changes as shown in Fig. 6.For large values of the mass and small inverse coupling, or equivalently small values of R 2 and β, the potential has a unique minimum.In this region, the Abelian Higgs model essentially corresponds to a pure U(1) gauge theory, describing nothing but a massless photon.Since pure gauge theories in their compact version show charge confinement, we refer to this part of phase diagram as the confining region.
Going into the opposite corner of the phase diagram in Fig. 6, characterized by small values of the mass and large inverse coupling, the potential V (|φ|) shows a "Mexican hat" structure with an infinite number of minima along the well of the Mexican hat.In this region, the photon acquires a mass as a result of the Brout-Englers-Higgs mechanism (see App. E for details) and we refer to this part as the Higgs region.The ground state in the Higgs region corresponds to one of the minima of the Mexican hat potential and, thus, spontaneously breaks the U(1) symmetry of the model.On a lattice with finite spacing, these two regions of the phase diagram are separated by a crossover that is disappearing for small values of the inverse coupling β.Taking the continuum limit, the crossover line ends in a Berezinskii-Kosterlitz-Thouless transition from a confining phase1 for large masses to the Higgs phase at small values of the mass [73].

Phase structure in the presence of a topological term
The phase structure of the model in the presence of a topological term has been investigated both theoretically [27] and numerically [21,22,74].Fig. 7 shows a sketch of the phase diagram as a function of ε 0 and R 2 for a fixed value of β.FIG.7: Sketch of the phase diagram in the R 2ε 0 plane at a fixed value of β.The dashed line indicates the critical line corresponding to a first-order phase transition which ends in a second-order quantum phase transition (indicated by the red dot).The confining region is indicated in orange, the Higgs region in blue.The inset shows the behavior of the EFD which, when crossing the phase transition, is showing a jump as a result of the first-order transition (see discussion in the main text).This behavior of the EFD serves as the signature of the phase transition for the here performed quantum simulation of the model.
Most notably, the physics is periodic in ε 0 with period 1.The origin of this periodicity can be understood intuitively by realizing that it is always possible to find an integer k, such that ε ′ 0 := k − ε 0 ∈ [0, 1).Inserting this into Eq.(5) we find for the electric-energy term while the other terms stay unchanged.The shifted electric field Ê′ n := Ên + k fulfills the same algebra as Ên in Eq. (3) and, thus, has the same spectrum as the original electric field operator (see Eq. (4a)).In other words, Ê′ n is unitarily equivalent to Ên .Moreover, since the spectrum of Ê′ n is given by the integer numbers Z, − Ê′ n is unitarily equivalent to Ê′ n .As a result, Ên + ε ′ 0 and Ên + ε 0 are related by a unitary transformation, thus showing that an integer shift in the background field does not change physics and we can restrict ourselves for the following discussion to ε 0 ∈ [0, 1).Moreover, the two points ε 0 = 0, 1/2 are special, because for these cases the Hamiltonian is symmetric under charge conjugation, meaning the exchange of particles and antiparticles (rigorous proofs for the periodicity and the charge conjugation symmetry are provided in App.F).
To get further insight into the phase diagram in Fig. 7, we can examine the Hamiltonian in Eq. (5) (or its equivalent version in Eq. ( 10)) in the limiting cases of large and vanishing mass of the Higgs field.Focusing first on the limit of small mass corresponding to R 2 ≫ 1, Eq. (10) which is a pure hopping Hamiltonian.In particular, we see that the expression above is independent of ε 0 and thus the physics does not depend on the topological term.Since we focus on the sector of vanishing total charge, the ground state in this limit is given by a superposition of all zero-total-charge states.
Looking at the opposite limit of large mass, or equivalently R 2 ≪ 1, the Hamiltonian in Eq. (5) reduces to The first term energetically favors a vanishing charge at every site, Q n = 0 ∀n.Inserting this into Gauss Law in Eq. (30) and taking into account that we focus on the sector of vanishing static charges, we see that for gauge invariant states the electric field has to take a constant value for all sites.For ε 0 < 1/2 the field configuration minimizing the electric energy in Eq. ( 23) is given by E n = 0 ∀n.In contrast, for ε 0 > 1/2 an electric field of E n = −1 ∀n minimizes the energy.Exactly at the point ε 0 = 1/2 both configurations yield the same electric energy contribution, and the ground state is doubly degenerate.These considerations show that in the limit of large masses the model undergoes a first-order quantum phase transition as we increase ε 0 from 0 to 1.This transition is accompanied by a discontinuity in the EFD (defined in Eq. ( 11)) which we expect for R 2 ≪ 1 to behave as as shown in the inset of Fig. 7.Moreover, the abrupt jump in the electric field configuration for ε 0 = 1/2 by one unit leads to a cusp in the electric-energy term.Hence, we expect the ground-state energy to show a non-smooth behavior at the transition point too.These expectations, which will be verified in Sec. 5, serve as the signature of the phase transition for the microwave-based quantum simulation of the model we perform.
Our discussion of the two limiting cases shows that the behavior of the model has to change going from large masses (corresponding to small R 2 ) to small masses (corresponding to large R 2 ).While for small masses a first-order quantum phase transition occurs for ε 0 = 1/2, this transition vanishes in the limit of large mass.Ref. [27] argued that the critical line at ε 0 = 1/2 ends in a second-order quantum phase transition belonging to the Ising universality class at a critical value R 2 c .The second-order transition at the endpoint of the critical line is accompanied by a spontaneous breaking of the charge conjugation symmetry, and has been observed in numerical MCMC simulations using a dual lattice formulation of the model [21,22].Moreover, despite the new features arising from the θ-term, the picture from the previous section remains true for large enough values of inverse coupling.For large mass, or equivalently small R 2 , there is a confining region whereas for large values of R 2 again the Higgs region occurs (see Fig. 7).
The microwave-based VQE approach detailed in Sec. 3 can be used to observe the first-order phase transition of the model in the presence of a topological term, as shown in Sec. 3. Since current quantum hardware is of small-scale, we will restrict ourselves to a small number of degrees of freedom and we have to work with lattices consisting only of a few sites.It can be shown that even for such a small system the considerations above hold true, but the location of the first phase transition shifts from ε 0 = 1/2 to some larger value of ε 0 .Focusing on the regime well below R 2 c and assuming that the mass of the Higgs field is large enough that the kinetic term can be neglected, we show in App.G that the first phase transition occurs at Equation (25) predicts the phase transition cannot occur before ε 0 = 1/2 for a finite lattice and only in the N → ∞ limit, the phase transition occurs exactly at ε 0 = 1/2.Moreover, the derivation shows that in this regime the the ground state before the phase transition is dominated by the state with vanishing charge at every site, |0⟩ ⊗N , whereas after the transition it holds a pair of charges with opposite sign and is dominated by |−1⟩ |0⟩ ⊗N −2 |+1⟩.

Spin Truncation and Matrix Product States
In this section we analyse the prospects for exploring the physics of the model using existing and nearterm small-scale quantum hardware.Since the number of modes available in a microwave cavity in the near future will be limited, we want to study in the following section if the limited size of the system that is studied allows to observe relevant physics.
To assess the feasibility of such an approach we first explore the effects of truncating the model to a small number of degrees of freedom numerically using matrix product states (MPS).For a system with N sites and open boundary conditions, the MPS ansatz reads (26) In the expression above, M i k k are complex square matrices of size χ for 1 < k < N , and are a basis for the d-dimensional local Hilbert space on site k.The parameter χ is called the bond dimension of the MPS and determines the number of variational parameters in the ansatz and limits the amount of entanglement that can be present in the state (see Refs. [40][41][42] for detailed reviews).The optimal set of tensors can be found variationally by iteratively minimizing the energy for each tensor while keeping the others fixed [75].
For numerical calculations with MPS, the dimension of the local Hilbert spaces has to be finite, which is in contrast to the infinite-dimensional degrees of freedom for each site in the Higgs Hamiltonian.Hence, they have to be truncated to a finite dimension.In Sec.5.1 we discuss a way of truncating the Hilbert spaces before numerically exploring finite-size effects on the phase structure in the presence of a θ-term in Sec.5.2.

Spin truncation
One possibility of truncating the theory is to replace the bosonic fields with integer spins where Ŝz and Ŝ+ are the z-component and raising operators for a particle with spin s, respectively, and Notice that the expression above approaches the correct commutation relation for bosonic field operators, [ φn , φ † n ′ ] = 0, in the limit s → ∞.The resulting spin Hamiltonian after applying the mapping in Eq. (27) to Eq. ( 10) can be addressed with MPS using standard methods, despite its long-range interactions [33,45,[76][77][78].

MPS results
In order to examine the effects of truncating the model to a small, finite number of degrees of freedom on the phase structure, we study the spin Hamiltonian at fixed coupling strength for a wide range of values of R 2 , s, and N .To estimate the error due to the finite size χ of the matrices in our numerical MPS simulations, we repeat the calculation for every combination of (R 2 , N, s) for a range of bond dimensions χ ∈ [10; 100].Afterwards we can extrapolate the results to the limit χ → ∞ following Ref.[33]. Figure 8 shows the MPS results obtained for the ground-state energy density E 0 /N and the EFD as a function of ε 0 for various system sizes and couplings.For all the results presented, we have chosen a penalty strength of ℓ = 3N which is sufficient to ensure that we are in the sector of vanishing total charge.
In general, we observe that truncation effects due to finite value of s are rather small compared to finite-size effects.In particular, for larger masses (corresponding to small R 2 ) even the simplest nontrivial truncation, s = 1, is sufficient to give the correct qualitative behavior.On the other hand, for the smallest masses we study, corresponding to R 2 = 1.0, results for s = 1 and s = 2 start to have a larger difference, in particular for larger values of ε 0 , and it is where the finite truncation shows the biggest effect.
In contrast to the spin truncation, finite-size effects are more pronounced.While we expect the physics to be periodic in ε 0 with period 1 as outlined in Sec.4.2, Fig. 8 shows that ground-state energy density as well as the EFD only show perfect periodicity throughout the entire range of ε 0 we study for our largest system size, N = 100.For smaller values of N the characteristic features still repeat, however the graphs for the energy density and the EFD show an overall increasing trend with ε 0 .Moreover, the period at which the characteristic features repeat is increasing with decreasing system size.
Focusing on our results for R 2 = 0.3 in Figs.8(a) and 8(b), corresponding to the largest value of the mass we consider, we clearly see the signatures of the first-order phase transition as discussed in Sec.4.2.The EFD shows sharp discontinuities accompanied by cusps in the ground-state energy for all system sizes we study.For our largest system size, N = 100, these transitions occur for ε 0 close to integer multiples of 1/2.Decreasing the system size, we observe that the location of the first transition gradually shifts to larger values of ε 0 , in agreement with Eq. (25).In particular, even for an extremely small system size of N = 4 the signatures of phase transition are still clearly visible in both the EFD and the ground-state energy density.
Going to a significantly smaller mass, the discontinuities and cusps in the EFD as well as the groundstate energy density vanish, as our data for R 2 = 1.0 in Figs.8(e) and 8(f) reveal.The smooth behavior of these observables is giving a clear indication that we are in the regime below the critical mass, or correspondingly R 2 > R 2 c , and the phase transition is gone.In general, finite-size effects seem to shift the critical mass at which the transitions vanishes towards smaller values (or equivalently R 2 c towards higher values), as our data for R 2 = 0.6 reveals.Looking at Figs. 8(c) and 8(d) we clearly observe a a first-order order transition signaled by a sharp discontinuity in the EFD for system sizes N ≥ 20.In contrast, the discontinuities are no longer present for the two smallest system sizes N = 4, 8 we study, and the EFD curves for these cases are smooth.
From our data in Figs.8(a) -8(d) we can extract the location ε ′ 0 of the first-order phase transition as the first discontinuity (cusp) in the EFD (energy density).Fig. 9 contains our results for R 2 = 0.3, 0.6 and various system sizes.Both panels show that finite-size effects shift the location of the first transition to values ε ′ 0 > 1/2.Comparing our numerical data to the prediction for the large mass limit R 2 ≪ 1 in Eq. (25), we observe good agreement for R 2 = 0.3.Going to a larger value of R 2 = 0.6, or equivalently to a smaller mass, the approximation in Eq. (25) breaks down and is no longer compatible with the numerical data.
In summary, our MPS data demonstrate that we can observe the relevant features of the model, even if we restrict ourselves to a small number of degrees of freedom.The truncation of the bosonic fields does not severely affect the physics in the presence of a θterm for the range of parameters we have considered.In particular, even for very small system sizes that are accessible with current quantum hardware, we can observe the characteristic features of the phase structure.For large masses, or equivalently small values of R 2 , the EFD and the ground-state energy density clearly show the signatures of the first-order phase transition which eventually vanish as we go to small masses or correspondingly large R 2 .We note that while for the parameters we have considered the MPS data is a useful benchmark for the proposed VQE, in general the VQE setup allows for going beyond these regimes, and, thus, for exploring physics that is not easily accessible with classical simulations.

Conclusions
We have designed a VQE protocol to be run on a superconducting microwave cavity.Here, the crucial advantage is the provision of bosonic degrees of freedom which allows to implement an infinitedimensional Hilbert space.In our case, truncation of different sizes are available for quantum simulation, offering the possibility to study and control the effect of the finite truncation.For a proposal to simulate gauge fields using the qudit space available in Rydberg atoms see Ref. [79].The readout of such a system is inherently different when compared to standard qubit-based simulations -it allows us to calculate different statistical moments from one single collection of measurement data.Moreover, it can be run continuously without the need for periodic, time-consuming recalibrations of the control electronics.Here we have shown that this platform is able to experimentally demonstrate an interesting feature of the lattice Higgs model with a topological term, namely its first-order quantum phase transition.In order to determine that indeed the physics of this model is probed, we used classical simulations based on matrix product states to study the effects of a finite lattice size and truncated operators.
The application of quantum simulations to lattice gauge theories has the potential to extend its reach beyond what has been possible so far with classical numerical methods, and has been a quickly developing field [11,12,50,[80][81][82][83][84][85].More specifically, studying lattice gauge theories that include topological terms is of particular interest, since they generate rich and interesting physics, and they are one of the settings that are, in general, intractable with conventional MCMC lattice methods.Our implementation of a bosonic VQE paves the way to extend this work in order to study topological terms in higher dimensions, with the ultimate goal to reach 3+1 dimensions.A work along these lines is [36].To that end, a VQE of a high-energy lattice gauge theory in three spatial dimensions is an important and ambitious milestone in the field of quantum simulation and a fundamental step towards answering some open questions about the universe.We note that in the case of higher number of spatial dimensions only part of the gauge fields can be eliminated through Gauss' law.This was done for example in Refs.[36,50,81].The remaining gauge degrees of freedom can be mapped appropriately to photonic modes similarly to what was done for the Higgs field in the present work, and there are no fundamental limitations to extending our work to higher spatial dimensions.An interesting extension of this work will be the study of real-time evolution of the lattice Higgs model in 1+1 dimensions.More specifically, the resource Hamiltonians presented in Sec.2.4 are extremely close to the terms appearing in the Hamiltonian, thus potentially offering the possibility of realizing a trotterized time evolution.Interesting dynamical effects include, for example, string breaking or the study of the phase diagram of the Higgs model with quenched systems [86].
For near term applications, there are many exciting opportunities to use the unique features of the microwave cavity platform (as explained in Sec.2.2) in upcoming VQE experiments.For example, the type of measurements allows for calculating higher-order correlation functions without requiring a change in the measurement protocol.This is in stark contrast to qubit-based VQEs, which are readout via the measurement of different Pauli strings.Their number quickly increases as soon as higherorder correlation functions should be accessed, which usually renders the experiment infeasible due to the, then necessary, exorbitant measurement budget.A further possibility that we have not yet exploited in the current work is the ability to perform simultaneous, two-tone pumping on the platform [15].This would allow for effectively combining the evolution generated by two resource Hamiltonians, without the need of approximating it with successive evolutions.Such technique could greatly expand the types of interactions available for future VQE implementations.Another promising path would be hybrid quantum platforms that include both qubits and bosons on the same hardware.In particular, these are an excellent fit for a lattice gauge theory simulation that involves fermionic matter fields.Given the highly tunable interactions available on this platform, it would also be suitable for VQEs of models outside of lattice gauge theories, including quantum chemistry and condensed matter physics.

A Introduction to VQEs
In the following we review the basic principles of variational quantum eigensolver (VQE) algorithms.For more details we refer to the review in Ref. [2] and the initial demonstrations in [87].In the work at hand, a VQE is employed to approximate the ground state of the target Hamiltonian ĤT , utilizing a closed feedback loop between the quantum processor and a classical optimizer.
The algorithm may be summarized as follows.First, the quantum device is initialized with an easily preparable state |ψ in ⟩.Next, a sequence of gates is applied to the initial state.The gates in this sequence are the unitaries exp −iθ k Ĥ(j) R , where Ĥ(j) R are the available resource Hamiltonians.Note that the set of available interactions does not need to be a universal gate set; a restricted gate set specific to the target Hamiltonian or the problem at hand is sufficient.The θ k are the variational parameters; their values are chosen by the classical optimizer and passed to the quantum device.In the circuit, the θ k manifest as the product of the interaction strength of Ĥ(j) R and the time for which the gate is applied.In a typical VQE circuit, a certain sequence of gates is often repeated with different variational parameters, where each elementary sequence of gates that is repeated is called a layer.
Once the VQE circuit has been applied to the initial state, the result is the VQE ansatz state Ψ( ⃗ θ) .For example, if the circuit employs r resource Hamiltonians in each layer and L layers in the circuit, then By tuning the the parameters ⃗ θ to their correct (but not necessary unique) values, this state approximates the the ground state of ĤT , provided the ansatz is expressive enough.. Once the ansatz state is prepared for a given set of variational parameters, the cost function Ψ( ⃗ θ) ĤT Ψ( ⃗ θ) is evaluated on the quantum device.Since quantum measurements are inherently stochastic, the preparation and measurement of Ψ( ⃗ θ) needs to be repeated in order to obtain a precise estimate of the cost function.In particular, if the state is measured M ≫ 1 times and σ 2 H = ⟨ Ĥ2 T ⟩ − ⟨ ĤT ⟩ 2 is the variance of the target Hamiltonian w.r.t.Ψ( ⃗ θ) , then the variance in the cost-function estimate is determined by σ 2 H /M .The estimate of the cost function is fed to a classical processor which employs an optimization algorithm to choose new values for the variational parameters, that are supposed to decrease the cost.These new parameters are fed back into the quantum device, where the process begins anew.This cycle is repeated until the convergence threshold for the classical optimizer is reached or the allocated measurement budget is exhausted.

B Elimination of the gauge fields
In the following, we express the Hamiltonian in Eq. (5) solely in terms of Higgs degrees of freedom.Generally, the complete elimination of the gauge field is always possible if only one spatial dimensions and open boundary conditions are considered.Such transformation results in an effective Hamiltonian that requires fewer modes to be simulated on the quantum hardware.
By considering the subspace of the theory that corresponds to a zero static charge on each vertex, Ĝn |Ψ physical ⟩ = 0, Eq. ( 7) can be rewritten as Ên = Ên−1 + Qn (30) for the physical states (see Sec. 2.

1).
Employing open boundary conditions, Eq. ( 30) can be solved recursively to obtain the electric field at every site resulting in This expression allows us to express the electric field operators entirely in terms of the charge operators.
In particular, the electric term of the Hamiltonian can be rewritten as In order to eliminate the gauge fields operators Ûn and Ên from Eq. (5), we can apply a residual gauge transformation analogously to the procedure presented in [5], which yields This induces the transformation, which in combination with Eq. (32) may be applied to Eq. ( 5) and subsequently results in the effective Hamiltonian given in Eq. (10).Here, the gauge degrees of freedom have been eliminated completely, which came at the cost of introducing long-range interactions between the Higgs fields.A method for obtaining the effective Hamiltonian in more than one spatial dimension is described in Refs.[50,81].

C Generation of ĤNN interactions in a parametric cavity
We now describe how the interaction term ĤNN in Eq. ( 16), which is one of the resource Hamiltonians used in the VQE, can be generated on the microwave platform.The derivation is analogous to the one presented in Ref. [51].Using a symmetric SQUID, the fourth-order term in the SQUID cosine potential is given by where âp (â n ) is the annihilation operator for the pump mode (cavity mode n) and g 0 is the intrinsic interaction strength between the pump and the cavity modes.The higher-order terms in the cosine potential are neglected.As a result, the full Hamiltonian of the superconducting cavity reads where ω p is the natural frequency of the pump and ω n is the natural frequency of mode n.We then transform Ĥ into the interaction picture w.r.t Ĥ′ 0 = ω p Np + n ω ′ n Nn , which yields the Hamiltonian The first term of Ĥint can be identified with the free rotation part of the unitary evolution given in Eq. ( 16), with Ω n = ω n − ω ′ n .For the remainder of this derivation we now focus on the second term in Eq. (37).The expansion of the fourth power factor yields a number of terms; as shown in Ref. [51], a suitable choice of the pump mode can selectively enhance the desired interaction among the terms available in the expansion.We are interested in generating interactions between photon number operators of different modes, which are terms that are time independent in the expansion.By choosing ω p much smaller than any other frequency that appears in the fourth power expansion, we can use the rotating wave approximation to neglect all terms that rotate at frequency larger than ω p .We assume that ω p is so small that it can be effectively considered zero for the duration of the experiment.Let us now calculate the effective Hamiltonian explicitly.There are two terms in the expansion that are time independent: a product of {â † i , âi , â † j , âj }, and a product of {â † i , âi , â † i , âi }.Using combinatorics and the bosonic commutation relations, the first type of terms sum to Combining Eqs.(38)(39) together gives Eq. (16) [the constant term is ignored since it contributes an overall phase to Eq. ( 16)].Finally, we assume that the pump has a strong coherent tone such that we can apply the parametric approximation and substitute âp with its classical amplitude |α|e iϕ .With this approximation the time evolution induced by Ĥint is exactly the one given in Eq. ( 16) with g′ = 2g 0 |α| cos ϕ.

D Truncation effects in the HOBM mapping
In this section, we consider the effects of using a truncated HOBM mapping (given in Sec.2.3), for our VQE simulations.In the following we omit the lattice site index as these expressions hold for all sites.Using the untruncated HOBM model, all of the commutation relations remain the same (see Sec. 2.1 for the commutation relations of the original model), except for The correct commutation relation is [ φ, φ † ] = 0, and it is recovered as N 0 → ∞.Let â(k) and N (k) be the truncated operators of size 2k+1.When we consider the truncated HOBM model (as is the case for our classical simulation) an error is introduced, and the following commutation relations can be derived To recover the correct commutation relation in Eq. (41a), we set k = N 0 , which gives Eq. (42a) is the same commutation relation as in the full Hilbert space [Eq.(1)] and, given that [ φ, φ † ] = 0 in the full Hilbert space, the truncation error in Eq. (42b) goes to zero as N 0 → ∞.In general, Eq. (42b) shows that the state with maximum photon number, |2N 0 ⟩ is most affected by truncation errors.

E Higgs mechanism and phase structure in the absence of a topological term
In this appendix we briefly review the Brout-Englert-Higgs mechanism for the continuum model and discuss its implications for the phase structure.
To this end let us start from the continuum Lagrangian where φ is a classical complex scalar field, D µ = ∂ µ + igA µ the covariant derivative with the gauge field A µ and the coupling strength g, F µν = ∂ µ A ν − ∂ ν A µ the field strength tensor and V (|φ|) = −m 2 |φ| 2 + λ 2 |φ| 4 with λ > 0. The first term of the action describes the kinetic energy of the scalar field and the coupling to the gauge field, the second term the kinetic energy of the gauge field, and the potential V (|φ|) contains the mass term and the self-interaction of the scalar field.It is straightforward to see that the action in Eq. ( 43) is invariant under U(1) gauge transformations given by φ where α(x) is a real, differentiable function.
To get an intuition about the physics of the model, it is instructive to derive a semiclassical picture.To this end, we assume that the potential term is dominant and in the ground state the field φ fluctuates only slightly around the vacuum expectation value φ 0 minimizing the potential V (|φ|).Rewriting the potential as we see that V (|φ|) is quadratic function of |φ| 2 and we can easily read off the minimum.Taking into account that |φ| 2 is a positive semi-definite quantity, we have to distinguish two cases depending on the sign of m 2 .(i) For m 2 ≤ 0 the potential is a parabola in |φ| 2 (see also lower left inset of Fig. 6) with a unique minimum at φ 0 = 0.For φ ≈ 0 Eq. ( 43) reduces to which is nothing but a pure gauge theory describing a massless photon and showing charge confinement [73].(ii) For m 2 > 0 the potential has the shape of a "Mexican hat" (see also upper right inset of Fig. 6) with the minima forming a level set given by In the expression above, we have defined the quantity v = 2m 2 /λ for convenience.To account for quantum fluctuations around φ 0 , we parameterize the field using polar representation where the real fields h(x) and ϕ(x) describe the fluctuations of the length of the field in radial direction and the phase.Combining this expression with Eq. (43) and considering only at the gauge part L gauge of the resulting Lagrangian we find showing that the photon is massive with mass m p = gv.Moreover, in that case the U(1) symmetry is spontaneously broken as the vacuum of the theory corresponds to a single one of the minima described by Eq. (44).Sending λ → ∞ while keeping the ratio m 2 /λ fixed, the radial fluctuations can be neglected and φ has a fixed length of v/ √ 2. Inserting this expression into Eq.(43) one finds the effective Lagrangian From Eq. ( 46) one can derive the Hamiltonian in Eq. (5) as shown in Ref. [39].
The simple semiclassical picture above shows that the model has two distinct regions.In the Higgs region the U(1) symmetry of the theory is spontaneously broken and the photon acquires a mass.In contrast, the confining region corresponds to a pure gauge theory describing a massless photon and the U(1) symmetry is intact.Going beyond this simple semiclassical picture and solving the the lattice discretization of Eq. (43) numerically using MCMC methods one finds that the intuition from the rather simple semiclassical picture also holds true more generally [38,73] and one obtains the phase diagram shown in Fig. 6 in the main text.

F Periodicity and Symmetries of the Hamiltonian
Here we briefly show that physics is periodic in ε 0 with period 1 and discuss the symmetries of the lattice Hamiltonian.For simplicity we work with the formulation in Eq. (5), in which the gauge field has not been integrated out yet.
To show that physics is periodic in ε 0 , we consider the transformation Let us first focus on k = 0.In that case Eq. (47) corresponds to a charge conjugation transformation Ĉn which exchanges particles and antiparticles.This transformation is unitary, and, in particular, we see that applying it twice we get the initial operators back, thus showing that Ĉ2 n = 1 and charge conjugation is a Z 2 symmetry.
To get further insight into the case k ̸ = 0, we look at the commutation relation in Eq. (3), from which follows that the unitary operator Ûn introduces integer shifts in Ên : Hence, by combining charge conjugation with Ûn to Ĉn ], we obtain an additional shift in the electric field by k positive (negative) units.Since both Ĉn and Ûn are unitary, Eq. (47) indeed describes a unitary transformation.Applying this unitary transformation to Gauss's law in Eq. (30) and the Hamiltonian in Eq. (5), we see that all terms are invariant except for the electric field energy which transforms according to Looking at that equation, we make the following observations.(i) For any value of ε 0 we can find a k such that the Hamiltonian is mapped to a unitarily equivalent one with k − ε 0 ∈ [0, 1).Consequently, physics is periodic in ε 0 with period 1, and we can restrict ourselves to ε 0 ∈ [0, 1) without loss of generality.(ii) For k = 0, ε 0 = 0 and k = 1, ε 0 = 1/2 the transformation from Eq. ( 47) is a symmetry of the Hamiltonian.In particular, for ε 0 = 1/2 there are two field configurations yielding the same electric energy.Thus, for R 2 → 0 the ground state of the Hamiltonian is doubly degenerate and has a Z 2 symmetry.This symmetry is preserved for nonvanishing R 2 along the critical line, before it is eventually spontaneously broken upon reaching the critical value R 2 c .

G Location of the first phase transition
In this appendix we derive Eq. (25) from Sec. 4.2.This formula gives the smallest, positive value of ε 0 where a phase transition occurs and the dominant term in the ground state after the phase transition.Assuming that R 2 is small such that the kinetic term in Eq. (10) can be ignored, the Hamiltonian (including the penalty term for vanishing total charge) be- Looking at Eq. (48), we see that the eigenstates of Ĥ in this approximation are given by charge eigenstates.Before the transition occurs, meaning in the regime of small ε 0 and R 2 , we expect the ground state to be dominated by the state with vanishing charge everywhere, |0⟩ ⊗N .The corresponding energy is given by (N −1)( Since we work in the subsector of vanishing total charge and for R 2 ≪ 1 the n c n Q2 n term in Eq. ( 48) strongly penalizes states with many nonzero charges, we expect the ground state after the phase transition to be of the form |ψ⟩ = |0⟩ ⊗n 1 −1 |k⟩ |0⟩ ⊗n 2 −n 1 −1 |−k⟩ |0⟩ ⊗N −n 2 with n 2 > n 1 and k ∈ Z.The energy of this state is given by where r = n 2 − n 1 .Since the ground-state energy is a continuous function of ε 0 , we can equate this expression to the energy before the transition which gives Since ε 0 > 0, that forces k < 0. So, for each value of k and r, this formula predicts where the first phase transition could occur.Of course, of all the possible choices of k and r, the one that predicts the smallest value for the phase transition point ε 0 gives the correct formula.Since k ̸ = 0 (no phase transition would occur in that case), we therefore choose k = −1 and r = N − 1.Thus, we predict the first phase transition to cause the ground state to transform from |0⟩ ⊗N to |−1⟩ |0⟩ ⊗N −2 |+1⟩ at the point

H General VQE protocol
In this appendix we provide the VQE strategy used to obtain Fig. 5. Approximating the ground state of the model in the vicinity of the critical point R 2 ≈ R 2 c requires a more complex and deeper variational circuit than the one used for detecting the first-order phase transition for smaller values of R 2 in Fig. 3. Furthermore, the number of variational parameters quickly becomes unfeasible for the type of Bayesian optimizer we considered [70,72], therefore we need to adopt a strategy where only one part of the circuit is optimized at a time.
The basic element of the circuit is a combination of beam splitters connecting all possible sites, alternated with the interactions in ĤNN [see Eq. (16)] where Ĉ(i,j) is a generalization of the operator given in Eq. ( 21), obtained by choosing N l = 1, and using the beam-splitter interaction acting on modes (i, j).The unitary Dk is then repeated for a suitable number of layers such that the number of variational parameters is feasible for the chosen optimization algorithm to obtain the block The optimization proceeds as follows.To the initial state |1223⟩ we apply the unitary B. For this block, both the rotations ⃗ Ω and the set of { ⃗ θ k } are considered as variational parameters.After the convergence of this first optimization we iteratively apply further unitaries B to the previously obtained state and optimize the parameters { ⃗ θ k } associated to each new block.The rotations ⃗ Ω are chosen by going into a suitable frame of reference and therefore have to be consistent throughout the variational circuit, hence they are only considered in the first block.To obtain the results show in Fig. 5 we found that the optimization in the case of the VQE without including statistical noise was successful by choosing N l = 2 for the first block and N l = 3 for each subsequent block.In total, we employed five blocks to obtain the results shown in the main text.

FIG. 1 :
FIG.1: Schematics of the VQE.The quantum processing unit (QPU) consists of a microwave cavity terminated by a SQUID.The variational parameters ⃗ θ determine how the SQUID is pumped and consequently which interactions are generated between the cavity modes to generate the ansatz state ψ( ⃗ θ) .A quantum-limited amplifier allows for performing heterodyne measurements on the cavity modes to obtain the cost function C( ⃗ θ).In our case the cost function is the expectation value of the Hamiltonian for which we want to find the ground state, as defined in Eq.(18).C( ⃗ θ) is used by the classical processing unit (CPU) to find the optimal variational parameters that minimizes it.

σ 2 HM
. The intrinsic variance of the indirect measurement σ 2

FIG. 4 :
FIG.4: Cost function (expectation value of the Hamiltonian) evaluations in an successful optimization for ε 0 = 1.9.As a guide for the eye, minimum values of the cost function are represented by bigger dots.The solid (dashed) horizontal line is the ground state (first excited state) energy.Inset: fidelity of the variational state with the ground state obtained from exact diagonalization.

FIG. 5 :
FIG.5: EFD of the ground state for different values of ε 0 , R 2 = 1 and β = 1.The ground state is calculated with a classical simulation of a VQE.In this parameter regime the quantum phase transition disappears, as discussed in more detail in Sec.4.2.Data points correspond to VQE simulations with increasing number of repeated measurements for one cost function evaluation.Blue triangles: M = 10 6 , orange squares: M = 10 9 , red circles: M → ∞ (i.e.statistical errors are neglected).

FIG. 6 :
FIG. 6: Phase diagram of the Abelian Higgs model with fixed length of the field in the R 2β plane.The Higgs region is indicated in blue whereas the confinement region is indicated in orange.The dashed line corresponds to a crossover between the Higgs and the confining region.The insets show a sketch of the form of the potential in the different regions.

FIG. 8 :
FIG.8: Ground state energy (left column) and EFD (right column) as a function of ε 0 calculated using MPS for β = 1.0 and R 2 = 0.3 (first row), 0.6 (second row) and 1.0 (third row).Dots indicate s = 1, triangles s = 2 and the different colors encode the different system sizes N = 4 (blue), 8 (orange), 20 (green) and 100 (red).As a guide for the eye the markers are connected with lines.The error bars represent the uncertainty from the extrapolation in χ.

FIG. 9 :
FIG.9: Location of the first discontinuity extracted from our numerical data as a function of system size for β = 1.0 and R 2 = 0.3 (a), and 0.6 (b).Dots indicate s = 1, triangles s = 2, and the solid line the prediction from Eq. (25) for the large mass limit.The error bars originate from our finite resolution in ε 0 .
).The resulting local Hilbert space is finite and has dimension d = 2s + 1.Using this mapping, the commutation relations of the original model stay intact (see Sec. 2.1 for more details on the commutation relations) except for