Quantum simulation of the Anderson Hamiltonian with an array of coupled nanoresonators: delocalization and thermalization effects

The possibility of using nanoelectromechanical systems as a simulation tool for quantum many-body effects is explored. It is demonstrated that an array of electrostatically coupled nanoresonators can effectively simulate the Bose-Hubbard model without interactions, corresponding in the single-phonon regime to the Anderson tight-binding model. Employing a density matrix formalism for the system coupled to a bosonic thermal bath, we study the interplay between disorder and thermalization, focusing on the delocalization process. It is found that the phonon population remains localized for a long time at low enough temperatures; with increasing temperatures the localization is rapidly lost due to thermal pumping of excitations into the array, producing in the equilibrium a fully thermalized system. Finally, we consider a possible experimental design to measure the phonon population in the array by means of a superconducting transmon qubit coupled to individual nanoresonators. We also consider the possibility of using the proposed quantum simulator for realizing continuous-time quantum walks.

Introduction. The achievement of the ground state of mechanical motion using nanoelectromechanical systems (NEMS), as demonstrated recently in remarkable experiments [1,2,3,4,5,6], opens up a new path for studying quantum behavior in macroscopic systems. Having been able also to coherently control and cleverly measure the state of the mechanical resonator [1,7,8,9,10,11], an immediate possibility to explore is the use of the NEMS as building blocks for fabricating analog quantum simulators to reproduce many-body quantum physics [12,13]. Analog quantum simulators are dedicated and controllable devices, which can imitate (within some accuracy) the evolution of certain types of Hamiltonians. Various quantum systems have already been investigated for quantum simulation. Previously, quantum simulators were experimentally implemented using ultracold quantum gases [14], trapped ions [15], photonic quantum systems [16] and superconducting circuits [17,18]. Recently, an array of optomechanical resonators has been suggested for simulating many-body nonlinear driven dissipative quantum dynamics [19].
One particularly important phenomena that emerges due to the wave-like nature of matter at the quantum regime is Anderson localization [20,21,22], a phenomena in which waves fail to propagate in disordered media due to interference. Anderson localization has been observed with many experimental setups, including microwaves [23,24,25], light waves [26,27,28,29,30,31,32], ultrasound [33,34] and matter waves [35,36,37,38,39,40]. Given the ability to achieve the quantum regime in mechanical resonators, it would be extremely appealing to observe mechanical localization in an array of NEMS, where thermal effects become relevant. Furthermore, arrays of NEMS could be used to simulate continuous-time quantum walk (CTQW) dynamics [41]. Quantum walks are the quantum version of random walks that have a crucial role in designing efficient quantum algorithms that outperform classical algorithms [42]. Implementation of CTQW has been already realized in a four-site circle using a two-qubit nuclear magnetic resonance quantum computer [43] and in a waveguide array [44]. A NEMS-based quantum simulator would provide a means to implement CTQW with phonons, opening up new opportunities for quantum algorithms and universal quantum computation [45].
In this paper, we propose a 1D array of coupled nanomechanical resonators for simulating the Anderson Hamiltonian, namely, a discrete tight-binding model without on-site interactions. The coupling between resonators is electrostatic (capacitive coupling), but could also be mechanic (elastic coupling) in order to be improved. The disorder can be induced in a controlled and predetermined manner through the appropriate variation in the design and nanofabrication of the geometrical dimensions of the nanoresonators in the array [46]. A qubit coupled to the chain can be used as a device for both, initializing and measuring the occupation probabilities of excitations. We also discuss the physical implementations of CTQWs using the proposed quantum simulator.

Model
An array of N capacitively coupled electromechanical resonators, as depicted in Fig. 1 (upper panel), can be described by the Hamiltonian where x j , p j , m j and ν j are the displacement from the equilibrium position, momentum, mass and frequency, respectively, associated with a single mechanical mode of each resonator. U j is the electrostatic interaction energy between the pair of nearest-neighbor resonators j and j + 1. The interaction between a resonator and its second nearest neighbor is small enough to be disregarded. Using a simple parallel plate capacitor model [see Fig. 1 (coupling)], the energy U j in Eq. (1) can be written as where C j and ∆V j = |V j+1 − V j | are the capacitance and the voltage difference between the resonators j and j + 1, respectively. The capacitance C j is expressed in terms of the vacuum permittivity 0 , the resonator area A and the equilibrium Figure 1 A chain of electrostatically coupled mechanical resonators. Every single resonator is charged by an individual source, to induce electric potential diferences, and is considered as non-ideal by being coupled to individual thermal reservoirs. Imperfections in the resonators fabrication naturally introduce a disorder in their respective principal oscillation mode, hence, Anderson localization of phonon excitations can be observed for a sufficiently large array and sufficiently low temperature. Depicted quantities are discussed within the text.
center-of-mass separation d. Note that the inevitable disorder appearing in d is here neglected because its effect would be to introduce a disorder in the coupling energies, which are too small to be taken into account. The oscillation amplitudes are close to the zero point fluctuations, hence, |δx j (t)| = |x j (t) − x j+1 (t)| d and the electrostatic potential (2) can be expanded in powers of δx j (t)/d, which up to the second order gives The linear terms will cancel out in the summation of U j , and only remaining linear terms are (x 1 (t) − x N (t))/d, which nonetheless will also be negligible in comparison to the remaining second order terms, due to the rapid oscillation of x j (t). We now quantize the system, by replacing the classical variables with the operators written in terms of the annihilation and creation operators of the resonator excitation, a j and a † j . In the interaction picture, the position and momentum operators are given byx respectively. Substituting them in Eq. 1, the rapidly oscillating terms vanish in the rotating wave approximation (RWA) and Hamiltonian (1) is converted to where the zero point fluctuation of the resonator j and ω j = ν j + 0A∆V 2 j 2d 3 mj νj is the rescaled mode frequency. ω j is randomly distributed over the range [ω − ∆, ω + ∆], where ω is the average frequency of the resonators and ∆ is called the disorder intensity. We assume that the difference between any pair of disordered mode frequencies is small, hence the resonance condition is almost satisfied. In other words ∆ ω. In fact for all simulations considered here we have assumed ∆ ≤ 10 −2 ω. Moreover in any practical situation ω g j , but ∆ ≥ g j (See Section 4 for discussion of engineering the system parameters, and Table 1 for experimental values), and therefore the validity of the RWA is always assured. Note that the coupling energies g j depend on several experimental parameters related to the resonator fabrication such as mass, frequency and geometry, as well as the chain parameters such as the distance between the resonators and the potential difference. Therefore they could also be random variables that would define an off-diagonal disorder. However, since ω j g j , the disorder in g j is negligible in comparison to the diagonal one introduced by ∆. To simplify the calculations we fix all the tunneling energies to g i = J.
The Hamiltonian (3) is the Bose-Hubbard Hamiltonian without inter-mode interactions, and a diagonal disorder in the on-site energies ω j . When the total number of phonons is restricted to 1, this Hamiltonian reduces to the Anderson tight-binding Hamiltonian (see below).

Closed System: Anderson Localization
Localization, from a general point of view, is a mesoscopic phenomena displayed by waves as they propagate through a disordered medium. It is built upon the interference of the many randomly scattered waves, which at sufficiently large distances, collude to produce a suppression in the amplitude of propagation. This behavior is characterized by the exponential decay of the wave-function, with a decay length ξ, known as the localization length. The celebrated scaling theory of localization [47] describes how the transition between the different diffusion regimes depend upon the size L, and dimension D, of the system, independent of the microscopic intricacies of the disorder. According to the theory, depending on L, three different transport regimes can be recognized: ballistic (ξ L), where the size of the system is too small and scattering events are rare; diffusive (ξ L), where some weak localization effects may take place; and strong localization (ξ L), for large systems [48]. Of course, the detailed form of the localization depends on the type of disorder potential considered and its energy spectrum.
In the Anderson model (which is the one of interest in this paper), disorder is modeled by a δ-correlated potential V , consisting of a series of spatially uncorrelated barriers [1] , with a finite maximum amplitude intensity ∆, called disorder strength. [1] For instance in 1D, this means that the correlation function of the potential One can as well picture a periodic lattice with randomly shuffled on-site energies. Strictly speaking, a proper phase transition from extended to localized states only exists for D = 3, at some critical disorder intensity or, interchangeably, some critical energy known as the mobility edge. For D < 2 it can be shown that all states are localized, but the localization length for D = 2 can easily exceed the size of the system for weak disorder or high enough particle energies, thus admitting the existence of a diffusive transport regime. In contrast, for D = 1, (almost [2] ) all single-particle eigenstates are localized even for a vanishingly small ∆ and there is no phase transition [3] .
For the description of the linear chain of resonators, we can disregard the free propagation between lattice sites, focusing only on the occupation amplitudes of each site [4] . In the single-phonon situation, Hamiltonian (3) turns into ( = 1): where the states |j are the occupation amplitudes of each resonator, J is the coupling rate between neighbouring sites, and the set of ω j are assumed to follow a uniform random distribution around the mean value ω. For this model, numerical simulations are easy to perform, for example using closed boundary conditions at both edges of the chain, and diagonalizing Hamiltonian (4) directly, varying the total number of sites N and disorder intensity ∆, for many realizations of the disordered potential [5] . To set the stage for further discussions, we start in fig. (2) by plotting the behavior of the averaged relative dispersion in the site population for the ground state of the resonators ψ 0 , given by ∆n/n av = [ N j j 2 |ψ (0) j | 2 −n 2 av ] 1/2 /n av , as a function of the system size (N ranging between 1 to 300 sites) and for different disorder strengths (∆/J = 2, 5, 10,15,20). Here n av = N j j|ψ (0) j | 2 , is simply the mean value of the site population, and ψ (0) j = j|ψ 0 are the different population amplitudes for the ground-state; the double angular brackets denote the average over the disorder realizations. It can be noted the sensibility of the model for small sizes, since although the eigenstates are always localized, the relative dispersion in population sites only converges for a large enough system, a situation that is more relevant as the disorder gets weaker as can be seen comparing the different curves for the region N 50. This dependence on the system size is clear in the curve in fig.(2) for ∆/J = 2. It would be required a much larger system for the curve to converge closer to the other ones.
The typical exponential profile for the population density is plotted in Fig. 3(a) for ∆/J = 15. A single excitation in the central site for a chain of 51 resonators is [2] Almost here means for almost every realization of the potential and for almost every energy, except perhaps for some pathological potentials. See [49] for example for more rigorous definitions of localization. [3] It is interesting to mention anyway that for some long-range correlated disorders an effective mobility edge can appear even in the 1D case [50]. [4] Discrete models bring essentially the same qualitative results as continuous ones regarding localization effects, and they are solvable for many disorder distributions. [5] In this work all the simulations were performed for 500 realizations of the disorder. assumed as the initial condition. In a Bose-Einstein Condensate (BEC), the measured quantity of interest is the density of states at a given position in the lattice, namely, |ψ(x)| 2 [35,36,37,38,39,40]. In our simulator, the quantity to be measured is the population of the first excited state of a given resonator j, which is ρ |ωj ωj | . Given the discrete nature of our system, we cannot expect a smooth Gaussian-toexponential transition in the population profile. However, we see that due to the presence of disorder in the diagonal terms of the Hamiltonian, ω j in Eq. (3), the spatial distribution of the excitations present in the chain remain always close to the initial spatial distribution. Note that in Fig. 3 the initial nonstationary regime is not shown. This short time which is disregarded corresponds to the transitory path to equilibrium. It is interesting to explore the physics that is difficult to be included in the Anderson model, and which is nonetheless present in a real implementation of the system. Specifically, the influence of thermal effects due to a surrounding reservoir for the resonators is considered in the following.

Open System
Here, we write a master equation for the Anderson Hamiltonian (3) to describe the effects of the environment on the chain and consequently to investigate the corresponding influence on the localization of the states. In this way, we need to take into account phonons and weak electromagnetic fields that are surrounding the nanoresonators. Assuming that each resonator in the chain is coupled to a bosonic bath, the corresponding interaction Hamiltonian takes the form where b j,l and b † j,l are the j-th resonator bath annihilation and creation operators and Γ jl is the coupling between resonator-bath modes. In order to derive the master  equation it is supposed that each resonator is weakly coupled to its respective bath, and that all bosonic baths are in thermal equilibrium, with mean thermal phonon number n. The master equation describing the dynamics of the chain is theṅ where ρ S is the system density operator, γ j is the effective resonator-bath interaction strength and n is the bath mean thermal phonon number, which is specified by using the Bose-Einstein statistics. We have considered γ j = γ = 1 MHz for all simulations in this work. Remark that since we are interested in considering the Anderson model for a single phonon population in the chain we must keep the temperature reasonably low. Considering actual temperatures reached in experiments for mechanical resonators in the quantum regime, this corresponds to very low n.
We have considered n varying from 10 −4 to 10 −1 [1], which nonetheless is sufficient to see a disturbance in the localization mechanism, without compromising our numerical simulation.
The effect of the reservoir on the system include loss of excitations due to imperfections in the medium. It simply means that the single phonons in the chain, eventually would leak into the medium. This phenomena can be understood like an incoherent spontaneous emission. The expected behavior of a system subject to an amplitude damping channel is the escape of all excitations in the chain. However, while such decoherence process is taking effect, we can argue the existence of enough quantum correlations as to assure the existence of localization phenomena. Working in a weak resonator-bath interaction regime (something like three orders of magnitude less than the resonator-resonator interaction), the correlations remain in the system until the leaking process has taken away all measurable probability of excitations. We can clearly see this behavior in Fig. 3(b) for γ = 1 MHz,n = 10 −2 , and normalized disorder ∆/J = 15. Similarly to Fig. 3(a) a single excitation in the central resonator is taken as the initial condition. We see after equilibration a behavior which similar to he one characterising localization in Fig. 3(a) but for a increasing uniform thermal excitation base. In Fig. 3(c) the same initial condition and dynamics is employed but for a smaller disorder ∆/J = 2. Now we see a typical thermal behavior of a non-localized density profile. Remark that for γ = 1 MHz there is enough time for equilibration.
In Fig. 4 it is depicted more clearly the effect of the thermal excitation over the phonon localization. There it is plotted the time evolution of the resonators phonon population dispersion, ∆n , as a function of the thermal phonon number,n, for two different disorder strengths: 4(a) ∆/J = 2 and 4(b) ∆/J = 15, and for a single excitation in the central resonator as initial condition. It can be seen that the higher the temperature (orn) gets, the stronger will be the delocalizing effect as expected, very much independent of the disorder intensity. However for the lower temperatures (green and orange curves), the phonon population remains localized around the initial excited resonator, even in the presence of leaking in the total population of excitations due to dissipative effects. Furthermore, for the higher temperatures (blue and red), the localization is lost due to thermal pumping of excitations into the array, producing a fully thermalized state. With an initial condition available experimentally, namely, exciting a central nanomechanical resonator in the chain (see Section 4), the fast dynamics creates correlations in nearby resonators. These correlations imply long-time entanglement between the resonators, which in turn give us the possibility to maintain localization until thermalization is reached. To investigate the correlations in a finite chain of nanoresonators we consider the concurrence which is a bipartite entanglement measure and defined as [51]: where e i are the square roots of the eigenvalues in decreasing order of the positive definite matrix ρρ, wherẽ in which ρ * is the complex conjugate of the density matrix ρ and σ y is the Y Pauli matrix. Concurrence belongs to the interval [0, 1] where C(ρ) = 0, 1 corresponds to separable states and maximally entangled states, respectively. we measure the concurrence between the initially excited nanomechanical resonator (at the centre of the chain) and each of the other nanomechanical resonators under different conditions of the disorder and coupling with the environment. The results are shown in Figs. 5 and 6.  6 show that the behavior of the entanglement is consistent with that of the population probability distribution from Fig. 3. The concurrence and the population both spread through the available sites or become localized due to the disorder. The presence of an environment does not introduce new dynamics. The only effect of the environment, as expected, is to destroy quantum coherences with a rate proportional to exp(−γt), regardless of the presence of localization. The figures also show the sudden death of the entanglement [52] in which the concurrence between the central nanomechanical resonator and all the others destroys periodically.

Figures 5 and
The inclusion of thermal effects does not immediately change the main feature of the localization (Fig. 6, bottom). The system is quickly localized as expected from the Anderson-like Hamiltonian with disordered diagonal terms. With dissipation, localization is still a quantum effect in the sense that correlations between the resonators remain until the dissipation takes away all possible dynamics. As time goes on, depending on the temperature, all the states become thermalized. Even for very low temperatures, before all excitations leave the chain, the equilibrium state will be a thermalized state. However, the thermal relaxation rate is slow enough that localized phonon populations could still be measured prior to decaying. Continuous-Time Quantum Walk. The Anderson Hamiltonian given in Eq. (3) also generates CTQW dynamics. For a small disorder, an initial phonon injected in the central nanoresonator spreads ballistically and the standard deviation of the corresponding probability distribution over the chain increases linearly with time. The measurement method which is described in the following section can be used to reconstruct the probability distribution, shown in Fig. 7. It is also possible to inject two or more phonons in the chain to investigate the multi-particle quantum walks. That, specifically, provides a means for simulating bosonic particles and is going to be addressed elsewhere.

Measurement
Quantum measurement is a crucial step in quantum simulation. For our simulator, we envision using a single superconducting transmon qubit [53,54] for the readout of the phonon density profile in the nanoresonator array. In this scenario, the nanoresonators would be metallized and coupled capacitively [54] to the transmon as shown in Fig. 8a. It is important that the location of the nanoresonators be determined. It can be done by mapping the nanoresonator frequencies in the array, which could be compiled either before or after the measurements using electron beam imaging (to determine with high precision the nanoresonators' geometries) and finite element simulations. Through a filtered circuit connection [55], a DC voltage V could be applied to establish the interaction between each individual nanoresonator and the transmon, which can be modeled using the Jaynes-Cummings Hamiltonian [1,56] -note the same applied voltage would serve to couple the nanoresonators to each other. The full measurement Hamiltonian, incorporating the entire array of nanoresonators, would thus be given bỹ where ν a is the transmon's lowest transition energy, λ j is qubit-resonator coupling strength, σ z is the Pauli Z matrix and σ + (σ − ) are the qubit raising (lowering) operator. The preparation and measurement protocol (Fig. 8b) would rely upon utilizing the resonant limit of this model. To proceed, the transmon would initially be detuned in energy from all of the modes in the array and prepared in its first excited state with a microwave pi-pulse applied through a coplanar waveguide (CPW) cavity [57].
Through the application of a flux pulse, the transmon would then be brought into resonance with one of the nanoresonators and allowed to interact for one-half a Rabi cycle (t Rabi = π/2λ j ), thus transferring its excitation to the nanoresonator mode. Next, after a predetermined delay, the location of the phonon in the nanoresonator array would be measured by scanning the transmon's transition energy ν a (via a flux ramp) through the range of nanoresonator energies. Upon achieving resonance with the populated nanoresonator mode, the transmon would transition through a Rabi transfer back to its first excited state, which could be measured through dispersive read-out of the transmon via the CPW cavity [57,58]. As me mentioned the precise location of the nanoresonator could be determined from a map of frequencies using electron beam imaging and finite element simulations. Also, it should be noted that the coherent exchange of excitations through the resonant interaction between a piezoelectric disk resonator and superconducting phase qubit has been demonstrated previously [1]; thus it is expected that a similar technique could be adapted to the case considered here.
For realization of the protocol, and to insure low thermal occupation (n < 10 −1 ) of the nanoresonator modes at dilution refrigerator temperatures, it would be necessary for the nanoresonator array to be composed of ultra-high frequency (UHF) modes in the GHz regime [59]. With proper design of the nanostructures' geometries, the third-order, in-plane flexural modes could be engineered to have frequencies varying over the range of 2 to 4 GHz -note that varying levels of disorder ∆ could be programmed into the array in a controlled manner by deliberating varying the dimensions of the nanostructures. We remark that the larger the disorder the smaller will the the required array for the observation of localization. However that cannot be freely varied, since a large disorder could mean that the resonators are far from resonance to each other. That would imply the necessity to consider the counterrotating terms in Hamiltonian (3). Therefore is good to keep the limit J ∆ ω j for the envisaged platform. We estimate that transmon-nanoresonator coupling strengths λ j ≥ 1 MHz can be achieved for this configuration by designing a tightly packed array that minimizes the spacing d between nanoresonators (Fig. 1). Table 1 provides estimates of λ j (and J) for geometrical parameters of the array that are realizable using standard nanolithographic techniques. With λ/2π ≥ 1 MHz, Rabi transfers between the transmon and nanoresonators would occur over a time-scale of at most 250 ns, which would set the minimum dwell time for each step in the flux scan to determine which nanoresonator mode is populated. For an array size of 50 nanoresonators, this would yield a maximum scan time of 25 micro-seconds, which in turn would set the minimum tolerable nanoresonator relaxation time. For GHz-range modes, this would require quality factors in excess of 100,000, which has been achieved with aluminum [10] and carbon-nanotube-based resonators [60] at milli-Kelvin temperatures, albeit at much lower resonance frequencies (10s MHz).

Tables 5 Conclusions and Perspectives
In this paper, we have devised a quantum simulator based on nanoelectromechanical systems, for analyzing the many-body effects in quantum systems. Actually, Figure 8 Illustration of the proposed layout and experimental protocol for implementing the simulation of the Anderson Hamiltonian. (a) Sample layout. The voltage-biased nanoresonator array would be capacitively-coupled to one pad of transmon's shunt capacitor C Q . The transmon would also be capacitively-coupled to a CPW cavity for read-out. Note that layout is not drawn to scale -an array of 40 nanoresonators would extend only several microns along the transmon pad. (b) Protocol for preparing an initial photon excitation in a nanoresonator (ω i ) and detecting the phonon after delay time τ d in a second nanoresonator (ω j ). The transmon transition energy νa would be controlled by application of flux ramp Φ(t) to the transmon. Note the ramp rates as shown are not drawn to scale. Table 1 Estimates of ω assuming the standard expression from thin-beam theory for the third in-plane flexural frequency with clamped-clamped boundary conditions. The nanoresonator is assumed to be aluminum, with the following dimensions: length L = 0.7 µm (0.6 µm) for the 2.5 GHz (3.5 GHz) mode; width w = 45nm; thickness t = 50nm; The transmon-NR coupling λ was calculated from circuit theory, assuming the systems to be in resonance, and is given by Here C = 20 aF, dC/dx = 6 × 10 −11 F/m are the coupling capacitance and its first derivative respectively, which are estimated from finite element simulations assuming a gap of d = 20 nm. C Q = 50 f F is the transmon shunt capacitance, which would yield a charging energy of E C /h = 400 MHz. Finally, m = 0.52 . ρwtL is the effective mass of the third mode, where ρ is the density of aluminum. The same parameter values are assumed for the calculation of J. a one-dimensional array of electrostatically coupled nanoresonators is suggested to simulate the Anderson Hamiltonian. A method is present for coupling the nanores-onator electrostatically, but other means could also be explored to establish larger coupling (such as mechanical links). By introducing a controlled source of disorder to the system, we studied the localization phenomena in an array of 50 resonators. Accordingly, the population of the first excited state of a given resonator was analyzed, however, due to the discrete nature of our system, one could not expect a smooth Gaussian-to-exponential transition in the population profile. We also studied the influence of thermal effects due to a surrounding environment, arising in real implementation of the system. By coupling the system to bosonic thermal baths, we studied the interplay between disorder and thermalization. For sufficiently low temperatures, so that the localization is experimentally detectable with the proposed simulator, the loss of phonons due to the dissipation does not immediately destroy the phonon population localization. For higher temperatures the localization is affected by the thermal pumping of excitations into the array, which generates a fully thermalized state.
Initializing the system in a known state and measuring the evolved system, effectively, are important steps in realizing a quantum simulator. Having coupled the chain of resonators to a superconducting transmon qubit, we have suggested detailed protocols to initialize and measure the system. By applying a flux pulse, the transmon qubit is brought into resonance with the desired nanoresonator which allows to interact with the resonator hence transferring its excitation to the nanoresonator. After system evolution, the location of the phonon in the nanoresonator array can be measured by scanning the transmon's transition energy through the range of nanoresonator energies. Having achieved resonance with the populated nanoresonator, the transmon would transfer back to its first excited state. The transmon can then be measured through dispersive read-out via the CPW cavity.
Besides simulating Anderson localization, we have also discussed the possibility of using the proposed quantum simulator for implementing the continuous-time quantum walk dynamics. The system allows to realize localization and decoherence in continuous-time quantum walks. The initialization and measurement protocols permit to inject several phonons in the chain to investigate the multi-particle quantum walks. A (non-trivial) two-dimensional version of the suggested simulator can be used to implement two-dimensional quantum walks. In such case, quantum algorithms can be implemented. It is worth to mention that to some extent the present discussion applies as well to optomechanical systems, such as in Refs. [61,62], whose architectures are worth to be explored for quantum simulation purposes. Moreover, given the ability to control individual resonator excitations by scanning the transmon frequency, it is possible to extend the present investigation to the situation including on-site interactions. When several realizations are taken into account the effective Hamiltonian averages out to (3) plus on-site interactions on all sites, in a similar way to the Bose-Hubbard model, allowing investigation of the allowed phases and respective quantum phase transitions when the parameters are varied [63,64]. In that way this would allow for investigation of nonequilibrium steady states of quantum many-body models in similar fashion to other proposed simulators [65,66].