Synthetic U(1) Gauge Invariance in a Spin-1 Bose Gas

Recent experimental realizations of the lattice Schwinger model [Nature 587, 392 (2020) and Science 367, 1128 (2020)] open a door for quantum simulation of elementary particles and their interactions using ultracold atoms, in which the matter and gauge fields are constrained by a local U(1) gauge invariance known as the Gauss's law. Stimulated by such exciting progress, we propose a new scenario in simulating the lattice Schwinger model in a spin-1 Bose-Einstein condensate. It is shown that our model naturally contains an interaction of the matter fields which respects the U(1) gauge symmetry but has no counterpart in the conventional Schwinger model. In addition to the Z2-ordered phase identified in the previous work, this additional interaction leads to a new Z3-ordered phase. We map out a rich phase diagram and identify that the continuous phase transitions from the disordered to the Z2-ordered and the Z3-ordered phases belong to the Ising and the 3-state Potts universality classes, respectively. Furthermore, the two ordered phases each possess a set of quantum scars which give rise to anomalous quantum dynamics when quenched to a special point in the phase diagram. Our proposal provides a novel platform for extracting emergent physics in cold-atom-based quantum simulators with gauge symmetries.

Introduction -Gauge invariance, which refers to the coordinated dynamics of matter and gauge fields being restricted by local symmetries at each spacetime location [1], has fundamentally shaped our understanding of interacting elementary particles in quantum electrodynamics (QED) [2] and quantum chromodynamics [3][4][5]. While a number of breakthroughs in synthesizing gauge fields in cold atoms have been made over the last decade [6], including the experimental realization of artificial electric [7] and magnetic fields [8], spinorbit coupling [9][10][11] and the density-dependent gauge field [12,13], none of them is essentially endowed with local symmetry. Very recently, two experimental simulations of the lattice Schwinger model in cold atoms have changed the situation [14,15]. In these experiments, the U(1) gauge symmetry is synthesized by locally tying the matter and gauge fields to each other via careful control of the tunneling and interactions of neural atoms. As a result, counterparts of such physical phenomena in particle physics as the spontaneous breaking of charge-parity symmetry [16], string inversion and meson formation [17,18] are expected to be observed. Very recently, it has been shown that a Rydberg chain with nearest-site Rydberg blockade [19] can also be mapped to the U(1) lattice Schwinger model [20].
Motivated by the recent experimental progress, we propose a new platform to simulate the U(1) lattice Schwinger model in a spin-1 Bose-Einstein condensate (BEC). The significance of quantum simulator lies not only in simulating existing models of interest, but also in the emergent new physics arising from the intrinsic properties of the simulators. Here we show that the particle collisions in the spinor BEC naturally lead to a term corresponding to the matter-field interaction that has no counterpart in the conventional Schwinger model [14][15][16][17][18]. Consequently, not only we recover the same phases (a disorder and a Z 2 ordered phase) observed in previous simulators, but we also identify a new ordered phase breaking the Z 3 translational symmetry. We prove that the second-order phase transitions from the disordered to the Z 2 -ordered and the Z 3 -ordered phases fall into the Ising and the 3-state Potts universality classes, respectively. The Potts criticality is intimately related to the anomalous quench dynamics and the quantum scars associated to the Z 3 -state, thus is of help in tracing the origin of Z 3 -related quantum scars. This new Z 3ordered phase exists in the experimentally realizable parameter regime of the commonly used atomic species such as 23 Na and 87 Rb [21,22], we therefore expect that these emergent physics can be experimentally observed in the near future.
Model -We consider a spin-1 BEC deeply confined in a one-dimensional optical lattice along the x-direction, as is schematically shown in Fig. 1(a). Under the tight-binding approximation, we label the lowest-band Wannier wave function of the site j by |j, σ with σ = {1, 0, −1} indicating the bare spin states. We construct the spin-dependent hopping using the technique of laser-assisted hopping [23]. To do so, we first introduce a biased magnetic field and a gradient potential, where the former defines a fixed quantization axis and provides the Zeeman shifts, and the latter provides a spinindependent tilt with strength ∆. Next, we shine the BEC with a traveling light on the z-direction which, together with the lattice beam, forms a Λ-type Raman transition that couples the states |j, σ = 1 and |j + 1, σ = −1 . Now, we write out the total Hamiltonian in the lab frame as (setting = 1) withĤ 0 the single-particle Hamiltonian given bŷ andĤ int the interaction Hamiltonian [24,25] given bŷ b j,σ are the local spin operators with S µ=x,y,z the generalized spin-1 matrices, andb j,σ is the bosonic field operator of spin σ at site j, accordinglŷ n j = σn j,σ is the local number operator withn j,σ = b † j,σb j,σ , andΞ j =n j,1 +n j,−1 . Here, p and q inĤ 0 represent the linear and quadrtic Zeeman shift, respectively; U 0 and U 2 inĤ int indicate the strengths of the spin-independent and the spin-dependent interaction, respectively. The last term inĤ 0 characterizes the Raman coupling with δω = ω 2 − ω 1 the frequency difference between the two Raman beams, and (−1) j λ 0 the staggered hopping, where the phase factor (−1) j can be achieved by making the net recoil momentum acquired by the atoms equal one-half the lattice wave vector [26].
We obtain the time-independent Hamiltonian in the rotating frame, i.e.,Ĥ 0 →Û (t)Ĥ 0Û is a unitary operator properly chosen to eliminate both the phase factor e iδωt in the hopping term and the gradient term ∆ [26]. Furthermore, an additional transformation b j,1 → e i(jπ)b j,1 andb j,0 → e i(jπ/2)b j,0 is applied, which removes the stagger phase factor (−1) j inĤ 0 while keepinĝ H int intact [26]. After these transformations, we havê with p = p − ∆/2 − δω/2. We will focus on the case under Raman resonance, i.e., p = 0. The resonant Raman coupling introduces two dressed states, denoted by |± and gapped by 2λ 0 as illustrated in Fig. 1(a), with the associated annihilation operators given bŷ Now we construct the U(1) lattice gauge model. First, we define the modeb j,0 on the spin-0 component as the gauge field, schematically denoted by the ovals in Fig. 1(b); and the lower-lying dressed modeâ j,− as the matter field, indicated by circles in Fig. 1(b). The spin-dependent interaction U 2 in Eq. (3) thus establishes the matter-gauge interaction. To prevent particles from being scattered into other modes (than the matter or the gauge modes defined above), we adopt the following two restrictions [26]: i) we require λ 0 q U 0 , in which case two neighboring lower-energy dressed states |− are resonantly coupled to the gauge field |j, 0 through the spin-exchange interaction, i.e., ∼â j−1,−b † j,0b † j,0â j,− , while the higher-energy dressed states |+ are far off-resonant; ii) we restrict that there are at most two particles on a gauge mode and at most one particle on a matter mode. Restriction ii) can be satisfied by a proper preparation of the initial state. These two restrictions help to further simplifyĤ and lead to the effective Hamiltonian [26] where we have definedN j =â † j,−â j,− the number operator of the matter field, m = q − λ 0 − U 0 /2 andŨ = (U 0 − U 2 )/2.
One can observe thatĤ eff possesses a global translational symmetry and a local U(1) gauge symmetry, where the latter is generated by the Gauss operatorĜ j : which is defined on a building block consisting of two neighboring gauge fields and one matter field, as illustrated in Fig. 1(b). Furthermore, to acquire the QED interpretation ofĤ eff , we perform the Jordan-Wigner transformation on the matter fieldsâ j,− and rewriteĤ eff as [26] Clearly, in the case ofŨ = 0,Ĥ f reproduces the quantum link expression of the lattice Schwinger model with the gauge fields being realized by spin-1/2 spinors [14,35]. Specifically, the first term ofĤ f characterizes the staggered mass of the charged ferminonsψ j and the last term denotes the matter-gauge interaction withσ ± j the raising/lowering operators of the photons (gauge bosons) [36].Ĥ f provides the following QED interpretation ofĤ eff . The occupation of even and odd matter sites respectively represent the electrons and the positrons (see Fig. 1(b)), and the spin-exchange interaction (U 2 term) describes the process that a pair of electron and positron annihilate with each other and in the mean time the electric field is flipped. In a building block, the local Gauss operatorĜ j ensures the total flux of the electric field being equal to the number of charged particles, representing the manifestation of the Gauss's Law. For finiteŨ , we additionally have a nearest-site matter-matter interaction (the second term inĤ f ) that has no counterpart in the conventional Schwinger model [16]. This term comes from the intrinsic interactions of the spin-1 BEC and will lead to a rich phase diagram as will be shown below.
Phase diagram-We discuss the equilibrium phases at 1/3 filling, i.e., there are totally L particles for a chain with L lattice sites, and focus on the gauge sector with no background charges, i.e., G j = 0 [16]. In this case, four occupation configurations, |0 1 0 , |2 0 0 , |0 0 2 and |1 0 1 , are allowed in a building block, as displayed in Fig. 2(a), where n j,0Nj n j+1,0 denotes the Fock basis. Since the state |1 0 1 is a dark state that are not coupled to the other three states through the U 2 interaction, we restrict our discussion within the subspace spanned by the remaining three states.
In the phase diagram, we identify phase transitions between the disordered phase and the ordered phases, D-Z 2 and D-Z 3 , to be of 2nd order, while the transition between the two ordered phases Z 2 -Z 3 of 1st order. The phase boundaries as well as the transition orders are determined by whether the 1st-or the 2nd-order derivatives of the ground-state energy with respect to the parameters (m or U 2 ) exhibit discontinuity or not [26]. Furthermore, the 2nd-order transitions D-Z 2 and D-Z 3 respectively belong to the Ising and the 3-state Potts universality classes, whose low-energy critical behaviors are described by the conformal field theory with different central charges c [39]. Practically, one can extract c through fitting the curve [40] where S(l A ) = −Tr(ρ A logρ A ) is the von Neumann entropy of the subsystem A with length l A , and s is a non-universal factor. In Fig. 2(c), we show the dependence of c as a function of the chain length L at two critical points (corresponding to the cross and the star in Fig. 2(b)), in which one can observe that the transitions D-Z 2 and D-Z 3 exhibit c = 0.5 and c = 0.8 in the thermodynamic limit 1/L → 0, clearly indicating the Ising and the Potts universality classes, respectively [39]. (c) The von Neumann entropy S plotted versus energy spectra i at the chiral (denoted by diamond in Fig. 2(b)), the middle (triangle) and the Potts critical (star) points on the line ofŨ = −2m in (c1)-(c3), where crosses and plus signs indicate the scar states with large overlap to the |Z2 and |Z3 states, respectively. In the calculation, we take L = 18 and S is calculated using lA = L/2.
interaction for 23 Na and 87 Rb is U 2 /U 0 ≈ 1% and −0.5% [21,22], respectively, and hence the emerged Z 3 phase, D-Z 3 and Z 2 -Z 3 phase transitions are expected to be directly observed in these two most commonly used atomic species. Quench dynamics and quantum scars -Since we fix the gauge sector (G j = 0), the matter field and the gauge field are no longer independent. The U(1) lattice gauge model can therefore be mapped to a spin-1/2 chain by eliminating the matter fields [20], i.e.,â j−1,−b † j,0b † j,0â j,− + h.c. ↔σ x j and N j ↔ −(σ z j +σ z j+1 )/2 using the Gauss's Law Eq. (7) witĥ n j,0 ↔ 1 +σ z j . Following this rule, the mapped spin Hamiltonian takes the form [26] withP j = (1 −σ z j )/2 the projection operator which projects out the cases of two neighboring spins being polarized up simultaneously. This projection is necessary to make sure that the system remains in the G j = 0 sector and the resulting states can be described by the three allowed configurations shown in Fig. 2(a). Particularly atŨ = m = 0 (denoted by the diamond in Fig. 2(b)),Ĥ s reproduces the PXP model [41] which was originally realized in a Rydberg chain [19]. The PXP Hamiltonian carries a symmetryχĤ sχ = −Ĥ s withχ = iσ z j . As a result, the energy spectrum is symmetric about = 0. This symmetry corresponds to the chiral symmetry of the original lattice gauge model [26]. The PXP model is well known to lead to dynamical revivals which refers to the phenomenon that the post-quench evolutions of the Rydberg |Z 2 and |Z 3 charge density waves (CDW) exhibit periodic recoveries and slow thermalization [19]. This revival can be attributed to the quantum many-body scar states [41][42][43], which are the low-entropy eigenstates of the PXP Hamiltonian that violate the eigenstate thermalization hypothesis.
The |Z 2 and |Z 3 ordered states in our current model correspond exactly to the Rydberg |Z 2 and |Z 3 CDW states, and hence our model would also exhibit the dynamical revivals by quenching these two states into the chiral point. We perform such numerics by exactly diagonalizing (ED) the effective HamiltonianĤ eff , and plot the evolution of the Loschmidt echo L(t) = | ψ(0)|ψ(t) | 2 and the occupation on one gauge site in Fig. 3 (a1) and (b1), where |ψ(0) is initialized by |Z 2 (solid line) and |Z 3 (dashed line) states, respectively. The periodic oscillating curves clearly demonstrate the dynamical revivals. In Fig. 3 (c1), we show the eigen-spectrum i with vertical axis S denoting the bipartite entropy of eigenstates | i , where the scar states responsible for the |Z 2 -and |Z 3 -revivals are marked by cross and plus signs, respectively. These scars are selected according to the projective probability | i |Z 2 | 2 and | i |Z 3 | 2 above the threshold ≥ 0.03. As one can see, the scar states possess equal energy intervals and relatively low entropy within the spectrum. The energy interval ∆ matches well with the revival period T of L(t) via the relation T = 2π/∆ . In comparison, the dynamics of the same quantities, quenched to the Potts critical point, are plotted in Fig. 3 (a2) and (b2), where the physical quantities exhibit fast thermalization without oscillation. To check the validity of these results, we also carry out numerical calculations based on the original Hamiltonian Eqs. (3) and (4) using the technique of matrix product states [44]. The results are in excellent agreement with the ED results when the condition λ 0 q U 0 is satisfied [26]. This also serves as a confirmation of the validity of the effective HamiltonianĤ eff .
The origin of the scars is of tremendous interest. Recently, Yao and co-workers observed that the |Z 2 -related quantum scars migrate from the low-energy low-entropy states of the Ising transition [45]. Considering the diagram Fig. 2 (b) possesses a Potts criticality, hence now we have an opportunity in tracing the origin of the scar states associated with the |Z 3 dynamics. We focus on the lineŨ = −2m (see the dotdashed line in diagram Fig. 2(b)) and show the |Z 2 -as well as the |Z 3 -related scar spectra at the chiral-symmetric (diamond), middle (triangle) and Potts critical (star) points in Figs. 3 (c1)-(c3), respectively. One may immediately observe that the spectra (c2) and (c3) are asymmetric about = 0 due to the chiral symmetry breaking induced by the finiteŨ interaction. Furthermore, as the Potts critical point is approached, |Z 3 -and |Z 2 -related scars respectively transfer to the lowand high-energy regimes, indicating that the scars associated with |Z 3 originate from the low-energy low-entropy states of the Potts transition.
Summary and discussion-We proposed a scheme to synthesize the U(1) gauge invariance in a spin-1 Bose gas. The effective model exhibits a matter-field interaction which gives rise to a new Z 3 -ordered phase. This ordered phase connects to the disordered phase by the Potts criticality whose lowenergy eigenstates are found to be the origin of quantum scar states responsible for the anomalous dynamical revivals of the |Z 3 states. However, several questions remain unclear at the current stage. For example, what is the interpretation of the emerged matter-field interaction in particle physics? Why do the ordered states, |Z 2 and |Z 3 , tend to be thermalized at quantum criticality? These questions will be addressed in the near future. Two very recent works [46,47] have explored the possibility of tuning the topological angle in the lattice Schwinger model. It will be also interesting to consider the similar possibility in our model and study the combined effect of topological angle and the matter-field interaction.  In this Supplemental Materials (SM), we provide additional information on this work. The SM is arranged as follows. In Sec. I, we show in detail the derivation of the effective Hamiltonian. In Sec. II, we discuss the relationship between our lattice model and the U(1) lattice Schwinger model, as well as how to transform our lattice model into a spin-1/2 chain. In Sec. III, we discuss more on quantum phase transitions in the equilibrium phase diagram. In Sec. IV, quench dynamics based on the time-evolution of matrix product states (tMPS) is included.

I. EFFECTIVE MODEL
Our model system concerns a spin-1 Bose-Einstein condensate (BEC) confined in a one-dimensional optical lattice along direction x. In Fig. S1 (a), we re-plot the architectures of our model with more details, which serves as a complement to Fig. 1(a) of the main text. The biased magnetic field B fixes the quantization axis long the z-direction. The lattice beam propagates with frequency ω 1 along the x-direction and is linearly polarized in the y-direction. To realize the laser-assisted hopping, we adopt another σ + -polarized traveling wave along the z-direction with frequency ω 2 , which together with the lattice beam forms the Λ-type Raman transitions [1][2][3], as schematically illustrated in Fig. S1 (b). The Raman beams couple the states |±1 while leaving the state |0 intact due to selection rule. Now, we have the single-particle Hamiltonian in the laboratory framê where the first two terms represent the kinetic energy and the lattice potential of the neural atoms, respectively, with k 0 = ω 1 /c the wave vector of the lattice beam. ∆ is a tilt field with a = k 0 /π the lattice constant, which can be realized using either the gravity [4] or another optical trap with intensity gradient [5]. The 4th term represents the Raman-assisted hopping, where M 0 is the two-photon Rabi frequency. The final two terms are the linear and the quadratic Zeeman splittings resulting from the biased field B.
Under the tight-binding approximation, we expand the above Hamiltonian in the lowest-band Wannier basis |j, σ , where j and σ indicate the site and the spin degrees of freedoms, respectively. Then we obtain where the natural (spin-independent) hopping has been ignored due to the deep lattice confinement. λ j = (−1) j λ 0 characterizes the amplitude of the nearest-neighboring Raman-assisted hopping, which couples the states |j, σ = 1 and |j + 1, σ = −1 , with λ 0 = M 0 dxw * (x) sin(k 0 x)w(x − a) the overlap integral and w(x) the Wannier function in coordinate representation. The staggered sign (−1) j can be attributed to the 2a period of the Raman coupling (see Eq. (S1)), and hence M 0 sin(k 0 x) is antisymmetric relative to the minima of the local potential well [2]. The off-resonant Raman coupling between the states |j, 1 and |j − 1, −1 has been neglected, as indicated in Fig. S1(b). Additionally, the Raman configuration would not result in on-site spin flips since the on-site integral naturally vanishes, i.e., dxw * (x) sin(k 0 x)w(x) = 0.
To discuss the effect of atomic interactions, we use the second quantization formalism. The second-quantizedĤ 0 is given in Eq. (2) of the main text. Under the unitary transformation we enter the rotating frame. Clearly,Û (t) contains two parts of contribution: the former related ton j is a spin-independent rotation which can eliminate the gradient term, while the latter with respect toF z j is to gauge out the off-diagonal phase factor e iδωt inĤ 0 and the remnant phase factor e i∆t due to former rotation. Moreover, another transformationb j,1 → e i(jπ)b j,1 and b j,0 → e i(jπ/2)b j,0 helps to eliminate the staggered sign in λ j and rewrite the single-particle Hamiltonian intô which is Eq. (4) of the main text. At p = 0, or namely δω = 2p − ∆, the Raman resonance occurs. Thus far, the total Hamiltonian is given byĤ =Ĥ 0 +Ĥ int , witĥ being the spin-1 interaction Hamiltonian also shown in Eq. (3) of the main text. The resonant Raman coupling opens up an energy gap 2λ 0 for the two dressed modes |j, ± (associated with annihilation operatorsâ j,± , respectively). These two dressed modes are related to the bare spin modes by Eq. (5) of the main text. We remark that the lower-lying dressed modes |j, − are treated as the matter fields in our lattice gauge model, and the spin-0 modes |j, 0 the gauge fields. In the parameter regime λ 0 q U 0 [namely the restriction i) indicated in the main text], only the matter fields, |j − 1, − and |j, − , are resonantly coupled to the gauge field |j, 0 satisfying E j−1,− + E j,− 2E j,0 ; whereas all the other processes involving the upper mode |j, + are far-off resonant, as clearly illustrated in Fig. S1(b). This condition suppresses the unwanted two-body scatterings (scattering out of the matter and gauge modes). Hence in our calculation, we expand the total HamiltonianĤ by the dressed modesâ j,± , exclude all theâ j,+ -related interactions, and then obtain withŨ = (U 0 − U 2 )/2 andN j =â † j,−â j,− , from which one can already observe thatĤ is locally symmetric with respect to the Gauss operatorĜ j =N j + (n j+1,0 +n j,0 )/2 − 1. Here, we emphasize that the U(1) gauge symmetry is protected by the energy gap 2λ 0 (or namely the restriction i)), without which the conservation ofĜ j will be broken. In Sec. IV of this SM, we perform a fully numerical calculation of the quench dynamics using the original Hamiltonian Eqs. (S4) and (S5) with all the modes in presence. There, the effects of the restriction i) onĜ j are demonstrated. To restrict our system within the gauge sector G j = 0, we adopt the restriction ii) as noted in the main text, which helps to further simplify the last two terms of Eq. (S6) in the following way: the former term vanishes if the matter fields are limited to at most single occupation at each site; the latter one can be recast to (U 0 /2) jn j,0 and merge with the single-particle term (q − λ 0 ) jN j as the gauge fields are restricted within the empty or the doubly occupied subspace, i.e., (n j,0 − 1) 2 = 1. Eventually, the effective HamiltonianĤ eff [Eq. (6) of the main text] can be obtained. From experimental point of view, various In the main text, we display a rich phase diagram where various phase transitions are identified. Here, we discuss the identification of phase transitions with further detail. According to the Ehrenfest classification [16], a n th -order phase transition would exhibit a discontinuity of the n th -order derivative of the free energy. At zero temperature, the free energy coincides with the internal energy which can be calculated by diagonalizing the HamiltonianĤ eff . In Fig. S2 (a) and (b), we respectively show the 1st-order and the 2nd-order derivatives of the ground-state energy as functions of m. In each subfigure, the line with squares denotes the case of U 2 = 0.05U 0 , while the line with circles denotes the case of U 2 = 0.3U 0 . The discontinuity of ∂ m 0 in Fig. S2 (a) clearly indicates the 1st-order Z 2 -Z 3 transition, while the discontinuities of ∂ 2 m 0 help to identify the 2nd-order D-Z 2 and D-Z 3 transitions, respectively.
On the 2nd-order critical points, we classify the transitions by calculating the central charge c, which can be realized through fitting the von Neumann entropy S (Eq. (9) of the main text) as a function of subsystem size l A . In Fig. S2(c), we show S at the D-Z 2 critical point (m 0.1U 0 , U 2 = 0.4U 0 ) and the D-Z 3 critical point (m −0.22U 0 , U 2 = 0.1U 0 ), respectively. These two  points are marked by the cross and the star on the diagram Fig. 2(b) in the main text. It is clearly shown that S behaves linearly with respect to s 0 = ln L π sin πl A L , and the line slopes yield the central charge c/3. Furthermore, to avoid the size effect of the system size L, we make the finite-scaling analysis with respect to c, which is shown in Fig. 2(c) of the main text, where the extrapolated c can be obtained in the thermodynamic limit 1/L → 0.

IV. QUENCH DYNAMICS USING TMPS
As mentioned in Sec. I, the restriction i) [λ 0 q U 0 ] helps to eliminate the off-resonantâ j,+ modes, which thus stands as the prerequisite for our effective HamiltonianĤ eff [Eq. (6) of the main text] and the gauge invarianceĜ j [Eq. (7) of the main text]. In this section, let us in detail discuss how this restriction would affect our lattice gauge model in the context of quench dynamics. To this end, we adopt the tMPS method [17] for the numerics based on the original Hamiltonian Eqs. (S4) and (S5). In our calculation, we set the lattice length L = 54 and quench the |Z 2 and |Z 3 states to the chiral pointŨ = m = 0. We focus on two operators that are experimentally measurable; one is the particle number of a single gauge site, e.g. n j,σ=0 , and the other is the averaged expectation of the Gauss operator, i.e., where Ĝ j , according to the definition [Eq. (7) of the main text], is simply the measurement on the local particle numbers within a building block. In the gauge sector G j = 0, the deviation ofḠ to zero reflects the breaking of the gauge invariance. In Figs. S3(a) and (b), we show the dynamics of these two observables with respect to the |Z 2 -and |Z 3 -state quench dynamics.
In each subfigure, the dot-dashed, dashed and dotted lines indicate the cases with U 0 /λ 0 = 1, U 0 /λ 0 = 0.1, and U 0 /λ 0 = 0.01, corresponding to the situations that the restriction i) being unsatisfied, just satisfied, and strictly satisfied, respectively. In comparison, the dynamics of n j,σ=0 obtained by exactly diagonalizing the effective HamiltonianĤ eff are also plotted by solid lines in Figs. S3(a1) and (b1). One can clearly observe that, for the two cases, U 0 /λ 0 = 0.1 and U 0 /λ 0 = 0.01 that the restriction i) is satisfied, the tMPS dynamics of n 1,0 are in good agreement with those obtained by the ED method, andḠ exhibits small oscillation around zero. In contrast, for the case of U 0 /λ 0 = 1 where the restriction i) is broken, the resulting n 1,0 exhibits a large discrepancy to the ED's results, accompanied by the large deviation of the local conservationḠ from zero. These results serve as a criterion for the validity of our effective modelĤ eff .