Ergodicity Breaking Under Confinement in Cold-Atom Quantum Simulators

The quantum simulation of gauge theories on synthetic quantum matter devices has gained a lot of traction in the last decade, making possible the observation of a range of exotic quantum many-body phenomena. In this work, we consider the spin-1 / 2 quantum link formulation of 1 + 1 D quantum electrodynamics with a topological θ -angle, which can be used to tune a confinement-deconfinement transition. Exactly mapping this system onto a PXP model with mass and stag-gered magnetization terms, we show an intriguing interplay between confinement and the ergodicity-breaking paradigms of quantum many-body scarring and Hilbert-space

The quantum simulation of gauge theories on synthetic quantum matter devices has gained a lot of traction in the last decade, making possible the observation of a range of exotic quantum manybody phenomena.In this work, we consider the spin-1/2 quantum link formulation of 1 + 1D quantum electrodynamics with a topological θ-angle, which can be used to tune a confinement-deconfinement transition.Exactly mapping this system onto a PXP model with mass and staggered magnetization terms, we show an intriguing interplay between confinement and the ergodicity-breaking paradigms of quantum many-body scarring and Hilbertspace fragmentation.We map out the rich dynamical phase diagram of this model, finding an ergodic phase at small values of the mass µ and confining potential χ, an emergent integrable phase for large µ, and a fragmented phase for large values of both parameters.We also show that the latter hosts resonances that lead to a vast array of effective models.We propose experimental probes of our findings, which can be directly accessed in current cold-atom setups.
Jad C. Halimeh: jad.halimeh@physik.lmu.de6 Conclusions and outlook 16 1 Introduction Gauge theories are quantum many-body models possessing gauge symmetries that dictate an intrinsic local relationship between matter and gauge fields [1][2][3].The quantum simulation of gauge theories has progressed tremendously in recent years across various platforms of synthetic quantum matter [4][5][6][7][8][9][10][11][12][13][14][15][16][17].This offers the opportunity to probe high-energy physics on low-energy tabletop platforms [18][19][20][21][22][23][24][25], establishing a valuable complementary venue to dedicated highenergy experiments, such as the Large Hadron Collider, and to classical computational methods, such as Quantum Monte Carlo and tensor networks.As an example, a cold-atom quantum simulator has recently been successfully employed in probing gauge invariance and thermalization dynamics in a U(1) gauge theory [13,14], with concrete proposals [26,27] of extending it to observe the confinement-deconfinement transition [28,29].In particular, Ref. [26] introduced a scheme that allows the realization of a tunable topological θ-angle term, which emerges in quantum electrodynamics [30][31][32] on account of the topological structure of the vacuum, and has profound consequences on quantum phases, dynamical behavior, and inherent symmetries.Tuning this angle can lead to confinement, an extremely active topic of research in gauge theories [17,33,34] and other many-body spin systems [35][36][37][38], and thus a tunable topological θ-angle term in a quantum simulator can allow the calculation of the time evolution of confined dynamics from first principles.The understanding of non-equilibrium properties of gauge theories sheds light on a variety of dynamical phenomena across high-energy particle physics, condensed matter and even the evolution of the early universe.It is thus important to study the ergodicity-breaking mechanisms that can appear in these theories.Gauge-theory quantum simulators offer the prospect of probing ergodicity-breaking phenomena that are relevant to condensed matter physics, such as disorderfree localization [40][41][42][43][44][45][46][47][48][49][50][51][52][53], quantum many-body scars (QMBS) [6,29,54,55], and Hilbert-space fragmentation (HSF) [56,57], the understanding of which can provide a deeper generic picture of far-from-equilibrium dynamics in closed quantum many-body systems [58][59][60].QMBS, a paradigm of weak ergodicity breaking, comprise angle, where charges are confined at θ=π in the limit of large negative mass, but otherwise deconfined.Also at θ = π there is Coleman's phase transition, which occurs at effective mass µ=0.3275κ, corresponding to the spontaneous breaking of a global Z 2 symmetry [39] that arises from the invariance of Hamiltonian (1) to transformations of parity and charge conjugation.(b) The staggered Bose-Hubbard model realizes the U(1) quantum link model with direct tunneling J, on-site interaction U , staggered potential δ, and a small linear tilting potential ∆.By further adding a period-4 potential χ, the system is tuned away from θ=π, leading to confinement.The U(1) quantum link model or equivalently the staggered Bose-Hubbard model can be further mapped onto the PXP model where spins reside on the gauge links, and doublons on the gauge links represent spin excitations (the mapping will be discussed in detail in Sec. 2).(c) Schematic representation of the rich phase diagram in terms of confinement potential χ and effective mass µ (see Appendix A for details).
states of anomalously low bipartite entanglement entropy that are usually equally spaced in energy across the entire spectrum of the Hamiltonian.
Preparing a system in their subspace and quenching it will lead to persistent oscillations in local observables, an anomalously slow growth of bipartite entanglement entropy, and a corresponding significant delay in thermalization beyond any relevant timescales [6,29,54,55].Indeed, Bernien et al. [6] was the first experiment that observed QMBS in Rydberg arrays, and Surace et al. [29] were able to show that the implemented quantum Ising-like model maps onto the spin-1/2 U(1) quantum link model, which is a regularization of the lattice Schwinger model where the infinite-dimensional gauge fields of the latter are represented by spin-1/2 operators [61][62][63].More recently, it has been shown that QMBS are ubiquitous in gauge theories, occurring over an extensive set of initial states [64,65], in 2+1D [66], and for spin-S representations of the gauge field when the system is prepared in highly excited vacua [67,68].This raises the question about the gauge-theoretic origin of QMBS and motivates their further exploration in realistic gaugetheory quantum simulators.
On the other hand, HSF is a mechanism of ergodicity breaking emerging in models where the Hilbert space is fragmented into exponentially many invariant subspaces due to a exponentially large commutant algebra [56,57,69].In the case of dipole-conserving models [56,57], HSF can be completely characterized by nonlocal integrals of motion [70].More recently, HSF has also been shown to emerge in models with strict confinement [71,72], including in gauge theories [73].Experimental investigations of HSF have been carried out [74,75], but exploring its presence in experimentally relevant parameter regimes of gauge theories would be very appealing, given recent large-scale cold-atom realizations of spin-1/2 U(1) quantum link models [13,14].In this vein, one wonders whether QMBS can occur in regimes of confinement in gauge theories where HSF also arises, or whether scars are a purely deconfined phenomenon.Indeed, while confinement has previously been linked to the appearance of non-thermal eigenstates in other theories with confinement [76,77] and to oscillations of local observables [78,79], a formal connection with QMBS has not been established.In the case of lattice gauge theories, Ref. [29] demonstrated a slowing down of the dynamics under confinement, while Ref. [73] showed the emergence of a symmetry in the limit of strong confinement at zero mass.However, a full analysis of the fate of scarring with both mass and confining potential at intermediate values of experimental interest is lacking.
In this work, we explore the regime of finite mass (denoted by µ) and confining potential (denoted by χ) using exact diagonalization and matrix product state (MPS) techniques [80,81].
We study the PXP model, which exactly corresponds to the spin-1/2 U(1) quantum link model (QLM) [29].Both have been realized experimentally on a Bose-Hubbard quantum simulator [14,64].The mapping between the different models is sketched in Fig. 1(a)-(b).We find prominent scarring regimes concomitantly with confinement, which we explain analytically and demonstrate numerically.We further map out the rich dynamical phase diagram where we find an ergodic phase for small values of the mass and confining potential.For larger values of these parameters, emergent conserved quantities at intermediate times appear, displaying a wide range of prethermal regimes with an emergent integrable phase for large mass in the deconfined regime, in addition to a fragmented phase for large values of both the mass and the confining potential, as sketched in Fig. 1(c) (see Appendix A for more details).We note that the emergent conserved quantities are only truly conserved for infinite values of µ and χ.However, for smaller values of these parameters, a long time can be needed to propagate between Hilbert space fragments, hence we used the prefix "prethermal" to reflect this.We note that in smaller systems this timescale can already be too long to probe in experiment, leading to strong signs of fragmentation at all achievable times.We will focus on the experimentally relevant timescales, and in that context we will study two regions of Fig. 1(c).We will first assess the effect of confinement on quantum many-body scarring when the mass is set to zero, at the bottom left of Fig. 1(c).We will then turn to the top right region where the mass and confining potential are of similar order and Hilbert space fragmentation initially dominates the dynamics.We will show that this region is not homogeneous but hosts special dynamics at different resonances between µ and κ.We will also show that the integrable and chaotic prethermal regimes at the top left can be understood as such resonances.
The remainder of this paper is organized as follows: In Sec. 2, we introduce the U(1) quantum link model with a topological θ-term and discuss its exact mapping onto the PXP model with a mass and staggered magnetization terms.In Secs. 3 and 4, we provide our analysis of scarring and fragmentation dynamics, respectively, in the presence of confinement.In Sec. 5, we discuss experimental probes in a state-of-the-art tilted Bose-Hubbard optical lattice that can facilitate the observation of our findings.We conclude and provide outlook in Sec. 6.Several Appendices and the Supplemental Material (SM) [82] contain supporting analytic and numerical results.

Models and mapping
For completeness, we start by introducing the models considered in this work and briefly outline the mappings between them.This section does not contain new results and the reader is referred to Refs.[26,29] for derivations of the mappings.

U(1) quantum link model
The Schwinger model on a lattice, which is 1+1D quantum electrodynamics, can be conveniently represented through the quantum link formulation, where the gauge and electric fields are represented by spin-S matrices [61][62][63], thereby facilitating experimental feasibility.Upon restricting to S = 1/2, and employing a Jordan-Wigner transformation to map the fermions onto Pauli operators and a particle-hole transformation on both the matter and gauge (electric) fields [26,83], the resulting U(1) quantum link model (QLM) Hamiltonian takes the form where now the Pauli operator σz ℓ represents the matter field on site ℓ, and the spin-1/2 operator ŝ+(z) ℓ,ℓ+1 represents the gauge (electric) field at the link between sites ℓ and ℓ + 1.We consider a lattice with L m sites and L m links under periodic boundary conditions (PBC).The deviation from the deconfined point θ=π is quantified by the confining potential χ=g 2 (θ−π)/(2π) [26].
At χ = 0, the U(1) QLM undergoes a secondorder quantum (Coleman) phase transition at µ/κ = 0.3275 related to the spontaneous breaking of a global Z 2 symmetry arising due to the invariance of Eq. (1) under a parity and charge (CP) transformation at χ = 0 [39,84,85].When χ ̸ = 0, the global Z 2 symmetry is explicitly broken by the last term in Eq. (1), and the Coleman phase transition is no longer present.The generator of the U(1) gauge symmetry of Hamiltonian which can be construed as a discretized version of Gauss's law.Throughout this work, we will work in the physical sector of Gauss's law, defined by the set of gauge-invariant states {|ψ⟩} satisfying Ĝℓ |ψ⟩ =0, ∀ℓ.
Despite being a regularization of the Schwinger model, the U(1) QLM (1) still captures a wealth of the physics of the Schwinger model, including Coleman's phase transitions at χ = 0 [39], and the confinement-deconfinement transition [86].

PXP Hamiltonian
The U(1) QLM with a topological θ-term, described by Hamiltonian (1), can be mapped onto the PXP model with a staggered-magnetization term [29,67,68], described by the Hamiltonian where we have assumed PBC, integrated out the degrees of freedom on the matter sites, and relabeled the gauge sites as ŝα l,l+1 → ŝα l for notational brevity. 1 The projector P annihilates any configuration with neighboring up-spins, in order to only allow configurations that respect Gauss's law.The Hamiltonian also admits a local formulation using the on-site projector Pl = 1  2 − ŝz that annihilates any component with s z l = +1/2.This Hamiltonian corresponds to the PXP model with detuning 2µ and staggered detuning χ.For χ = 0, the PXP model has been extensively studied in the context of QMBS both theoretically and experimentally [6,29,55,64].It has also been studied in the context of HSF in the limit of χ/κ ≫ 1 [73].
In the following sections, we will explore both the ergodic and fragmented regimes of the PXP model with a staggered magnetization term.The approximate values of µ and χ to which these correspond for L m ≈ 30 are schematically illustrated in Fig. 1 (see Appendix A for actual data).
3 Quantum many-body scarring in the presence of confinement For simplicity, we shall set κ = 1 as the overall energy scale and focus on the PXP model (3) in this section.The PXP model is known to host QMBS associated with the Néel state [6,55].Scarring was also observed experimentally from the polarized state |•• . . .••⟩ for µ ≈ ±0.4 [64].Here, for historical purposes [6], a filled (empty) circle denotes a Rydberg atom in the excited (ground) state, but one can equivalently think of it as a spin-1/2 particle in the up (down) eigenstate.The Néel and anti-Néel states correspond to the two vacua of the PXP model, i.e., its degenerate ground states at µ → ∞, while the polarized state corresponds to the charge-proliferated state, which is the non-degenerate ground state of the PXP model at µ → −∞.These two realizations of QMBS occur in the ergodic regime at the bottom left of Fig. 1(c).This means that other initial states display the expected thermalizing behavior, while conventional probes of chaos, such as level spacing statistics, are consistent with those of a chaotic system.In this section, we focus on this region and discuss the impact of confinement on QMBS.One would generally expect confinement to strongly modify the dynamics and lead to the disappearance of revivals from scarred initial states.While this indeed happens for the polarized state at intermediate values of χ, we will show that this is not true for the Néel state, where the confining term has an interesting nonmonotonic effect and leads to a pronounced enhancement of scarring at interme-diate values.In this section, we first present the numerical results illustrating the effect of confinement on the QMBS dynamics and entanglement growth.We then provide analytical explanation of the observed enhancement of QMBS revivals.Finally, we discuss the effect of confinement on the dynamics of defects in QMBS initial states.

Effect of confinement on scarred revivals and entanglement growth
In order to assess the effect of the confining potential on QMBS dynamics, we compute the selffidelity F(t)=| ⟨ψ 0 | e −i ĤPXP t |ψ 0 ⟩ | 2 , where |ψ 0 ⟩ is the initial state.In a generic many-body system, the fidelity starts out at F(t = 0) = 1 and then decays to a value exponentially small in system size, F ∼ exp(−L m ), after a time ∼ ℏ/κ.In a QMBS system, following the initial decay, the fidelity periodically rises again to a value ∼ O(1), corresponding to the revivals of the wave function.We will denote the first revival peak of the self-fidelity by F 1 , and the self-fidelity at half that time by F 1/2 .The signature of QMBS therefore is F 1 ≈ 1 and F 1/2 ≈ 0. Hence, we will use the quantity F 1 − F 1/2 as a probe of QMBS, where it expected to take values close to 1.This quantity distinguishes QMBS from the trivial situation where the initial state is close to an eigenstate, in which case F(t) ≈ 1 at all times and the fidelity difference is strongly suppressed, In Fig. 2(a), we plot both F 1 and F 1 −F 1/2 as a function of staggered magnetization for the PXP model with µ = 0 and the Néel initial state.We observe interesting nonmonotonic behavior as a function of χ, whose mechanism will be explained in detail in the subsequent section.In Fig. 2(b) we see that this nonomonotic behavior of the fidelity is also reflected in the overlap between the Néel state and the energy eigenstates, plotted for several values of χ.The enhancement of scarring at χ ≈ 0.3 can be related to a larger separation between the top band of eigenstates and the bulk of the spectrum, as well as due to the more even energy separation between eigenstates belonging to the top band, as shown in the middle panel of Fig. 2(b).At large values of χ, fragmentation is clearly visible by eigenstates forming bands, separated in energy by multiples of χ.This regime will be addressed in detail in Sec. 4.
Focusing on the Néel state, we find that the growth of entanglement entropy is also highly However, as χ becomes large, the Néel state effectively becomes an eigenstate, leading to an increase of F 1/2 and therefore a decrease of F 1 − F 1/2 .The revival period T is in good agreement with the prediction based on spin precession, 2π/ χ 2 + 1/2, and first-order perturbation theory denoted by E ′ s0 , see Sec. 3.2 for details.The latter predicts the optimal value χ = 1/ √ 8, shown by the vertical black dashed line, which practically coincides with the optimal revival point from the numerics.(b) The nonmonotonic behavior of F as a function of χ can be related to changes in the structure of the overlap between the Néel state and the energy eigenstates |E⟩, here plotted for χ = 0, χ = 0.3 and χ = 1.52.For the anti-Néel initial state or for χ < 0, the dynamics is the same but the spectrum is flipped, E ↔ −E.
suppressed at χ = 0.35, as shown in Fig. 3(a).This nonmonotonicity in the post-quench behavior as a function of χ can be witnessed in local observables, measurable in experiment.For a local observable, we choose the excitation density ⟨n e ⟩ on the even sublattice.In order to avoid finitesize effects, we consider only the times where the dynamics is converged in system size L m = 32.As the sublattice density oscillates in time, we first identify the local maxima and then fit its envelope.The details of this procedure are explained in the SM [82], and the results are shown in Fig. 3  crease around χ = 0.3 compared to case without confining potential (see legend).Finally, we see that for χ = 0.6 the peaks themselves show an oscillatory behavior and are no longer well approximated by a decaying exponential.We emphasize that this is not a finite-size effect but rather a sign that we enter the low-energy regime where only a handful of eigenstates participate in the dynamics.The oscillation in the maxima is then a beating frequency linked to the mismatch in the energy spacing of these eigenstates.

Mechanism of revival enhancement due to confinement
The effect of the staggered magnetization on its own in the PXP model has been studied in Ref. [73] in the limit of χ ≫ κ.In this regime, the Hilbert space fractures and an emergent symmetry appears due to an approximate su(2) algebra.The same spectrum-generating algebra has been linked to QMBS in the PXP model [87][88][89][90][91].In the algebraic picture, the PXP Hamilto-nian can be thought of as the global Ĵx operator, while the Néel state is the highest-weight state of the corresponding global Ĵz operator.Thus, the scarred dynamics is interpreted as precession of a collective spin in a transverse magnetic field.If the algebra were exact, the Néel state would undergo perfect revivals.The PXP model preserves this algebra but only approximately, hence the revivals are present but with a decaying fidelity amplitude.Here we will generalize this picture and identify the mechanism that explains the observed revival enhancement in Fig. 2(a).At µ = 0, the operators defining the mentioned su(2) algebra can be explicitly constructed [92] by dimerizing the chain into 2-site unit cells and interpreting the states of each cell as basis states of a local spin-1, i.e., |↑↓⟩ ≡ |−⟩, |↓↓⟩ ≡ |0⟩, and |↓↑⟩ ≡ |+⟩, where site j of the spin-1 model maps to sites 2j − 1 and 2j of PXP 2 .The spin-1 representation automatically satisfies the constraint within each cell, as no |↑↑⟩ state is admitted, however it can still violate the constraint between the unit cells (e.g., . . .+ − . ..).Hence, one needs to also add operators that annihilate the forbidden matrix elements.To obtain the algebra generators, one then simply maps the global spin-1 operators, Ŝα , into the spin-1/2 space.The PXP Hamiltonian can then be written as (for µ = 0 and κ = 1) where Ĥ⊥ is the Hamiltonian acting on the part of the spin-1 Hilbert space that is annihilated by the PXP constraint [92].The perturbations that take care of the constraint are Without Ĥ1 and Ĥ2 , the Hamiltonian in Eq. ( 4) would be a free spin-1 model.There would then be perfect revivals from the Néel state (|+ + • • • + +⟩ in the spin-1 language) with a period 2 One can alternatively take a mapping where site j of the spin-1 model maps to sites 2j and 2j + 1.In that case the results are the same except that the minus sign in front of the P Ŝz P term in Eq. 4 is replaced by a plus sign.This simple expression matches well the exact result for PXP in Fig. 2(a).The agreement improves with χ, which is due to the fact that, as Ŝz becomes larger, the influence of Ĥ1 becomes weaker.
As χ becomes large, we saw in Fig. 2(a) that F 1 − F 1/2 starts to decrease as the minimum fidelity between the revivals increases.This can be understood easily in the spin-precession picture with the help of the Bloch sphere, illustrated in Fig. 4. The Néel state is at the North Pole, along the Z axis.For χ = 0, the precession axis is the X axis, lying in the equatorial plane.As χ is increased, the precession axis is tilted out of plane and moved closer to the Z axis.Hence, the antipodal point on the trajectory is no longer the anti-Néel state (at the South Pole), but another state placed on the opposite side of the precession axis.As χ is further increased, this opposite point gets closer and closer to the Néel state.In the limit χ → ∞, the precession axis becomes the Z axis and the Néel state an exact eigenstate of the system.The effect of this "axis tilting" can also be seen directly in the overlap between the Néel state and the corresponding PXP eigenstates at various values of χ, as shown in Fig. 2(b).The top band of scarred states is always visible.However, as χ is increased, the eigenstates with the highest overlap shift from the middle of the spectrum towards its left edge, where the ground state resides.
While the precession picture above gives a relatively good approximation of the revival period, it also predicts the revivals would increasingly get better as χ becomes larger.This is in stark contrast with the observed nonmonotonic behavior, in particular with the strong revival enhancement around χ ≈ 0.35.In order to understand the nonmonotonic behavior, we construct the states |S n ⟩, which are the eigenstates of − Ŝx / √ 2−χ Ŝz projected into the constrained PXP space, with n = −L m /2, −L m /2 + 1, . . ., L m /2.This construction provides a good approximation of the scarred eigenstates in the case χ = 0 [92].By computing the energy variance of these approximate eigenstates, one observes nonmonotonic behavior, with local minima at χ = 0 and χ → ±∞ (see Appendix B).This suggest that the revivals will display better fidelity at these points and get worse as χ moves away from them.While this matches well with the numerical result in Fig. 2(a), it does not explain the behavior near χ = 0.35.
To understand why fidelity is maximal around that point, a finer analysis of the effects of the perturbations Ĥ1 and Ĥ2 , Eq. ( 5), on the |S n ⟩ states is needed.It is particularly important to understand how the energy spacing between eigenstates is modified by these terms, as a more regular spacing can lead to significantly better revivals.We first note that as P Ĥ2 = 0 but P Ĥ1 = Ĥ1 , only Ĥ1 will contribute.Hence we can focus on computing ∆E n , the first order correction of the energy of the |S n ⟩ states due to H 1 .As the |S n ⟩ are (projected) eigenstates of a noninteracting Hamiltonian, their structure can be written down and ∆E n can be computed analytically for any n and χ.In particular, we find that the in the thermodynamic limit the energy spacing between two consecutive scarred eigenstates is given by (see Appendix B) where s ≡ (2k)/L m ∈ [−1, 1] can be thought of as energy density.From this, we can immediately see that for χ = 1/ √ 8 ≈ 0.354 the energy spacing is independent of s.At this special point, all eigenstates are exactly equally spaced at first order.This leads to largely improved revivals, accounting for the observed fidelity peak.
A few remarks are in order.The identification of the optimal χ from Eq. ( 7) in the thermody-namic limit also applies equally well to finite systems.In finite systems, the Néel state is still far from the edges of the spectrum at χ = 1/ √ 8 and the level statistics indicate a chaotic system, see Appendix A. The revivals at this optimal χ converge very quickly in system size, and they are not visible for other initial states, which generically thermalize relatively fast [82].Moreover, the result in Eq. ( 7) also allows to get a much better estimate of the revival frequency.To achieve that, we only need to find the value of s which maximizes the overlap with the Néel state.This can also be done analytically (see Appendix B) and yields in the thermodynamic limit.The corresponding period 2π/E ′ s 0 was shown in Fig. 2(a) and gives a very good approximation of the exact data for finite systems.Finally, in SM [82], we show that the optimal revival point χ ≈ 0.35 can also be witness within the semiclassical limit obtained by projecting the dynamics onto a variational manifold [93].Interestingly, in that case the classical orbit linked to the scarred dynamics is not exactly periodic, in contrast with the vast majority of example of QMBS.

Defects and confinement
While adding a nonzero χ affects the algebra and the revivals from the Néel state, from a lattice gauge theory point of view it also causes confinement [26][27][28][29].Notably, each particle-antiparticle pair experiences an energy cost ∝ χd, where d is the distance between the particle and the antiparticle.For χ = 0, there is no cost and so the particles can spread ballistically in opposite directions, i.e., they are deconfined.As soon as χ > 0, any separation distance is penalized, leading to their confinement.
In the the PXP model, a configuration with neighboring unexcited sites |••⟩ implies the presence of a particle or antiparticle between them (depending on whether the sites are even-odd or odd-even).Similarly, an excited site next to an unexcited site (|••⟩ or |••⟩) means that there is no particle/antiparticle between them.Thus, a single particle-antiparticle pair on top of the vacuum takes the form of a single "defect" on top of the Néel (or anti-Néel) state.Specifically, we will use the state which corresponds to the Néel state with an excitation removed on site M near the middle of the chain.This can also be understood as creating two domain-walls between Néel domains •••• and anti-Néel domains ••••.We quench this state and monitor the excitation occupancy on each site ⟨n j ⟩.We also monitor the presence of domain-walls using the quantity ⟨(1 − nj )(1 − nj+1 )⟩ = 1−⟨n j ⟩−⟨n j+1 ⟩, where we have utilized that ⟨n j nj+1 ⟩ = 0. Finally, we track the ZZ connected correlator The corresponding results at χ = 0 and χ = 0.35 are shown in Fig. 5, which have been obtained by matrix product state (MPS) techniques with OBC.MPS allows us to probe large system sizes where the boundary is not reached even at late times (more details about this can be found in SM [82]).The case χ = 0 was previously studied in Ref. [29] for a smaller system, and we find similar results here.There are oscillations both between the two domain walls (i.e., the anti-Néel domain in the middle of chain) and outside of the domain walls (the Néel domain).In the two regions, the dynamics returns to the Néel state with the same period but with a phase difference of π/2 [29].Meanwhile, for χ = 0.35 we find that the dynamics in the two domains synchronizes, with the phase difference quickly going close to 0. At later times it is difficult to see any clear trace of the original defects.While it is not certain if this is directly linked to the enhancement of revivals without any defect, this kind of "selfcorrection" of defects is also experimentally desirable.Finally, the ZZ correlator shows a strong difference between the two values of χ, with clear confinement at χ = 0.35.
In order to quantify confinement more precisely, we compute the "root mean square spread" of the ZZ correlator as We plot this quantity at all times for χ = 0 and χ = 0.35 in Fig 6 .There, we can see a clear difference between the two, as for χ = 0.35 the spreading reaches a plateau around t ≈ 30, while for χ = 0 the correlations continue to spread until they reach the boundary of the system.We thus conclude that at the optimal point for QMBS χ = 0.35, there are clear signs of confinement even in finite systems.

Hilbert space fragmentation induced by confinement
The mass term in Eq. (3) does not fit into the approximate algebra discussed above, and it generically destroys scarring and revivals from the Néel state.However, the combination of nonzero mass and confining potential allows a more complex picture to emerge due to the interplay between them.
In the regime χ, µ ≫ κ at the top right corner of Fig. 1(c), different configurations can have large energy differences compared the dynamical term strength κ.Transitions between them is then heavily suppressed and resonant processes (between configurations with the same energy) dominate the dynamics.The effective Hamiltonian can be obtained by performing a Schrieffer-Wolff (SW) transformation.We will show that in many cases the leading order of that effective Hamiltonian shows exact Hilbert space fragmentation.While in general the higher-order terms can connect different sectors, they are suppressed by factors of approximately κ/(µ ± χ) and so will not matter until long-times.This means that at shorter times the wavefunction will be trapped in the sector in which it was initially, leading to a prethermal fragmentation regime.Our focus is on regimes that can be probed in quantum simulators.In these systems, only shorter times can be probed before decoherence starts to appear.So the impact of higher-order terms is minimal and we mostly focus on the leading order that will dominate the dynamics.
For general values of µ and χ, we extend the approach of Ref. [73], which studied the case µ = 0. We show the main results in this section, while all details are relegated to the SM [82].We find that the odd-order terms in the SW Hamiltonian are identically zero.The leading term is then the second-order one, which is diagonal and reads As this term is diagonal, we need to go to the next nonzero order for dynamics to occur in the Z basis.This happens at fourth order with an effective Hamiltonian given by These results are similar to the µ = 0 case explored in Ref. [73], in which the diagonal secondorder term was understood as stemming from the approximate su(2) algebra.While the addition of the mass term "breaks" it, the additional terms due to a nonzero µ are non-resonant at low order except in the special case of χ = 0.So in general, the cases µ ̸ = 0 and µ = 0 will only show different dynamical terms at higher orders.

χ = ±2µ Resonance
One clear feature of both Hamiltonians (10) and (11) is that some of their terms diverge if χ = ±2µ due to a resonance condition.Let us focus on the case χ = −2µ.Then we have that 2µ + (−1) j χ = 2µ 1 − (−1) j is equal to 4µ for odd j and to 0 for even j.Thus, only spinflips on odd sites lead to an energy change.In the limit χ, µ ≫ κ, the odd sites are effectively frozen as the energy cost to change their state is too high.Hence, the Hilbert space fractures into 2 Lm/2 sectors corresponding to all the possible combinations of the spins on odd sites.The sites that are frozen in the excited position also freeze their neighboring sites due to the PXP constraint.The remaining spins on even sites that have no up-spin neighbors can be excited freely and independently, and thus the effective model at leading order is that of an free spin-1/2 paramagnet.
We can extend our analysis to other values of µ and χ around the resonance such that χ = −2µ + γ and µ, χ ≫ κ, γ.The leading Schrieffer-Wolff term is at first order and reads where K is the set of all even sites with both neighbors frozen in the down position and E 0 is the (constant) contribution of all frozen even sites (so with at least one neighbor in the up position).This Hamiltonian is then clearly integrable for all values of κ and γ.
We show the effects of this resonance in quenching from the Néel and polarized states in Fig. 7.While this fragmentation affects all initial states, unlike QMBS discussed in Sec. 3, we only show these two initial states for brevity as they are the most relevant one for experimental preparation.For both of them, strong revivals of the wavefunction can be seen around the resonance already at relatively small values of χ and µ.For the Néel state, only χ = −2µ leads to revival and not χ = 2µ.Indeed, the Néel state only has excitations on even sites.So for χ = −2µ these can be freely changed between the excited and unexcited states.Meanwhile, for χ = 2µ they cannot be deexcited, and thus also freeze the odd sites due to the constraints.
For the polarized state the pattern is more complicated.This might seem surprising as this state is in the same Hilbert space sector as the Néel state.However, the terms in the Hamiltonian connecting the sectors (changing the state of odd sites) are much less suppressed for the polarized state than for the Néel state, as for the latter two even sites must first be deexcited before such a move can be done.Hence, even for large values of µ, the odd sites cannot be considered as totally frozen and we have some weak couplings between the Hilbert space sectors.This leads to further restriction for revivals in order to make all sectors aligned in energy, as shown in Fig. 8.The additional resonance condition can be cast as with n an integer (see Appendix C for details).
The case γ = 0 yields the condition for χ to be half-integer or integer.Alternatively, we can get rid of γ and write the relation between µ and χ directly as We recognize here the equation for a hyperbola, which was plotted in Fig. 7(b) for various values of n.Fig. 8 illustrates how all eigenstates are approximately equally spaced as long as these equations are followed.

Other resonances
Beyond χ = ±2µ, there are other resonant ratios between µ and χ.Exciting an even site leads to a change in energy of ∆E 1 = −2µ − χ while on an odd site it costs ∆E 2 = −2µ + χ.In order to have resonances, we want ∆E 1 and ∆E 2 to be commensurate.The two simplest cases are if (a) Self-fidelity after a quench.As γ is varied the period changes but the revivals remain close to perfect (b) Overlap between polarized state and the eigenstates.The color indicates the occupation density of the odd sites, and and so the various sectors.The red dashed lines are all spaced in energy by 1 + γ 2 , while the gray dashed line is placed with energy difference 2χ from the highest overlap eigenstate.In all cases of γ we see that all eigenstates show close to equal spacing.∆E 1 is an integer multiple of ∆E 2 or vice versa.This leads to the resonance condition with n an integer (which can also be negative or equal to 0).These resonances are shown in white dashed lines in Fig. 7(b).The only notable ones that lead to changes in the effective Hamiltonian at order 4 or below, are χ = 0 (n = 1) and µ = 3 2 χ (n = 2).Note that unlike the resonance at χ = ±2µ, in these two cases we do not expect to see revivals in Fig. 7 at large values of χ and µ, as the polarized and Néel states become eigenstates and the fidelity stays close to 1 at all times.Nonetheless, revivals can be seen at intermediate values along these resonances for the polarized state.
The effective Hamiltonian around χ = 0 also leads to interesting properties.By studying this special case we can understand the appearance of the effective chaotic and integrable regions at the top right of Fig. 1(c).If µ ≫ χ, the total number of excitations becomes conserved.The effective Hamiltonian at second order then gains a nearest-neighbor XY type term, and a small χ will add some additional diagonal term at first order.The effective Hamiltonian then becomes There is no fragmentation for this Hamiltonian, as each U(1) sector is fully connected.The case χ = 0 is known to be integrable [94].It is interesting to note that the effective Hamiltonian resembles the lattice fermions M 1 supersymmetric model introduced in Ref. [95].The two terms that compose the Hamiltonian are the same, but their prefactors are different.Nonetheless, the same mapping to the an XXZ-type model can be performed, albeit with a different value of ∆.
Consequently, in each of the U(1) sectors the effective Hamiltonian is integrable [94].However, as soon as χ ̸ = 0, this Hamiltonian no longer maps to an XXZ-type model, and when χ and κ 2 8µ are of comparable strength we find level spacing statistics close to Wigner-Dyson (see Appendix A).We emphasize that this does not mean that the full system is ergodic.It still effectively splits into various sectors, but in each of them the effective Hamiltonian is ergodic and there is no further fragmentation.As they are well separated in energy, the level spacing statistics is dominated by the spacings within each sector, which follow a Wigner-Dyson distribution.
Overall, the interplay of the mass and confining terms allows for a vast array of regimes.It should be noted that while our Schrieffer-Wolff transformations are performed assuming that µ or χ are much larger than κ, we observe the onset of fragmentation already for relatively small values of these parameters such as µ = 0, χ = κ/2.This is likely a finite-size effect as larger values of χ and µ are needed to escape the ergodic region as L m is increased.5 Experimental probes in a Bose-Hubbard simulator The U(1) QLM has been experimentally implemented on a Bose-Hubbard quantum simulator for the deconfined case of χ = 0, where gauge invariance was directly observed [13], and thermalization dynamics was probed [14].More recently, an experimental proposal outlines how these experiments can be feasibly updated to study confinement in this model [26,27].Using this proposal, we will show that our results in Sections 3 and 4 can be experimentally observed.
We begin by reviewing the Hamiltonian of the tilted Bose-Hubbard model employed in Refs.[26,27], which is where J is the tunneling strength, U is the on-site interaction strength, bj and b † j are bosonic ladder operators satisfying the canonical commutation relations bj , b † r =δ j,r , nj = b † j bj is the bosonic number operator at site j, and ∆ is an overall tilt.The staggering potential δ distinguishes between matter sites (odd j) and gauge links (even j).Connecting this to Eq. (1) means that ℓ cor-  responds to odd bosonic sites j, while the link between ℓ and ℓ+1 corresponds to even bosonic sites j, where the bosonic model (17) hosts a total of L=2L m − 1 sites, with L m the number of matter sites.The second staggering potential, is related to the topological θ-term, and in the bosonic lattice distinguishes between odd and even gauge links, but has no effect on the matter sites.
The mapping between the bosonic and QLM representations is such that on an odd bosonic site, which represents a site of the QLM, the presence of a single boson represents matter occupation, while no bosons means matter is absent.On an even bosonic site, which represents a link of the QLM, zero (two) bosons indicate that the local electric field points down (up).As such, we need to enforce these local occupations in an experiment.This is achieved in the regime of U, δ ≫ J, µ, where Eq. (1) derives from Eq. ( 17) up to second order in perturbation theory [13,26].
The dominating terms of Eq. ( 17) are then diagonal in the on-site bosonic number operator, and can be collected as where we have resorted to the QLM indexing, which relates to that of Eq. (17) as ℓ corresponding to an odd bosonic site j, while the link between sites ℓ and ℓ+1 corresponds to j+1, the even (i.e., gauge) site between the odd (i.e., matter) sites j and j+1.We can now define a "proto-Gauss's law" with the generator Relating the configurations allowed by Gauss's law in both the bosonic and QLM representations, we find that µ = δ − U/2, which we can insert in Eq. ( 19) and utilize Eq. (20) to get up to an inconsequential energy constant, with c ℓ =2∆(−1) ℓ ℓ.This formulation with the "Stark" coefficients c ℓ has recently been shown to suppress couplings between different gauge sectors, and hence stabilize gauge invariance, up to all numerically accessible times [96] based on the concept of linear gauge protection [97].Looking at Eq. ( 21), we see that a large on-site potential constraints the local bosonic configurations on sites to {|0⟩ ℓ , |1⟩ ℓ } and on links to {|0⟩ ℓ,ℓ+1 , |2⟩ ℓ,ℓ+1 }, as desired.The latter correspond to the local eigenstates of the Pauli operator σz ℓ and ŝz ℓ,ℓ+1 of Eq. (1).In this regime, Eqs.(20) and (2) become equivalent.Using degenerate perturbation theory, the parameters of Eqs. ( 17) and (1) can be related through where, as mentioned, µ = δ−U/2 is the fermionic mass.
The extended BHM Hamiltonian (17) can be realized by a three-period optical superlattice formed by standing waves of lasers.The main lattice laser with wavelength λ forms a lattice with period λ/2, and two additional lattices with period λ and 2λ, respectively.These lasers can be phase stabilized with respect to each other to generate the desired superlattice potential as described in Ref. [26].The two vacuum states in the gauge theory correspond to states |000200020 . ..⟩ and |020002000 . ..⟩ in the BHM.They can be prepared with site-selective addressing techniques using the spin-dependent superlattice, also described in Ref. [26].Now that we have reviewed the experimental implementation of the QLM with confinement in the tilted Bose-Hubbard model as proposed in Refs.[26,27], we can show how our findings have clear signatures in that model.In keeping with experimental relevance, we set the microscopic parameters of Hamiltonian (17) to U = 1368 Hz, J = 58 Hz, and ∆ = 57 Hz, which are close to the values of these parameters as employed in the experiment of Ref. [14].We first investigate the crossover from deconfined scarring to fragmentation.In the deconfined case (χ = 0), when quenching the vacuum state to µ = 0, the system displays scarred dynamics in form of "state transfer" between the two vacuum states, see Fig. 9(b).We obtain the numerical results for the BHM using the TenPy toolkit [98] (see [82] for details).
The system shows persistent many-body oscillations between the two vacuum states, as can be seen in Fig. 9(d).
When the confinement potential is increased to χ = 0.35κ, the system remains ergodic, while scarring from the Néel state persists.This regime of confined scarring is discussed in Sec. 3. The transfer to the opposite vacuum state is suppressed due to confinement, but dynamics between the unit cells are still present, as shown on Fig. 9(e).At late times, we actually see an occupation of the gauge sites closer to 2 than in the χ = 0 case, indicating that the state returns nearer to the Néel state.
The system becomes fragmented when confinement potential χ is above 0.5κ.Here, we use the resonant case with χ = κ and µ = −0.5κas an example.As shown in Fig. 9(c), unit cells marked in gray shades are tuned on resonance with µ = −χ/2, while the adjacent unit cells are out-of-resonance.Therefore dynamics are restricted within the unit cells, see the numerical results in Fig. 9(f).
We further probe confinement and fragmentation by investigating the spread of a particleantiparticle excitation in the vacuum background [26,29,99].
In the Bose-Hubbard quantum simulator, this state corresponds to a |. . .0101 . ..⟩ impurity at the center of the |. . .0020 . ..⟩ background, see Fig. 10(a).This state can be prepared by a local addressing operation in the optical superlattice, as described in [26].In the deconfined case, the excitation couples with neighboring unit cells, see Fig. 10(b).The light-cone spreading can be seen in the numerical results in Fig. 10(e) and (g), calculated numerically using MPS.Additional to the mean density ⟨n j ⟩ in panel (e), in panel (g) we calculate the density-density correlation between site 43 (a gauge site between the two excitations) and the rest of the gauge sites, ⟨n M nj ⟩ c .This quantity exactly maps to the ZZ correlator used in the PXP model.
In the fragmented case, however, the dynamics are restricted within each unit cell, as we demonstrate for the case χ = κ and µ = −0.5κ in Fig. 10(c).The excitation is disconnected from the rest of the system and goes on its own dynamics, see Fig. 10(d), therefore, no spreading can be observed in the density-density correlation, see Fig. 10(f).
We note that the PXP model with a confining potential can also be implemented in Rydberg atoms [6,29] and in a tilted Bose-Hubbard simulators [64,99].While in this implementations only the physical sector with no background charges can be probed, it presents other advantages compared to the mapping discussed in this section.Notably, even with the confining potential only two lattices (with period λ and 2λ) are needed.As the matter degrees of freedom have been integrated out, with the same number of bosonic sites one can probe twice the number of gauge sites.

Conclusions and outlook
Using a combination of analytic and numerical methods, we have calculated the dynamical phase diagram of the U(1) quantum link model with a topological θ-term, or equivalently, the PXP model with mass and staggered magnetization terms [29].By tuning the topological θ-angle on a quantum simulator, as proposed in Refs.[26,29] and realized in Ref. [99], it is possible to controllably induce confinement in this model.We map out various ergodicity-breaking phases that can be observed experimentally by tuning the confining potential χ and the effective mass µ in a coldatom setup.Starting from the ergodic phase at small values of the mass and confining potential, we study the crossover to the prethermal fragmented phase at large χ.Additionally, we have identified various resonant processes in the (χ, µ) dynamical phase diagram, which we have analyzed in detail.We further uncovered regimes of robust quantum many-body scarring in the presence of confinement.
Our results are readily accessible in modern cold-atom quantum simulators [13,14] by adding a staggering potential in the associated tilted Bose-Hubbard optical superlattice [26,27].Our findings highlight a very interesting phenomenological and technological aspect.They show that confinement and quantum many-body scarring, two paradigmatic phenomena of gauge theories, can arise concomitantly.Importantly, our results also show that this behavior can be captured on current quantum simulators, which is encouraging given the current drive of probing high-energy physics on table-top quantum devices.Our work also further sheds light on the question of the gauge-theoretic origin of quantum many-body scarring.Previously, scarring and confinement were thought to arise in different regimes, but here we show that they can coexist and give rise to interesting dynamics.Confinement is a fundamentally gauge-theoretic phenomenon, so it is interesting to further investigate the nature of confined scarring in more generic gauge theories such as Z 2 lattice gauge theories [100] or quantum link models in higher spatial dimensions [66].

A Level statistics with µ and χ
In order to identify non-ergodic regimes, we compute the level-spacing statistics characterized by the ⟨r⟩ value [101].This is only meaningful once symmetries have been resolved, and we consider this quantity only in the most symmetric sector.The relevant symmetries and the sector will change depending whether µ and χ are zero or not.Let us denote the translation operator by T and the spatial inversion operator by L. When χ = 0, the Hamiltonian commutes with both T and L, and we consider the symmetry sector with momentum k = 0 and spatial inversion eigenvalue p = +1.When χ ̸ = 0, the Hamiltonian no longer commutes with either T nor L (as we consider even system sizes).However, it does commute with T 2 and T L. Just as T and L commute in the zero-momentum sector, T 2 and T L commute when the momentum k 2 (with respect to T 2 ) is 0. So in these cases we focus on the symmetry sector k 2 = 0 and p T = +1.Finally, as shown in SM [82], when µ = 0 the Hamiltonian has a particle-hole symmetry and consequently a large number of zero modes.We always truncate the spectrum to remove these modes and we also exclude the edges of the spectrum with the low density of states.
Once the symmetries have been resolved, for small values of µ and χ we expect the Hamiltonian to consist of a single connected component which is non-integrable.However, when fragmentation sets in, we have several disconnected sectors that do not have level repulsion between them.We have also seen that away from resonances the effective Hamiltonian at second order is non-interacting and so trivially integrable.As such we expect a quick departure from the Wigner-Dyson level statistics as soon as fragmentation starts to occur.Figure 11 shows that a value of χ ∼ 0.5κ is sufficient to break ergodicity for L m = 28 sites.The departure from ergodicity around χ/κ ≈ 0.5 is confirmed by Fig. 12 (a) that shows the scaling of the r ratio for various values of χ and L m .As expected, a larger perturbation is needed to induce signs of fragmentation in the level statistics as L m is increased.Nonetheless, spectrum statistics essentially probe infinite-time properties of the system.At the values where we see fragmentation in Figure 11, at finite times after a quench we still expect to see signatures of fragmentation in larger systems.Hence for all system sizes we can denote this regime as "prethermal fragmented".This allows us to use Figure 11 to sketch the phase diagram of the system in Figure 1(c).
Moving away from µ = 0, we see a nonmonotonic behavior of ⟨r⟩ as χ is increased.This is manifested by the lighter "rays" emanating from the ergodic region in Fig. 11.These seemingly special ergodic lines are also visible in Fig. 12(b), which shows the influence of χ for various fixed µ.In order to show that these regions are not fully ergodic, we will concentrate on the case with µ = 0.8 and χ = 0.39, which corresponds to a local maximum of ⟨r⟩.We compute the integrated density of states which simply counts the fraction of energy levels below energy ϵ.For a fully chaotic system in which symmetries have been resolved, we expect G(ϵ) to be a smooth function.However, plotting G for µ = 0.8 and χ = 0.39 in Fig. 13(a) shows visible plateaus, indicating gaps in the spectrum where no levels are found.If we now consider each set of eigenstates between the plateaus as independent spectra, we can unfold them and plot the distribution of level spacings as in Fig. 13(b).Individually, they show good agreement with the Wigner surmise, and consequently the ⟨r⟩ value is close to 0.53.However, due to this splitting  into bands, the full system is clearly not ergodic.Thus, we understand the non-uniformity of the "prethermal fragmented" region in Fig. 11 as the result of interplay between disconnected components.Indeed, the first few orders of the effective Hamiltonian for generic µ and χ do not connect the entire Hilbert space.Thus, for large enough values of these parameters, we end up with an extensive number of disconnected components that can be diagonalized independently.Each of these will have energies centered around some energy E, which depends on µ and χ.Around reso- nances, several disconnected components can end up close in energy.This can happen without them becoming connected by the low-order terms of the effective Hamiltonian.In that case, their energy levels will overlap, thus impacting the energy spacing and the average ⟨r⟩ value.This will lead to lower ⟨r⟩ value, as in the case of unresolved symmetries.
For µ the picture is different, as even for µ = 1.2 and χ = 0.1 the ⟨r⟩ value is close to 0.53.However, further investigation shows that this does not indicate the Hamiltonian is fully ergodic.Indeed, for a small χ and a large µ the dominant Hamiltonian is in the one in Eq. ( 16), which is non-integrable for any nonzero χ.We still have L m /2 + 1 U(1) sectors (corresponding to the total excitation number), but they are very far apart energy due to µ being large.So their energies do not overlap and as such do not create any accidental degeneracies (or near-degeneracies).The level spacing statistics is dominated by the spacings within each sector, which all obey Wigner-Dyson distribution.Hence starting in a sector will lead to prethermalization within that sector for an exponential long time depending on µ/κ, before then starting to explore the other sectors.This means that, at short and intermediate times, the dynamics will only explore a part of the system and so it is non-ergodic.If instead χ = 0, then the Hamiltonian in each U(1) sector is integrable and we see much faster convergence towards Poisson statistics with µ, as shown in Fig. 12(c).Nonetheless, as we have noted above, for µ ≫ χ, κ there is an emergent U(1) conservation law at short and intermediate time but no fragmentation within each U(1) sector.

B Spin-1 parent Hamiltonian
In this section we use the spin-1 picture of the PXP model [92] to derive the optimal χ that enhances the QMBS revivals from the Néel initial state.We consider the eigenstates | Sk ⟩ of the integrable part of the spin-1 parent model of PXP, − 1 The variance is shown in Fig. 14(a), where we observe very different behavior depending on the value of k.In order to pick out the relevant states, we need to understand their overlap with the Néel state.At χ = 0, the Néel state can be written as a superposition of the | Sk ⟩ states as [92] |Z For arbitrary χ, a similar expansion holds as A simple argument as to why this is possible is that the Néel state is the highest-weight eigenstate of S z in the spin-1 picture.This implies that it also has maximal total spin, and the | Sk ⟩ always form a basis of maximal total spin states.
We can then use this to compute the average of the σ 2 E k weighted by it.This is shown in Fig. 14(b), where we recover the anticipated nonmonotonic behavior in χ.The local minima at χ = 0 and χ → ∞ are visible, and this explain why as we go away from these points the fidelity of the revivals starts to decay.However, the peak at χ = 1/ √ 8 is not visible here.To understand its appearance we need to take into account the spacing of the eigenvalues.
For that we need to compute the perturbation of the energies of the | Sk ⟩ states due to the presence of Ĥ1 in the full Hamiltonian in Eq. ( 4).Thanks to the simple structure of the | Sk ⟩, it is possible to compute analytically the first-order correction as ∆E k = ⟨ Sk | Ĥ1 | Sk ⟩.The derivation is similar to Appendix D of [92] and can be found in SM [82].The final result is In general, this cannot be factorized further.However, in the special case χ = 1/ √ 8, we find it reduces to −k/(5 √ 10).This is a remarkable result: for any state | Sk ⟩, at first order its energy will be given by As a consequence, the energy spacing between two consecutive eigenstates will simply be 23   10 √ 10 , which is independent of both k and N b .Thus, at first order, all scarred states are exactly equidistant for any system size.This means that any superposition of these states will have good revivals, which explains the peak in fidelity around this value of χ.
This special property of the spacing at χ = 1/ √ 8 is linked to a change in the behavior of ∆E k .Namely, for |χ| < 1/ √ 8, ∆E k is larger near the edges of the spectrum, while for |χ| > 1/ √ 8, ∆E k is at its largest in the middle of the spectrum.This becomes much clearer after taking the thermodynamic limit, N b → ∞.We can the define s = k/N b , where s ∈ [−1, 1], to get Adding the nominal reference energy of sN b 1/2 + χ 2 , we get As the period of revivals is determined by the energy difference between the eigenstates, we need to compute dEs ds , which gives The change of sign of the term linked to s 2 at χ = 1/ √ 8 is now fully apparent.In order to predict the revival frequency, we need to target the value of s which maximizes the overlap with the Néel state.Eq. (25) gives us the overlap with | Sk ⟩; this is not the same as |S k ⟩, however the difference between the two is polynomial in N b .As we will show later, this means we can neglect this as a subleading correction in the limit N b → ∞.In the same limit, we can use the Stirling approximation where is the Shannon entropy function.We can then rewrite the d 2 k using s = k/N b to get Let us first focus on the case χ > 0. In that regime, √ 2χ − 1 + 2χ 2 goes to 0 as χ becomes large while This is expected as the Néel state gets closer to the ground state for larger positive χ.A consequence is that for positive χ the maximum value will always be in the range [−1, 0].We can then restrict to this range and drop the signum function and absolute values.In the end, this allows us to rewrite with We recognize here that in the thermodynamic limit the most important term is the exponential e −N b f (s) and we use a saddle point approximation.This means that the maximum of d 2 s will simply be the minimum of f (s), which yields s 0 that was given in Eq. ( 8) in the main text.This justifies why we can neglect the difference in norm between the | Sk ⟩ and |S k ⟩: as it is polynomial in N b and s, it will irrelevant in the thermodynamic limit and lead to the same s 0 .
A quick sanity check shows that Eq. ( 8) correctly yields s 0 = 0 at χ = 0 and s 0 → −1 as χ → ∞.A similar derivation can be done for χ < 0, yielding the same formula but without the minus sign in front.This allows us to estimate the frequency of oscillations for a given χ after a quench from the Néel state by plugging s 0 into Eq.(7).While this derivation is done in the thermodynamic limit, Fig. 2 in the main text shows a good agreement in a finite system with N b = 13.

C Polarized state at χ ≈ ±2µ
At the resonance point χ ≈ −2µ, the Hilbert space fractures when χ, µ ≫ κ as creating an excitation on an odd site is extremely costly in energy.However, when 1 < χ/κ < 10 such moves are still possible from the polarized state.As a consequence, the dynamics is not restricted to the Hilbert space sector with no excitations on the odd sites, but the sector with one excitation of these sites also has some nonzero contribution.For a detuning γ = 2µ + χ, the Hamiltonian in each fragment is given by Eq. (12) in the main  (a) Self-fidelity after a quench.When χ is integer or halfinteger, the first peak is significantly higher.(b) Overlap between polarized state and the eigenstates.The color indicates the number of excitations on odd sites and so differentiates states in different sectors.The red dashedlines are all equally spaced by one unit of energy.When χ is integer or half-integer, the spacing between sectors is a multiple of the energy spacing within each of them and all eigenstates are are approximately equally spaced.
text.It is straightforward to see that the energy spacing between eigenstates in that system is 1 + γ 2 for κ = 1.Meanwhile, the energy spacing between the "center" of fragments with respectively 0 and 1 excitations on odd sites is largely unaffected by γ and is approximately 2χ.As a consequence, in order for all energies to be spaced by a regular amount we require one of these quantities to be a multiple of the other.As χ ≫ γ, it must be that 2χ > 1 + γ 2 and so the condition can be expressed as with n an integer.Alternatively, one can combine both equations to remove γ and get 0.0 0.5 1.0 which parameterizes a hyperbola.There is such a hyperbola for any n, although the agreement gets better as n increases as the condition µ, χ ≫ γ, κ is better satisfied.At the same time, points between hyperbolas also get better revivals between as µ and χ increase.This is best seen on the resonance line γ = 0.In that case, the additional condition for energies to align is simply to require χ to be half-integer or integer multiples of κ, as the spacing within each sector is simply κ.The effect of this can be seen on Fig. 15 for κ = 1.
When the condition is met, the polarized state revives with period 2π.
One can monitor that period as well as the fidelity at the first revival as χ = 2µ is increased.Figure 16 shows this for the Néel and polarized states.As expected, the fidelity gets better as χ is increased because other Hilbert space sectors are more suppressed.
In this Supplementary Material, we provide further analysis and background calculations to support the results in the main text.In Sec.S1, we provide additional details on scarring from the Néel state at low values of χ.In Sec.S2, we derive the first-order energy shift from the free spin-1 parent model.In Sec.S3, we discuss the particle-hole symmetry at µ = 0 and its consequences.In Sec.S4, we derive the effective Hamiltonians in the non-ergodic regimes using a Schrieffer-Wolf transformation.In Sec.S5, we provide additional details of the MPS simulation of the spreading of defects in the PXP and Bose-Hubbard models.In Sec.S6, we show that the nonmonotonicity of the revivals from the Néel state as χ is increased can also be seen in the semi-classical limit using the time-dependent variational principle.

S1
Additional results on scarring from the Néel state with µ = 0 In this section we provide additional results on quenches from the Néel state with µ = 0. Fig. S1 (a) shows that after quenches from the Néel state at the optimal point χ = 0.342, the fidelity density log 10 [F(t)/] /L m is well converged already at small system sizes L m ≈ 14.Fig. S1 (b) and (c) show that the nonmonotonic fidelity behavior witnessed for this state is not visible when quenching from other initial product states, and that at the optimal point for the Néel state (around χ = 0.342) other states display relatively quick thermalization in large systems.These results confirm that the revivals observed around χ = 0.342 for the Néel state are genuine scarring and are not due to fragmentation (which would affect most initial states) or to finite-size effects.
We now show more details on how Fig. 3(b) in the main text is obtained.We perform quenches from the Néel state and monitor the excitation density on the even sublattice.We do this for several system sizes and identify the local maxima for each of them.We can then see which of these are converged in system size and discard the other ones.We finally fit these maxima with a de-caying exponential of the form f (t) = a + be −t/τ .Fig. S2 shows the raw data for various system sizes, the scaling of the height of each peak with system sizes, and finally the plot of the converged peaks for χ = 0. We reproduce the same procedure for χ = 0.15, χ = 0.3 and χ = 0.6.For brevity we do not show the full data here but only the converged peaks and their exponential fit (see Fig. 3(b)in the main text).We note that as χ gets larger the finite-size effects occur at increasingly longer times.Overall we find a clear difference between χ = 0.3 and the other values already after 3-4 revivals.The decay time τ also shows a much larger value around χ = 0.3.

S2 Correction to the energy spacing in the spin-1 parent Hamiltonian
In this section we compute the first-order perturbation to the energy eigenvalue of the | Sn ⟩ states due to Ĥ1 .Our derivation follows closely the one in appendix D of Ref. [92].We start from the spin-1 Hamiltonian with spin 1 and N b = L m /2 sites and consider its eigenstates with maximal total spin.We denote these states by | Sn ⟩, with n = 0 to 2N b = N (where we use N = 2N b = L m to keep the notation of Ref. [92]).As this Hamiltonian is the same as for χ = 0 up to a rescaling and a rotation around Ŝy , the | Sn ⟩ will also be identical up to a rotation.The eigenstate n will have energy Due to the rescaling, Ĥ1 is also modified and becomes − 1 2 (|+, 0⟩ + |0, −⟩) ⟨+, −|.Apart from that, the construction from Ref. [92] is not fundamentally modified, meaning that the state | SN−n ⟩ can still be computed by applying a raising operator on | SN ⟩.We can then also target a pair of spin-1 sites and decompose it in the basis of the local − 1 ) A non-zero value of χ will break that as one of the state goes towards the ground state and the other one towards the ceiling state.Nonetheless, when µ = 0 they still show identical dynamics for any χ due to a symmetry of the spectrum.Let us start with the PXP model with µ = 0 and an even number of sites L m = 2M and rewrite it as For µ = χ = 0, it is well known that the spectrum has a particle-hole symmetry with respect to the operator as Ĥ0 , Ẑ = 0.At the same time, the Hamiltonian is also symmetric under the spatial inversion operator L that takes site j to L m + 1 − j and so Ĥ0 , L = 0. Interestingly, the staggered potential term Ĥχ has the exact inverse relations with these two operators and Ĥχ , Ẑ =0. while Ĥχ , L = 0.Both these operators are hermitian and square to the identity: L2 = Ẑ2 = 1.They also commute with each other.It is then straightforward to prove that Ĥ0 + Ĥχ anti-commutes with L Ẑ. Indeed, we have that This has a few direct consequences, such as the presence of the same number of zero-modes as for χ = 0 and the fact that the spectrum is symmetric around E = 0.It also has less obvious consequences for the Néel and anti-Néel states, as both L and Ẑ have special relations with these two states.They are eigenstates of Ẑ with eigenvalue (−1) M , while L takes one to the other as |Z 2 ⟩ = L |Z 2 ⟩.These properties allow us to directly show that for any value of χ, the return fidelity of both states will be exactly identical as long as µ is 0. Let us first consider the energy expectation value of both states: We indeed find that the degeneracy between them is broken.The exact same procedure can be done for any arbitrary power of the Hamiltonian to get This allows us to easily compute the return fidelity by doing the power expansion of the exponential of the Hamiltonian as (S32) From there it is trivial to see that (S33) Similar equivalences can be proven for the expectation of an operator Ô after quenches from the Néel and anti-Néel states as long as Ô, L Ẑ = 0.So while the spatial symmetries of the Hamiltonian do not guarantee any equivalence between the Néel and anti-Néel states after a quench, the particle-hole symmetry does.
This immediately breaks down when a nonzero µ is added.The simplest way to see this is by checking the relation between the chemical potential term Ĥµ and L and Ẑ.A quick calculation yields Ĥµ , L = Ĥµ , Ẑ = 0.So in the procedure used in this section neither L nor Ẑ can be used to change the sign of Ĥµ and it will not be possible to get a term equal to − ĤP XP .

General case, µ/χ irrational
If χ/µ is irrational, each sector is characterized by the number of excitations on each sublattice, n o and n e .Using Eq. (S38), it is straightforward to derive the next order terms in the SW transformation.The leading order is then the second order one which is diagonal, giving The next non-zero term is at 4th order and reads ) These terms are all similar to the ones studied in Ref. [73] for µ = 0, and match them exactly in that limit.As the off-diagonal terms are the same (up to an overall factor), the connectivity of the Hilbert space is unaltered.Thus, we will not focus more on this case as it has already been studied extensively, but instead, we will show that adding the mass term µ allows to engineer different types of fragmentation due to resonances between the χ and µ terms.

Resonances
There is technically an infinite number of resonances between µ and χ.The simplest case in which they happen is when the energy difference on one sublattice is equal to n times (n being integer) the energy change on the other sublattice.Hence two configurations can have a different number of excitations on their sublattices but still the same h 0 .The resonance condition can be expressed as and setting these values can lead to a resonance at order |n| + 1 at the minimum.Indeed, the allowed processes that match the resonance require adding n excitations on one sublattice and removing one on the other, so they require n + 1 moves.We can also chose n to be a negative integer.In that case we need to do the same operation (add or remove) on both sublattices.In general, as we have a dynamical term at 4th order for any value of µ and χ, anything beyond |n| = 3 will never be the dominant off-diagonal term.
More complex resonances can also occur, in which the ratio between the energy cost of an excitation on each sublattice is not an integer but a rational number.In this case we have the same condition as in Eq.S46 but with n a rational number.Let us denote by a and b the numerator and denominator of the irreducible form of n.
In that case the resonance can only lead to new terms in the effective Hamiltonian at order a + b.As integers can be represented with b = 1, this includes the integer-n case.Let us also note that taking a fraction with a = 1 will lead to the same result as the integer case with n = −b.With that in mind, the simplest case that does not fit in the integer limit is a = 2 and b = 3 or vice-versa.As this can only lead to new effective terms at order 5 or more, we will not consider these and thus we focus on integer n. n = 0: χ = ±2µ For this particular relation between the two parameters, we have that µ and χ add up on one sublattice and cancel on the other one.For brevity, we will discuss the case χ = −2µ, for which h 0 = 2χ(n o − M 2 ).Here we see that the sectors are much larger than in the trivial case.
The approach used to compute the SW transformation terms in Eq.S39 fails as V now has terms within individual sectors.This results in an operator Ŝ(1) with some infinite matrix elements (see Eq. S38).However, it also means that the first order effective Hamiltonian, Eq. (S40), is non-zero.It is the leading order Hamiltonian and is given by Thus, all the odd sites are frozen while even sites can flip freely if their neighbors are frozen in the unexcited position.
As the odd sites are frozen, two configurations with the same n o but different arrangements of them will be disconnected.Therefore, in each sector with a fixed h 0 , the Hilbert space shows further fragmentation.The number of subsectors is given by the number of unique configurations on the odd sites: The size of each sector is then 2 k , where k denotes the number of even sites with unexcited neighbors.This is not uniquely determined by L m and n 0 , as even with these two parameters fixed different configurations on the odd sites can lead to different values k.
Note that all sites are either completely frozen (odd sites and the neighbors of excited odd sites) or allowed to flip freely without any interaction (even sites with both neighbors frozen in the unexcited positions).So the effective model at first order is that of a free paramagnet.Let us denote the by K the set of all non-frozen (even) sites, then we can recast the Hamiltonian as The same transformation can be done is the vicinity of the resonance χ = −2µ + γ as long as γ is not much larger than κ.The expression of h 0 becomes h 0 = (2χ − γ)(n o − M 2 ) and the odd sites are still frozen, but now the first order Hamiltonian gains a diagonal contribution as ) where E 0 is the (constant) contribution of all frozen even sites.This Hamiltonian is still integrable, but the energy of the eigenstates and the revival period from the initial Fock states is altered.
If χ = +2µ, the same analysis can be done but the even sites are frozen and the odd ones can flip freely if their neighbors are unexcited.n = 1: χ = 0 If n = 1, we get µ = ±∞χ, meaning that µ can take any value but that χ must be 0.That regime is quite special, as it does not have confinement and the only thing slowing down the dynamics is the large mass of the fermionic matter.Nonetheless, we can treat it like the non-resonant case, and we find that h 0 = 2µ(n o + n e − M ).The Ŝ (1) operator and its commutators with respect to V are the same, but the projector P0 is changed to reflect the new h 0 .It now leaves Pj−1 σ+ j σ− j+1 Pj+2 unaffected as even and odd sites are on an equal footing.So moving an excitation between neighbors preserves h 0 , which only depends on the total number of excitations, and there is an emergent U(1) symmetry.The effective Hamiltonian at 2nd order is then We note that in this regime there is no Hilbert space fragmentation.Indeed, in each U(1) sector all states are connected.
It is interesting to note that the effective Hamiltonian in Eq.S51 resembles the M 1 supersymmetric model of lattice fermions introduced in Ref. [95], which is an integrable model that can be mapped to an XXZ-type chain.The two terms that appear in the Hamiltonian are the same, but their prefactors are different.However, the same mapping to the XXZ model can be done, albeit with a different value of ∆.Consequently, in each of the U(1) sectors the effective Hamiltonian is integrable, as was already derived in Ref. [94].
We can also study what happens if we add a small χ such that κ, χ ≪ µ.In this case we have The PXP term will lead to the same 2nd order term, while the (−1) l ŝz l term will add a first order term.We end up with    do not lead to additional terms in the effective Hamiltonian up to 4th order.This can be seen in the expressions in Eqs.(S41)-S43, as none of the terms that get projected out by P0 in the incommensurate case survive here.For n = −2, we would need terms with 2 ŝ+ on one sublattice and one ŝ+ on the other sublattice.For n = −3, we would need terms with 3 ŝ+ on one sublattice and one ŝ+ on the other sublattice.For n = 3, we would need terms with 3 ŝ+ on one sublattice and one ŝ− on the other sublattice.It is straightforward to see so none of these terms show up, and as such there is no meaningful difference with the incommensurate case up to 4th order.

Resonances at intermediate values of µ and χ
For values of χ and µ which are equal to only a few time κ, the system already starts to form fragments with a definite value of h 0 .However, they are still weakly connected and multiple of them contribute to the dynamics.At the resonances discussed, the sectors have an energy spacing that is always commensurate and so this can lead to revivals.This occurs for example for the Néel state around χ = 0.6 and µ = −0.9.

S5 MPS methods
In order to reach large systems to assess confinement in the PXP model and Bose-Hubbard model, we use matrix product state (MPS) techniques as implemented in the TeNPy package [98].

PXP model
For PXP, we use the time-dependent variational principle (TDVP).We use two-site TDVP at first to grow the bond-dimension until it saturates to a desired value.At that point, we switch to onesite TDVP as it is faster and exactly conserves the energy.We also use OBC instead of PBC to avoid correlations spanning the whole length of the chain.The initial state chosen is then as it leads to suppressed boundary effects.The results were compared to exact ones for 29 sites, and different bond-dimensions were used to assess convergence.Fig. S4 shows the data for the RMS ZZ correlator, introduced in the main text, for two different bond dimensions d max , as well as their difference.The results show that convergence has been achieved for the larger d max used.The results in Fig. 5 in the main text can also be compared with their equivalent using exact time-evolution for L m = 32 and PBC in Fig. S3.Times as long as 70 units cannot be reached as the particle and antiparticle have scattered with each other already at time ≈ 50.Nonetheless, we qualitatively see the same as in the larger system with OBC.Not only do we see confinement, we also see the same synchronization of oscillations in the case with χ = 0.35.

Bose-Hubbard model
The simulations for the Bose-Hubbard model were done using time-evolving block decimation (TEBD) with two-site unit cells.Fig. S5 shows a comparison of the occupation of the sites between two different maximum bond dimensions d max .Even at late times, the difference are small compared to the occupation values.

S6 Symmetric subspace and TDVP semi-classical limit
In order to study the PXP model for large systems, on can use the symmetric subspace approximation [104].The 'symmetric subspace' is defined as a set of states that are fully symmetric under the exchange of atoms on the even (or odd) sublattice, but not necessarily symmetric under the permutations that exchange excitations between the two sublattices.This is a form of "mean-field" theory for the dynamics of the Néel state in the PXP model, as it retains only the information about the total number of excitations on each sublattice (as opposed to the detailed positions of the excitations).The basis states defining the symmetric subspace approximation are, by construction, independent of χ and µ.
The fidelity after the quench from the Néel state within the symmetric subspace is computed in Fig. S6.The nonmonotonic behavior is well captured within this approximation.For generic values of χ, the fidelity decays as L m increases.The only values of χ for which this does not happen are χ = 0, χ = ∞ and χ = 0.41, hinting that these are all special points.This is expected for χ = 0 and χ = ∞, as in these two cases one of the two terms in the Hamiltonian is completely removed.However it is quite a strong result for χ = 0.41 as it implies that the fidelity improvement is not a mere finite-size effect but a more fundamental change.
For an infinite-system, the dynamics in the symmetric subspace admits a classical description using the time-dependent variational principle (TDVP) Ansatz [93,103,104].The quantum state is defined using four classical variational parameters: two θ angles that relates to the excitation density on each sublattice (θ o and θ e ) and two ϕ angles that describe phases (ϕ o and ϕ e ).This corresponds to two spin-coherent states (one on each sublattice), on which a projector was applied to eliminate all states that violate the PXP constraint.This Ansatz was developed to give a semi-classical limit for the Néel state revivals at χ = µ = 0, but it also holds for χ ̸ = 0.The equations of motion are the same as in Eqs.(9a) and (9b) of Ref. [103] for C = 2, except that now for µ z = +χ for φ1 and µ z = −χ for φ2 .We integrate numerically the equations of motion starting from the Néel state, and compute the excitation density on the even (⟨n e ⟩) and odd (⟨n o ⟩) sublattices from them.Fig. S7 compares the TDVP and exact diagonalization (ED) results in a system with L m = 20.The TDVP Ansatz assumes a system of infinite size, so we do not expect it to match ED result at longer times where finite-size effects kick in.Nonetheless, we see a good agreement at shorter times, for all values of χ.We note that there is a major difference between the cases χ = 0 and χ ̸ = 0.For the former, the two ϕ angles are dynamically frozen when starting from the Néel state and there are only two independent degrees of freedom.As a consequence, chaos cannot develop and every trajectory is either periodic (like the Néel and anti-Néel one) or it goes asymptotically towards a fixed point.For χ ̸ = 0, the ϕ angles are no longer frozen and the system has 3 independent degrees of freedom (taking into account energy conservation) and thus chaos can develop.We find that for χ ̸ = 0 the trajectory is no longer exactly periodic, highlighting the instability of the Néel periodic orbit at χ = 0. Despite this, we still find that the excitation density shows persistent oscillations.The expectation of ⟨n e ⟩ after one period (meaning at its first local maximum) shows the same nonmonotonic behavior as in the quantum system, at least qualitatively, as shown in Fig. S8 (a).
In order to have a better understanding of the

Figure 1 :
Figure 1: Tuning the topological θ-angle in the Bose-Hubbard quantum simulator.(a) The U(1) quantum link model Hamiltonian (1) involves a topological θ-angle, where charges are confined at θ=π in the limit of large negative mass, but otherwise deconfined.Also at θ = π there is Coleman's phase transition, which occurs at effective mass µ=0.3275κ, corresponding to the spontaneous breaking of a global Z 2 symmetry[39] that arises from the invariance of Hamiltonian (1) to transformations of parity and charge conjugation.(b) The staggered Bose-Hubbard model realizes the U(1) quantum link model with direct tunneling J, on-site interaction U , staggered potential δ, and a small linear tilting potential ∆.By further adding a period-4 potential χ, the system is tuned away from θ=π, leading to confinement.The U(1) quantum link model or equivalently the staggered Bose-Hubbard model can be further mapped onto the PXP model where spins reside on the gauge links, and doublons on the gauge links represent spin excitations (the mapping will be discussed in detail in Sec. 2).(c) Schematic representation of the rich phase diagram in terms of confinement potential χ and effective mass µ (see Appendix A for details).

Figure 2 :
Figure 2: (a) Revival period and amplitude from the Néel initial state for the PXP model with staggered detuning and L m = 26.The first fidelity peak F 1 is large for all values of χ, indicating the presence of revivals.However, as χ becomes large, the Néel state effectively becomes an eigenstate, leading to an increase of F 1/2 and therefore a decrease of F 1 − F 1/2 .The revival period T is in good agreement with the prediction based on spin precession, 2π/ χ 2 + 1/2, and first-order perturbation theory denoted by E ′ s0 , see Sec. 3.2 for details.The latter predicts the optimal value χ = 1/ √ 8, shown by the vertical black dashed line, which practically coincides with the optimal revival point from the numerics.(b) The nonmonotonic behavior of F as a function of χ can be related to changes in the structure of the overlap between the Néel state and the energy eigenstates |E⟩, here plotted for χ = 0, χ = 0.3 and χ = 1.52.For the anti-Néel initial state or for χ < 0, the dynamics is the same but the spectrum is flipped, E ↔ −E.
(b) for several values of χ.Overall, we find a clear enhancement of QMBS revivals around χ = 0.3 compared to other values.The decay time τ shows an approximately 3-fold in-

Figure 3 :
Figure 3: (a) Bipartite von Neumann entanglement entropy after a quench from the Néel state for L m = 26 and various values of χ.The entropy growth is strongly suppressed around χ = 0.35.(b) Excitation density on even sites for L m = 32 and various values of χ and their exponential fit.Only times for which the expectation values are converged in system size are used.Once χ ≥ 0.5 oscillations are visible in the peaks and an exponential decay no longer describes their behavior with time.

Figure 4 :
Figure 4: Schematic picture of the effect of the confining potential on the dynamics after a quench form the Néel (or vacuum) state.The X, Y and Z labels refer to the the approximate su(2) algebra, which matches with the spin-1 description.

Figure 5 :
Figure 5: Observable dynamics after a quench from the Néel state with a defect for L m = 61 and OBC.The data is truncated to only focus on sites near the middle of the chain where the defect was initially localized.(a)-(c) χ = 0 (d)-(f) χ = 0.35.The effect of confinement is visible for the ZZ correlator in panels (c) and (f).At late times, the Néel and anti-Néel domains are still clearly distinguishable in panels (a) and (b) but not in panels (d) and (e).

Figure 6 :
Figure6: Root mean square of the spread of the connected correlator between Pauli Z matrices at the middle site and other sites.The exact and MPS data coincide well at shorter times, when the influence of the chain boundaries are small.In both cases, the confinement induced by χ is clear.

Figure 7 :
Figure 7: Self-fidelity from (a) the Néel state and (b) polarized state in the PXP model with L m = 20.The black dash-dotted line shows the resonance line χ = −2µ, while the black hyperbolas show the optimal revival around it.The white dashed lines show the other resonances µ = − n+1 n−1 χ.

Figure 8 :
Figure 8: Quenches from the polarized state in the PXP model with L m = 24 around the resonance line as χ = 5 1 + γ 2 and µ = γ−χ 2 .(a)Self-fidelity after a quench.As γ is varied the period changes but the revivals remain close to perfect (b) Overlap between polarized state and the eigenstates.The color indicates the occupation density of the odd sites, and and so the various sectors.The red dashed lines are all spaced in energy by 1 + γ 2 , while the gray dashed line is placed with energy difference 2χ from the highest overlap eigenstate.In all cases of γ we see that all eigenstates show close to equal spacing.

Figure 9 :
Figure 9: (Color online) Proposed experimental realization of the confinement-induced scarring and fragmentation in a Bose-Hubbard quantum simulator.(a) By tuning confinement potential χ from 0, regimes of the deconfined scarring (χ = 0), confinement-enhanced scarring (χ = 0.35κ), and confinement-induced fragmentation (using the resonant case χ = κ, µ = −0.5κas an example) can be realized.(b) and (d): Schematic and MPS simulation of the deconfined scarring, the "state transfer" between the two vacuums.(c) and (f): In the fragmented regime, at resonance, the Schwinger pair creation mechanism still happens, but dynamics are localized within each building block due to confinement.(e) MPS simulation of the confined scarring at χ = 0.35κ.Simulations are for L = 87, with data near the edges truncated for clarity.

Figure 10 :
Figure 10: The experimental probe of a particle-antiparticle excitation in the vacuum background.(a) Preparation of the particle-antiparticle excitation in the Bose-Hubbard quantum simulator.(b) Schematic of the coupling between the excitation and the background in the deconfined case.Coupling between the unit cells leads to the spreading of the excitation.(e) and (g): MPS simulation of the light-cone spreading.(e) The particle density, (g) the density-density correlation between a center gauge site with the particle and the rest of the gauge sites.(c) In the fragmented case, at resonance χ = κ, µ = −0.5κ, the defect is unable to couple to the neighboring unit cells due to the staggering χ, and is hence unable to spread, resulting in decoupled dynamics between the excitation and the background.(d) and (f): Same simulation as for (e) and (g), for the fragmented case.The excitation is localized due to confinement.Simulations are for L = 87, with data near the edges truncated for clarity.

Figure 11 :
Figure 11: Level spacing statistics in the PXP model with PBC and L m = 28.All values below 0.39 are set to this value and all values above 0.53 are set to 0.53.The black lines represent the approximate boundaries between different regimes.

Figure 12 :
Figure 12: Level spacing statistics in the PXP model with PBC for various values of µ, χ and L m .The black and red dashed lines indicate the expected values for Poisson and Wigner-Dyson distributions, respectively.(a)Results for µ = 0, with each colored curve corresponding to a different system size.The inset shoes the system size scaling at the optimal value of χ that maximizes revivals from the Néel state.(b) Results for L m = 30, with each curve corresponding to a fixed value of µ.The values for χ = 0 are not included as they have different symmetries.As µ is increased, a lower value of χ is needed to deviate from 0.53.For µ ≥ 0.4, one can see peaks in the ⟨r⟩ value.However, further inspection shows that these regimes are not truly ergodic (see the text).(c) Results for L m = 30 and two different values of χ.A small χ is enough to break the emergent integrability that appears at large µ.

Figure 13 :
Figure 13: Level spacing statistics in the PXP model with µ = 0.8 and χ = 0.39 in a system size L m = 28 with PBC.The staircase-like behavior of G(E) in (a) is emblematic of the presence of sectors with many energy levels separated by empty intervals where G(E) does not grow.The red dashed lines in the inset indicate the fraction of eigenstates chosen as limits between sectors.In each sector the spectrum was then unfolded and the probability distribution of level spacings s was plotted in (b).

√ 2 Ŝx
− χ Ŝz .If we denote by N b = L m /2 the number of sites in the spin-1 model, the index k runs between −N b and N b .The states | Sk ⟩ are all eigenvalue of angular momentum with total spin and energy E k = 1/2 + χ 2 .We can define their projected counterparts |S k ⟩, and compute their energy variance in the PXP model,

Figure 14 :
Figure 14: Energy variance of the |S k ⟩ states in the PXP model with L m = 26.(a) Data for individual states.(b) Data weighted by the overlap with the Néel state, which accounts for their relevance in the dynamics.The nonmonotonic behavior in χ with local minima at χ = 0 and χ → ∞ is clearly visible.This explains why these points are local maxima of revival fidelity.

Figure 15 :
Figure 15: Quenches from the polarized state in the PXP model with L m = 24 along the resonance line χ = −2µ.(a)Self-fidelity after a quench.When χ is integer or halfinteger, the first peak is significantly higher.(b) Overlap between polarized state and the eigenstates.The color indicates the number of excitations on odd sites and so differentiates states in different sectors.The red dashedlines are all equally spaced by one unit of energy.When χ is integer or half-integer, the spacing between sectors is a multiple of the energy spacing within each of them and all eigenstates are are approximately equally spaced.

Figure 16 :
Figure 16: Self-fidelity and period of revival in the PXP model with L m = 22 along the resonance line χ = −2µ for the polarized and Néel states.As χ gets large, contributions from sectors with excitations on odd sites become smaller and the further resonance condition for χ to be integer or half-integer matters less.

Figure
Figure S1: (a) Fidelity after quenches from the Néel state χ = 0.342 and various system sizes.(b) fidelity dynamics after quenches from various product states with χ = 0.342 and L m = 26.(c) Fidelity at the first revival from various product states with χ = 0.342 and L m = 24.

√ 2 Ŝx
− χ Ŝz operator with total spin-2.Let us denote these local eigenstates as | T2 , m⟩, where m is the magnetization quantum

2 ⟨
Figure S2: Occupation density on the even sublattice after quenches from the Néel state for χ = 0 and various system sizes.(a) data for all system sizes.The squares indicate the local maxima.(b) Local maxima as a function of system size.Each line correspond to a different peak and each color to a different system size.The red rectangle highlights the peaks that are converged in system size at L m = 32 and that will be used to extract a decay timescale.(c) Exponential fit of the converged peaks for L m = 32 to get the decay exponent of the envelope.

)
This Hamiltonian no longer maps to the XXZ model, and when χ and κ 2 8µ are of comparable strength we find level spacing statistics close to Wigner-Dyson in each U(1) sector.

Figure S3 :
Figure S3: Observables dynamics after a quench from the Néel states with a defect with L m = 32.(a)-(c) χ = 0 (d)-(f) χ = 0.35.The effect of confinement is visible for the ZZ correlator in panels (c) and (f).At late times, the Néel and anti-Néel domains are still clearly distinguishable in panels (a) and (b) but not in panels (c) and (d).

35 Figure
Figure S4: (a) Root mean square of the spread of the connected correlator between Pauli Z matrices at the middle site and other sites in the chains for various values of χ for two different bond dimensions.(b) Absolute value of the difference of the ZZ correlator between d max = 250 and d max = 200.

35 Figure
Figure S5: TEBD simulation of quenches from the Néel state in the Bose-Hubbard model with L = (a) Occupation of the middle (gauge) site for various values of χ for two different bond dimensions.No difference is visible upon visual inspection, showing the convergence with maximum bond dimension.(b) Absolute value of the difference in the occupation of each site between d max = 250 and d max = 200.

20 Figure S6 :
Figure S6: Effect of the confining potential in the symmetric subspace approximation.(a) Revival fidelity from the Néel state for various system sizes.(b) Standard deviation of the energy spacings of the scarred eigenstates.Both the best revival fidelity and the most equal energy spacing show optimal points between 0.3 and 0.5.The only values of χ for which the fidelity does not decay with system size are χ = 0, χ = ∞ and χ = 0.41.

Figure S7 :
Figure S7: Sublattice occupation after a quench from the Néel state as obtained with ED and L m = 20 and with TDVP.The data agrees well up to the first revival for all values of χ.