Realizing the 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

Here, we show how to realize the two-dimensional (2D) hard-core Bose-Hubbard model (HCB) illustrated in Fig. 1 using an array of transmon qubits [17], 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 [18]. Outside of one dimension (1D), this system has no known analytical solution, and its study has been mostly through numerical methods limited in their scope. The most successful approach has been the use of tensor network methods, which focus on finding the ground state energy [19,20]. 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 * yariv@lps.umd.edu 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 many-body properties of its excited states. Previous experiments have realized the HCB in 1D [21], where the model can be solved by analytical meth-ods [22,23] and has the dynamics of a free fermion gas. 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.
The remainder of this article is structured as follows. In Section II we give an overview of the operating regimes of a superconducting quantum many-body simulator (QMBS) and the phase diagram of the simulated systems. In Section III we narrow our focus to the hardcore Bose-Hubbard model and give a theoretical description of the model and some of its interesting many-body properties. In Section IV we outline the physical realization of the HCB using a superconducting circuit based on transmon qubits. Finally, in Section V we describe the process of preparing states that allow us to observe many-body properties.

II. SUPERCONDUCTING QUANTUM MANY-BODY PHYSICS SIMULATOR
We consider the implementation of a QMBS with a superconducting quantum circuit made up of multiple repetitions of small basic circuits implementing qubits [24] and coupling elements. We describe the system with the HamiltonianĤ where summation is over all qubits i and over all coupled pairs i, j . The termsĤ Q i andĤ J i,j describe the basic qubit and coupling circuits, respectively. The system in Fig. 1 is one example of the Hamiltonian of Eq. (1), with circles (qubits) representingĤ Q i and diamonds (coupling elements) representingĤ J i .

Energy scale Description
Frequency variance 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 many-body physics to appear. Where ∆ω J ωq, the behavior is particle-like and we expect to see a version of the Bose-Hubbard model with A playing the role of on-site interaction. Where ωq J |A| 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.
We observe that the QMBS can be characterized by four energy scales derived from these Hamiltonians, outlined in Table I. The qubit frequency ω q and hopping strength J are the typical energy scales of the qubit and coupling, respectively, and the anharmonicity A describes the deviation of the qubits from harmonic level spacing. The frequency mismatch ∆ω is the scale of nonuniformity across the system, including, e.g., variation introduced during fabrication. We neglect deviations in the coupling strength, and assume that the deviations in the first level spacing are typical of the rest of the spectrum.
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 many-body dynamics to appear in the lattice. The ratios J/|A|, J/ω q then determine which states are effectively coupled, and so which theoretical models are accessible [25]. These features are collected in the form of a phase diagram in Fig. 2 and discussed in further detail below.

A. Hopping dominant
In the regime at the top right corner of Fig. 2, the dominant energy scale is J, the coupling energy. 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 lowtemperature 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 [26,27].
We note also that in other cases 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 elementŝ

B. 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 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â † i is the creation operator for site i andn i =â † iâ i is its energy level. Here, the anharmonicity plays the role of the on-site interaction strength, while the coupling term allows for hopping or nearest neighbor interaction. The sub-regime where J < |A| ω q -the working point of the transmon qubit [17] -is the most experimentally accessible parameter regime and is widely adopted by the superconducting circuit community in a multitude of experiments (see e.g. [16]), including recent implementations of 1D Bose-Hubbard lattices [21,28]. In this manuscript, we focus on this regime from Section III onwards.
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.
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 excitationconserving 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 the ultra-strong or deep-strong coupling regime [29], and it is more challenging to realize experimentally. However, superconducting artificial atoms are more suitable for its realization than real atoms coupled to an electromagnetic cavity, as their coupling strength to a harmonic oscillator mode is not necessarily limited by the fine structure constant [30,31]. In general, physical couplings in the deep-strong coupling regime can be achieved with strongly non-linear qubits and highimpedance circuits [31]. A promising qubit modality to reach such high couplings is the flux qubit, where ω q |A|, as demonstrated experimentally [32,33]. The fluxonium qubit [34], an extension of the flux qubit, has recently been demonstrated to preserve long coherence times while in the high anharmonicity regime [35].

III. THE HARD-CORE BOSE-HUBBARD MODEL
For the remainder of this article, we focus our attention on the regime ∆ω J ω q , |A|.
This combines the two constraints mentioned in Sections II B and II C: 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  13)] as a function of energy for each eigenstate. We plot the sectors n = 6, 7, 8, near n ≈ N/2, because we expect these sectors to be strongly dominated by many-body effects. We see a transition from area-law behavior at the band edge to volume-law behavior at the center. 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 (see Section V).
whereσ z i ,σ ± i =σ x i ± iσ y i are the Pauli z and 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 |n, where where E G is the ground state energy. This spectrum is sketched out, in the rotating frame, Fig. 3(a). The HCB is difficult to solve except in some specific cases. The 1D chain can be solved through fermionization [22,23], and the case of n/N 1 (1 − n/N 1) can be understood analytically by perturbative corrections to the free particle (free hole) problem [36]; both these 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 [19,20]. 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 tensor network methods in large systems. Beyond this, as we discuss below, an experimental realization can access to the system's entire spectrum.
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. 3, 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.
a. Correlation length: The typical scale beyond which different sites are no longer correlated serves as an order parameter for phases with long-range order [37]. The correlation length is a limiting factor for the applicability of tensor-network methods, which can be used only when correlations are finite [38]. 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 This correlation is straighforward to measure directly in the experimental setup described in Section IV B. We then extract the correlation length of a state |n, by fitting C z i,j (|ψ ) to the form over all pairs of nearest neighbors and next-nearest neigh-bors. Here |r i − r j | = |x i − x j | + |y i − y j | is the Manhattan distance between the sites i, j. We plot the correlation length ξ as a function of eigenstate energy in Fig. 3(b) for a 4×4 lattice at n = N/2. 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.
b. Entanglement entropy: For a state with density matrixρ, the entanglement entropy of some subset X of the lattice is given by the entropy of the state when the rest of the system is traced out, whereρ X = Tr ∀i / ∈Xρ is the reduced density matrix of the subsystem X. 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") [18]. Volume-law states are also harder to approximate using tensor-network methods.
To describe the growth law for an eigenstate |n, , we extract the parameters s V , 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. 3(c). 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 in Section V.

IV. TRANSMON IMPLEMENTATION OF THE HARD-CORE BOSE-HUBBARD MODEL
The transmon qubit [17] is a natural building block for the implementation of the HCB with a superconducting circuit. It 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 [17]. The self-Kerr non-linearity of the transmon Hamiltonian maps directly onto the on-site interaction term in the Bose-Hubbard model [21]. Since the hard-core Bose-Hubbard model operates in a regime where J/|A| 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 ∝ N J ∼ |A|. Generally, this effect can be treated as a perturbative correction to the hard-core 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 [21], 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 µs to 100 µs [16], corresponding to Γ/2π 10 kHz.
Fabrication variations translate to variations in transmon transitions which may exceed 200 MHz [39], yielding disorder in the emulated model on the order of ∆ω ∼ |A|.
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 [21]). Additionally, dynamic flux control allows rapid detuning of the qubits, which freezes the system dynamics and allows for state preparation and readout of effectively decoupled qubits (see Section IV B).
As discussed above, an experimental implementation operating in the regime J/|A| C P /C P = 0 C P /C P = 0.2 C P /C P = 0.5 C P /C P = 0.8 The relevant capacitances are the direct shunting capacitance C sh between the qubit electrode pads, the capacitance CG of each pad to the ground, and the coupling capacitance CJ between electrodes of adjacent qubits. By maintaining this coupling between left electrodes (blue) and right electrodes (green), the phase of the mutual couplings around a four-qubit plaquette sums to zero. A different choice could be made to create an effective gauge field with flux π. The two Josephson junctions form a dc-SQUID with total effective critical current Ic, which allows the qubit frequency to be tuned. We also show the dispersive readout and control scheme. Here, three control lines lines ( = 1, 2, 3) are coupled to groups of qubits (marked S1, S2, S3) via detuned resonators. The different resonators on each line have different frequencies, and thus detunings from the uniform qubit frequency. Each resonator is characterized by its frequency ωq + ∆i, position along the input line τi×speed of light, linewidth κi, and qubit coupling strength gi.   [40][41][42][43], in contrast to recent realizations where one of the electrodes is grounded (e.g. Xmon qubits [44]). In circuit designs with an increasing number of qubits, circuit elements can come close in the layout or even overlap when using crossover fabrication techniques [45]. This results in unwanted spurious coupling, referred to as cross-talk. 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 [46,47]. 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 of a floating and a grounded transmon. We assume a (parasitic) capacitive coupling C P , C P between the resonator and the electrodes of floating trans-mon (circuit diagram in Fig. 4(b)), or C P to the electrode of a grounded transmon (Fig. 4(c)). 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 P as well as the capacitance to the ground C G . Assuming without loss of generality C P ≤ C P , circuit analysis (see Appendix B) yields an effective coupling capacitance We note C eff if C P ≈ C P . This corresponds to an effective confinement of electric fields, which is advantageous in larger and more complex circuits. The relation of the effective coupling capacitances C eff is plotted in Fig. 4(d) for typical parameters C P , C 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 , 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 cross-talk is to physically separate circuit elements by introducing a multi-layer chip layout (3D integration) [45]. This approach is particularly beneficial in the implementation of a 2D grid of qubits, since the circuit topology prevents in-plane access to interior qubits. In a planar layout, this can be resolved by using airbridges to cross over signal lines [48], but these are naturally prone to unwanted cross-talk. 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 lowloss superconducting trenches etched inside the silicon substrate [45].

B. Qubit readout
Individual qubit readout and control in devices with only few qubits is 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 [49,50]. We sketch out an example of this setup in Fig. 4(a).
As implied by the color coding in Fig. 4(a), 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 [45].
Measuring the qubit state requires integrating out the background shot noise. For a homodyne measurement, this results in a typical measurement time of [51] T meas 4κ −1 + 4κnχ 2 where κ is the resonator linewidth, χ its effective dispersive coupling andn is the mean number of photons in the resonator during readout. Faster readout can be achieved by increasing the resonator linewidth κ or cavity occupationn, but these are limited by induced loss through the readout resonator and by the dispersive regime requirement. If we keep fixed where γ P is the Purcell decay rate [52], we find the shortest measurement time in the hard-core regime to be When we compare this measurement time to the system dynamics, which generally evolve on a time scale characterized by 1/J, we find that the qubit anharmonicity determines our ability to read out information: • If J 0.1|A|, we can choose system parameters such that T meas 1/J. In this case, we can measure the state of the system faster than it evolves, essentially getting a snapshot of the state at a specific time.
• On the other hand, if J 0.1|A|, we necessarily have T meas 1/J. This is the weak continuous measurement regime [53], 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.
These regimes are shown in Fig. 5, where we see how the homodyne measurement captures either the momentary state or the average over many cycles.
The strong measurement regime can be reached by reducing the hopping energy so that J 0.1|A|. This is limited by two factors. First, as discussed in Section II, the hopping energy must exceed the frequency disorder across different qubits. Second, in order to observe the system's full many-body dynamics, we require J LΓ, where Γ is the qubit decoherence rate and L the system length. The strong-measurement regime thus requires With typical parameters of a transmon implementation described above, |A|/2π ∼ 250 MHz, Γ/2π ∼ 10 kHz, and ∆ω/2π ∼ 100 kHz, the energy scales are two orders of magnitude too close to satisfy these constraints.
An alternate approach makes use of frequency-tunable qubits. In this case, we can effectively freeze out the interactions between the qubits by 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, for instance, qubits can be detuned in a checkerboard pattern, where all "white" qubits remain at the original frequency and all "black" qubits are shifted away from resonance by ∆ω q . Experimentally, this could be implemented with a shared flux bias line for all "black" qubits, reducing the necessary number of signal lines connecting to the chip.
To effectively implement such a freeze-out, we require the detuning ∆ω q to be large enough to suppress hopping between coupled qubits, Similarly, the detuning must be picked so that it does not bring the system out of the hard core regime, i.e.
Finally, the rate of change in the frequency (i.e. the time to vary the magnetic field) must be fast enough that the state transfer is diabatic, leading to a variant of the Landau-Zener formula, The application of the freeze-out is shown in Fig. 5(d) and compared to the direct measurement approach discussed earlier, showing how a measurement of the momentary state can be made even for J 0.1|A|.

C. Microwave control
While 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 qubits -which 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 (see Fig. 4(a)),α 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 Appendix A for the derivation of this operator and the values ofg,α i . While 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. This in turn would allow the measurement of quantities such as the entropy of a subsystem [54,55]. 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.

V. COHERENT-LIKE STATES
As we have seen in Section III, the most interesting behavior of the HCB manifests 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. As described in Section IV C, 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 outlined in Section III. Instead of preparing a specific known eigenstate, we apply a weak drive using the operators of Eq. (24) 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 manybody properties is robust with regards to experimental control limitations on chip.

A. Preparing coherent-like states
To understand this process, we begin by rewriting the Hamiltonian of Eq. (7) in its eigenmode basis, where |n, 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. (24) in the same basis, ρ n+1 ρ n n + 1, |Σ † |n, |n + 1, n, |. (26) Consider first the perturbative limit, where the driving is very weak compared with the energy spacing, ∀n, , : g n + 1, |Σ † |n, 2 ρ n+1 ρ n 1.
In this case, the driving operator will couple only eigenstates differing exactly by the detuning, and we can approximate it as a combination of definedenergy raising operators where n = + n × δ. Observe that eachÂ ε couples a subset of eigenstates of the form |n, + n × δ . In the spectrum outlined in Fig. 3(a), 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 ,Â † 0 . Inserting the operators Eq. (29) into the Hamiltonian of Eq. (23), 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 × δ , i.e.
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(a). This form can be observed in Fig. 6(a) for a numerical simulation of a system with very weak driving. In practice, the approximation of Eqs. (27) and (29) 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. (27). For weak driving, the qualitative picture remains similar but the prepared state seen in Eq. (32) acquires a finite width in energy space, proportional to the driving strength. These features are seen in Fig. 6. FIG. 6. Results of the state preparation method described in Section V. We show here numerical results for the same 4 × 4 system shown in Fig. 3, ∆ω = 0.2J, with a driving term described byĤ dr =g e i(ωq+δ)tΣ + e −i(ωq+δ)tΣ † , whereΣ is as described in Eqs. (A12) and (A13), taking realistic experimental parameters. We plot the overlap of the state with different eigenvalues, | n, | ψ(t) | 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) = |0, 0 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 , as discussed around Eq. (32). (b) we plot the prepared state |Ψ = |ψ(t = 8 × J/g 2 ) at varyingg. For stronger driving, the energy width of the prepared state grows as ∆E ∝g.

B. Observing many-body properties
We've 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 thê σ z basis we can effectively project the state into a definite n sector by post-selection. This is useful for measuring, e.g., the correlation length shown in Fig. 3(b).
Second, we have observed in Fig. 3(c) that the manybody properties that we are interested in behave similarly in different n sectors of the spectrum. For these, we expect the state in Eq. (32) 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. 3(b), approximates very well the eigenmode correlation length for states in similar energy show in Fig. 7(b). Similarly, the entanglement entropy measured as shown in Fig. 7(c) exhibits the same behaviors we pointed out in Fig. 3(c).

VI. OUTLOOOK
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 many-body wavefunction behaviors such as a crossover from volume-law 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 propose in Section IV could be realized in the near term.
In the body of this paper we've 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 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, as mentioned in Section II C, to realize spin systems and explore their rich physics, including probing phase transitions and understanding spin liquids.
whereξ L ,ξ R are Gaussian white noise operators describing the vacuum fluctuations of the left-moving and rightmoving 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, We review the Hamiltonian of a floating transmon qubit coupled to a harmonic oscillator mode, depicted in Fig. 8(a). 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. 8(b). In Section IV 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. 57 we can write down the Lagrangian for the circuit in Fig. 8(a) 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Č the capacitance matrix defined by Eq. (B1). In order to recover the relevant transmon degree of freedom, we perform a variable transformation in the transmon subspace to 'plus-minus' variables Φ ± = Φ 1 ± Φ 2 . With the transformation matrix Here Φ = S · Φ = Φ + Φ − Φ 3 T , and the capacitance matrix in the transformed basis becomesČ = S −1Č S −1 . A Legendre transformation yields the circuit Hamiltonian where q = q + q − q 3 . 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 is the matrixČ −1 with the column and row corresponding to the + mode removed.
The effective Hamiltonian of Eq. (B7) has the same form as the Hamiltonian resulting from analysis of the circuit Fig. 8(b), with the substitutions Φ 1 → Φ − , q 1 → q + . We can then find the effective pa-rameters 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. 4(b), we find the parasitic coupling between the floating transmon qubit and the resonator