Two-dimensional hard-core Bose–Hubbard model with superconducting qubits

The pursuit of superconducting-based quantum computers has advanced the fabrication of and experimentation with custom lattices of qubits and resonators. Here, we describe a roadmap to use present experimental capabilities to simulate an interacting many-body system of bosons and measure quantities that are exponentially difficult to calculate numerically. We focus on the two-dimensional hard-core Bose–Hubbard model implemented as an array of floating transmon qubits. We describe a control scheme for such a lattice that can perform individual qubit readout and show how the scheme enables the preparation of a highly excited many-body state, in contrast with atomic implementations restricted to the ground state or thermal equilibrium. We discuss what observables could be accessed and how they could be used to better understand the properties of many-body systems, including the observation of the transition of eigenstate entanglement entropy scaling from area-law behavior to volume-law behavior.


INTRODUCTION
Analog quantum simulators have evolved in the last two decades from a theoretical concept to an experimental reality (see e.g., refs [1][2][3]. Initial experimental success was predominantly achieved with atomic systems, including neutral gases and trapped ions [4][5][6][7][8] . More recently, superconducting circuits have emerged as a viable quantum simulation platform [9][10][11][12][13] . This modality-based on "artificial atoms"-features a high degree of experimental controllability and stability 14 . The flexibility of the superconducting platform has enabled several successful quantum simulation experiments [15][16][17][18][19][20] . Here, we show how to realize the two-dimensional (2D) hardcore Bose-Hubbard model (HCB) illustrated in Fig. 1 using an array of transmon qubits 21 , the current workhorse qubit design in superconducting circuits. The HCB is a strongly interacting system that displays some of the critical properties of interacting quantum systems, including the area-law to volume-law transition of the entanglement spectrum that has been extensively studied in many-body systems 22 . Outside of one dimension (1D), this system has no known analytical solution, and its study has been conducted mostly through numerical methods limited in their scope. The most successful approach has been the use of tensornetwork methods, which focus on finding the ground state energy 23,24 . An experimental realization of a 2D HCB could offer new and complementary insights about the eigenstates and dynamics of many-body systems. It could also be used to validate the results of tensor-network methods in large systems, and test their underlying assumptions on the nature of many-body wavefunctions. An experimental realization also offers access to the system's entire spectrum, allowing one to measure the manybody properties of its excited states.
Previous experiments have realized the HCB in 1D 15,25 , where the model can be solved by analytical methods 26,27 and has the dynamics of a free fermion gas. Recent realizations have also explored entanglement propagation in ladders and a 3 × 7 array 20 . Here, we propose the implementation of the 2D HCB with state-of-the-art transmon qubits. We calculate the requirements on qubit uniformity and lifetime, and describe the control systems required to measure the array's many-body properties. Finally, we propose a technique to generate highly excited states that enable one to more completely explore the system's spectrum and observe its many-body properties.
Superconducting quantum many-body physics simulator We consider the implementation of a quantum many-body physics simulator (QMBS) with a superconducting quantum circuit made up of multiple repetitions of small basic circuits implementing qubit and coupling elements. Note that while here and throughout the paper we take each site to be a qubit, i.e., a twolevel system, the discussion here applies equally to systems of spins or particles on a lattice. We describe the system with the Hamiltonian The behavior of the QMBS depends on the ratio of J to the other three scales. At J ≪ Δω, exchange of energy between different qubits is suppressed, and the system will behave as a collection of uncoupled circuits. In this case, there is no many-body physics to speak of, and the system decomposes into multiple systems with a single degree of freedom each. Thus J ≳ Δω is required for manybody dynamics to appear in the lattice. The ratios J= A j j and J/ω q then determine which states are effectively coupled and thereby which theoretical models are accessible 28 . These features are collected in the form of a phase diagram in Fig. 2 and discussed in further detail below. This describes systems such as quantum rotor models in the paramagnetic phase. We note that generally, for a large number of sites, these have a large density of states, and it may be difficult to prepare the system in a low-temperature quantum state. In that case, the system can be understood by a semiclassical description, and it is hard to observe uniquely quantum dynamics. Such experiments have been performed for large numbers of Josephson junctions 29,30 .
We note that this regime can be avoided by choosing a different basis of states to describe the Hamiltonian, i.e., by switching the choice of which circuits describe the qubitsĤ Q i and the coupling elementsĤ J i;j .
• Particle-like models: In the central portion of the phase diagram, the hierarchy of scales is Here, the rotating wave approximation is valid, and the coupling elements can move an excitation between sites but will not change the total number of excitations. This regime is equivalent to models of bosonic particles, and we may describe the system with the Bose-Hubbard Hamiltonian, whereâ y i is the creation operator for site i andn i ¼â y iâ i is its energy level. Here, the anharmonicity plays the role of the onsite interaction strength, while inter-qubit coupling generates transverse hopping terms (J xx ) and longitudinal interaction terms (J zz ).
The sub-regime where J < A j j ( ω q -the working point of the transmon qubit 21 -is the most experimentally accessible parameter regime and is widely adopted by the superconducting circuit community in a multitude of experiments (see e.g., ref. 17 ), including recent implementations of 1D Bose-Hubbard lattices 25,31 . In this manuscript, we focus on this regime.
We note that a subset of the particle-like regime, where Δω J, can be used to simulate disordered systems. This can be achieved either by intentionally varying the qubit frequency across the lattice, or by decreasing the hopping energy at a constant residual disorder.  Fig. 1 The two-dimensional hard-core Bose-Hubbard model (HCB). a Sktech of a sample 4 × 4 HCB lattice. Each circle represents a qubit, constrained to two energy levels with energy difference ω q . The diamonds represent coupling between each pair of nearest neighbors at strength J. One magnified version of each element shows these energies. b The spectrum of the same system, here with J = ω q /10. On the left we show the entire spectrum; for a system with N qubits, it is composed of N + 1 sectors defined by the total excitation number n. On the right, we show a close up of a particular sector, with an energy bandwidth ΔE / 8n 1 À n=N ð Þ J.
Frequency variance Here, e j i i ( f j i i ) is the state with qubit i in its first (second) excited state and all others in the ground state. X i (X i;j ) denote the average of X i (X i,j ) over all qubits (all coupled pairs). We take ħ = 1 and the ground state energy to be zero. Fig. 2 Accessible models with a QMBS based on a superconducting circuit. Which models can be realized depending on the ratio of the coupling strength J to the qubit frequency ω q and anharmonicity A. The diagonal cyan shading highlights the area accessible with transmon qubits, which are our focus here. The coupling must be greater than the frequency spread, J ≳ Δω, for any kind of manybody physics to appear. Where Δω ≲ J ≪ ω q , the behavior is particlelike and we expect to see a version of the Bose-Hubbard model with A playing the role of on-site interaction. Where ω q tJ ( A j j we have a spin-like model, where each unit acts as a two-level qubit while the coupling does not conserve excitation number. When J dominates all other scales, we expect semiclassical behavior.
Here, the anharmonicity dominates the coupling term, ensuring that each unit cell remains within the qubit manifold. However, the coupling elements are strong enough to change the qubit state in a non-excitation-conserving way. The rotating wave approximation then breaks down, and the system is best understood by a spin-like model, whereσ μ i are the Pauli operators on site i. This regime, where the coupling strength becomes similar to the transition frequencies of the coupled systems, is known as the ultra-strong or deep-strong coupling regime 32,33 and it is more challenging to realize experimentally. However, superconducting artificial atoms are more suitable for its realization than natural atoms coupled to an electromagnetic cavity, as their coupling strength to a harmonic oscillator mode is not necessarily limited by the fine structure constant 34,35 . In general, physical couplings in the deepstrong coupling regime can be achieved with strongly nonlinear qubits and high-impedance circuits 35 . A promising qubit modality to reach such high couplings is the flux qubit, where ω q ( A j j, as demonstrated experimentally 36,37 . The fluxonium qubit 38 , an extension of the flux qubit, has recently been demonstrated to preserve long coherence times while in the high anharmonicity regime 39 .

RESULTS
For the remainder of this article, we focus our attention on the regime This combines the two constraints mentioned in our analysis of the possible working regimes: the system operated with these parameters both conserves the number of excitations and remains within the qubit manifold. This is a bosonic model, where each site can be either empty or occupied by a single particle. The system is then described by the effective Hamiltonian whereσ z i ;σ ± i ¼σ x i ± iσ y i are the Pauli z, raising and lowering operators on site i.
As the Hamiltonian is number preserving, its spectrum decomposes into N + 1 distinct sectors defined by the total excitation number n. Each sector is composed of N n levels, defined by their rotating-frame energy ϵ, with bandwidth proportional to J. The eigenstates of Eq. (7) are then given by where E G is the ground state energy. This spectrum is sketched out in Fig. 3. The HCB is difficult to solve except in some specific cases. The 1D chain can be solved through fermionization 26,27 , and the case of n/N ≪ 1 (n/N ≈ 1) can be understood analytically by perturbative corrections to the free particle (free hole) problem 40 ; both regimes exhibit noninteracting behavior that is much simpler than what we describe below. In addition, small systems can be exactly diagonalized, as we do here for a 4 × 4 lattice. Beyond these limits, research into the model has generally used tensor-network methods and focused on the ground state energy 23,24 . An experimental realization of a 2D version of Eq. (7) can therefore contribute significantly to our understanding of the eigenstates and dynamics of many-body systems, and also the validity and limits of tensor-network methods in large systems. Beyond this, as we discuss below, an experimental realization can access the system's entire spectrum.
Many-body physics in the hard-core Bose-Hubbard model We consider two particular measures of the system's many-body spectrum: the correlation length and the behavior of entanglement entropy for each eigenstate. In Fig. 4, we show these quantities exhibit transitions along the spectra within each sector: as we go from the edges of the band to the center, the correlation length grows from finite to infinite, and the entanglement entropy of subsystems evolves from obeying an area-law dependence on the subsystem's size to a volume-law dependence.
• Correlation length: The typical scale beyond which different sites are no longer correlated serves as an order parameter for phases with long-range order 41 . The correlation length is a limiting factor for the applicability of tensor-network methods, which can be used only when correlations are finite 42 . Having experimental access to the correlation length therefore provides significant insight into the many-body properties of the system.
For our purpose, we define the correlation length in terms of the correlation function We then extract the correlation length of a state n; ϵ j i by fitting C x i;j ψ j i ð Þ to the form jC x i;j ðjn; ϵiÞj 2 ' Aðn; ϵÞ exp½Àj r i ! À r ! j j=ξðn; ϵÞ; over all pairs i, j of nearest neighbors and next-nearest neighbors. Here j r i ! À r ! j j ¼ jx i À x j j þ jy i À y j j is the Manhattan distance between the sites i, j.
We plot the correlation length ξ as a function of eigenstate energy in Fig. 4a for a 4 × 4 lattice near halffilling, where we expect many-body effects to dominate. As discussed above, we observe it goes from finite and short for states at the edge of the band to effectively infinite for states at its center. Fig. 3 The spectrum of the hard-core Bose-Hubbard model in the rotating frame. Here, calculated for a single realization of a 4 × 4 square lattice with nearest-neighbor hopping ( Fig. 1) at Δω = 0.2J. The full spectrum comprises 17 distinct sectors with fixed n, separated by ω q with width proportional to J. ϵ is the rotating-frame energy, as defined in Eq. (9).
• Entanglement entropy: For a state with density matrixρ, the entanglement entropy of some subset X of the lattice is the entropy generated when it is severed from the rest of the system, where Sσ ð Þ is the entropy ofσ, andρ X ¼ Tr 8i= 2Xρ ,ρ X ¼ Tr 8i2Xρ are the reduced density matrices of the subsystem X and the remainder of the lattice, respectively. Note that if the initial density matrixρ ¼ ψ j i ψ h j is a pure state, Sρ ð Þ vanishes while the entropy of both subsystems must be identical, so that Throughout this paper, for the purpose of numerical calculations, we use the second Rényi entropy, The entanglement entropy is a measure of entanglement between different parts of the lattice, and has been an important tool in the study of many-body systems. In particular, there has been significant study of the difference between states where it is proportional to the size of the subsystem X ("volume-law") and where it is proportional to the size of its boundary ("area-law") 22 . Volume-law states are also harder to approximate using tensornetwork methods.
To describe the growth law for an eigenstate n; ϵ j i, we extract the parameters s V and s A by fitting the entanglement entropy to the form over different lattice subsets X. Here V X is the number of sites in X (its "volume") and A X is the number of coupling terms between sites in X and the rest of the lattice (its "area"). The fit parameters can then be understood as s V Bulk entanglement entropy per site; s A Boundary entanglement entropy per bond: Thus, the ratio s V /s A determines whether the entanglement entropy obeys an area-law-like or volume-law-like behavior. We plot this quantity for a 4 × 4 lattice in Fig. 4b. We see the transition from area-law behavior for states at the edges of the band to volume-law behavior at its center. We also see little variation in this behavior between different sectors with similar n/ N. This allows us to explore the behavior of the entanglement entropy by preparing coherent-like superposition states across multiple sectors, as described below.
Global measures such as the entanglement entropy are key to understanding many-body properties, but observing them in the lab poses experimental challenges. Naively, the entropy of a state is derived from the density matrixρ and extracting it requires full state tomography. The challenge here is two-fold: first, the number of measurements scales exponentially as 2 2N 43 ; and second, one must have sufficient control to apply any combination of rotationsσ ± i to all sites concurrently. The situation, however, is not quite so dire. Multiple recent proposals have suggested alternative approaches for measuring non-local observables such as n-time correlation functions 44 and the second Rényi entropy [45][46][47][48] . These proposals substitute random unitaries for the full set of rotations mentioned above, easing the control requirements. They also require fewer unitaries than does full state tomography, though the number of measurements needed still scales exponentially with system size. We note, though, that even as the total size of the system increases, the scaling coefficients s V , s A can be determined from the entanglement entropy of fixed-size subsystems (e.g., a block of sites of size 3 × 3 and all its subsystems), leaving the required number of measurements constant even if we increase N.
As noted above, the emergence of many-body behavior requires relatively uniform qubit frequency, Δω ≲ J. In Fig. 4c, we quantify the tolerable amount of variation for the metrics discussed here. We do so by calculating the behavior of the entanglement entropy at the center of the band and at its edge at varying disorder strength, averaged over multiple realizations of the lattice. We find that up to Δω ≈ 0.5J, one can observe distinctly different physics in different parts of the spectrum. At larger frequency disorder, the variation between lattice realizations dominates this effect.
Prospects for transmon implementation The transmon qubit 21 is a natural building block for the implementation of the HCB with a superconducting circuit. It Numerical evidence for many-body behavior in the HCB. a, b Calculated for a single realization of a 4 × 4 square lattice with nearestneighbor hopping and Δω = 0.2J. We expect the physics of the sectors n = 6, 7, 8, near n ≈ N/2, to be dominated by many-body effects. For each eigenstate in these sectors, we calculate and plot a the C x correlation length [Eq. (11)] and b the ratio s V /s A between the volume coefficient and area coefficient of the entanglement entropy [Eq. (15)]. We observe a clear variation in physics along the spectrum, going from a finite correlation length and area-law behavior of the entanglement entropy at the edges of the band to diverging correlation length and volume-law behavior at the center of the band (Note that in a 4 × 4 system the largest separation between qubits is L = 6, and so ξ ≳ 6 hints at long-distance order). We also observe that the behavior is similar for the three sectors with similar n/N. The robustness of this signature across filling number enables us to probe entanglement entropy using coherent-like states. c For the sector n = 8, we examine the effects of increased frequency variation, Δω. We plot the ratio between s V /s A at the center of the band (ϵ = 0) and its edge (ϵ = 10J), averaged over 10 realization of the disorder. The shaded area gives the range of results over one standard deviation. For Δω ≤ 0.5J, we can clearly observe the change in physics over the spectrum; for Δω > 0.5J, the variance due to different realizations dominates.
behaves as a weakly non-linear oscillator with a fundamental transition frequency in the range of ω q /2π~5 GHz. Each lattice site is represented by a single transmon qubit, with the local site energy corresponding to the qubit transition frequency ω q .
The anharmonicity of the transmon qubit is negative, typically in the range of A/2π~−250 MHz or~5% of its frequency 21 . The self-Kerr non-linearity of the transmon Hamiltonian maps directly onto the on-site interaction term in the Bose-Hubbard model 25 .
As the hard-core Bose-Hubbard model operates in a regime, where J= A j j ( 1 (Mott insulator phase), the population of the same lattice site with two or more particles is strongly suppressed due to the presence of the self-Kerr term, irrespective of its sign. One may note that for large enough lattices, the kinetic energy may reach the scale of the anharmonicity ϵ / NJ $ A j j. Generally, this effect can be treated as a perturbative correction to the hardcore approximation, as we expect to see only a small number of sites out of a large occupation ∝ N in the forbidden state.
It is straightforward to connect transmon qubits via capacitive coupling 25 , leading to the hopping term in the Bose-Hubbard model with nearest-neighbor coupling energy J. Typical achievable coupling strengths are tens of megahertz, rendering the qubit-qubit interaction well within the strong coupling regime J > Γ, where Γ denotes the qubit decoherence rate. Contemporary transmon qubits feature reproducible coherence times in the range of 20 to 100 μs 17 , corresponding to Γ/2π ≲ 10 kHz.
As discussed above, an experimental implementation operating in the regime J= A j j ( 1 suppresses transitions to the second and higher levels, and implements the HCB. In order to observe manybody physics, we generally require the qubit lifetime to be much longer than the characteristic time scale for information to traverse the system, 1/Γ ≫ L/J, where L is the number of hops to go across the system (its length). With five orders-of-magnitude in separation, A j j \ 10 5 Γ, this is easily achievable with transmon lattices of 100 qubits or more. Generally the sweet spot in this case MHz. Qubit coherence in the proposed transmon HCB lattice is expected to be at the level of individual state-of-the-art transmon qubits 17 , limited by a combination of material defects 49 and parasitic coupling to stray modes in the sample package 50 . Scaling to a larger number of qubits typically requires a chip and sample package of larger dimensions, with the risk of introducing additional parasitic modes at frequencies at or close to the qubit frequencies and therefore impairing qubit performance. In previous implementations of arrays with 24 and 53 qubits, energy relaxation times averaged at around 10 μs 18 and 15 μs 19 .
Fabrication variations translate to variations in transmon transitions which may exceed 200 MHz 51 , yielding disorder in the emulated model on the order of Δω $ A j j. To compensate for such variation, we consider a lattice of frequency-tunable transmon qubits. This is achieved by replacing the single Josephson junction of the qubit with a dc-SQUID, facilitating a frequency tunability of several GHz. In an experiment, this enables one to tune the individual qubit frequencies mutually on resonance (to within their spectral linewidth 25 ).
Individual frequency control requires N slow (dc) control lines for flux biasing each of N qubits. Such low-frequency wiring can be straightforwardly routed in dilution refrigerators and connected to the sample package in large numbers, as the necessary connectors are compact and bulky attenuation at multiple temperature stages is not required.
Frequency variation in the lattice is mitigated experimentally by calibrating the (dc) flux crosstalk matrix, containing information about the frequency shift of qubit i responding to a flux bias applied to bias line j (1 ≤ i, j ≤ N). In large lattice implementations, qubits are physically located far away from flux bias lines of other qubits. By taking into account only nearest-neighbor and nextnearest neighbor parasitic flux coupling, the resulting flux crosstalk matrix is sparse, reducing the number of matrix elements from O N 2 À Á to O N ð Þ. In general, flux crosstalk calibration requires the measurement of sections of all N qubit spectra while consecutively biasing each of the N flux control lines. As the spectra can be measured simultaneously with multiplexed readout, this requires O N ð Þ individual measurement scans and therefore scales linearly with lattice size.
In addition, dynamic (ac) flux control allows for rapid frequency tuning of the qubits. By detuning a qubit away from its neighbors, we can effectively decouple it from the lattice. For example, in a square lattice, system dynamics can be entirely frozen out, enabling state preparation and readout, by detuning every other qubit in a checkerboard pattern, where all "white" qubits remain at the original frequency and all "black" qubits are shifted.
This scenario requires N/2 qubits to be equipped with fast flux lines, such that even a large lattice of size 10 × 10 requires only 50 flux control lines. Assuming individual bias lines used, enabling full control on each qubit, the number of required coaxial lines is still moderate compared with recent implementations using 50 and 200 coaxial control lines for a 24-qubit and 50-qubit chip, respectively 18,19 .
Implementation with floating transmon qubits Figure 5a shows a possible circuit implementation of the 2D HCB based on a grid of transmon qubits each consisting of two floating electrodes [52][53][54][55] , in contrast to recent realizations where one of the electrodes is grounded (e.g., Xmon qubits 56 ). In circuit designs with an increasing number of qubits, circuit elements can be proximal or even overlap when using crossover fabrication techniques 57 . This may result in unwanted spurious coupling, referred to as crosstalk. Such spurious couplings can exist between signal lines, readout resonators, and qubits. In order to minimize this effect, it is advantageous to confine electric fields by decreasing the mode volume; this however comes at the expense of an increased electric field strength, leading to enhanced surface defect loss 49,58 . The floating layout is advantageous since it suppresses parasitic couplings in the circuit.
To see this, we compare the parasitic coupling of a resonator mode to a floating and a grounded transmon. We assume a (parasitic) capacitive coupling C P ; C 0 P between the resonator and the electrodes of the floating transmon (circuit diagram in Fig. 5b), or C P to the electrode of the grounded transmon (Fig. 5c). While the coupling capacitance for the grounded transmon is simply the effective coupling capacitance between resonator and the floating transmon depends on the parasitic capacitances C P ; C 0 P as well as the capacitance to the ground C G . Assuming without loss of generality C 0 P C P , circuit analysis (see the Methods section) yields an effective coupling capacitance We note C  Fig. 5d for typical parameters C P ; C 0 P ; C G . The argument remains valid if the capacitances of the two transmon electrodes to ground are not identical.
Another potential benefit of the floating transmon design is that it provides a tuning knob for the strength of long-range interactions within the lattice. In particular, the coupling range between non-adjacent qubits can be adjusted by controlling C G / C sh (Braumüller et al., unpublished), the ratio of the qubit shunt capacitance to the capacitance of the two pads to the ground. For the implementation of the HCB Hamiltonian, Eq. (7), next-nearest neighbor couplings must be suppressed, which can be achieved in the limit where C G ≫ C sh , but the use of floating transmons opens the possibility of exploring models with non-local interactions in the future.
Another strategy to mitigate unwanted crosstalk is to physically separate circuit elements by introducing a multi-layer chip layout (3D integration) 57 . This approach is particularly beneficial in the implementation of a 2D grid of qubits, as the circuit topology prevents in-plane access to interior qubits. In a planar layout, this can be resolved by using airbridges to crossover signal lines 59 , but these are naturally prone to unwanted crosstalk. The 3D integration approach allows for a separation of coherent elements (qubits) on one layer and signal lines on another layer, with their respective electric fields well separated. Couplings between qubit and control or readout lines are achieved via a flip-chip approach and connectivity to the other substrate surface is facilitated by through-silicon vias (TSV), which are low-loss superconducting trenches etched inside the silicon substrate 57,60 .
Qubit readout and control Individual qubit readout and control in devices with only few qubits can be achieved by connecting a separate signal line to each qubit. For a QMBS-style device with a large number of qubits, this approach is limited by the available number of signal lines as well as by geometric constraints. Instead, efficient multiplexed readout can be performed by coupling multiple qubits to a single signal line through individual dispersive readout resonators with frequencies spaced at intervals large compared to their linewidths 61,62 . We sketch out an example of this setup in Fig. 5a.
As implied by the color coding in Fig. 5a, signal lines must cross qubit pads or qubit coupling elements in a planar circuit implementation in order to reach qubits inside the lattice, leading to experimental challenges. A possible strategy to address this issue is the use of 3D integration techniques 57 .
The particulars of dispersive readout for individual qubits are well established. The challenge in reading out a large, degenerate array of qubits is the interplay between measurement and the ongoing dynamics. To get a snapshot of the system at a particular time, we must generally measure the qubits on a time scale T meas ≪ 1/J, or else freeze the dynamics.
For a homodyne measurement, typical measurement time scales as 63 where κ is the resonator linewidth, χ its dispersive shift between qubit states g, e; and n is the mean number of photons in the cavity during readout. There are two limiting factors for this readout speed: cavity occupation must remain below the critical photon number in order to ensure to operate in the linear dispersive regime, and the induced Purcell decay of the qubit, γ P , must remain small 64,65 , Here η PF accounts for protection resulting from a Purcell filter 65 , and L is the maximum distance between any two qubits. We have taken the anharmonicity A j j to be much smaller than the qubitresonator detuning. If we keep fixed then the measurement time is We find that the operating regime for measurement depends on the ratio J= A j j, favoring a smaller J. The design is limited, however by the factors we have previously discussed. The hopping energy must exceed the frequency disorder across different qubits, J ≳ Δω; with individual (dc) flux bias lines, one can reasonably expect to achieve Δω/2π ≈ 100 kHz, independent of lattice size. The rate must also be fast enough to allow information to travel across the entire system before decoherence kicks in; at a conservative qubit lifetime of T 1 ≈ 10 μs, this translates into the requirement J/2π ≫ L × 15 kHz. Thus, a larger 10 × 10 lattice, with L = 20, requires J/2π ≳ 3 MHz, J= A j j \ 0:01. We find that one of several experimental approaches can be taken: • If J t 0:03 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ε 1 ε 2 =η PF L p A j j, we can choose system parameters, and in particular κ > J, such that JT meas ≲ π/10. In this case, we can easily readout the state of the system faster than it evolves. Taking a typical ε 1 = ε 2 = 0.2, this regime can be reached by using narrow Purcell filters, η PF ≤ 0.01, and large bandwidth cavities, κ/2π ≳ 20MHz. Multiplexed nondemolition qubit readout with similar parameters was demonstrated in less than T meas = 50 ns 65 .
The experimental overhead for this approach is large in bigger lattices. As the Purcell filters must be spectrally very narrow, only one cavity can be brought in direct resonance with each filter, and so each qubit needs a readout resonator and a separate Purcell filter. Cavity frequencies must be spaced sufficiently far apart, at intervals of Δ ≳ 100 MHz. Another downside to this approach is that narrow filters, while increasing the qubit lifetime, make it hard to drive the qubits through the readout line. As we discuss below, we find that this is a useful tool in preparing states that explore the system's many-body properties, and if the readout cannot be used in this way separate drive lines would be necessary.
j j, measurement speed is limited. However, if the Stark shift generated by driving the cavity detunes the measured qubit away from the lattice, δω q $ nχ $ ε 1 A j j ) J, its state is frozen and we can once again readout a snapshot at a given time.
A design of this form would call for χ % Àκ ð Þ =2π % 5 MHz. If neighboring qubits are to be measured simultaneously, the driving pulses must be carefully calibrated to maintain a frequency detuning, which means the protocol may not be robust.
• Finally, if J \ ε 1 A j j, we necessarily have T meas ≳ 1/J. This is the weak continuous measurement regime 66 , and the amount of information that can be extracted about the system is reduced: we would not, for example, be able to obtain the probability statistics required to measure the entropy of a state. In this regime, full readout can be enabled by turning off interactions, either directly 67 , or by making use of frequencytunable qubits.
We can effectively freeze out the interactions between the qubits by mutually detuning their frequencies, essentially shifting the system into the individual particle regime. Note that we do not need an infinite array of frequencies, as only coupled qubits must be detuned from each other. In a square lattice, qubits can be detuned in a checkerboard pattern, as described above. This freeze out can be achieved by attaching fast flux lines to N/2 qubits, requiring comparable or reduced overhead to the use of individual Purcell filters, similar to previously realized setups 18,19 . This method also allows for more flexibility in measuring observables other thanσ z , as rotation pulses can be applied to the qubit between detuning and measurement, possibly through the cavity array, as described below.
Although the resonator configuration discussed above enables selective readout of specific or all qubits, it does not facilitate individual qubit control with microwave drives when all the qubits in the lattice are degenerate.
This issue can be overcome in several ways. Most directly, it may be useful to couple control lines to a single or few specific qubits to allow for direct microwave control, e.g., to prepare a certain initial state in the lattice. Alternately, the use of tunable qubitswhich we have suggested above for the purpose of a freeze-out prior to qubit readout-allows one to address an individual qubit or a subset of qubits if they are detuned in frequency away from the otherwise degenerate lattice.
In addition, the readout layout described above can be used to effect a specific form of system-wide driving. We note, when a signal line ℓ is driven at near resonance, at ω d ≈ ω q , the effective Hamiltonian becomeŝ where the driving operator is given bŷ Here, the summation is over the set of qubits S ' coupled to the signal line ℓ (Fig. 5a),α i is the effective relative coupling to that qubit, determined by the resonator's parameters, andg is a coupling energy proportional to the driving strength. See the Methods section for the derivation of this operator and the values ofg;α i . Although the set of operatorsΣ ' does not allow us full control of the system, driving at different strengths or for different lengths of time allows us access to a set of defined unitary transformations. As mentioned above, this would allow the measurement of quantities such as the entropy of a subsystem 45,46 . As we discuss below, it also enables the preparation of many-body states whose nature is determined by the detuning of the drive from the qubit frequency and can be used to probe the spectrum of the system.

Coherent-like states
As we have seen, the most interesting behavior of the HCB is manifest in the finite-excitation density sectors where 0 < n/N < 1. Within these sectors, energy eigenmodes vary in their behavior between the edges of the band and its center, exhibiting manybody properties such as different entanglement entropy laws. To study these properties, we must be able to prepare such states, which is challenging. In our proposed implementation, state preparation can be performed by applying drive pulses that reach the qubits via the readout resonators. To prepare a specific eigenstate, we would not only have to tailor a series of specific pulses, but also know the wavefunction of the prepared state, negating the premise of a quantum simulator to access states which are not understood theoretically.
Here, we propose an alternate route to observing the spectral properties of the HCB. Instead of preparing a specific known eigenstate, we apply a weak drive using the operators of Eq. (23) at some detuning from the joint qubit frequency. This prepares the lattice in a coherent state-like superposition of eigenstates in multiple n sectors, but with definite kinetic energy within each sector. This strategy of extracting many-body properties is robust with regards to experimental control limitations on chip.
To understand this process, we begin by rewriting the Hamiltonian of Eq. (7) in its eigenmode basis, where n; ϵ j i are the eigenstates of Eq. (8) and ρ n ϵ is the density of states for the sector with n excitations. Then, we rewrite the driving operator of Eq. (23) in the same basis, Consider first the perturbative limit, where the driving is very weak compared with the energy spacing, In this case, the driving operator will couple only eigenstates differing exactly by the detuning, and we can approximate it as a combination of defined-energy raising operators where ϵ n = ϵ + n × δ. Observe that eachÂ ε couples a subset of eigenstates of the form n; ϵ þ n δ j i . In the spectrum outlined in Fig. 3, these can be identified as the states sitting on a line with slope δ and intersecting n = 0 at ϵ.
Thus, if we initialize the system in the ground state, it is affected only byÂ 0 ;Â y 0 . Inserting the operators of Eq. (28) into the Hamiltonian of Eq. (22), we find at later times it has a form reminiscent of a coherent state, While this wavefunction is difficult to evaluate theoretically, it is composed only of states of a defined energy, n; n δ j i , i.e., ψ t ð Þ j i% X n e Àiω d nt c n t ð Þ n; n δ j i (31) for some time-dependent functions c n t ð Þ. As described above, these eigenstates lie along a line with slope δ in the spectrum shown in Fig. 3. This form can be observed in Fig. 6a for a numerical simulation of a system with very weak driving.
In practice, the approximation of Eqs. (26) and (28) are insufficient to describe the dynamics. For any fixed n/N, the energy spacing between states shrinks exponentially with N as we increase the size of the lattice, violating the assumption of Eq. (26). For weak driving, the qualitative picture remains similar but the prepared state seen in Eq. (31) acquires a finite width in energy space, proportional to the driving strength. These features are seen in Fig. 6.
We have discussed above how to prepare the HCB system in a coherent-like state. This state has a defined kinetic energy per excitation, but it does not have a definite excitation number. We argue that this is not an impediment to measuring the many-body properties described above.
First, we note that for any measurements purely in theσ z basis we can effectively project the state into a definite n sector by postselection.
Second, we have observed in Fig. 4b that the many-body properties that we are interested in behave similarly in different n sectors of the spectrum. For these, we expect the state in Eq. (31) to exhibit the same behavior as a function of its kinetic energy.
As such, preparing these coherent-like states may allow us to measure many-body properties of the spectrum by varying the detuning δ. We verify this numerically in Fig. 7, where we described a state prepared this way and measure its many-body properties. We find that the correlation length of the state, shown in Fig. 4a, approximates very well the eigenmode correlation length for states in similar energy show in Fig. 7b. Similarly, the entanglement entropy measured as shown in Fig. 7c exhibits the same behaviors we pointed out in Fig. 4b.

DISCUSSION
We have offered here a roadmap for the realization of a quantum many-body simulator of the 2D hard-core Bose-Hubbard model using a superconducting circuit made up of transmon qubits. An experimental realization of this setup would allow the exploration of this analytically hard-to-solve model in regimes where it has not been realized before. In particular, we have shown how such a realization could access non-equilibrium states that exhibit manybody wavefunction behaviors such as a crossover from volumelaw to area-law entanglement. As discussed throughout, the experimental parameters we consider in this article are within reach of current fabrication and control systems. The system we have proposed could be realized in the near term.
In the body of this paper, we have presented numerical results for a 4 × 4 HCB lattice, which can be diagonalized on a moderately powerful computer. However, the difficulty of this task grows exponentially, and a system of 6 × 6 or 7 × 7 sites is beyond numerical reach for any reasonable resource expenditure. An experimental realization would thus provide an example of Fig. 6 Coherent-like state preparation. We show here numerical results for the same 4 × 4 system shown in Fig. 4, Δω = 0.2J, with a driving term described byĤ dr ¼g e i ωqþδ ð ÞtΣ þ e Ài ωqþδ ð ÞtΣy , whereΣ is as described in Eqs. (43)- (46), taking realistic experimental parameters. We plot the overlap of the state with different eigenvalues, n; ϵjψ t ð Þ h i j j 2 at different times, and with different driving strength (arbitrary scale, darker colors denote greater overlap). Here, we drive the system at δ = −J. a We plot the evolution of the state from the initial ψ 0 ð Þ j i¼ 0; 0 j i for very weak driving,g ¼ 0:5J. We see that at any time the state can be described by a superposition of eigenstates n; δ n j i , as discussed around Eq. (31). b We plot the prepared state Ψ j i ¼ ψðt ¼ 8 J=g 2 Þ at varyingg. For stronger driving, the energy width of the prepared state grows as ΔE /g. quantum simulation beyond our theoretical and numerical abilities.
Beyond the model presented here, the paradigm of the QMBS can be used to explore a variety of other systems. Two immediate extensions of the model include changing the lattice topology or varying individual qubit frequencies to understand the role of disorder in this many-body system. In the longer term, it would be interesting to explore other parts of the phase diagram in Fig. 2. In particular, a reliable and long-lived qubit with large anharmonicity would allow us to realize spin systems and explore their rich physics, including probing phase transitions and understanding spin liquids.

Driving through a line coupled to multiple qubits
Here, we give the derivation for the driving operator of Eq. (23).
The system used, schematically shown in Fig. 5a, is described by the Hamiltonian whereĤ HCB , given in Eq. (7), describes the qubits,Ĥ L ' the signal line ℓ, andĤ ' i the resonator coupling line ℓ to qubit i, Here, ℓ sums over the different signal lines; for each line,R ' ν (L ' ν ) are the annihilation operators for its right (left) moving modes with energy ν, and S ' is the set of qubits coupled to it through resonators. For each resonator coupled to qubit i,ĉ i is the annihilation operator for a photon in the resonator, and Δ i , g i , κ i , τ i are that resonator's detuning, its coupling to the qubit i, its linewidth, and its distance from the termination of the signal line (divided by the speed of light), respectively. This setup is outlined in Fig. 5(a).
Using standard input-output theory 68 , the Heisenberg-Langevin equations of motion for the operatorsĉ i are are Gaussian white noise operators describing the vacuum fluctuations of the left-moving and right-moving modes, respectively, on the line ℓ coupled to i.
If we drive the line at frequency ω d , we have where Ω is the driving field. We find, in steady state, In a dispersive readout scheme the linewidths of the cavities are narrow 61 , If we then drive near the qubit frequency, we can approximate Now, from Eqs. (32) and (34), we have that the driving Hamiltonian can be described bŷ and for ω d~ωq , we can take the rotating wave approximation and combined with Eq. (40) we find Fig. 7 Probing many-body properties of the HCB with coherent-like states. We show here numerical results for the same 4×4 system shown in Fig. 4, Δω = 0.2J, with the drivingĤ dr ¼gðe i ωqþδ ð ÞtΣ þ e Ài ωqþδ ð ÞtΣy Þ, as in Fig. 6, applied for time t ¼ 8 J=g 2 to prepare the state. Here, we maintain the driving strengthg ¼ J and vary over the detuning δ. Ψ j i is prepared withΣ is as described in Eqs. (43)- (46), as in Fig. 6, while Ψ ϕ is prepared withΣ ¼ P e iϕ iσ À i for uniformly distributed, random ϕ i . a We plot the overlap of the prepared state with different eigenvalues, n; ϵjΨ h i j j 2 (arbitrary scale, darker colors denote greater overlap). We see that different values of δ access different parts of the many-body spectrum. b, c We compare the many-body properties of the two prepared wavefunctions at various values of the detuning (black and red lines) to those of the equivalent eigenmodes we expect it to be composed of (colorful dots). These are reproduced from Fig. 4 with the energy axis rescaled for comparison. We find remarkable agreement both for b the correlation length [Eq. (11)] and c the ratio s V /s A between the volume coefficient and area coefficient of the entanglement entropy [Eq. (15)] for the prepared states. wherê as in Eq. (23), and Circuit analysis of the floating transmon qubit We review the Hamiltonian of a floating transmon qubit coupled to a harmonic oscillator mode, depicted in Fig. 8a. This allows us to extract effective values for the qubit capacitance C q and the coupling capacitance C c comparable to those of a grounded transmon qubit, shown in Fig. 8b. In the Results section, we utilize this result to compare unwanted crosstalk in an architecture with floating transmon qubits versus an architecture that makes use of grounded transmons. Following the node flux representation described in ref. 69 we can write down the Lagrangian for the circuit in Fig. 8a as where Φ i are node fluxes and ϕ i = 2πΦ i /Φ 0 node phases, with Φ 0 the magnetic flux quantum. E J = Φ 0 I c /2π is the Josephson energy of the Josephson junction with critical current I c . The kinetic part of the Lagrangian can also be written as where Φ ! ¼ ðΦ 1 Φ 2 Φ 3 Þ T and C the capacitance matrix defined by Eq. (48). In order to recover the relevant transmon degree of freedom, we perform a variable transformation in the transmon subspace to 'plusminus' variables Φ ± = Φ 1 ± Φ 2 . With the transformation matrix we can rewrite the capacitive part of the Lagrangian as Here Φ ! 0 ¼ S Á Φ ! ¼ ðΦ þ Φ À Φ 3 Þ T , and the capacitance matrix in the transformed basis becomes C ¼ S À1 CS À1 .
A Legendre transformation yields the circuit Hamiltonian where q ! 0 ¼ ðq þ q À q 3 Þ T . Since the '+'-mode of the transmon does not have an inductive term in the Hamiltonian, its frequency is not relevant for qubit operation. Conversely, the Josephson energy of the transmon enters via the '−'-mode. We can therefore trace over the q + degree of freedom, to find where is the matrix C À1 with the column and row corresponding to the + mode removed.
The effective Hamiltonian of Eq. (53) has the same form as the Hamiltonian resulting from analysis of the circuit Fig. 8b, with the substitutions Φ 1 → Φ − , q 1 → q + . We can then find the effective parameters of the reduced circuit by identifying We therefore find the effective transmon capacitance, including coupling capacitances to the resonator, from the diagonal entry, and from the off-diagonal entries we can extract the effective coupling capacitance between the floating transmon qubit and the resonator Applied to the circuit in Fig. 5b, we find the parasitic coupling between the floating transmon qubit and the resonator C ðfÞ eff ¼ taking C g1 = C g2 = C G .

DATA AVAILABILITY
The results of the simulations generated during the study are available from the corresponding author on reasonable request and with the approval of our US Government sponsor. Circuit diagrams for a resonator coupled to a a floating transmon qubit and b a grounded transmon. Independent nodes are labeled with their respective node phase ϕ i = 2πΦ i /Φ 0 , relating to node fluxes Φ i .