Kinetically Constrained Quantum Dynamics in Superconducting Circuits

We study the dynamical properties of the bosonic quantum East model at low temperature. We show that a naive generalization of the corresponding spin-1/2 quantum East model does not posses analogous slow dynamical properties. In particular, conversely to the spin case, the bosonic ground state turns out to be not localized. We restore localization by introducing a repulsive interaction term. The bosonic nature of the model allows us to construct rich families of many-body localized states, including coherent, squeezed and cat states. We formalize this finding by introducing a set of superbosonic creation-annihilation operators which satisfy the bosonic commutation relations and, when acting on the vacuum, create excitations exponentially localized around a certain site of the lattice. Given the constrained nature of the model, these states retain memory of their initial conditions for long times. Even in the presence of dissipation, we show that quantum information remains localized within decoherence times tunable with the parameters of the system. We propose an implementation of the bosonic quantum East model based on state-of-the-art superconducting circuits, which could be used in the near future to explore dynamical properties of kinetically constrained models in modern platforms.

In this work, we explore the low-temperature dynamical properties of the bosonic quantum East model, a gen-eralization of the spin-1/2 model studied in Refs. [26,59], in which spin excitations can only be created on sites to the "east" of a previously occupied one. Our contributions can be summarized as follows. (i) We show that repulsive density-density interactions are necessary to entail localization in the ground state, in contrast to East models with a finite-dimensional local Hilbert space. (ii) We exploit the properties of the localized phase and the bosonic nature of the model, to construct families of non-Gaussian many-body states that are useful for quantuminformation processing. (iii) We illustrate how localization enhances the robustness of these states against decoherence. (iv) Finally, we propose an implementation of the bosonic quantum East model based on chains of superconducting qubits.
In the spin-1/2 case, evidence has been provided in support of a dynamical transition from a fast thermalizing regime to a slow, non-ergodic one [26,59]. In particular, in Ref. [59], it has been argued that the slow dynamics is a byproduct of the localized nature of the low-energy eigenstates of the model. Namely, the corresponding wavefunctions contain nontrivial excitations only on a small compact region of the lattice and they are in the vacuum state everywhere else. This has direct consequences for the dynamical properties of the system, as the localized states can be used as building blocks to construct exponentially many "slow" states in the size of the system.
The dynamical transition observed in Ref. [59] is not guaranteed to survive in the bosonic case. In fact, we provide strong numerical evidence that this is not the case for the most naive bosonic generalization of the spin-1/2 model. In order to restore localization at low temperature, we consider a modified model in which densitydensity interactions -absent in the bare spin case -play a crucial role. More precisely, we show that the ground state remains localized as we increase the finite cutoff of the local Fock-space dimension only in the presence of arXiv:2112.08387v2 [cond-mat.stat-mech] 8 Jun 2022 repulsive interactions. We support our findings by combining numerical and analytical approaches. Within the localized phase, the ground state is well approximated by a product state for any value of interaction. It is therefore well approximated by a matrix product state, making large system size and local Fock space dimension numerically accessible (cf. Secs. II and III).
The bosonic generalization of the spin-1/2 East model opens up a number of directions including the construction of many-body versions of archetypal states that are relevant for quantum information applications such as coherent states, squeezed states, and cat states [60]. These states possess the same properties as their singlemode counterparts, although they are supported on a few neighboring sites. We provide a formal description of these objects by proposing a simple adiabatic protocol that defines a set of superbosonic creation-annihilation operators (Sec. IV). These operators fulfill the canonical bosonic commutation relations and they are exponentially localized in the neighborhood of a given site on the lattice. This allows us to construct an effective, noninteracting, theory at low temperature in terms of these operators, in which the Hamiltonian is reminiscent of the l-bit construction in many-body localization (MBL) [61][62][63][64].
In Sec. V, we couple the system to different noise sources and, via a detailed numerical analysis, we show that localized states retain some memory of their initial condition even in the presence of strong dissipation (see Fig. 1). First, we consider the effects of dephasing noise coupled to bosonic occupations, which preserves the "East symmetry" (see the definition in Sec. II). In this scenario, the localized states are barely altered by the environment. We show that the fidelity between the time-evolved state and the initial state decays exponentially with a long decoherence time, controlled by the parameters of the Hamiltonian, the initial state, and the strength of the noise. Second, we consider the effects of particle losses that break the "East symmetry." As expected in this situation, the magnitude of the fidelity decays exponentially fast in time, with a decoherence time that is parametrically small in the loss rate. It is important to stress that as the localized states have non-trivial structure only on a small support, any external noise that does not act in their immediate vicinity leaves them essentially invariant. This set of noise-resilient properties renders the many-body states studied in this work qualitatively different from localization induced by disorder, which is inherently fragile to decoherence (for studies on MBL systems coupled to a bath or external noise see Refs. [65][66][67][68][69][70]). In particular, in Sec. VI we argue that our localized states can be manipulated on timescales shorter than the characteristic relaxation and decoherence times of superconducting qubit wires.
In fact, our proposal for an implementation of the bosonic quantum East model based on superconducting qubits is one of the key findings of this work. In recent years, unprecedented quantum control of interacting superconducting qubits with microwave photons has In the red box we write the lowenergy effective interaction between the j-th and (j + 1)-th superconducting qubits. (b): A sketch of a localized state subject to external noise (arrows). The visibility of the initial peak with respect to the rest of the system (measured by the imbalance I(t)) decays exponentially with a time τ much larger than the characteristic operational timescales of stateof-the-art superconducting circuits.
been reached in circuit-QED platforms [71][72][73][74][75][76][77][78][79]. These circuits allow quantum-information-processing tasks and the quantum simulation of paradigmatic light-matter interfaces. Superconducting Josephson junctions allow us to introduce nonlinearity in quantum electrical circuits, which is a key factor in protecting quantum resources, by making these platforms resilient to noise and errors. This is a key factor of merit for any superconducting qubit, ranging from the established transmon to, for instance, the more recently developed superconducting nonlinear asymmetric inductive element (SNAIL) [80,81]. Here, we consider a chain of superconducting qubits (see Refs. [80,[82][83][84][85][86][87][88][89]), which can be described as anharmonic oscillators, coupled via a hopping term (cf. Fig. 1). In the limit of weak coupling and low anharmonicity, we find an effective description of such superconducting qubits array in terms of the bosonic quantum East chain.
The paper is organized as follows. In Sec. II, we introduce the Hamiltonian of the model, enumerate its symmetries, and compare it to previous works on similar models. In Sec. III, we explore the localization properties of the ground state of the model. In particular, we show when the transition point is independent of the size of the cutoff of the local Fock-space dimension and how the localization length behaves in the proximity of the transition. On the localized side of the transition, we quantitatively compare results extracted with tensor-network methods and mean field, and we show that they are in excellent agreement. In Sec. IV, we introduce a description in terms of superbosonic operators, which allows us to generalize coherent, squeezed, and cat states. In Sec. V, we study the robustness of these localized states against noise source. In Sec. VI, we present the implementation of the Hamiltonian for the bosonic quantum East model, based on a chain of superconducting qubits.

II. BOSONIC QUANTUM EAST MODEL
We investigate the following Hamiltonian with open boundary conditions whereâ j andâ † j are bosonic annihilation and creation operators acting on site j respectively; e −s controls the constrained creation and annihilation of bosons; is the on-site density-density interaction; and U is the nearestneighbor density-density interaction.
As discussed in Sec. I, Eq. (1) is a kinetically constrained "East" model. The unidirectional constrained feature has consequences for the accessible portion of the Hilbert space by the dynamics. Namely, any initial state with a product of vacua from the left edge up to a given site in the bulk will exhibit nontrivial dynamics only on the right side of the lattice after the first occupied site. For sake of concreteness, let us consider the state |00100 . . . 0 . Via subsequent application of the Hamiltonian given in Eq. (1) we have,   where → represents the action of the constrained creation and annihilation of bosons at each step of perturbation theory. The occupation of the first nonvacant site and of those at its left cannot change as a consequence of the "East" constraint. More formally, the Hamiltonian commutes with the projectors where P s,j = |s j j s| is the projector on the Fock state with s particles on site j, 1 j is the identity acting on site j, and k and n 0 are, respectively, the position and occupation of the first nonvacant site. We can split the Hilbert space into dynamically disconnected sectors H n0,k , such that the action of P (n 0 , k) is equivalent to the identity, while the action of the other projectors gives zero. For example, the state |00100 . . . 0 ∈ H 1,2 (note that the first site index is 0). Furthermore, since L k=0 ∞ n0=1 P (n 0 , k) = 1 these sectors {H k,n0 } constitute a complete and orthogonal basis of the whole In the following, we focus on a certain block specified by k, n 0 , and the number of "active" sites L right next to the k-th one. Since the action of H on sites to the left of the k-th one is trivial, the index k is physically irrelevant for our purpose and we therefore choose k = 0 without any loss generality. Exploiting this property, we write the Hamiltonian given in Eq. (1) withĥ 1 ≡ − 1 2 n 0 e −s â 1 +â † 1 − n 0 − Un 1 − 1 and n 0 ∈ N + . Furthermore, since H L+1 (n 0 ) commutes with the operators acting on the (L + 1)-th site, we can represent it as the sum of an infinite number of commuting terms H L+1 (n 0 ) = βr H βr L (n 0 ) ⊗ Π βr L+1 , where Π β L is the projector over the eigenstate |β r with eigenvalue β r = rU − e −2s /U of the operator In Sec. III, we focus on the properties of the ground state of the Hamiltonian given in Eq. (5) within a certain symmetry sector.
The Hamiltonian given in Eq. (1) can be linked to its spin-1/2 version [59] by setting U = = 0 and replacing the bosons with hard-core ones. Since the Hilbert space of each spin is finite, the "East" symmetry is largely reduced with respect to the bosonic case. Each symmetry sector H k,n0=1 is specified only by the position of the first excitation, since n 0 is bound to be zero or one. The ground state properties within a symmetry sector H k,n0=1 , where the position k of the first nonempty is again irrelevant, have been investigated in Ref. [59]. It has been observed that the probability of finding an occupied site in the ground state decays exponentially fast around the first occupied site when s > 0, namely where the expectation value is taken on the ground state and we introduce the localization length ξ > 0. The localization length ξ is the typical distance from the first occupied site such that the state becomes a trivial product state that is well approximated by the vacuum. In Sec. III, we investigate the conditions for localization for different values of s at fixed nearest-neighbor densitydensity interaction U = 1. We fix L = 15, a cutoff Λ = 30 to the maximal occupation number, and n0 = 1. In the plot, we do not display the occupation n0 of the zeroth site that fixes the "East symmetry" sector. The dashed lines are the exponential fit, the slope of which is −1/ξ, where ξ is the localization length (cf. Eq. (6)).
of the ground state at finite values of s upon trading spins (hard-core bosons) for bosons. Such generalization is not granted. The amplitude for "eastern" particle creation can now be enhanced by the prefactor n 0 , suggesting that the transition may be qualitatively established when (n 0 e −s ) ∼ 1. This would imply a critical value s c ∝ log n 0 , which is parametrically large in n 0 , pushing the extension of the localized phase up to s → ∞. Nonetheless, we show in Sec. III that a localized phase still occurs for s > 0 whenever repulsive interactions are included in Eq. (1).

III. LOCALIZATION TRANSITION
In this section we show that the Hamiltonian in Eq. (5) displays a localization-delocalization transition at finite s and U > 0. We give numerical evidence corroborated by analytical observations that repulsive interactions are necessary to observe such a transition at finite s. We use the inverse localization length ξ −1 controlling the decay of the average occupation number in space (cf. Eq. (6)), as proxy for the transition.
In the following, we fix = 0 and the symmetry sector β r=0 in Eq. (5), unless mentioned otherwise. The additional nonlinear term proportional to would complicate the analysis from a technical standpoint without altering the main contents of the paper. For the sake of clarity, Appendix A shows that, for U = 0 and > 0, the localization properties of the ground state remain qualitatively similar to those discussed in the main text.
In order to investigate the properties of the ground . The cutoff is saturated over many sites. The staggered feature is due to the repulsive nearest-neighbor interaction. In the right panel, we consider a typical localized ground state (s = 0.05). Along each site j, the probability of having k bosons, P k,j , drops exponentially fast with k. The light color means that the value is smaller than 10 −12 .
state, we resort to a combination of mean-field arguments, exact diagonalization (ED), and density matrix renormalization group (DMRG) methods [90]. Since we aim to explore large system sizes, we mainly resort to the DMRG and we use ED as a benchmark when both methods can be used. Interestingly, we find that mean field is able to analytically predict the location of the transition point obtained via the DMRG. We compute the ground state |ψ 0 (n 0 ) at fixed n 0 , s, and U . We fix the system size at L = 15. This value is sufficiently large to capture the localized tail of the ground state, without relevant finite-size effects. Although the local Fock space is infinite, in order to treat the model numerically, we need to fix a finite cutoff Λ. We work with Fock states |0 through |Λ , such that the spin-1/2 case of Ref. [59] is recovered at Λ = 1. In Appendix B, we show how localization is only mildly dependent on the sector selected by the occupation n 0 of the zeroth site. Accordingly, in the following, we set n 0 = 1.
The Hamiltonian is one dimensional, local, and gapped at finite Λ; therefore, its ground state can be efficiently accessed via a matrix product state (MPS) formulation of the DMRG [90]. The main source of error is given by the finite cutoff Λ. Indeed, the properties of |ψ 0 (n 0 ) can change nontrivially as a function of Λ. More precisely, for any finite cutoff Λ, the model falls into the class of localized systems studied in Ref. [59]. As a result, |ψ 0 (n 0 ) is always localized for a large enough s at finite Λ but this does not imply localization for Λ → ∞. Indeed, although U > 0 makes the spectrum of the Hamiltonian in Eq. (1) bounded from below, it does not ensure that its ground state is still localized in space when s is finite. In the following, we extract the Λ → ∞ limit via a scaling analysis. In Fig. 2, we show the average occupation number n j as a function of site j for some values of s at fixed U = 1. For s not large enough, the average occupation does not change smoothly with the site j and it saturates the cutoff Λ, meaning that there are strong finite-cutoff effects. In contrast, for s large enough, the occupation decays exponentially in j, matches Eq. (6) well, and does not change upon increasing the cutoff Λ. The value of s at which this change of behavior occurs depends on U , as we discuss in more detail in this section.
In order to check the effects of a finite Λ cutoff, we compute the probability of having k bosons on site j, namely the expectation value of the projector P k,j = |k j j k|, where |k j is the Fock state with k particles on site j. In Fig. 3, we show P k,j as a function of k and j for typical localized and delocalized ground states, respectively. The results in the delocalized phase are not reliable, since the observable suffers finite-cutoff effects. Instead, in the localized phase, with ξ F,j > 0 for any site j. The exponential decay in the localized phase sheds additional light on the fact that the system is well described by a finite effective cutoff (for additional details, see Appendix C).
For each value of U and Λ, the inverse of the localization length goes from values smaller than or equal to zero to positive values as s increases. We identify the region where 1/ξ ≤ 0 as the delocalized phase, while the region where 1/ξ > 0 is identified as the localized phase. In the delocalized phase, strong finite cutoff effects can lead to a positive localization length ξ. In order not to mistakenly identify these points as belonging to the localized phase, we fix a threshold λ > 0 and for each Λ and U we identify the transition point s c (U, Λ) as the value of s such that 1/ξ ≤ λ and 1/ξ > λ for s smaller and greater than s c (U, Λ), respectively. We choose λ ≈ 10 −1 . The results are weakly affected by this choice of λ. Furthermore, the precise location of the transition point s c (U, Λ) is beyond the scope of this work, since we are interested in engineering states deep in the localized phase, as we discuss extensively in Sec. IV.
As discussed above, in the delocalized phase, results are strongly dependent on the cutoff, since the average occupations always saturate their artificial upper bound. This circumstance allows us to draw only qualitative conclusions on the physics at s < s c in the case of the bosonic East model (Λ → ∞).
In Fig. 4, we show the inverse of the localization length ξ swiping s for different values of Λ at fixed U . For U = 0, the transition point s c (U = 0, Λ) always increases with Λ. Instead, when U > 0, the transition point converges to a finite value independent of Λ for Λ → ∞. In Fig. 4.(a), we show the numerically extracted transition point s c (U, Λ) as a function of Λ and U . For U > 0, it is possible to extract a finite value of s c (U ) ≡ lim Λ→∞ s c (U, Λ). Instead, for U = 0, the transition point scales as s c (U = 0, Λ) ∝ log(Λ), suggesting that in the actual bosonic system we have s c (U = 0) = ∞, meaning that there is no transition. Therefore, whenever U > 0, the system undergoes a delocalized-localized transition at finite s c (U ). In Fig. 5, we show the inverse of the localization length ξ as a function of s for different values of U at fixed Λ. The transition point s c depends on the competition between the dynamical term, controlled by e −s , and the nearest-neighbor density term, proportional to U . The former favors the delocalization of the state, while the latter favors its localization. Indeed, in the U → 0 limit, we provide evidence that the bosonic system is always delocalized if s < ∞. Instead, in the large U limit, the Hamiltonian is approximated by U in ini+1 +n i , the ground state of which in a specific symmetry sector at given total particle number is simply |n 0 |00 . . . 0 .
The role of the interaction term U in the localization of the bosonic system can be appreciated in a mean-field treatment. We project the Hamiltonian into the manifold of coherent product states |φ = L j=1 |α j j , witĥ a j |α j j = α j |α j j . We evaluate the Hamiltonian given in Eq. (4) in this basis: where |α j | 2 is the average number of particles in the coherent state at site j.
From unidirectionality of the interaction, we can write φ|H( For energetic stability the effective field h j (α j+1 , s, U ) on site j should be negative: Since the system does not conserve the number of particles there can be an unbounded number of excitations in the ground state within a fixed symmetry sector. Therefore, in order to have localization at a mean-field level it is necessary that Eq. (9) holds for any value of α j+1 ∈ [0, ∞), namely s > max αj+1 s c (α j+1 ), and for all sites. For U > 0, such condition is satisfied if , which turns to be in very good agreement with the DMRG numerical findings (see Fig. 5). Instead, for U ≤ 0, there is no finite value of s that fulfills Eq. (9) for all α j+1 .
The excellent agreement between the DMRG and the mean-field analysis can be explained by observing that the ground state |ψ 0 (excluding the zeroth site, which fixes the symmetry sector) obtained via the DMRG is well approximated via a product state, namely |ψ 0 ≈ L j=1 |φ j . To further investigate the nature of the state |ψ 0 , we consider the correlator ∆ j ≡ ( n jnj+1 − n j n j+1 ). We use this operator as a proxy for non-Gaussian correlations. We compare ∆ j computed on the ground state obtained via the DMRG and the one computed assuming that the same state is Gaussian in the operators {â ( †) j } L j=1 , using Wick's theorem. As shown in Appendix D, the closer we are to the transition point s c , the more the state develops non-Gaussian features at distances j ξ. On the contrary, deep in the localized phase, the Gaussian ansatz captures the actual correlations at all sites well. Indeed, in the large s limit, the Hamiltonian turns out to be diagonal in the number basis, namely H(s 1) ∼ j (n jnj+1 +n j ), the ground state of which is |n 0 |00 . . . 0 , which is a product state of Gaussian states (excluding the zeroth site, which fixes the symmetry sector).
The localized tail can be explained in a more intuitive way via the adiabatic theorem. Indeed, the Hamiltonian is gapped in the localized phase when U > 0; therefore, we can adiabatically connect two ground states within it. In particular, we can link any localized ground state to the one at s = ∞. This choice is particularly convenient since the Hamiltonian is diagonal in the number basis at s = ∞, H(s → ∞) = j=1 (Un jnj+1 +n j )/2 and its ground state at fixed symmetry sector is simply |n 0 L j=1 |0 j . Then, the evolution with the adiabatically changing Hamiltonian will dress the initial site with an exponentially localized tail. In Sec. IV, we further exploit the adiabatic theorem to design the many-body version of a variety of states that are relevant in quantuminformation setups, such as coherent states, cat states, and squeezed states.

IV. LOCALIZED STATES ENGINEERING
In Sec. III, we have discussed the localization properties of the ground state of the bosonic quantum East model within each symmetry sector specified by the occupation n 0 of the first nonvacant site. In this section, we show that the ground states of different symmetry sectors are connected via bosonic creation and annihilation operators. We use this infinite set of localized states to construct the localized versions of cat, coherent, and squeezed states that are relevant for quantuminformation purposes. These states share the same prop-erties as their single-mode counterparts, although they are supported on a few neighboring sites toward the East as the ground states.
Starting with a given symmetry sector fixed by n 0 , our aim is to find operators A and A † that obey the bosonic canonical commutation relations A, A † = 1, with the defining property where |ψ 0 (n 0 ) is the localized tail of the ground state at fixed symmetry sector n 0 and N is a constant. In other words, by acting n 0 times on the bosonic vacuum state with the operator A † , we aim to retrieve the localized ground state of the Hamiltonian in Eq. (1) in the symmetry sector with n 0 particles on the first nonvacant site. From now on, we refer to these operators as superbosonic creation and annihilation operators since, in contrast to single site annihilation and creation operators, they act on a localized region of the system, by creating or destroying a bosonic localized tail along the chain. Likewise, we refer to the localized ground states | n 0 as superbosons.
In order to find an explicit form for such operators, we employ the adiabatic theorem. From numerical evidence our Hamiltonian, is gapped within the whole localized phase (see Fig. 6). Therefore, there exists a slow tuning of s that enables us to connect two localized ground states at fixed values of U and n 0 . We consider such a unitary transformation U(s, U ) linking the ground state for s = ∞ with the target one at s > s c (U ) in a fixed symmetry sector specified by the occupation n 0 of the first nonvacant site. We fix s = ∞ as our starting point since the Hamiltonian is diagonal in the number operator when s → ∞ and its ground state is simply the tensor product |n 0 ⊗ j≥1 |0 j . By the adiabatic theorem, the unitary operator takes the following form [91,92]: where T indicates the time-ordering operator and s(t) is a function that interpolates from s(t = 0) = ∞ and s(t = T ) = s. The function s(t) has to be chosen such that it satisfies [91,92], at all times t. In Eq. (12), the state |Ψ n (t) is the n-th excited eigenstate of the Hamiltonian computed at time t;Ḣ(t) is the time derivative of the Hamiltonian, which encodes the information about the specific protocol; finally, The inset (a) shows the maximum matrix element maxn Vn(s)/n0 ≡ maxn ψn(s)|V |ψ0(s) /n0 of the perturbation V = jn j (âj+1 +â † j+1 ) between the n-th excited state and the ground state at fixed s. We fix a system size L = 6, cutoff Λ = 3 and nearest-neighbor density-density interaction U = 1. The transition point is at sc(U = 1) ≈ 0. The results are weakly affected (of the order of few percent) by the finite cutoff Λ for s 2.
J(t) = −e −s(t) /2. The time derivative of the Hamiltonian then readsḢ(t) =J(t)V . Let us focus on the perturbation V and the gap ∆ at first and then on the specific protocol J(t). In Fig. 6, we show the gap of the Hamiltonian and the maximum matrix element max n V n (s) ≡ ψ n (s)|V |ψ 0 (s) connecting the ground to the n-th excited state as a function of s at fixed U . Within the localized phase, the gap is O(1) and the maximum matrix element max n V n (s) ∼ n 0 , where n 0 is the occupation of the first nonempty site fixing the symmetry sector. Due to the kinetic constraint, the largest matrix element max n V n (s) is between the localized ground state and the second localized state perturbatively close to the product states |n 0 100 . . . (note that this is not necessarily the first excited state). Therefore, the leading contribution comes from the first few sites, since the other terms are exponentially suppressed in the localization length of |Ψ 0 . Let us consider, as a possible adiabatic protocol, the linear ramping J(t) = −e −s t/(2T ), where t ∈ [0, T ], with T as the total duration time. From Eq. (12), the total time T has to satisfy T n 0 e −s . Recall that we set the on-site bare frequency of the bosons as our energy scale and therefore the time T is expressed in that unit as well. In Sec. VI, we propose a possible experimental implementation of the bosonic quantum East model based on superconducting qubits. The typical onsite bare frequency of superconducting qubits is O(GHz), leading to T (n 0 e −s )ns ∼ 1ns, which is within the typical coherence time of O(1µs) of state-of-the-art superconducting qubits [71]. For s(t) that satisfies Eq. (12), we obtain where θ is a phase acquired during the adiabatic time evolution [91,92]. where We can straightforwardly generalize Eq. (14) taking into account the position j starting from which we want to embed the state | n 0 . We define j . Therefore, they are bosonic operators. As anticipated, we call the operators A j (s, U ) ( †) superbosonic annihilation(creation) operators.
Since the transition point s c is essentially independent of the value of n 0 (see Appendix B), we can design a protocol that obeys the adiabatic theorem for any initial state |n 0 ⊗ |0 . . . 0 .
Furthermore, since these states belong to dynamically disconnected symmetry sectors, H k=0,n0 , for any values of s and U , it is possible to adiabatically evolve them independently of each other. Therefore, any linear combination of initial states turns under the adiabatic protocol into where θ(n 0 , s, U, T ) is the phase acquired during the adiabatic time evolution. As discussed in Appendix B, deep in the localized phase the spectrum depends linearly on n 0 , with small corrections. Since the phase acquired during the adiabatic evolution depends on the energy of the given state during the protocol, we have θ(n 0 , s, U, T ) ∼ n 0 f (s, U, T ), where f (s, U, T ) is a function that is dependent on the specific protocol. This has important consequences for the state engineering we discuss in the following. As an example, let us consider as initial state of the adiabatic preparation the coherent state |α ≡ |α 0 j≥1 |0 j , where Using Eq. (15), the state |α turns into where α = αe if (s,U,T ) . In Fig. 7, we compute the overlap between U(s(t), U )|α and the superbosons We expect that when α is large, the fidelity achieved by the protocol becomes small, since corrections to the linear dependence of θ(n, s, U, T ) from n become important. We call the localized version of a coherent state | α ≡ U(s, U )|α a supercoherent state. Analogously, we perform the same analysis considering as initial state a cat state |C on site j = 0. Indeed, since the phase factor e if (s,U,T ) does not depend on α, given a cat state where N is a normalization constant, its localized version is where | C ≡ U(s, U )|C , and α = αe if (s,U,T ) . We call | C a supercat state. We can extend Eq. (17) to states of the form where ρ n ∈ R and β, θ ∈ C. Indeed, if we apply the adiabatic protocol to the state defined in Eq. (20), the phase acquired can be absorbed into β. Coherent states, cat states, and squeezed states all fall into the class described in Eq. (20). In other words, using the adiabatic protocol, not only can we engineer the localized versions of states such as coherent and squeezed states but we can do so preserving their single-mode properties. For instance, the localized versions of coherent and squeezed states can be implemented either via the adiabatic time evolution or the application of an operator M that is linear or quadratic in the superbosonic operators A. The operator M can be obtained applying the adiabatic protocol to its single-site counterpart M , namely M = U(s, U )M U(s, U ) † . For instance, we define the dressed displacement operator, where α ∈ C is the displacement parameter, and the dressed squeezed operator, where ζ ∈ C is the squeezing parameter, the action of which on the vacuum creates a supercoherent and supersqueezed state, respectively. However, the most natural way to prepare such states is by starting from their single-mode version and then adiabatically turning on the off-diagonal term ∝ e −s in the Hamiltonian. Note that these states are Gaussian with respect to the superbosonic operators A ( †) and not with respect to the bare operatorsâ ( †) . We call these states super-Gaussian.
We find that superbosons | n 0 , with different n 0 and the same position j of the first nonvacant site, are connected via the operators A ( †) j . We see that their localized feature makes their energies approximately evenly spaced as a function of n 0 (cf. Appendix B). The evenly spaced energies of different ground states and the fact that the different ground states are connected via a bosonic operator A j (s, U ) ( †) resemble the features of a quadratic Hamiltonian, such as the one-dimensional harmonic oscillator. Adding up these properties, the action of the interacting Hamiltonian H(s, U ) in Eq. (1) in the manifold of the ground states is approximatively equivalent to a free theory in the superbosonic operators A j (s, U ) ( †) , namely the eigenstates of which are  23) captures the equilibrium properties in the localized phase and the dynamical features of states such as the supercat state and supersqueezed state well when the interacting part bewteen superbosons can be neglected. In this regard, the properties of the localized phase of quantum East models are reminiscent of the l-bits construction in MBL [61][62][63][64]. Let us consider a supercat state |ψ(t = 0) = | C defined in Eq. (19) as initial state in order to test the effective quadratic theory in Eq. (23). We evolve it and compute the fidelity As shown in Fig. 8, the fidelity displays almost perfect oscillations at short times, followed by a drop and almost perfect revivals. The short-time behavior is compatible with a rotation of the supercat state in the dressed phase- | C generated by Eq. (23) as where α(t) = α(t = 0)e −i 0t . The state in Eq. (25) is a rotating supercat state in the dressed space. From Eq. (25) we can estimate the expected fidelity. In Fig. 8, we compare the expected value and the numerical results. The former matches the numerical results up to times parametrically large in s and 1/α. On the one hand, nonlinear corrections are suppressed the more the system is localized. On the other, corrections to the linear dependence of the energies n|Ĥ| n become important the larger n is or, equivalently, α, leading to dephasing processes [93]. The revivals can be explained considering nonlinear effects; indeed, perfect revivals are observed for single-mode cat states with self-Kerr interaction [94] (for a circuit-QED implementation, see Ref. [95]). Differently from the latter case, we have an extended state and nearest-neighbor density-density interactions. As a consequence, pushing the simulations to longer times we observe no perfect revivals as in the case of single bosonic modes with Kerr nonlinearities. Such behavior might be captured by improving the effective theory introduced in Eq. (23), adding nonlinearities in the basis of superbosonic degrees of freedom. This is beyond our current scope and therefore left as a potential interesting follow-up. We can extend these dynamical properties to any state prepared via the adiabatic protocol starting from a state of the form given in Eq. (20). Indeed, these states evolve analogously to the supercat state under the effective quadratic theory defined in Eq. (23). The super-Gaussian states fall into this class. Once again, we highlight that these states are Gaussian with respect to the superbosonic operators A ( †) but not with respect to the bare operatorsâ ( †) .
We have discussed the application of the adiabatic protocol to a single-site state embedded in the vacuum; however, this extends directly to more general initial states. For instance, we could have started from a product state made of single-body states separated by a large number of empty sites, with respect to the localization length ξ, or from a superposition of those. At the end of the protocol, each one will be dressed independently from the others. Therefore, the final state will be made of localized states concatenated one after the other.

V. EFFECTS OF DEPHASING AND LOSSES
In this section, we investigate the dynamical properties of the localized states introduced in Sec. IV when coupled to the environment. Here, we study the effects of two different couplings with an external bath, namely a global dephasing due to a noise coupled to the local densities, which commutes with the "East" symmetry, and global losses, which break the "East" symmetry. Both of these couplings are experimentally relevant in superconducting circuits setups [71], which are at the core of the experimental implementation we propose in Sec. VI. We provide numerical evidence that local information is erased very slowly when the environment is coupled via densities to the system. We show how the characteristic time scales depend on the parameters of the Hamiltonian, the initial state, and the strength of the coupling to the environment. On the contrary, we show that losses are highly disruptive and that the time scales are dependent on the strength of the coupling to the environment and the initial state, while the underlying coherent dynamics does not play a substantial role. At the end of the section, we show that the typical couplings to the environment currently achieved in superconducting circuits are small enough to make the effects of the coherent dynamics appreciable and observable in the presence of losses.
We consider the following Linbland master equation: whereρ is the state of the system,Ĥ is the Hamiltonian in Eq. (1) with = 0,L j is the quantum jump operator acting on site j and γ is the corresponding rate. In order to efficiently simulate the Lindbland master equation in Eq. (26), we resort to the quantum trajectories algorithm, which is based on defining the effective non-Hermitian HamiltonianĤ and alternating the action of the Hamiltonian given in Eq. (27) with the jump operators {L j } based on a stochastic process (for the details, we refer to Refs. [96,97]). The dynamics of any observableÔ result from averaging over N different uncorrelated stochastic trajectories labeled by η ∈ [1, N ], where |ψ(t) η is the state for a given stochastic trajectory η ∈ [1, N ] at time t and · η denotes the average over the different trajectories. We resort to tensor-network methods for performing the simulations (see Appendix E). We consider two different jump operators, namelyL j =n j andL j =â j . The former corresponds to dephasing, while the latter corresponds to losses. We choose such jump operators in order to investigate the effects of the environment when it preserves the "East" symmetry, as for the dephasing process, or when it does not, as for the global losses. Both situations are relevant in superconducingcircuit setups [71]. We compute the observables averaging over 1000 to 3000 stochastic realizations depending on the value of γ and the jump operator.
We study the dynamical properties of superbosons | n defined in Eq. (10), since they constitute the building blocks of any localized state that we can engineer. Then, we turn our attention to a paradigmatic superposition of superbosons, namely the supercat state, providing arguments to extend our findings to a class of states to which supersqueezed and supercoherent states belong. We consider as initial state |ψ k (t = 0) = k−1 j=−∞ |0 j ⊗ | n , where the subscript k in |ψ k (t = 0) refers to the position of the first site of the embedded superboson. Since | n is localized with localization length ξ (cf. Eq. (6)), we can truncate its support to L ξ sites. Thus, our initial state is where L is the size of the superboson support. In a generic non-integrable system, we expect information about initial states encoded in local observables to be washed out fast. Here, we want to study how localization and slow dynamics instead protect the information encoded in local quantities. We compute the fidelity and the imbalance. The fidelity (cf. Eq. (24)) provides global information about the state and sets an upper bound on the time dependence of the expectation value of any local observable. Nonetheless, the fidelity is highly sensitive to any local perturbation of the state. Indeed, it is enough to have even a single occupied site far from the superbosons | n to make Eq. (24) negligibly small. Among all the possible local observables, we want to investigate if the initial localized peak remains well resolved. We therefore compute the imbalance between the occupation of the initial peak and the second highest peak in the system, namely where k is the position of the first site of the embedded state (cf. Eq. (29)). The imbalance I ∈ [−1, 1] and for I > 0 the initial peak is the largest one in the system. When dissipation enters in the form of a dephasing noise coupled to the bosonic densities, the Lindbland equation respects the "East" symmetry. The jump operators commute with the operator in Eq. (3). Thus, the n excitations on the first site of the superbosons | n and the empty sites to its left are conserved. Furthermore, since the the jump operators are not able to generate excitations out of the vacuum and the state is exponentially localized, we can keep only a few empty sites to the left of | n L without introducing relevant size effects. For the set of parameters that we choose, restricting the superboson support to L ≈ 10 sites and keeping only one empty site to its right turns out to be sufficient. Thus, our initial state is In Fig. 9 we show the dynamics of the fidelity and imbalance for different values of s and noise strength γ keeping U = 1, starting from the state in Eq. (31) with n 0 = 1. The imbalance displays an exponential decay I(t) ∼ I(0)e −t/τ , with τ dependent on the initial state, the parameters of the Hamiltonian, and the coupling strength γ with the external bath. The decay time τ increases the more the system is in the localized phase and the larger is the initial occupation n 0 , while it decreases with the noise strength γ as τ ∝ 1/γ. Therefore, the time decay τ can be enhanced by either tuning the parameters of the Hamiltonian or embedding a superboson with n 0 large (cf. Eq. (31)). On the one hand, increasing s or U helps to protect the local memory at longer times, at the cost of making the initial state less entangled. Indeed, in the s, U → ∞ limit, the Hamiltonian tends to ∝ i (U n i n i+1 + n i ), the ground state of which is a product state of eigenstates of number operators. On the other hand, we can exploit the bosonic nature of the system and embed a superboson with a larger initial n 0 , keeping s small and enhancing the initial state entanglement. It is important to stress that despite the exponential feature of the decay, the time scale τ is generally very large with respect to the time scales of the coherent dynamics of the system. From Eq. (30), and inspecting the late times average occupation number, the initial peak remains still well resolved and so does the information encoded within it. The fidelity decays exponentially fast in time F(t) ∼ e −t/τ , with a decoherence time τ dependent on the parameters of the Hamiltonian, the initial state, and the strength of the noise. Analogously to the decay time τ of the imbalance, the decoherence time τ increases the more the system is in the localized phase and decreases with the noise strength γ as τ ∝ 1/γ. Contrary to the imbalance, the fidelity drops faster the larger is n 0 . Indeed, the conserved initial occupation n 0 pumps excitations on the next site, reducing the typical coherent time scales by approximatively 1/(n 0 e −s ) and effectively enhancing the effects of the environment. Under the action of single-body losses, the dynamics no longer preserve the "East" symmetry. Indeed, losses can deplete the occupation of the first site, which fixes the "East" symmetry sector. Since the vacuum is invariant under the action of losses and coherent dynamics cannot create excitations to the West of the initial embedded superboson, we can still consider Eq. (31) as our initial state. In Fig. 10, we show the dynamics of the fidelity and imbalance for different values of n 0 , keeping U = 1, s = 1.5 and γ = 0.1 fixed. Losses turn out to be detrimental to the initial state independent of the parameters of the Hamiltonian. Instead, the height of the initial peak plays a substantial role in enhancing the conservation of the imbalance. Intuitively, if the first site n 0 is highly occupied at time t = 0, it will require longer times to drain all the particles. This leads to an initial plateaux in the imbalance, followed by an exponential decay toward the minimum value I(t → ∞) = −1.
The decay is well fitted by I(t) = Ae −t/τ − 1 at long times, where τ ∝ 1/γ is the relaxation time and A is a constant. The insensitivity of the time decay with respect to the parameters of the Hamiltonian indicates that the slow dynamics do not provide additional protection against this type of coupling to the environment. Indeed, the decay of the imbalance is due to the emission of particles from the first occupied site, which fixes the symmetry sector, and since the coherent dynamics cannot create excitations on top of it the initial peak is depleted in time ∝ 1/γ. The fidelity drops to zero exponentially fast, as expected, with a decay time that is parametrically small in the occupation of the initial peak. Indeed, the higher the peak is, the larger is the probability that the emission occurs, which immediately produces a state orthogonal to the initial one. Despite losses being more detrimental with respect to dephasing, we show at the end of the section that the coherent dynamics takes place on time scales that are small with respect to the relaxation time in typical superconducting circuits (cf. Sec. VI for the experimental implementation of the bosonic quantum East model).
Note that we can immediately extend our analysis to a large variety of states. For instance, we can consider states given by the superposition of superbosons embedded in different regions of the systems, namely where |ψ k (t = 0) is defined in Eq. (29), θ is a phase, and |s − k| ξ. These two states are weakly coupled by the coherent and dissipative dynamics. In a first approximation, we can apply our analysis to each of them separately, and therefore predict their dynamics easily.
The extension of these results to superposition of superbosons embedded in the same support (cf. Eq. (15)) is less trivial and depends on the specific coupling to the environment. For instance, a coupling that does not preserve the "East" symmetry makes the different states dynamically connected, likely leading to different results from the ones observed for the single superbosons. On the other hand, a coupling which preserves the "East" symmetry can also lead to additional phenomena such as dephasing processes between the superimposed states.
Indeed, we observe that coupling to the densities is also detrimental. We give further details in Sec. V A, exploring the effects of local dephasing in the system.

A. Local dephasing
We now investigate the effects of local dephasing in the dynamical properties of a state given by the superposition of superbosons embedded in the same support. Among the possible choices, we consider a paradigmatic super-Gaussian state, namely the supercat state, and then we generalize.
We consider local dephasing due to noise coupled to the densities (see e.g. [98]). In the case of local dephasing acting on a compact support S, the effective theory in Eq. (27) turns intô where the summation is along the support S. We con-siderL j =n j as jump operator.
We study the impact of the dephasing as a function of the strength γ and the extension of its support S. Since the dephasing preserves the "East" symmetry, we can once again focus on system comprising a few sites without introducing relevant finite-size effects. We initialize our system in the state where | C L is a supercat state (cf. Eq. (19)) with support L and average number of particles |α| 2 . A support of L = 10 turns out to be large enough for the parameters explored (α = 1.50, s = 1.5 and U = 1). In Fig. 11 we show the dynamics of the fidelity as a function of the coupling strength γ and support S. The supercat state is still localized in space for any γ and S. Nonetheless, the coherence of the state is highly dependent on γ and S. Indeed, local dephasing is highly disruptive in an exponentially localized region around the peak, where the state is mostly located. If, instead, the local dephasing acts on a region far from the localized peak it does not produce any appreciable effect. More precisely, we estimate that the typical time τ at which the embedded state is appreciably affected by the noise scales as τ ∼ min |k−j|∈S 1/(γ n j ) ∼ min |k−j|∈S e |k−j|/ξ /γ, where k is the site where the peak is located. We numerically verify the polynomial dependence of τ on γ. On the contrary, it is not possible to extract the dependence on the support S with high enough accuracy from the times explored, because of the slowness of the decay. Our findings can be extended to other channels that do not necessarily preserve the "East" symmetry. For instance, losses acting far from the localized peak will not affect local information encoded in the localized state. Furthermore, we expect that the observed dynamical properties can be easily extended to any state prepared via the adiabatic protocol from a state of the form given in Eq. (20), to which super-Gaussian states belong.
In this section we have discussed the effects of dephasing and losses, without much emphasis on the actual value of the coupling strength γ to the environment in typical superconducting circuits (cf. Sec. VI for the implementation). As previously mentioned, we set the on-site bare frequency of bosons as our energy scale, which is O(GHz) in typical superconducting circuits [71]. The typical strength of the coupling to the environment γ is O(MHz) [71]. Therefore, γ ≈ 10 −3 in our nondimensional units. As a consequence, coherent dynamics take place on smaller scales with respect to the operational times of typical superconducting platforms of O(1µs), hinting that the physics of localized states is potentially observable in state-of-the-art experiments. Corroboration of this statement with more quantitative calculations would require an ab initio study of the dynamics of the architecture introduced in Sec. VI, which constitutes an interesting follow-up project per se.

VI. SUPERCONDUCTING CIRCUIT IMPLEMENTATION
In this section, we propose an experimental imple-mentation of the Hamiltonian in Eq. (1) in terms of a simple superconducting-circuit setup. We consider a chain of driven superconducting qubits. A superconducting qubit is basically a quantized LC oscillator with capacitance C and nonlinear inductance L [71]. This nonlinear dependence can be achieved via a Josephson junction working in the superconducting regime without introducing undesired dissipative effects [71,99,100]. In particular, we consider here the SNAIL introduced in Ref. [81] as our building block. We consider specifically the SNAIL parameters in Ref. [101], where kinetically constrained terms (at just two sites) are obtained using the second-order nonlinearity ∝ (â †â † a + h.c.) of the SNAILs. Differently from Ref. [101], we do not use the second-order nonlinearity of SNAILs. Indeed, any superconducting qubit that can be approximated as an anharmonic oscillator with positive anharmonicity could be a suitable candidate for our setup (e.g., the C-shunt flux qubit [102]).
We consider an array of L driven superconducting (SC) qubits coupled via an exchange interaction as our starting point. We retain all the energy levels of each SC qubit. The Hamiltonian can be decomposed as a sum of three terms, whereâ † j (â j ) creates (destroys) an excitation in the j-th SC qubit; H 0 is the bare Hamiltonian of the SC qubits with qubit frequencies {ω j } L j=1 , and anharmonicity E C > 0 [71]; H drive describes the action of classical drive fields on the bare SC qubits; and V describes hopping processes that can be engineered by a common bus resonator [103] or a direct capacitance [104]. An illustration of the scheme of Eq. (35) is given in Fig. 1(a). We work in the weak-coupling regime g |ω j − ω j+1 | and in the low-anharmonicity limit E C |ω j − ω j+1 | for all j. The former condition is necessary in order to have far-detuned processes connected by V , and therefore to treat V in perturbation theory [105]. The low-anharmonicity limit is necessary to retrieve a bosonic model in the effective perturbative Hamiltonian achieved after treating V with a Schrieffer-Wolf (SW) transformation in the small g limit. Each SC qubit j ∈ [1, L − 1] is driven by a classical drive field of amplitude Ω j and frequency α j . These classical drive fields give rise to the desired interaction together with undesired single-site fields in the low-energy effective Hamiltonian [106]. In order to get rid of them, we add another drive field on each SC qubit j ∈ [2, L] of amplitude j and frequency α j−1 [107,108].
We are interested in exploiting the multilevel (bosonic) structure of SC qubits. We do not reduce each component of the system to a qubit. We therefore introduce the ladder operatorŝ whereĉ ,j is the ladder operator which destroys an excitation in the ( + 1)-th level and creates an excitation in the -th level on the j-th SC qubit. Analogously, we can define its Hermitian conjugate,ĉ † ,j .
We work in the dispersive regime, g ∆ j,j+1 , where ∆ i,j = ω i −ω j . We perturbatively diagonalize the Hamiltonian H 0 + V to second order in g via a SW transformation S [109]. The drive field terms in H drive are modified by the same SW transformation. From now on, we neglect terms of order O(g 2 Ω j /∆ 2 j,j+1 ) and higher. We move to the frame that rotates at the frequencies of the drives and we neglect the fast oscillating terms by employing the rotating-wave approximation (RWA). Before detailing the calculations, we discuss the physics of each term in the Hamiltonian defined in Eq. (35). The bare Hamiltonian H 0 provides the necessary anharmonicity that we desire. The perturbation V gives rise to the nearest-neighbor interaction, a renormalization of the bare energies of the SC qubits, and some additional two-excitation processes. The drive field yields the constrained termsn j (â j±1 +â † j±1 ) toward "East" and "West". The time dependence of the drive fields in the laboratory frame enables us to get rid of the undesired processes, such as the two-excitation processes and the "West" terms, passing in the rotating frame of the drive fields and employing the RWA.
In order to find the explicit form of the SW transformation, we follow the prescription in Ref. [110]. First, we compute η = [H 0 , V ]; we consider η with arbitrary coefficients as an ansatz for S. Finally, we fix these coefficients, imposing the condition [S, H 0 ] = −V . We obtain (cf. Appendix F 1) where∆ ,j = (ω j +E C ), the first summation is along the system, while the second summation is along all the levels of the SC qubits. Using the Baker-Campbell-Hausdorff expansion, the Hamiltonian in Eq. (1) after the SW transformation reads After lengthy yet standard calculations, we obtainH explicitly dependent on the ladder operatorsĉ ,j introduced in Eq. (36) and with coefficients dependent on the site and internal levels (see Appendix F 2). Our aim is to writeH as a function of the bosonic operatorsâ ( †) j . We need to find a regime in which the coefficients inH are approximately independent of the specific level, so that we can use Eq. (36). These coefficients are similar to the one appearing in Eq. (37). In order to make them level independent, we need which holds if | − s| |∆ j+1,j |/E C . Since the SC qubit can have an infinite number of excitations, we have ( − s) ∈ (−∞, +∞). This means that Eq. (39) cannot be satisfied for all possible and s if E C = 0. Nonetheless, it can be achieved up to a certain value N of and s, such that N |∆ j+1,j /E C |. Therefore, the coefficients inH satisfy Eq. (39) up to the N -th energy level, leading to a bosonic Hamiltonian that approximates the action of the full Hamiltonian to states with an occupation that is small with respect to N (cf. Appendix F 3). The bosonic H still displays undesired processes, such as hopping and local fields. We move to a rotating frame of reference via the unitary transformation and we neglect all the oscillating terms by employing the RWA (cf. Appendix F 4). In doing so, we get rid of almost all the undesired processes except for some local fields at the sites j ≥ 2. These fields can be eliminated via the additional drive fields of amplitudes { j }, analogously to what has been done in similar scenarios (see, e.g., Refs. [107,108]). We tune their amplitudes such that they cancel the undesired local terms. We obtain the matching condition j = gΩ j−1 /∆ j−1,j , with j ≥ 2. This leads to the effective Hamiltoniañ We now evaluate the couplings in Eq. (41), considering the SNAIL as our SC qubit and using the parameters of Ref. [101]. We work in the parameter regime in which the SNAILs Hamiltonian is given by H 0 in Eq. (35). We fix E C ≈ 150MHz, g = 75MHz and ω j ≈ 3GHz. We consider the classical drive fields with amplitude Ω j = −100MHz (the amplitude has to be negative to have the correct sign for the constrained hopping), which can be achieved by adding a π phase to the external drive fields. Any real system is inevitably coupled to the environment and SC circuits are no exception. In the context of SC circuits, two different time scales are defined, namely T 1 j = 1 j = 2 j = 3 j = 4 j = 5 αj(GHz) 0. 75  and T 2 [71]. The time scale T 1 is the typical time at which the coupling with the environment leads excited states to decay to lower-energy states. The time scale T 2 quantifies the coherence time of the system. For consistency with the chosen parameters (taken from Ref. [101]), we also consider, as T 1 and T 2 , the values from Ref. [101], which are T 1,2 ≈ 1µs. We fix the qubit frequencies ω j and the drive field frequencies α j in order to satisfy: (i) the dispersive regime, valid for g/∆ j,j+1 1; (ii) the low-anharmonicity limit, E C ∆ j,j+1 ; (iii) the validity of the RWA, namely |α j | Ω j , |α j+1 − α j | Ω j and |α j+2 − α j | > gΩ j+2 /∆ j+1,j+2 ; (iv)ω j ≈ ω j − α j−1 > 0 for j > 1, necessary in order to have localization; (v) 1/T 1,2 small with respect to the typical energies in the effective Hamiltonian in Eq. (41); and (vi) the system is in the localized phase. The more stringent conditions are given by (ii) and (v). A good trade-off between (ii) and (v) is obtained at |∆ j,j+1 | ≡ ∆ ≈ 5E C ≈ 750MHz, for which the typical time scale of the kinetically constrained term is approximately T 1,2 /2. We have g/∆ j,j+1 ≈ 0.1, meaning that (i) is reasonably satisfied. Condition (iii) is satisfied by a staggered configuration of the drive field frequencies with an additional dishomogeneity between next-neighbor drive field frequencies, for instance: α j = α j−1 + (−1) j (δ + (j − 1)ζ) for j ∈ [2, 4] and boundary condition α 1 Ω (for larger systems, it is enough to periodically repeat the configuration of the frequencies), with δ Ω, α 1 Ω, and ζ gΩ/∆ ≈ 10MHz. Condition (iv) is satisfied by a staggered configuration of the qubit frequencies as well: ω j+1 = ω j + (−1) j ∆ for j ≥ 2, ω 2 = ω 1 + ∆, with boundary condition ω 1 > α 1 . For instance, we can consider α 1 =750MHz, δ = 750MHz, ζ = 100MHz, and ω 1 = 3GHz. These conditions lead to Eq. (41) being almost translationally invariant (except for dishomogeneities in the frequenciesω j of the order of approximately 5%, which can be eliminated via a more fine-tuned choice of {ω j }). Moreover, condition (vi) is satisfied for these set of parameters. In Table I, we summarize a possible set of parameters available in stateof-the-art superconducting circuits for implementing the bosonic quantum East model.

VII. PERSPECTIVES
The implementation of a kinetically constrained East model using superconducting circuits represents a bridge between the two communities of circuit-QED and nonergodic quantum dynamics. It has the potential to attract the former toward fundamental questions regarding dynamical phase transitions and to stimulate the latter toward the search for quantum-information and metrological applications of constrained dynamics. Our explicit construction of localized analogs of squeezed and cat states relying on the East constraint represents a first stepping stone in this direction.
A fruitful prosecution of this work is the study of an analog of the mobility edge separating localized from delocalized states in the spectrum of East models (for the mobility edge in MBL see Refs. [20,21]). An understanding of how such a mobility edge scales with Λ, is essential for predicting the onset of dynamical transitions in platforms with unidirectional constraints, as well as of practical interest. For instance, a mobility edge at finite energy density is a feature of direct relevance for experimental realizations, since it would yield the conditions for performing efficient quantum manipulations deep in the localized phase when finite-temperature or heating effects are present. A related interesting question is the survival of the effective integrable description of the localized phase discussed in Sec. IV upon increasing the density of energy above the ground state. This would have implications for heat and particle transport features of the East model in the nonergodic phase, which would be governed by the effective integrable description in (23), as it happens for MBL systems [111].
The insensitivity to noise acting away from localized peaks could open up a path toward the study of the pro-tection of spatially separated macroscopic superpositions of superbosonic states. Given the slow decay of localized wave packets in the presence of noise, one could conceive the storage and noise resilience of long-lived many-body entangled states in faraway regions, with applications to quantum communication.
To conclude, we observe that the implementation discussed in Sec. IV may be easily adapted to retain kinetic terms with both East and West symmetries. This could, for instance, lead to the formation of localized modes at edges of the wire, with exciting perspectives for novel forms of topological states in kinetically constrained models that are realizable with circuit QED. We are currently focusing our research efforts in this direction. the role of on-site density-density interactions, focusing on the localization properties of the ground state and comparing with the statements in the main text resulting from numerics performed at U > 0 and = 0. Starting from the Hamiltonian in Eq. (5), we consider U = 0 and ≥ 0. For = 0, the model does not display localization at finite s in the bosonic limit, as extensively discussed in Sec. III. On the other hand, for > 0, the ground state is localized for s > s c in the bosonic limit, with s c being parametrically small in . We perform the same scaling analysis as a function of the cutoff Λ discussed in Sec. III. In Fig. 12, we show the inverse of the localization length ξ swiping s for different values of Λ at fixed . The scaling analysis suggests that the transition point s c (Λ, ) converges to a finite value independent of Λ for Λ → ∞. The overall qualitative picture is therefore unaffected if one considers on-site or nearestneighbor nonlinearities. A nonzero value of introduces, however, anharmonic spacings between ground states with different values of n 0 . Indeed, we have, for the energy of the ground state, E(n 0 ) ≈ n 0 /2 + n 2 0 /2. This additional anharmonicity has an impact on the adiabatic protocol discussed in Sec. IV, since each adiabatically evolved state U|n 0 0 ⊗ j>0 |0 j in Eq. (15) would acquire a phase with a nonlinear dependence in n 0 , which technically complicates state preparation without altering the main physical message. Nonetheless, it is still possible to tame the effect of this nonlinearity by considering a small enough , at the cost of having a smaller e −s (larger s) and therefore working effectively deeper in the localized phase. In this appendix, we discuss the properties of the ground state upon changing the symmetry sector specified by the occupation n 0 of the first nonempty site. We show that the transition point and the exponential decaying tail of the ground state occupation is weakly dependent on n 0 . We discuss the dependence of the ground state energy on n 0 , which is relevant in the state preparation via the adiabatic protocol discussed in Sec. IV.
We perform the same scaling analysis as a function of the cutoff Λ discussed in Sec. III (see Fig. 13). We extract the transition point s c for different values of n 0 from the inverse of the localization length ξ. The existence of a finite critical point s c in the Λ → ∞ limit turns out to be weakly dependent on the specific symmetry sector n 0 at fixed U . We investigate the dependence of the localized tail of the ground state |ψ 0 (n 0 ) as a function of n 0 (we exclude the first site, which fixes the symmetry). To this end, we compute | ψ 0 (n 0 = 1)|ψ 0 (n 0 ) | 2 , with n 0 ≥ 1 (see Fig. 14). We fix n 0 = 1 as a reference as we want to see whether or not the tail is weakly dependent on n 0 . All the ground states are computed by fixing Λ = 30. The overlap | ψ 0 (n 0 )|ψ 0 (n 0 = 1) | strongly depends on s and U . Indeed, the more the system is in the localized phase, the more the exponentially localized  tail is weakly dependent on n 0 . Therefore, deep in the localized phase, |ψ 0 (n 0 ) is approximately independent on the specific sector n 0 and we can write | n 0 ≡ |n 0 ⊗ |ψ 0 (n 0 ) ≈ |n 0 ⊗ |ψ 0 , where |ψ 0 is explicitly independent of n 0 . The weak dependence of |ψ 0 (n 0 ) with respect to n 0 has consequences on the ground state energy. Indeed, the expectation value of the Hamiltonian on Eq. (B1) is E 0 (n 0 ) ≡ n 0 |Ĥ| n 0 ≈ 1 2 n 0 + O(n 0 e −1/ξ(n0) ), (B2) where n j ∼ e −j/ξ(n0) since we are in the localized phase.
In Fig. 15, we give a numerical evidence of Eq. (B2).

Appendix C: Scaling analysis in Λ
In the main text, we show that the bosonic system displays a delocalized-localized transition at finite s if U > 0. Here, we show that the ground state is not only localized but it is weakly dependent on the physical cutoff Λ. This provides quantitative proof that we can investigate the bosonic system with a finite Λ in the localized phase.
We fix the symmetry sector n 0 and (s > s c (U ), U > 0) in the localized phase. We compute |ψ 0 (Λ) for different values of Λ. We calculate 1 − | ψ 0 (Λ)|ψ 0 (Λ + 1) | 2 as a function of Λ (see Fig. 16). The fidelity | ψ 0 (Λ)|ψ 0 (Λ + 1) | 2 approaches 1 exponentially fast in Λ. The more the system is in the localized phase and n 0 is small, the faster is the convergence. This gives the first evidence that the ground state of the actual bosonic system is well described with small effective cutoffs.
We compute the variance of the Hamiltonian given in Eq. (1) over the ground state |n 0 0 ⊗ |ψ 0 (Λ) , taking into account the bosonic nature of the original Hamiltonian in Eq. (1). This quantity is exactly zero if the state |n 0 0 ⊗ |ψ 0 (Λ) is an eigenstate of H. We aim to see how this quantity goes to zero as a function of Λ. In order to do so, we write the Hamiltonian given in Eq.  label the sectors on which H ± acts nontrivially as the H ± sectors, respectively. We apply the same procedure to the number operator and the annihilation(creation) operator: The commutator [n − ,n + ] = 0, while [â − ,â + ] = Λ(Λ + 1)|Λ − 1 Λ + 1| = 0. This is because the operatorsâ ( †) ± connect the two sectors H ± . From Eq. (C1), we straightforwardly obtain the expressions for H ± : In our numerical scheme we fix a finite cutoff Λ. Therefore we are computing the ground state |ψ 0 (Λ) of H − . Sinceâ ± are noncommuting operators, the two Hamiltonians H − and H + do not commute as well. Therefore, it is not ensured that |ψ 0 (Λ) is an eigenstate of the full Hamiltonian H. We compute the variance ∆H over |ψ 0 (Λ) of the Hamiltonian H = H − + H + , to check whether |ψ 0 (Λ) is an eigenstate of H. The terms in H ± that preserve the sectors H ± give a zero contribution in Eq. (C3). Indeed, the ones that keep the system in the H − sector give a zero contribution since |ψ 0 (Λ) is an eigenstate within this sector by definition. Instead, the ones that keep the system in the H + sector trivially give zero since we do not have any occupation larger than Λ. The only contribution comes from the operatorsâ ( †) ± or, more precisely, the term √ Λ + 1|Λ + 1 Λ| + h.c. which connects the two sectors. Using Eq. (C2), we straightforwardly obtain ∆H = Λ e −2s 4 where P j,k = |k j j k| is the projector on the Fock state with occupation k on site j. The first term of the sum (j = 0) encodes the information about the fixed symmetry sector, since n 2 0 = n 2 0 . The variance given in Eq. (C4) depends on the mean occupation number and on the projector over the Fock space on Λ. In the main text, we show that the system displays a localized phase in the bosonic limit, Λ → ∞, if U > 0. This enables us to estimate Eq. (C4) in the localized phase. In the localized phase, the average occupation number of the ground state is n j ∼ e −j/ξ (cf. Eq. (6)). The exponential decay of the occupation number along the chain reflects on the behavior of the expectation value of P k,j , which decays exponentially fast in k (cf. Eq. (7)). Therefore, the series in Eq. (C4) is finite for Λ → ∞ and L → ∞, since each term is exponentially suppressed. In Fig. 17, we numerically compute the variance ∆H over |ψ 0 (n 0 , Λ) for different values of Λ and n 0 . Rigorously, the cutoff Λ limits the accessible n 0 , since n i ≤ Λ. Nevertheless, because n 0 appears as a constant in the Hamiltonian, we can also compute the ground state |ψ 0 (n 0 , Λ) for n 0 > Λ. The numerical results match Eq. (C4) perfectly. The variance goes exponentially fast to zero. Therefore, an eigenstate  c † ,j c s,j+1 + h.c , where ω ,j = (ω j −E C /2) j+E C j 2 /2 and we introduce p ,j ≡ | , j , j| for convenience. We compute the generator where∆ ,j = ω +1,j −ω ,j = (ω j +E C ). Following Ref. [110], the ansatz for the generator S of the SW transformation is S =