Deterministic Many-Resonator W Entanglement of Nearly Arbitrary Microwave States via Attractive Bose-Hubbard Simulation

Multipartite entanglement of large numbers of physically distinct linear resonators is of both fundamental and applied interest, but there have been no feasible proposals to date for achieving it. At the same time, the Bose-Hubbard model with attractive interactions (ABH) is theoretically known to have a phase transition from the superfluid phase to a highly entangled nonlocal superposition, but observation of this phase transition has remained out of experimental reach. In this theoretical work, we jointly address these two problems by (1) proposing an experimentally accessible quantum simulation of the ABH phase transition in an array of tunably coupled superconducting circuit microwave resonators and (2) incorporating the simulation into a highly scalable protocol that takes as input any microwave resonator state with negligible occupation of number states |0>and |1>and nonlocally superposes it across the whole array of resonators. The large-scale multipartite entanglement produced by the protocol is of the W-type, which is well-known for its robustness. The protocol utilizes the ABH phase transition to generate the multipartite entanglement of all of the resonators in parallel, and is therefore deterministic and permits an increase in resonator number without increase in protocol complexity; the number of resonators is limited instead by system characteristics such as resonator frequency disorder and inter-resonator coupling strength. Only one local and two global controls are required for the protocol. We numerically demonstrate the protocol with realistic system parameters, and estimate that current experimental capabilities can realize the protocol with high fidelity for greater than 40 resonators.


I. INTRODUCTION
Entanglement is an essential resource for a wide range of fundamental and applied uses. The generation of entanglement among linear resonators, however, is a difficult problem, for the fundamental reason that nonlinear processes are required to generate nonclassical states. Entanglement generation in linear resonators therefore requires the assistance of other nonlinear systems such as atoms or qubits, which makes the entanglement of physically distinct resonators a significant challenge. In 2011, Wang et al. [1] met this challenge by demonstrating deterministic entanglement of photons across two separate on-chip superconducting resonators with the assistance of tunable phase qubits and an auxiliary resonator. The protocol that they employed is extendible, with an increase in complexity, to more than two resonators, but the increase in complexity makes the protocol unfeasible for the many-resonator regime. Feasible protocols to achieve the entanglement of a large number of physically distinct resonators remain unformulated to date.
In a separate vein, the standard Bose-Hubbard model [2] consists of repulsive on-site interactions that compete with intersite hopping to give rise to the well-known superfluid-to-Mott-insulator quantum phase transition. In the superfluid phase, all of the bosons occupy a single reciprocal mode of the lattice and are consequently delocalized, but the phase is coherent across the whole lattice. For integer ratios of boson number to total number of lattice sites, the insulating phase arises when the on-site interaction energy is sufficiently larger than the hopping energy, so that the bosons are localized to single sites and number fluctuations on each site are suppressed to zero. The superfluid-to-Mott-insulator phase transition was observed in the landmark experiment of Ref. [3]. In contrast, the attractive Bose-Hubbard model (ABH) [4][5][6][7][8] has attractive on-site interactions and supports a (quasi-) [9] quantum phase transition wherein, as the attractive interactions increasingly dominate over the hopping, the ground state jÉ gs i changes from the superfluid phase to a phase in which the bosons are collectively confined to the same site but nonlocally superposed over all sites and number fluctuations at each site are amplified: jÉ gs i% 1 ffiffiffi ffi M p Â ðjNi 1 j0i 2 j0i 3 ...j0i M þj0i 1 jNi 2 j0i 3 ...j0i M þj0i 1 j0i 2 jNi 3 ... j0i M þÁÁÁþj0i 1 j0i 2 j0i 3 ...jNi M Þ, where M is the number of sites and N is the number of bosons. As the ratio of hopping energy to interaction energy goes to zero, the ground state becomes exactly a W state. The ABH phase transition has been theoretically considered within the context of cold atoms [4,5], trapped ions [8], and polaritons in cavity arrays [10], but experimental demonstration in these platforms remains a difficult and unmet challenge. Meanwhile, on-chip superconducting-circuit systems have emerged as a very effective platform for quantum electrodynamics [11][12][13]. Theoretical and experimental activity regarding such circuit QED systems has recently begun to move toward the many-body regime, where quantum coherence can potentially be achieved over 1000 or more interconnected microwave resonators [13]. Further, this platform allows the unique capability of in situ tunable coupling between resonators [14][15][16][17]. Access to the many-body regime and tunable coupling makes circuit QED an attractive option for quantum simulation of many-body Hamiltonians [13,18], but a proposal for the simulation of the ABH phase transition has not yet been formulated in this platform.
In this theoretical work, we show that recent experimental advances in superconducting circuits can be used to realize a circuit QED system wherein the ABH phase transition may be simulated in an array of tunably coupled superconducting microwave resonators. Further, we present a protocol that uses the ABH phase transition to convert almost any input state jc in i of a single microwave resonator into a W-type state 1 ffiffiffi ffi . . jc in i M Þ that spans all M resonators of the array, thereby deterministically generating discretevariable multipartite entanglement of many resonators. The only fundamental restriction on jc in i in our protocol is that it must have negligible occupation of Fock states j0i and j1i. The scalability of our protocol to arrays with large M is due to the fact that the ABH phase transition entangles all of the resonators simultaneously, rather than one by one, so that only one local and two global controls are required for the protocol, regardless of the number of resonators. We estimate that existing technology makes the protocol feasible for up to M % 40 resonators. Our protocol complements well the experimental capability demonstrated in Ref. [19] of producing arbitrary microwave states jc in i in single on-chip microwave resonators.
The ability to deterministically create many-resonator W entanglement of nearly arbitrary microwave states in superconducting circuits is of significance for many reasons. For instance, the qubit W state 1 ffiffiffi ffi . . j1i M Þ is known for the robustness of its entanglement under loss [20,21]. In particular, the global entanglement decay of such a state under both phase and amplitude damping is known to be independent of the number of qubits [21]. The W-type states that are output by our protocol (given above) may also be considered as qubit W states because the jc in i are always orthogonal to j0i because of the restriction on jc in i mentioned above, and they therefore have the same robustness of entanglement. Therefore, although the mapping of the single-resonator state to a W state does not increase the lifetime of the state itself, the entanglement that is generated by the process is well suited to the many-resonator regime. Further, it is of particular use for tests of nonlocality that utilize large-scale W states [22,23]. Another significant aspect of our protocol is that when jc in i are coherent states (with a sufficiently large amplitude), the protocol outputs entangled coherent states (ECSs) [24,25], which are of both fundamental and applied significance in their own right [25][26][27][28][29]. For example, coherent states of resonators are quasiclassical states, and many-resonator ECSs therefore constitute large-scale Schrödinger cats. To our knowledge, manypartite ECS generation has not been feasibly considered in any platform. Finally, microwave resonators in superconducting circuits may interface with other types of systems [30], such as various types of circuit and noncircuit qubits [30], mechanical resonators [31,32], and (potentially) optical fields [33][34][35][36][37][38][39][40]. The protocol presented here may therefore provide an indirect route toward large-scale W entanglement in those systems that may be more feasible than other, more direct approaches.

II. SYSTEM AND MODEL
As schematically illustrated in Fig. 1, we consider a onedimensional lattice of microwave-frequency superconducting coplanar waveguide (CPW) resonators, each embedded with a SQUID that intersects the center conductor line, as conceptually introduced in Ref. [41]. The resonators have nearest-neighbor coupling via the tunable coupler demonstrated in Refs. [16,17]. As theoretically analyzed in Refs. [41,42], the SQUIDs can induce a negative Kerr nonlinearity in the microwave modes of each CPW such that, looking at only the fundamental mode c of each resonator and considering uniform parameters across the lattice, the Hamiltonian of the system in the interaction picture is where is positive, may be positive or negative [16,17], M is the total number of resonators in the lattice, j is the resonator index, and periodic boundary conditions are assumed. Both and may be tuned in situ via flux biases, and we designate their ranges as 0 max and À max max . (We note that a similar circuit QED system with open boundary conditions and no tunable coupling is theoretically studied in the driven-dissipative regime in Ref. [43]. We also note that in such a setup may not be made strictly zero, but its minimal value may be roughly 2 Â 10 kHz [41,42]. This residual value of may be canceled in each resonator through dispersive coupling of a qubit, as in Ref. [44], which demonstrated a positive Kerr nonlinearity of roughly 2 Â 1 MHz.) Experimental considerations related to damping, higherorder terms of the nonlinearity, and resonator-frequency disorder are discussed in Sec. V. Equation (1) is precisely the ABH model that is theoretically studied in Refs. [4][5][6][7][8]. For a fixed total number of quanta N and N > 1, the ABH phase changes qualitatively as a function of the parameter ¼ jj ðNÀ1Þ , as depicted in Fig. 2. 1 % 0:25 is a characteristic constant of the ABH that is largely independent of lattice size [5,7]. When < 1 , the attractive on-site energy dominates the hopping energy and the quanta in the eigenstates are collectively localized to single sites but superposed across all sites to form the nearly degenerate W states where k is an integer in the range 0 k M À 1. These W states become the exact eigenstates as ! 0, and our proposal for ABH simulation is the only one to date that is able to access this perfect W-state regime, as it is the only one in which is tunable to zero so that ¼ 0. 2 is another characteristic constant of the ABH, such that for > 2 , the ground state is a superfluid state. It is always the case that 1 2 . When 1 2 , the ABH phase is intermediate in nature between the W-state phase and the superfluid phase. 2 increases monotonically with M: for M ! 5 [5], which is approximately linear in M and is plotted in Fig. 3. The approximate phase diagram for the ABH is represented in Fig. 4 with the regime 1 > > 2 omitted. We also consider the case when only the first site in the lattice may have a nonzero attractive interaction with strength 1 and range 0 1 max : For N > 1, the two parameter regimes of interest are 1 ðN À 1Þ=2 ) jj, in which the quanta in the lowestenergy eigenstate are localized only at the first site, and 1 ðN À 1Þ=2 ( jj, in which the quanta in the lowestenergy eigenstate form a superfluid state across the whole lattice for any finite value of . The approximate phase diagram for this model is represented in Fig. 5.   Each resonator has a SQUID (green shaded boxes) connected to its center conductor and is coupled to its nearest neighbors with flux-tunable couplers (line-hatched boxes). The SQUIDs induce a flux-tunable nonlinearity into the resonator modes. Two bias lines are used to globally tune the SQUIDs and couplers separately. The first resonator (not shown) has a separate bias line for its SQUID.
Through local control of the nonlinearity 1 of the first resonator and global control of jÞ1 and , our proposed system can simulate both Eqs. (1) and (4). An initial state of the proposed system with N ! 2 in the lowest-energy eigenstate and zero quanta in the other eigenstates may therefore adiabatically transition between the single-site localized phase, the delocalized phase, and the W state. If instead the system starts in a superposition of different N in the lowest-energy eigenstate, the linearity of quantum mechanics dictates that each component of the superposition with N ! 2 may independently undergo the phase transitions. The entanglement protocol revolves around these phase transitions and therefore necessarily excludes initial states with number distributions that have a significant occupation of N < 2.

III. ENTANGLEMENT PROTOCOL
Mention is made in the literature of using adiabatic transitions of Hamiltonians to generate entangled Fock states [10,45,46]. In particular, Hartmann et al. [10] suggest using an atom-cavity realization of the ABH to create a polaritonic (photon-atom hybrid) approximate W state of single Fock states jni (with n > 1) in multiple optical cavities via adiabatic transition from the superfluid regime ( > 2 ) to the W-state regime ( < 1 ) by tuning . Here, we show that the ABH implemented in our proposed superconducting circuit can be used to create perfect W-type entanglement of arbitrary superpositions of photonic Fock states jni (n > 1).
With the system in the vacuum state and ¼ ¼ 0, the first resonator of the system is initialized in a state jc in i ¼ P n C n jni, where jni denotes a Fock state with n quanta and C n are complex amplitudes, so that the initial total system state is The probability distribution of jc in i in the number basis is therefore the probability distribution of the total number of quanta N in the system. jc in i therefore has the constraint that jC n j % 0 for n < 2. As mentioned at the end of Sec. II, each total quantum number n from the distribution of N may be treated independently. We may therefore express the total system state at all times as where jÉ n ðtÞi denotes the system state at time t with definite total quantum number n, and the evolution of each jÉ n ðtÞi may be considered separately. The entanglement protocol employs manipulations of 1 , , and that simultaneously evolve each jÉ n ðtÞi from jÉ n ð0Þi ¼ jni j¼1 (a Fock state of the first resonator) into spanning all of the resonators) for all n > 1. The total system state at the end of the protocol is therefore a W-type distribution of the input state jc in i: Because the entanglement protocol maintains each jÉ n ðtÞi in a single eigenstate of the system, we may express the time evolution as jÉ n ðtÞi ¼ e Ài n ðtÞ jnðtÞi eig ; where n ðtÞ ¼ 1 @ R t 0 E n ðt 0 Þdt 0 , and E n ðtÞ denotes the energy of the eigenstate jnðtÞi eig that is occupied by the n quanta at time t. We first describe the trajectory of jnðtÞi eig during the entanglement protocol to explain how it evolves from Q rÞj j0i r , and then explain how e Ài n ðTÞ ¼ 1 is achieved so that In general, n ðtÞ Þ m ðtÞ when n Þ m, which can result in a distortion of the phase information contained in the initial state.
The steps of the entanglement protocol are illustrated with numbered arrows in the phase diagrams of Figs. 4 and 5 and in the pictorial representation of Fig. 6, which depicts the evolution of the full system state jÉðtÞi [see Eq. (6)].
Let t x denote the time in the entanglement protocol after step number x. The first three steps of the protocol (see Figs. 5 and 6) use the model of Eq. (4) to convert each jni localized in the first resonator into a delocalized superfluid state across the whole lattice. In the first step, 1 is rapidly increased to max in a time Át 1 , which does not change the occupied eigenstate: jnðt 1 Þi eig ¼ jnð0Þi eig . The attractive interaction thus introduced into the first resonator provides an energy barrier against which may be tuned adiabatically in step 2 from 0 to À max in a time Át 2 ) 4 max =½ max ðn À 1Þ 2 . In step 3, 1 is tuned to zero (adiabatically with respect to the hopping dynamics that redistributes the quanta in real space) in a time Át 3 ) max ðn À 1Þ=2 2 max . Because of the negative sign of and FIG. 6. Illustration of the entanglement protocol for the case of four resonators. Square boxes denote single resonators, solid black lines shared between square boxes denote inter-resonator coupling () off, dashed black lines denote coupling on, green shading denotes resonator nonlinearity () on, white shading denotes nonlinearity off, periodic boundary conditions are implied, and dashed blue arrows denote steps of the entanglement protocol. jc Ã in i ¼ P n e Ài n ðtÞ C n jni denotes a version of jc in i with modifications to the phases of each number component. A detailed account of the process in terms of the individual jni is provided in the main text. In the first step, the nonlinearity is turned on for the first site only, which creates an energy barrier against which may be turned on adiabatically in step 2. jc Ã in i may then adiabatically transition to the superfluid state by adiabatically turning off 1 (step 3). At this point, jc Ã in i is in the highest excited superfluid state because of the lattice asymmetry introduced in step 1 and the sign of . After switching the sign of in step 4, jc Ã in i is in the superfluid ground state.
Step 5 ( is adiabatically turned on across the whole lattice) and step 6 ( is adiabatically turned off) induce the ABH phase transition and make ¼ 0, so that the perfect W state is achieved, with jc Ã in i nonlocally superposed in each resonator. Finally, turning off the resonator nonlinearities (step 7) restores the lattice to an array of uncoupled linear resonators that now contains the entangled state. With correct timing of the last step, n ðTÞ ¼ 0 (T is the total time of the protocol) and jc Ã in i is restored to jc in i.
the single-site nonlinearity during the first three steps, jnðt 3 Þi eig is actually the highest-energy eigenstate of the system. In step 4, is rapidly tuned from À max to þ max in a time Át 4 so that jnðt 4 Þi eig is the ground state of the system. jnðt 4 Þi eig can therefore adiabatically connect to the lowest-energy W state jÉ ðk¼0Þ W i of Eq. (2). Steps 5 and 6 of the protocol (see Figs. 4 and 6) achieve this adiabatic connection by using the ABH [Eq. (1)]. The adiabaticity of these steps is with respect to two different physical processes: In the region = 2 > 1, Ref. [5] shows that the distribution of the quanta among the reciprocal modes is roughly constant with , and the adiabaticity must therefore be with respect to the hopping dynamics that changes the intersite phase relationships; in the region 0 < = 2 < 1, Ref. [5] shows that the reciprocal-space quantum distribution changes approximately linearly with , and the adiabaticity must therefore be with respect to the on-site interaction, which redistributes the quanta in the reciprocal basis. (Analytical details are given in Sec. V B.) In step 5, the global parameter is adiabatically tuned from zero to max in a time Át 5 , so that ¼ max = max ðn À 1Þ.
Step 5 may place jnðt 5 Þi eig in any of the three phases of the ABH (see Fig. 2), depending on the exact values of max , max , and n. In step 6, is adiabatically tuned from max to zero in a time Át 6 , so that ¼ 0 and jnðt 6 Þi eig is the W state jÉ ðk¼0Þ W i. In this step, although the eigenstates become degenerate as ! 0, the on-site attractive interaction energy proportional to serves as an energy barrier that protects the system from transitioning out of the ground state as long as the tuning is adiabatic with respect to . As is now off, the superposition is locked into place and may be rapidly tuned in step 7 from max to zero in a time Át 7 without altering the occupied eigenstate: As the nonlinearity is now off, the CPWs are not hybridized with the SQUIDs and the state of the system is purely photonic.
The effects of the entanglement protocol on n ðTÞ may be understood by separately considering the contributions of the on-site nonlinearity terms and the nearest-neighbor hopping terms of Eqs. (1) and (4). In the case of a single site, the nonlinearity contributes a number-dependent phase n R t 0 j ðt 0 Þ 2 dt 0 [47], so that appropriate timing can make R T 0 j ðt 0 Þ 2 dt 0 equal to an integer multiple of 2 and the effective contribution to n ðTÞ is zero for all n. It may be conjectured that with appropriate timing, the same type of cancellation for all n could occur in the multisite case. However, this conjecture is far less trivial due to the changing distribution of the quanta with time.
Remarkably, numerical investigations show the conjecture to be correct. For the hopping terms, it can also be conjectured that the contribution to n ðTÞ can be approximately canceled by simply making Át 2 ¼ Át 6 and Át 3 ¼ Át 5 , so that the phase evolution due to the hopping terms when < 0 (steps 2 and 3) may cancel the phase evolution due to the hopping terms when > 0 (steps 5 and 6) because of the opposite signs of the hopping energy and because step 5 is in some sense the reverse of step 3 superposed on each site. Remarkably, numerical investigations reveal that this cancellation is also exact, even in the absence of adiabaticity. It is therefore possible to achieve n ðTÞ ¼ 0 through appropriate timing of the steps of the entanglement protocol, as we demonstrate in the numerical simulations below. As an alternative to the cancellation of n ðTÞ through timing, a calibration may be done whereby the phase n ðTÞ is measured after test runs of the protocol with jc in i ¼ jni 1 for different n, so that the arbitrary input states generated by, for example, the method in Ref. [19] may be prepared with appropriate offset phases for each number component, so that the n ðTÞ accumulated through the protocol are all canceled.

IV. DETECTION
Verification of W-state creation may be done by employing bipartite Wigner tomography [1] between different pairs of resonators to reconstruct their joint density matrix. This tomography is the technique that was used in Ref. [1] to show the creation of a NOON state of two superconducting-circuit resonators, and it should therefore be readily applicable to our proposed system. It does require, however, coherent drive access to each individual resonator that is to be measured as well as qubits coupled to each such resonator. Individual drive access should be possible, however, for proof-of-principle experiments in the few-resonator regime. The tomography entails a series of identical state preparations of the system, each one followed by destructive measurements of the two selected resonators by their corresponding qubits. A sufficient number of such measurements yields enough information to approximate the joint density matrix of the two resonators. By performing this process with different pairs of resonators after the entanglement protocol is run, the nature of the full system state after the entanglement protocol runs may be inferred. The technical details of the tomography process may be found in the Supplemental Material of Ref. [1].
Alternatively, if after the entangled state is prepared the system is made into a linear network by quenching from zero to max , the protocol of Tufarelli et al. [48] allows for a single qubit tunably coupled to any resonator to be used to reconstruct the state of the entire array. This protocol may be more suitable for systems with larger numbers of resonators.

V. EXPERIMENTAL CONSTRAINTS
With present-day devices [16,17,42], max =2 and max =2 may reach as high as hundreds of MHz and may be tuned on a time scale of a few nanoseconds. However, should be limited to about max =2 ¼ 30 MHz in the interest of preserving high experimental fidelity of the hopping Hamiltonian with the tunable coupler [49]. Modes c j have flux-dependent T 1 and T because of the flux-dependent hybridization of the CPWs with the SQUIDs. Q factors of over 1 Â 10 6 have been demonstrated for on-chip CPWs [50], which gives T 1 > 20 s for the 2 Â 7:5 GHz CPWs that we assume here. As mentioned in Ref. [42], the SQUIDs embedded in the CPW resonators should be able to achieve very long coherence times, considering the recent T 1 and T 2 measurements of the Josephson junction qubits in Refs. [51][52][53]. The 2D Xmon qubit demonstrated in Ref. [53], for example, shows a maximum T 1 of approximately 44 s and T 1 % T =2 % 20 s at the fluxinsensitive point. Further, considering an asymmetric SQUID as in Ref. [42] allows for T * 1 ms at the flux value, where is maximum [54]. For our purposes here, we therefore approximate a constant T 1 value of 20 s and a flux-dependent T value that varies linearly from 1 s to 300 s as is increased linearly from 0 to max . In this section, we show that these and other experimental capabilities make our protocol feasible for array sizes of greater than 40 resonators. In particular, for the choices of max =2 ¼ 30 MHz and max =2 ¼ 14 MHz, the ranges 2 n 40 and 2 < M 42 become available, which can accommodate high-amplitude many-partite ECSs.

A. Higher-order terms of CPW nonlinearity
The full expression for the nonlinearity (H NL ) introduced into the CPWs by the embedded SQUIDs is given in Eq. (27) of Ref. [42]. Looking only at the fundamental mode and using the relations È 0 ¼ h=2e, E J ¼ ð È 0 2 Þ 2 =L J , L J ¼ L 0 = l , and L 0 ¼ 2E 0 C =ð! 2 c;0 e 2 Þ, we find is the inductive participation ratio dependent upon the flux through the SQUID loop and ! c;m is the frequency of the CPW mode m in the absence of the SQUID. After the rotating-wave approximation, to be able to neglect nonlinear terms higher than ðc y cÞ 2 , we require Larger values of max come at a cost of reduced n max but allow for quicker adiabatic tuning of , as per the discussion below. Smaller values of max enable larger array sizes, also as per below. Larger values of n max may be accommodated for a fixed max by using higher-order CPW resonator modes. If we assume a fundamental mode frequency of ! c;1 =2 ¼ 7:5 GHz, selecting max =2 ¼ 14 MHz allows n max ¼ 40 when the fundamental mode is used, selecting max =2 ¼ 25 MHz allows n max ¼ 23 when the fundamental mode is used, and selecting max =2 ¼ 100 MHz allows n max ¼ 6 when the fundamental mode is used. dt ( max and d= 2 dt ( 2 max in the respective regions. We denote the value of at the end of step 5 as Ã ¼ max max ðnÀ1Þ . The consideration in the next subsection shows that it is of interest to minimize the time taken to tune from 2 to 1 . We find that the minimization occurs when Ã % 2 . (Although the minimization is not simultaneously possible for all n, it is sufficient to assume so as a rough estimate.) In this case, the adiabatic constraint to tune from 0 to max in a time Át 5 is d dt & 1 5 2 2 ðn À 1Þ. Since the constraint only becomes applicable as approaches 2 , the approximation ! 00 is made on the right-hand side to yield Át 5 * 5= max :

B. Adiabaticity
For step 6, we have the requirement d dt & 1 10 2 max ðn À 1Þ 2 . Tuning from max to ð 1 Þ in a time Át 6a requires Át 6a * 10 2 À 1 max 2 , and tuning from ð 1 Þ to zero in a time Át 6b requires Át 6b * 10 1 max 2 . The constraint on step 6 is therefore For max =2 and max =2 on the scale of tens of MHz, this constraint puts steps 5 and 6 on the 10-100 ns time scale, which is well below the estimated T 1 and T of the resonators.

C. Disorder
Disorder in inter-resonator coupling frequencies can be neglected when resonator frequencies greatly exceed the coupling frequencies [55], which is the case here. We focus therefore upon disorder in resonator frequencies.
The inevitable frequency spread Á! of the resonators relates to the discussion in Ref. [4] concerning the fact that the nondegeneracy of the lattice sites induces nonuniformity in jðn; tÞi eig in the site basis at a rate 1=Á!. This undesirable dynamics is relevant in the intermediate regime 1 < < 2 , which is traversed in step 6, given the assumptions above. If this intermediate regime is traversed on a time scale comparable to or longer than 1=Á!, the quanta have time to localize in the lowest-energy site. In order to achieve symmetric W entanglement with high fidelity, we therefore require Át int ( 1=Á!, where Át int is the time it takes to traverse the intermediate regime. From the previous subsection, Át int ¼ Át 6a . Letting Át int ¼ 1=10Á!, we find 1 Á! * 100 max ð1 À 1 = 2 Þ, where it is assumed that 2 ¼ max max ðn min À1Þ . This inequality gives which gives a lower bound on max . Smaller values of max (or n min ) enable larger 2 (and therefore larger M) for a fixed n min (or max ) but also require smaller values of Á! and larger Át 6 . Current capabilities demonstrate Á! * 2 Â 1 MHz for GHz-frequency CPW resonators [55]. Choosing n min ¼ 2 and max =2 ¼ 30 MHz, we thereby find the constraint max =2 ! 14 MHz. Using Eqs. (3) and (11), this constraint translates to M 42 and n max 40, which is sufficient for many-partite ECS creation. If Á! is reduced by a factor of 2, we find max =2 ! 7:5 MHz and M 80 for n min ¼ 2. Reducing Á! to 2 Â 0:1 MHz gives max =2 ! 1:6 MHz and M 360 for n min ¼ 2.

VI. SIMULATIONS
We demonstrate the entanglement protocol by numerically integrating the master equation where is the density matrix for the CPW chain, H is the system Hamiltonian, D½c j ¼ c j c y j À c y j c j =2 À c y j c j =2 is the amplitude-damping operator, G½c j ¼ c y j c j c y j c j À ðc y j c j Þ 2 =2 À ðc y j c j Þ 2 =2 is the phasedamping operator, 1=T 1 is the amplitude-damping rate, and 1=T is the phase-damping rate. The integration is done using the fourth-order Runge-Kutta method with a time step of size 10 À12 s. As discussed in Sec. V, we model the damping with a constant T 1 value of 20 s and a fluxdependent T value that varies linearly from 1 s to 300 s as is increased linearly from 0 to max .
In the first simulation, a system of three resonators (M ¼ 3) is considered with an input state jc in i ¼ 1 ffiffi 2 p j2i þ 1 ffiffi 2 p j3i. The system parameters max =2 ¼ 100 MHz and max =2 ¼ 30 MHz are used, which satisfies the constraints discussed in Sec. V, and the timings of the protocol steps are chosen to respect the adiabaticity constraints. The Hilbert space at each site is truncated at n ¼ 3. The fidelity jhÉðtÞjÉ W ij 2 of the system state jÉðtÞi with the target state j¼1 jc in i j Q rÞj j0i r during the final step (step 7) of the protocol is shown in Fig. 7. The total physical time of the simulation is 0:1064 s, at the end of which the fidelity is about 97:5%. The same simulation is also performed without damping (not shown), and a peak fidelity of 98:6% is found. Further numerical investigations (not shown) reveal that in the case of no damping, the peak fidelity is limited only by imperfect adiabaticity. Finally, the same simulation is also performed (not shown) without damping but with added disorder, such that the first site has a frequency that is 2 Â 0:5 MHz (2 Â 1 MHz) higher than the second site, and the third site has a frequency that is 2 Â 0:5 MHz (2 Â 1 MHz) lower than the second site. In this case, the fidelity peak corresponding to the first peak in Fig. 7 drops to about 96:8% (92:5%) and the second peak drops to about 94:7% (85:9%). As there is no damping, this discrepancy between the first and second fidelity peaks can be understood as being due to the disorder, which causes the phases to evolve at different rates on each side. This effect of the disorder indicates that disorder may place limits on the total time length of the protocol for cases where intersite phase differences would be undesirable.
In the second simulation, we use M ¼ 3, max =2 ¼ 120 MHz, and max =2 ¼ 30 MHz. The input state is jc in i ¼ 1 ffiffi j¼1 jc in i j Q rÞj j0i r during the final step of the entanglement protocol with system parameters max =2 ¼ 100 MHz, max =2 ¼ 30 MHz, and M ¼ 3 and damping as explained in the text. The input state is jc in i ¼ 1 ffiffi 2 p j2i þ 1 ffiffi 2 p j3i. The fidelity oscillations are due to the phases 2 ðtÞ and 3 ðtÞ of the respective number components evolving at different rates because of the number-dependent frequency induced by the Kerr nonlinearity at each site. The frequency of the oscillations decays to zero as is tuned to zero, and the timing Át 7 of the step is chosen such that the oscillations cease at a fidelity peak of about 97:5%.
at each site is truncated at n ¼ 4. The fidelity jhÉðtÞjÉ W ij 2 during step 7 is shown in Fig. 8. The timings of the protocol steps are based on trial simulations that determine the level of adiabaticity needed to achieve the high peak fidelity shown. The total physical time of the second simulation is 0:3556 s.

VII. CONCLUSIONS AND OUTLOOK
We have proposed a circuit QED simulation of the attractive Bose-Hubbard model in which the as-yet experimentally unobserved superfluid-to-W-state quantum phase transition may be realized with existing experimental capability. A unique aspect of this proposal is the capability of the in situ tuning of the hopping energy to enable access to the perfect W-state regime. We have further presented a protocol built around the attractive Bose-Hubbard simulation that deterministically produces W-type entanglement of nearly arbitrary single-resonator states over a large number of microwave resonators in parallel. We have numerically demonstrated our protocol with complex input states in an array of three resonators using realistic parameters and have shown the attainability of high outputstate fidelity with the target state. The highly entangled large-scale states that the protocol is able to produce have both fundamental and applied significance.
Looking ahead, considering the equilibrium physics of the attractive Bose-Hubbard model in a circuit QED setup opens new prospects due to the unique flexibility of the platform. Studies of the ABH involving different types of lattice geometries, couplings, spatial modulations of parameters, defects, and controlled disorder that were not feasible in other platforms are a possibility in circuit QED. Further, it is interesting to consider the equilibrium and nonequilibrium phenomenologies that would result from interweaving the ABH and Jaynes-Cummings-Hubbard models in the same circuit. Also, the manyresonator entanglement protocol presented here for microwave resonators may offer a path toward creating large-scale W-type entanglement in mechanical and optical degrees of freedom because of the potential that superconducting circuits hold for interfacing with those platforms. Finally, it is worthwhile considering how lattice geometries and couplings different than the onedimensional, nearest-neighbor coupling case considered here could allow the protocol to scale to larger numbers of resonators.
FIG. 8. Fidelity of the system state jÉðtÞi with the target state j¼1 jc in i j Q rÞj j0i r during the final step of the entanglement protocol with system parameters max =2 ¼ 120 MHz, max =2 ¼ 30 MHz, and M ¼ 3 and damping as explained in the text. The input state is jc in i ¼ 1 ffiffi 6 p j2i þ 2 ffiffi 6 p j3i þ e i=9 1 ffiffi 6 p j4i. The fidelity oscillations are due to the phases of the number components evolving at different rates because of the number-dependent frequency induced by the Kerr nonlinearity at each site. Each period of the oscillations has three peaks because of the three frequency differences of the three number components of jc in i. The frequency of the oscillations decays to zero as is tuned to zero, and the timing Át 7 of the step is chosen such that the oscillations cease at a fidelity peak of about 95%.