# Quantum jamming brings quantum mechanics to macroscopic scales

Maurizio Fagotti

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maurizio Fagotti. Quantum jamming brings quantum mechanics to macroscopic scales. Physical Review X, 2024, 14 (2), pp.021015. 10.1103/PhysRevX.14.021015. hal-04179314

HAL Id: hal-04179314

https://hal.science/hal-04179314v1

Submitted on 5 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantum Jamming Brings Quantum Mechanics to Macroscopic Scales

Maurizio Fagotti

Université Paris-Saclay, CNRS, LPTMS 91405, Orsay, France

<figure><img src="image_3.png" /></figure>

(Received 1 September 2023; revised 20 January 2024; accepted 11 March 2024; published 23 April 2024)

A quantum spin-1/2 chain with an axial symmetry is normally described by quasiparticles associated with the spins oriented along the axis of rotation. Kinetic constraints can enrich such a description by setting apart different species of quasiparticles, which can get stuck at high enough density, realizing the quantum analog of jamming. We identify a family of interactions satisfying simple kinetic constraints and consider generic translationally invariant models built up from them. We study dynamics following a local unjamming perturbation in a jammed state. We show that they can be mapped into dynamics of ordinary unconstrained systems, but the nonlocality of the mapping changes the scales at which the phenomena manifest themselves. Scattering of quasiparticles, formation of bound states, and eigenstate localization become all visible at macroscopic scales. Depending on whether a symmetry is present or not, the microscopic details of the jammed state turn out to have either a marginal or a strong effect. In the former case or when the initial state is almost homogeneous, we show that even a product state is turned into a macroscopic quantum state.

DOI: 10.1103/PhysRevX.14.021015

Subject Areas: Condensed Matter Physics, Quantum Physics

## I. INTRODUCTION

Low-dimensional quantum many-body systems easily capture the interest of theoretical physicists as a result of their unusual properties, e.g., strong correlations [1], integrability in the presence of interactions [2] in 1D, and topological order [3] in 2D. Since nowadays such systems have become experimentally accessible [4,5]— turning the concept of quantum simulator [6–8] into a reality—the interest in them has sparked far beyond the academic level. One-dimensional systems, in particular, have been intensively investigated in nonequilibrium set- tings [9–12]. On the one hand, they can exhibit exotic behaviors such as relaxation to nonthermal states [13–15]; on the other hand, there are interacting 1D models that allow for nonperturbative solutions, either exact [16] or asymptotic [17,18], and, especially in quantum spin chains, numerical approaches based on tensor networks [19] are particularly efficient.

Among the most remarkable nonequilibrium phenomena observed in 1D systems, we point up prethermalization [20–28] and prerelaxation [29,30]. Such slow dynamics appear, for example, in strong coupling limits as a result of emergent symmetries [29,31]. One can then develop

effective descriptions that remain valid for possibly large but still intermediate times; we mention, for example, the strong-weak duality that allows one to replace strong- coupling Hamiltonians with weak-coupling ones [32]. In a quantum spin chain, if the strong-coupling model has an axial symmetry, the interplay between the conservation of a component of the total spin and the emergent conservation of the operator that is strongly coupled in the Hamiltonian could result in constrained dynamics— see, e.g., Ref. [33]—which have been recently studied especially in the context of Hilbert space fragmentation [34]. This is an example of a physical situation in which (quasi)particles, at high enough density, can get stuck. When that happens, in analogy with classical physics, we shall say that the system exhibits “quantum jamming.” To the best of our knowledge, the quantum version of jamming was originally identified by mapping classical dissipative systems at finite temperature into quantum systems at zero temperature [35]—see also Refs. [36,37]. The interest in it had come from the idea that jamming might be linked to a new kind of phase transition [38], akin to the more famous glass transition, which could be approached by varying a thermodynamic quantity such as the density.

The key question we address here is which phenom- enology should be expected close to quantum jamming when the system is perturbed out of unstable equilibrium rather than driven from or to thermodynamic equilibrium. In the systems we consider, indeed, the space of jammed states is spanned by product states, which could be reproduced, for example, in optical lattices with ultracold

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
atoms [39]. This opens the door to investigations into the dynamical effects of local perturbations after having prepared the system in a jammed state.

Aiming at both generality and simplicity, we focus on a class of quantum spin-1/2 chains with local Hamiltonians that exhibit quantum jamming due to the fact that spins up are not allowed to jump over other spins up and, if there are none of them along the way, only jumps by multiples of a given number y are allowed. Specifically, the Hamiltonians are constructed (through sums and products) with the building blocks

$$T_{\ell, \ell+y}^{s, -s} = \sigma_{\ell}^{s} \prod_{j=1}^{y-1} \frac{1 - \sigma_{\ell+j}^{z}}{2} \sigma_{\ell+y}^{-s}, \quad \sigma_{\ell}^{z}, \quad (1)$$

where $s = \pm 1$, $\sigma_{\ell}^{\pm} = (\sigma_{\ell}^{x} \pm i\sigma_{\ell}^{y})/2$, and $\sigma_{\ell}^{\alpha}$ act like Pauli matrices on site $\ell$ and like the identity elsewhere. Although the building blocks could seem unnatural, the underlying constraints echo the blockade effect in chains of Rydberg atoms [8,40], and indeed similar terms can arise in effective descriptions of conventional systems, such as in the strong anisotropy limit of the Heisenberg XXZ chain [33]; interactions of that form have also been recently shown to trigger coexistence of phases in fermionic systems [41]. Beyond question is that constrained systems of that kind are attracting more and more attention [42,43].

The model with Hamiltonian $\sum_{\ell} T_{\ell-1, \ell+1}^{+} + T_{\ell-1, \ell+1}^{-}$, known as “dual folded XXZ” [33] or, following the terminology of Refs. [43,44], “hard-rod deformation of XX with length 2,” is a special case of an integrable system introduced by Bariev in 1991 [45]. Recently, Ref. [46] pointed out that, in that model, a local unjamming perturbation in a jammed product state has macroscopic effects. Reference [47] exhibited a similar result in a simplified setting where it was possible to map the dynamics of a single impurity over a jammed state into the dynamics of a single spin up in the XX model ($\sum_{\ell} \sigma_{\ell}^{x} \sigma_{\ell+1}^{y} + \sigma_{\ell}^{y} \sigma_{\ell+1}^{x}$). Such a transformation is an instance of a duality mapping, described in the following, that has the practical effect of removing the kinetic constraints. In the specific case of the dual folded XXZ, it translates trivial noninteracting dynamics in XX into unusual dynamics that incorporate the inhomogeneity of the underlying jammed state into effective deformations of the space; but not only that. The effect of the perturbation reaches macroscopic scales; i.e., the perturbation indelibly affects local observables within a region with extent proportional to the time.

In that respect, Refs. [46,47] have provided a first proof of concept: Macroscopic effects from local perturbations, which we compactly call MELPs, can be observed in systems described by local Hamiltonians with kinetic constraints even if they are prepared in product states. Later, Ref. [48] discovered MELPs in generic systems

(crucially) without axial symmetry from a completely separable excited state, and numerical evidence was found that MELPs could be associated with the emergence of macroscopic entanglement. Which are the key properties of the systems exhibiting MELPs? How general is the connection between MELPs and macroscopic quantumness? These questions have remained unanswered.

We reconsider here nonequilibrium time evolution in systems with quantum jamming from a more general perspective and show that the MELP appears also in the presence of more impurities, and it survives both integrable and nonintegrable interactions, at least as long as the Hamiltonian has a symmetry under “transmutation of particle species” (to be compared with the “inertness” property of Ref. [42]). We show that the state that builds up from a local unjamming perturbation is quantum at a macroscopic scale; this sets quantum jamming apart from its classical counterpart and includes a first proof of macroscopic entanglement associated with MELPs in a system with an axial symmetry. The case without symmetry under transmutation of species requires a separate treatment, and the phenomenology turns out to strongly depend on the inhomogeneity of the jammed state. Specifically, MELPs can still be observed whenever the jammed state is (almost) homogeneous; the more disordered the sequence of species, however, the more hindered the spreading of the perturbation.

## II. A SURVEY OF THE MODEL

The beauty of the constrained model that we propose is that a significant part of its phenomenology can be understood without specifying the details of the interactions. This puts us in the privileged position to discern the effects of quantum jamming from the complex dynamical properties of the underlying systems, which can be described by essentially any spin-1/2 Hamiltonian with an axial symmetry.

### A. From spins to particles, from positions to macropositions

The building blocks in Eq. (1) commute with the total spin in the z direction $S^z = \frac{1}{2} \sum_{\ell} \sigma_{\ell}^{z}$; hence, spins up can be reinterpreted as particles that are conserved in number. In addition, the building blocks preserve the sequence of positions of spins up modulo y; therefore, it is reasonable to distinguish y species of particles, characterized by their position modulo y. We can then separate the information about the species from the position $\ell$ by defining the macroposition $\ell'$ as

$$\ell = y\ell' - b, \quad (2)$$

where $b \in \{0, 1, ..., y-1\}$ labels the particle species. The macroposition plays the role of the position for the
particles; hence, there can be up to y particles in the same macrosite.

The symmetry under transmutation of species that we mentioned before is present when the interactions are blind to the species. This happens, in particular, when the Hamiltonian can be written in terms of the local averages over y sites of the building blocks: $(1/y) \sum_{n=0}^{y-1} T_{\ell-n, \ell-n+y}$ and $(1/y) \sum_{n=0}^{y-1} \sigma_{\ell+n}^z$.

We note that a family of jammed excited states can be readily identified independently of the Hamiltonian: Any state in which all consecutive particles have macrodistance not larger than 1 (i.e., the spins up are at a distance not larger than y) is jammed.

## B. Multispecies particle-hole duality transformation

Since the sequence of species $b$ is conserved, it is convenient to split the representation of the state in two parts, one of which stores just the information about the sequence $b$ and, hence, is stationary. By removing the information about the species, particles become indistinguishable; let us associate each of them with a pseudospin down. A pseudospin up represents instead the absence of a particle, so the number of pseudospins up between two pseudospins down should match the macrodistance (i.e., the difference between ordered macropositions) between the corresponding particles minus one. We can interpret the mapping of spins into pseudospins for given sequence of species as a multispecies particle-hole duality transformation.

The prescription that we describe could look plain at first sight but belies an important complication: Macropositions, pseudopositions (i.e., the position of the pseudospins), and the relative position of a particle in the sequence of species $b$ are not trivially related to one other. That is to say, there is arbitrariness in assigning the positions to a given spin, a given pseudospin, and a given particle, but, once the choice is made in a state, there is no freedom anymore and a link between the different positions is induced in any state.

For the purposes of this work, we can restrict ourselves to a sector with a fixed number $\nu$ of impurities (pseudospins up) over a jammed state. Let us then agree to assign macroposition 0 to the particle with relative position 0 in the state in which the impurities are at pseudopositions 0, 1, ..., $\nu - 1$; we also agree that, in that state, the particle with relative position 0 corresponds to the pseudospin down at pseudoposition $-1$. The macroposition of the particle with relative position $\tilde{j}$ is then given by

$$ \ell'(\tilde{j}) = \tilde{j} + \delta\tilde{j} + \sum_{k=1}^{\nu} \theta(n_k + 1 - k < \tilde{j}), \quad (3) $$

where $n_1, ..., n_\nu$ are the pseudopositions of the impurities and $\delta\tilde{j}$ parametrizes the deformation of the space due to

the different distance that particles of different species have:

$$ \delta\tilde{j} = \sum_{\tilde{n}=0}^{[\tilde{j}]-1} \begin{cases} -\theta_H(b_{\tilde{n}+1} - b_{\tilde{n}} < 0) & \tilde{j} \ge 0 \\ \theta_H(b_{-\tilde{n}} - b_{-\tilde{n}-1} < 0) & \text{otherwise.} \end{cases} \quad (4) $$

An instance of such a deformation has been already pointed out and discussed in Ref. [47] for a single-impurity problem with $y = 2$. Since, to our present purposes, the deformation plays only a marginal role, we refer the reader to Ref. [47] for additional comments.

In the following, we denote the state with the sequence of species $b$ and impurities at $n_1, ..., n_\nu$ by $|n_1, ..., n_\nu; b\rangle$ (the absence of impurities is denoted by $|0; b\rangle$). An example with $y = 2$ follows

$$ \begin{gathered} |\underbrace{\cdots\downarrow}_{\text{jammed}}\downarrow\downarrow\downarrow\downarrow\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow\cdots\rangle \\ \ell'=-2 \quad \underbrace{\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}_{-1} \quad \underbrace{\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}_{0} \quad \underbrace{\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}_{1} \quad \underbrace{\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}_{2} \quad \underbrace{\downarrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow}_{3} \quad \underbrace{\cdots\cdots\cdots}_{\text{jammed}} \end{gathered} \\ \equiv |\downarrow\uparrow\downarrow^{(1)}\uparrow\downarrow^{(0)}\downarrow_{j=0}^{(0)}\uparrow\downarrow^{(1)}\downarrow^{(0)}\downarrow\rangle \\ \equiv |-4, -2, 1; b\rangle, \quad (5) \end{gathered} $$

where $b = ..., 1, 0, 0, 1_{\tilde{j}=0}, 0, ..., \downarrow$ stands for a semi-infinite string of pseudospins down, and positions, macropositions, pseudopositions, and relative positions are denoted by $\ell, \ell', j$, and $\tilde{j}$, respectively; the first representation in the second line shows the pseudospins with the species as superscripts; the last representation displays the state by indicating the pseudopositions of the pseudospins up followed by the sequence of species. Figure 1 shows how some allowed interactions act on Eq. (5).

For given $b$, the dynamics can be projected onto the space of pseudospins by applying the duality transformation to the Hamiltonian. A procedure is described in Sec. V B. We find that symmetric translationally invariant operators with local densities are represented in the pseudospace by operators with densities that are still local and, in addition, are independent of the sequence of species. Nonsymmetric operators, instead, have a sequence-dependent semilocal representation. Here, we are using the term “semilocal” to indicate representations of local observables that do not commute with all local operators acting on the pseudospins at arbitrarily large pseudodistance; an example follows with Eq. (19).

## III. RESULTS

A qualitative understanding of MELPs in the class of models we consider can be gained without specifying the details of the model. We do it focusing on the local magnetization $\sigma_\ell^z$. This is an observable with diagonal form factors in the jammed sector; hence, it is arguably the simplest operator we could study that witnesses MELPs in a transparent way. We then corroborate the qualitative analysis considering specific examples.
<figure><img src="image_4.png" /><figcaption>FIG. 1. Cartoon of some actions allowed by the model with y = 2, where orange ovals represent the (impenetrable) particles. The green (blue) arrows on the bottom (top) show the movement of the particles (impurities) in the exemplar state (5). The spins that are flipped [with respect to Eq. (5)] are highlighted in orange.</figcaption></figure>

A. Local magnetization

Since the total spin in the z direction is conserved, the local magnetization is a charge density. As such, it satisfies a continuity equation, which typically admits of hydro-dynamic descriptions. As pointed out in Ref. [47], this kind of large-scale approach is rather ineffective in the present case. Thus, we work out time evolution employing standard Lehmann representations of the correlations functions. We start here with the one-point function—the local

magnetization—whereas a two-point function is discussed
in Sec. III D. We introduce the notation

$$
:O: = O - \langle 0; \underline{b} | O | 0; \underline{b} \rangle. \quad (6)
$$

In the basis $|n_1, \dots, n_\nu; \underline{b}\rangle$, the local magnetization is a
diagonal operator. Using Eq. (3), we derive in Appendix A
the following matrix elements:

$$
\begin{equation}
\begin{split}
\langle n; \underline{b} | :\sigma_{y\ell'-j}^z : | n; \underline{b} \rangle ={}& \langle 0; \underline{b} | \sigma_{y(\ell'-1)-j}^z | 0; \underline{b} \rangle \theta_H(n + \delta n + 1 < \ell') \\
& - \langle 0; \underline{b} | \sigma_{y\ell'-j}^z | 0; \underline{b} \rangle \theta_H(n + \delta n - \theta_H(b_{n+1} < b_n) < \ell') \\
& + \delta_{\ell',n+\delta n+1} [\theta_H(b_{n+1} < b_n)\delta_{b_{n+1},j-1}] \\
& + \delta_{\ell',n+\delta n} [\theta_H(b_{n+1} < b_n < b_{n-1})\delta_{b_{n-1},j} \\
& \qquad - \theta_H(b_{n+1} < b_n)(1 - \delta_{b_{n,j}})]
\end{split}
\tag{7}
\end{equation}
$$

and

$$
p_{n,t}^{(j)} = \sum_{n'} |\langle \dots, n'_{j-1}, n, n'_{j+1}, \dots | e^{-i\hat{H}_t} | \tilde{\Psi}_0 \rangle|^2. \quad (10)
$$

$$
\begin{align}
& \langle n_1, n_2, \dots, n_\nu; \underline{b} | :\sigma_\ell^z : | n_1, n_2, \dots, n_\nu; \underline{b} \rangle \nonumber \\
& = \sum_{k=1}^\nu \langle n_k + 1 - k; \underline{b} | :\sigma_{\ell'-y(k-1)}^z : | n_k + 1 - k; \underline{b} \rangle. \tag{8}
\end{align}
$$

Roughly speaking, the local magnetization in a product
state with a single impurity—Eq. (7)—is identified with
either one local magnetization or another (y spins apart) in
the underlying jammed state, depending on whether the
impurity is on the left or on the right. The matrix elements
with more impurities—Eq. (8)—have a similar interpreta-
tion and can be reduced to those with a single impurity.

Time evolution is readily obtained:

$$
\langle :\sigma_{\ell}^{z}: \rangle = \sum_{k=1}^{\nu} \sum_{n} p_{n+k-1,t}^{(k)} \langle n; \underline{b} | :\sigma_{\ell'-y(k-1)}^{z}: | n; \underline{b} \rangle, \quad (9)
$$

where $p_{n,t}^{(j)}$ is the probability that the $j$th pseudospin up is at pseudoposition $n$ at time $t$, and it is given by

Here, $|\tilde{\Psi}_0\rangle$ is the projection of the initial state onto the pseudospin space, $\tilde{H}$ is the projected Hamiltonian, and the sum is over all integers $n'_1, \dots, n'_{j-1}, n'_{j+1}, \dots, n'_\nu$ satisfying $n'_1 < \dots < n'_{j-1} < n$ and $n < n'_{j+1} < \dots < n'_\nu$.

For finite $\nu$, the qualitative behavior of the magnetization
can be readily understood if $\tilde{H}$ is local (i.e., the original
Hamiltonian is symmetric) and does not exhibit exceptional
localization properties in the eigenstates. It is then reason-
able to expect that $p_{n,t}^{(j)}$ approaches zero in time.
Consequently, the terms in the last two lines in Eq. (7)
give a subleading contribution to Eq. (9), and the argument
of the step functions on the first line in Eq. (7) can be
shifted without affecting the asymptotic behavior of the
local magnetization. We finally end up with the following
simplified expression:

$$
\langle :\sigma_{\ell}^{z}: \rangle \sim \sum_{k=1}^{\nu} \langle 0; \underline{b} | \sigma_{\ell'-yk}^{z} - \sigma_{\ell'-y(k-1)}^{z} | 0; \underline{b} \rangle F_{t}^{(k)}(x_{\ell}), \quad (11)
$$
where $F_t^{(k)}(x) = \sum_{n \le x} p_{n,t}^{(k)}$ is the cumulative distribution and $x_\ell$ is the effective pseudoposition associated with site $\ell$, i.e.,

$$ \frac{y}{\ell}(x_{\ell} + \delta x_{\ell}) \sim 1 + O(\ell^{-1}). \quad (12) $$

For its generality and transparency, Eq. (11) is the most important formula of this paper.

(i) First, by modifying the underlying jammed state, the coefficients of $F_t^{(k)}(x_\ell)$ in Eq. (11) (which belong to $\{-2, 0, 2\}$) can be changed in an arbitrary way without even affecting the asymptotics of $x_\ell$. This allows us to ignore the coefficients altogether and focus on the cumulative distributions.

(ii) Second, as clarified in the following examples, $p_{n,t}^{(k)}$ typically decays as $1/t$ (at least for small enough $\nu$), and, hence, the typical variation of the cumulative distribution is over ballistic scales $x \sim t$. Thus, we predict a general nontrivial profile for the local magnetization as a function of $x/t$, which is the unusual light cone behavior pointed out in Refs. [46,47] in a specific model with a single impurity.

(iii) Corrections should be expected if there is a finite number of localized eigenstates, whose more explicit effect would be to create effective discontinuities in the relevant scaling limit of $F_t^{(k)}(x)$.

(iv) On the contrary, in order for the impurities not to change the local magnetization in a region that grows in time for a generic jammed product state, the probabilities $p_{n,t}^{(k)}$ should remain localized around the original positions of the impurities, which is possible only if the sector with $\nu$ impurities is spanned by almost only localized excited states. This is practically impossible when the Hamiltonian is symmetric and the interactions are local, but we show in the examples below that it becomes an option in the presence of nonsymmetric interactions.

We emphasize that these qualitative rules apply to the entire class of constrained models constructed with the building blocks in Eq. (1), which include noninteracting, interacting integrable, and nonintegrable systems.

## B. Example: Symmetric Hamiltonian

Reference [47] has proven MELPs in the special case in which the Hamiltonian is local and noninteracting in the pseudospace. In order to prove the robustness of that phenomenology under interacting and integrability breaking terms, we consider the following Hamiltonian with $y = 2$:

$$ H = J \sum_{\ell} (\sigma_{\ell-1}^{x} \sigma_{\ell+1}^{x} + \sigma_{\ell-1}^{y} \sigma_{\ell+1}^{y} + \Delta \sigma_{\ell-1}^{z} \sigma_{\ell+1}^{z}) \frac{1 - \sigma_{\ell}^{z}}{2} \\ + g(\sigma_{\ell-1}^{x} \sigma_{\ell+3}^{x} + \sigma_{\ell-1}^{y} \sigma_{\ell+3}^{y}) \frac{1 - \sigma_{\ell}^{z}}{2} \frac{1 - \sigma_{\ell+1}^{z}}{2} \frac{1 - \sigma_{\ell+2}^{z}}{2}. \quad (13) $$

The first line is a hard-rod deformation of the XXZ model and is interacting integrable. The second line is a symmetric interaction that breaks integrability. [49] The Hamiltonian in Ref. [47] is recovered by setting $\Delta = g = 0$. The projection of $H$ onto the pseudospace reads

$$ \tilde{H} \sim J \sum_j (\tau_j^x \tau_{j+1}^x + \tau_j^y \tau_{j+1}^y + \Delta \tau_j^z \tau_{j+1}^z \\ + g(\tau_{j-1}^x \tau_{j+1}^x + \tau_{j-1}^y \tau_{j+1}^y)) \frac{1 + \tau_j^z}{2}, \quad (14) $$

where we indicate the pseudospin operators by $\frac{1}{2}\vec{\tau}_\ell$ to avoid confusion with the spin operators.

### 1. Single impurity: Noninteracting representation

Quite generally, the physics of a single impurity can be described by an effective noninteracting Hamiltonian. Our example is no exception, and, in fact, we find that the coupling constants $\Delta$ and $g$ that set this example apart from the problem studied in Ref. [47] are irrelevant. Specifically, the additional interactions results in just a shift of the excitation energy, given by $E(p) = 4J(\cos p - \Delta)$, that does not have any impact on the local magnetization.

### 2. Two impurities: Integrable representation

Much more interesting is the case in which there are two impurities in the initial state. In this sector, interactions start playing a role; however, integrability is effectively preserved, in the sense that the Bethe ansatz [2,50] provides sufficient degrees of freedom to solve the problem. We obtain the scattering phase

$$ S(p_1, p_2) = - \frac{1 - 2[\Delta + g \cos(p_1 + p_2)]e^{ip_2} + e^{i(p_1 + p_2)}}{1 - 2[\Delta + g \cos(p_1 + p_2)]e^{ip_1} + e^{i(p_1 + p_2)}}. \quad (15) $$

Importantly, the corresponding integrable model exhibits also a bound state in the two-particle sector, provided that the momentum $p$ satisfies $1 + \cos p < 2(\Delta + g \cos p)^2$. The bound state consists of two impurities with total momentum $p$ traveling at a typical pseudodistance $1/\log|\Delta + g \cos p/\cos(p/2)|$ with energy $E_2(p) = 4J[\cos^2(p/2) + g^2 \cos^2 p - \Delta^2]/(\Delta + g \cos p)$. The large time limit of $F_t^{(k)}(x)$ takes contributions from both magnons ($m$) and bound states ($b$), $F_t^{(k)}(x) = F_{t,m}^{(k)}(x) + F_{t,b}^{(k)}(x)$. They are worked out in Appendix C and are given by
$$F_{t,m}^{(j)}(x) \sim \int_{-\pi}^{\pi} \frac{d^2 p}{8\pi^2} \theta_H \left(\frac{x}{t} - v^{(j)}(p_1, p_2)\right) \times |e^{i(n_1^{(0)} - n_2^{(0)})(p_1 - p_2)} + S(p_1, p_2)|^2, \quad (16)$$

$$F_{t,b}^{(1)}(x) \sim F_{t,b}^{(2)}(x) \\ \sim \int_{-\pi}^{\pi} \frac{dp}{2\pi} \theta_H \left(\frac{x}{t} - E_2'(p)\right) \\ \times \max \left[ \left( \frac{\Delta + g \cos p}{\cos \frac{p}{2}} \right)^2 - 1, 0 \right] \\ \times \left( \frac{\cos^2 \frac{p}{2}}{(\Delta + g \cos p)^2} \right)^{n_2^{(0)} - n_1^{(0)}}, \quad (17)$$

where

$$v^{(j)}(p_1, p_2) = \begin{cases} \min[E'(p_1), E'(p_2)] & j = 1, \\ \max[E'(p_1), E'(p_2)] & j = 2. \end{cases}$$

In the equations, $n_1^{(0)} < n_2^{(0)}$ denote the pseudopositions of the impurities at the initial time (note that we took the limit $x, t \to \infty$ for given $n_1^{(0)}$ and $n_2^{(0)}$, so the impurities in the initial state correspond to $x \sim 0$).

Figure 2 highlights the effect of the interactions by comparing the cumulative distributions with those obtained

<figure><img src="image_5.png" /><figcaption>FIG. 2. The cumulative distributions for the (symmetric) example in Eq. (13) with $\Delta = 0.5$ and $g = -1$. Each solid curve corresponds to a different value of the pseudodistance $n_2^{(0)} - n_1^{(0)}$ between the impurities at the initial time. The associated dashed lines are obtained by trivializing the scattering phase. The region delimited by each curve represents the contribution from the bound states.</figcaption></figure>

by trivializing the scattering phase, i.e., $S(p_1, p_2) \to -1$. The bound state contributes significantly only when the impurities at the initial times are close enough. The dissimilarity between the maximal velocities of magnons and bound states is manifested in the appearance of cusps within the light cone. For the readers familiar with Bethe ansatz, we specify that those can be interpreted as bare velocities; indeed, the state differs from the reference state by a finite number of excitations. Note, however, that there is an implicit dressing coming from the deformation of the space, i.e., from the relation [Eq. (12)] between $x$ and $\ell$, which can be read as a simplified instance of the effects pointed out in Ref. [51].

### C. Example: Nonsymmetric Hamiltonian

The symmetry of the Hamiltonian shown in Eq. (13) is broken by the interaction

$$H_I = J V \sum_{\ell} \sigma_{\ell}^{z} \sigma_{\ell+1}^{z}, \quad (18)$$

which has the following semilocal representation in the pseudospace:

$$\tilde{H}_I = 4JV \sum_j \sum_{\tilde{j}} \frac{1 - (-1)^{b_j - b_{j+1}}}{2} \int_{-\pi}^{\pi} \frac{d\omega}{2\pi} \\ \times e^{i\omega \left( j - j - 1 + \sum_{j'= -\infty}^{j} \frac{1 + \tau_j^{z}}{2} \right)} \frac{1 - \tau_j^{z}}{2} \frac{1 - \tau_{j+1}^{z}}{2}. \quad (19)$$

Here, the integral over $\omega$ represents a Kronecker delta. In the following, we consider time evolution under $H + H_I$. In the specific case of a single impurity, the total Hamiltonian acting on the pseudospace can be represented by the following noninteracting operator:

$$\tilde{H}_1 \sim J \sum_j \tau_j^x \tau_{j+1}^x + \tau_j^y \tau_{j+1}^y - V[1 - (-1)^{b_j - b_{j+1}}] \tau_j^z, \quad (20)$$

from which we removed an irrelevant term proportional to the $z$ component of the total pseudospin.

#### 1. Quasihomogeneous jammed state

The semilocality of the interaction is masked if the sequence of particles consists of a repeating pattern. In that case, by wrapping the chain around a cylinder in such a way that a complete revolution around its axis corresponds to the periodicity of the initial jammed state, the system can be reinterpreted as a spin ladder with helical interactions and the jammed state is homogeneous along the legs.

In order to emphasize the differences with respect to the symmetric case, we study the same example considered in Ref. [47]: The sequence of species of the underlying jammed state consists of the repeating pattern $\{1, 0, 0, 1\}$ (corresponding to a three-leg helical ladder) and the
perturbation creates an impurity by removing a particle. We can fix the notations in such a way that the initial state is given by $|\cdots\uparrow\uparrow\downarrow\uparrow\downarrow\downarrow\uparrow\uparrow\cdots\rangle \sim |0; b\rangle$, where the sequence of species reads

$$b = \dots, (1, 0, 0, 1), (1, 0, 0, 1_0), (0, 0, 1, 1), (0, 0, 1, 1), \dots \tag{21}$$

Here, we display relative position 0 and group the species to highlight the repeating pattern. The almost periodicity of the species translates into the nearly two-site shift invariance of $\tilde{H}_1$ in Eq. (20), a symmetry that is broken around $j=0$, where the effective staggered spin coupled by V changes sign. Generally, step potentials can give rise to localized excited states. This is confirmed by our numerical investigation, which shows that the initial state has a significant overlap with a localized excited state. Nevertheless, there is a delocalized contribution that triggers MELPs. This is shown in the first row in Fig. 3, where we can also see that the localized excited state creates a discontinuity in the cumulative distribution function (the sharp change of color around $x=0$) that is more pronounced as V is increased.

The problem becomes considerably more complicated when the initial state has two impurities, but it is still reasonable to expect an integrable description almost everywhere and, hence, MELPs. We surmise, however, that the effects could become less pronounced due to the localized inhomogeneities of the interactions.

## 2. Disordered jammed state: Localization

The picture presented in the previous example breaks down if we add a disordered region in the sequence of

species. This creates a region in the pseudospace in which the term coupled by V in Eq. (20) is disordered, which is known to trigger Anderson's localization [52]. As shown in the second row in Fig. 3, for large enough V the information about the impurity is not able to cross the barrier. Thus, the jammed state on one side of the disorder barrier is protected against a localized perturbation on the opposite side, provided that the width of the barrier is large enough (with respect to $1/V^2$).

## D. Macroscopic quantum states

The problem under investigation shares some similarities with the MELP studied in Ref. [48], which considered local perturbations in low-entangled excited states that cannot be interpreted as ground states of local conservation laws. Following Ref. [48] (see also Ref. [46]), let us then imagine to prepare the system in a jammed product state $|\Psi(0^-)\rangle$ (with spins along $\hat{z}$). After identifying a site $\ell_0$ with a spin up whose flip would unjam the state, we perform a projective measurement of $\hat{n} \cdot \hat{\sigma}_{\ell_0}$, with $\hat{n}$ different from $\hat{z}$. The state after the measurement $|\Psi(0^+)\rangle$ will be in a superposition of the original jammed state and an almost identical quasi-jammed state with an impurity, $|\Psi(0^+)\rangle = \sqrt{p}|0; b\rangle + e^{i\phi}\sqrt{1-p}|n_0; b(-)\rangle$ ($\phi$ and $p$ depend on $\hat{n}$), with probability $p$. By the Lieb-Robinson bound [53], the speed at which the information about the impurity spreads is bounded from above; therefore, the state remains separable outside a light cone growing linearly in time. We restrict to the corresponding subsystem $\Omega_t$ (with extent proportional to the time, $|\Omega_t| \propto t$). Let us then consider an extensive operator of the form $S_s^z = \frac{1}{2} \sum_l s_l \sigma_l^z$, with $s_l = \pm 1$, which has the property not to fluctuate outside

<figure><img src="image_6.png" /><figcaption>FIG. 3. The cumulative distribution $F_t^{(1)}(x)$ for a single impurity as a function of pseudoposition (x axis) and time $Jt$ (y axis) in the (nonsymmetric) example described by Eq. (20) starting from $|0; b\rangle$ with the sequence in Eq. (21). Each column corresponds to a different value of V = {0.1, 0.3, 0.5} (the values of $\Delta$ and $g$ are irrelevant). The first row is for the configuration without disorder. The second row is for a configuration with a small disordered region delimited by the two orange vertical lines.</figcaption></figure>
the light cone. Because of that, its variance matches the variance of its restriction to $\Omega_t$.

The variance of such an operator per unit $|\Omega_t|$ provides a lower bound to the minimal system size for which, to use the words of Ref. [54], “one has to assume validity of quantum mechanics in order to explain a nonclassical phenomenon.” Let us then compute the variance of $S_s^z$. Denoting its expectation value and its variance at time $t$ by $\langle S_s^z \rangle_t$ and $\langle S_s^z \rangle_t$, respectively, we have

$$ \langle S_s^z, S_s^z \rangle_t \sim \frac{p}{1-p} (\langle S_s^z \rangle_t - \langle S_s^z \rangle_{0-})^2 + (1-p) \langle n_0; \underline{b}^{(-)} | e^{iHt} S_s^z, S_s^z e^{-iHt} | n_0; \underline{b}^{(-)} \rangle. \quad (22) $$

Assuming Eq. (11), the right-hand side of this equation grows as $t^2$, because the difference in the expectation values grows linearly in time. In fact, we are about to show that, at least in the presence of the symmetry under transmutation of species, the same conclusion applies also to the marginal case with $p=0$, which is the basic problem considered outside this section. The variance in the second line in Eq. (22) can indeed be expressed as follows

$$ \begin{aligned} \langle n_0 | S_s^z(t), S_s^z(t) | n_0 \rangle &= \frac{1}{4} \sum_{\ell} [1 - \langle n_0 | \sigma_{\ell}^z(t) | n_0 \rangle^2] \\ &\quad + \frac{1}{4} \sum_{\ell_1 \neq \ell_2} s_{\ell_1} s_{\ell_2} \left[ \sum_n p_{n,t} (1 - p_{n,t}) \langle n | : \sigma_{\ell_1}^z : : \sigma_{\ell_2}^z : | n \rangle \right. \\ &\quad \left. - \sum_{n \neq n'} p_{n,t} p_{n',t} \langle n | : \sigma_{\ell_1}^z : | n \rangle \langle n' | : \sigma_{\ell_2}^z : | n' \rangle \right], \end{aligned} \quad (23) $$

where the dependency on $\underline{b}^{(-)}$ is understood and the last two lines depend on $n_0$ through $p_{n,t}$. The first term on the right-hand side describes a standard extensive contribution to the variance; the interesting behavior comes from the remaining terms. As discussed before, when the Hamiltonian is symmetric we can generally use the simplified expression

$$ \langle n; \underline{b} | : \sigma_{\ell}^{z} : | n; \underline{b} \rangle \sim \langle 0; \underline{b} | \sigma_{\ell-y}^{z} - \sigma_{\ell}^{z} | 0; \underline{b} \rangle \theta_H(n < x_{\ell}), \quad (24) $$

which results in

$$ \begin{aligned} \langle S_s^z, S_s^z \rangle_t \sim & \frac{1}{4} \sum_{\ell_1 \neq \ell_2} \{ F_t[\min(x_{\ell_1}, x_{\ell_2})] - F_t(x_{\ell_1}) F_t(x_{\ell_2}) \} \\ & \times s_{\ell_1} s_{\ell_2} \langle 0; \underline{b} | \sigma_{\ell_1-y}^z - \sigma_{\ell_1}^z | 0; \underline{b} \rangle \\ & \times \langle 0; \underline{b} | \sigma_{\ell_2-y}^z - \sigma_{\ell_2}^z | 0; \underline{b} \rangle. \end{aligned} \quad (25) $$

Let us then choose $s_{\ell} = \langle 0; \underline{b} | \sigma_{\ell}^z | 0; \underline{b} \rangle$; we find

$$ \begin{aligned} \langle S_s^z, S_s^z \rangle_t \sim & \sum_{\ell_1 \neq \ell_2} \{ F_t[\min(x_{\ell_1}, x_{\ell_2})] - F_t(x_{\ell_1}) F_t(x_{\ell_2}) \} \\ & \times \langle 0; \underline{b} | \frac{1 - \sigma_{\ell_1-y}^z \sigma_{\ell_1}^z}{2} | 0; \underline{b} \rangle \\ & \times \langle 0; \underline{b} | \frac{1 - \sigma_{\ell_2-y}^z \sigma_{\ell_2}^z}{2} | 0; \underline{b} \rangle. \end{aligned} \quad (26) $$

Generally, this grows as $t^2$, i.e., $|\Omega_t|^2$; indeed, $F_t(x)$ approaches a function of $x/t$ and $\langle 0; \underline{b} | 1 - \sigma_{\ell-y}^z \sigma_{\ell}^z / 2 | 0; \underline{b} \rangle$ has a strictly positive local average. This implies that also the second line in Eq. (22) grows quadratically in time.

We can, therefore, conclude that the local perturbation brings $\Omega_t$ to a macroscopic quantum state. That is to say, the MELP under investigation goes beyond classical physics and should be regarded as a signature of quantum jamming.

## IV. DISCUSSION

Macroscopic effects from local perturbations have been recently observed in several settings: ground states in ordered phases with spontaneous symmetry breaking [55,56], symmetric initial states in the presence of semi-local conservation laws [57,58], separable quantum scars [48], and jammed states in constrained systems [46,47]. This work falls within the last category and has lifted the observations of Refs. [46,47] to a new level of generality: We have analytically studied the robustness of MELPs under the inclusion of interactions, no matter if integrable or not; we have identified the key role played by a symmetry and the change in the phenomenology due to its explicit breaking; we have shown that the described phenomenon is genuinely quantum. This was discussed in full generality within a class of spin-chain models that allow for a particle description and exhibit a sector of eigenstates spanned by jammed product states. To explicitly show the most important features identified, we have worked out the dynamics in simple examples.

Quantum jamming appears as a phenomenon that mixes up microscopic with macroscopic, classical with quantum. We discussed a duality transformation that, when looked over in the opposite direction, maps conventional local U(1) symmetric models into constrained local models with jammed sectors, and it brings, in turn, aspects that are typically microscopic, such as scattering phases and formation of bound states, to macroscopic scales. We have seen, in particular, that, in a constrained system, the expectation value of the local magnetization—which is linearly related to its probability distribution—becomes, in the dual unconstrained system, the cumulative distribution.

A question that remains open is about the emergence of a diffusive MELP, in which the effect spreads out as $\sqrt{t}$. This was reported in Ref. [46] starting from a quasijammed product state with a single impurity ($\nu = 1$) but a fluctuating sequence of species. A related question that we have not
addressed is about the first signatures of integrability breaking, which we expect in the sector with three impurities, although some preliminary numerical checks have not revealed any qualitative difference in the behavior of the local magnetization in the presence of more impurities [59].

In the integrable case, it is instead still unclear the relationship between the light cone dynamics emerging because of quantum jamming and those emerging because of the existence of infinitely many conservation laws, which have become the signature of the theory of generalized hydrodynamics.

Finally, we have studied spin chains, but the kinetic constraints that we considered can be easily generalized to higher dimensions. Such systems would still exhibit jammed product states, but what would be the effect of a local perturbation?

# V. METHODS

## A. Explicit construction of the Hamiltonian in the example

We clarify here why the Hamiltonians in the examples belong to the class of models characterized by the building blocks in Eq. (1). Specifically, it turns out that those Hamiltonians have both a linear and a quadratic dependency on the building blocks: Up to a constant and a term proportional to $S^z$, which is conserved and irrelevant, the Hamiltonian in Eq. (13) reads

$$H \sim 2J \sum_{\ell} \bar{T}_{\ell, \ell+2}^{+-} + \bar{T}_{\ell, \ell+2}^{-+} + 2\Delta \{ \bar{T}_{\ell, \ell+3}^{+-}, \bar{T}_{\ell, \ell+3}^{-+} \} \\ + g(T_{\ell, \ell+2}^{-+} T_{\ell-2, \ell}^{-+} + T_{\ell-2, \ell}^{+-} T_{\ell, \ell+2}^{-+}),$$

where we express, in particular, the interaction coupled by $\Delta$ in terms of the averaged building blocks, $\bar{T}_{\ell-1, \ell+2}^{+-} = (T_{\ell-1, \ell+1}^{+-} + T_{\ell, \ell+2}^{+-})/2$, so as to make it explicit the symmetry under transmutation of particles. The symmetry of the last term can be understood as it describes a double hopping in the same direction, completely independent of the parity of the spin up that hops.

The additional interactions that we added to break the symmetry under transmutation of species do not need additional comments, because they were already expressed in terms of the building blocks $\sigma_{\ell}^{z}$.

## B. Duality mapping of shift invariant operators

A key ingredient to solve the dynamics in the class of constrained models that we considered is the Hamiltonian acting on the space of pseudospins. Using translational invariance, an effective procedure to determine the representation of one-site shift invariant operators in the pseudospace is to inspect their action when applied to fragments of a product state that are sufficient to reconstruct the state and that are mapped into elementary fragments in the

pseudospace. For example, the action of the identity (which is trivial but not simply represented), of $\sum_{\ell} \sigma_{\ell}^{z}$ (which is as trivial as the identity because it is conserved), and of $\sum_{\ell} \sigma_{\ell}^{z} \prod_{j=1}^{m-1} (1 - \sigma_{\ell+j/2}) \sigma_{\ell+m}^{z}$ for $m \ge 1$ can by inferred by applying them to a sequence of spins starting with a spin up and ending before the next spin up, $|...\uparrow_j...\downarrow...\downarrow(\uparrow_{j+n+1})...\rangle$. Such a fragment is

indeed mapped into the elementary fragment $|...\downarrow^b \uparrow...\uparrow(\downarrow^{b-n-1} (\text{mod } y))...\rangle$, where arrows represent $[n/y]$

pseudospins. It is then easy to identify the pseudospin operator that acts on the dual fragment in the same way as the original spin operator does on the original fragment. This is explicitly done in Appendix B.

## C. Bethe ansatz with two impurities

In order to solve the two-impurity sector of the Hamiltonian in Eq. (13), we have used the Bethe ansatz

$$\sum_{n_1<n_2} [e^{in_1p_1+in_2p_2} + S(p_1, p_2)e^{in_1p_2+in_2p_1}]|n_1, n_2; \underline{b}\rangle.$$

In the symmetric example, we obtain the scattering phase in Eq. (15) and energy $E(p) = 4J(\cos p - \Delta)$. By applying an artificial finite-volume regularization with periodic boundary conditions, we find that the effective integrable model with scattering phase $S(p_1, p_2)$ exhibits also bound states in the two-particle sector. The corresponding excited states can be parametrized as follows:

$$\sum_{n_1<n_2} e^{i[(n_1+n_2)/2]} p \left( \frac{\cos \frac{p}{2}}{\Delta + g \cos p} \right)^{n_2-n_1} |n_1, n_2; \underline{b}\rangle$$

provided that $1 + \cos p < 2(\Delta + g \cos p)^2$. Their energy reads $E_2(p) = 4J[\cos^2(p/2) + g^2 \cos^2 p - \Delta^2]/(\Delta + g \cos p)$. The reader can find more details in Appendix C, which also includes the derivation of Eqs. (16) and (17).

# ACKNOWLEDGMENTS

I am grateful to Saverio Bocini for sharing some data of time-evolving block decimation simulations (not included here); I also thank him for useful discussions. This work was supported by the European Research Council under Starting Grant No. 805252 LoCoMacro.

# APPENDIX A: EXPECTATION VALUES OF THE LOCAL MAGNETIZATION

Equation (3) links the macroposition of the spins up with their relative position in the sequence. Since, in a product state of spins aligned along $\hat{z}$, $\langle 1 + \sigma_{\ell}^{z} / 2 \rangle$ is nonzero and equal to 1 only if there is $\tilde{j}$ such that $\ell = 2\ell'(\tilde{j}) - b_{\tilde{j}}$, we have
$$ \begin{aligned} \sigma_{\ell}^{z} &= \sigma_{2[\ell]}^{z} - [(-\ell)_{\text{(mod } y\text{)]}} \\ &\rightarrow -1 + 2\sum_{\tilde{j}} \delta_{b_j, (-\ell)_{\text{(mod } y\text{)}}} \int \frac{d\tilde{\omega}}{2\pi} e^{i\tilde{\omega}(\ell'(\tilde{j}) - [\ell/y])}. \end{aligned} \quad (\text{A1}) $$

This expression can be simplified at a fixed number of impurities, as is done in the following.

### 1. Zero impurities: Jammed states

In a jammed state (which is the vacuum $|0\rangle$ in the pseudospace) we have

$$ \ell'(\tilde{j}) = \tilde{j} + \delta\tilde{j}. \quad (\text{A2}) $$

Thus, the local magnetization takes the following value:

$$ \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle = -1 + 2 \sum_{\tilde{j}} \delta_{\ell', \tilde{j}+\delta\tilde{j}} \delta_{b_j, j}. \quad (\text{A3}) $$

### 2. One impurity

We denote by $|n\rangle$ the state with an impurity at position n in the pseudospace. We have

$$ \ell'(\tilde{j}) = \tilde{j} + \delta\tilde{j} + \theta_H(\tilde{j} > n) \quad (\text{A4}) $$

and, hence,

$$ \langle n | \sigma_{y\ell'-j}^{z} | n \rangle = -1 + 2 \sum_{\tilde{j} \in B} \delta_{b_j, j} \delta_{j+\delta\tilde{j}+\theta_H(\tilde{j}>n), \ell'}. \quad (\text{A5}) $$

We are going to express these matrix elements in terms of those in the underlying jammed state using Eq. (A3). We start from the identities

$$ \begin{aligned} \sum_{\tilde{j} > n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j} + 1, \ell'} &= \frac{1 + \langle 0 | \sigma_{y(\ell'-1)-j}^{z} | 0 \rangle}{2} \\ &\quad - \sum_{\tilde{j} \le n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j} + 1, \ell'}. \end{aligned} \quad (\text{A6}) $$

$$ \sum_{\tilde{j} \le n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j}, \ell'} = \frac{1 + \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle}{2} - \sum_{\tilde{j} > n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j}, \ell'}. \quad (\text{A7}) $$

Since $\tilde{j} + \delta\tilde{j}$ is a nondecreasing function of $\tilde{j}$, the sum on the right-hand side of the first equation vanishes if $n + \delta n + 1 < \ell'$. On the other hand, the sum on the left-hand side of the first equation vanishes when $n + 1 + \delta(n + 1) + 1 > \ell'$. An analogous argument applies to the second equation. Using

$$ \delta(n+1) - \delta n = -\theta_H(b_{n+1} < b_n), $$

we finally get

$$ \sum_{\tilde{j} > n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j} + 1, \ell'} = \begin{cases} \frac{1 + \langle 0 | \sigma_{y(\ell'-1)-j}^{z} | 0 \rangle}{2} & n + \delta n + 1 < \ell', \\ 0 & n + \delta n + 1 - \theta_H(b_{n+1} < b_n) \ge \ell', \\ \delta_{b_{n+1}, j} & b_{n+1} < b_n, \ell' = n + \delta n + 1, \end{cases} \quad (\text{A8}) $$

$$ \sum_{\tilde{j} \le n} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j}, \ell'} = \begin{cases} \frac{1 + \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle}{2} & n + \delta n - \theta_H(b_{n+1} < b_n) \ge \ell', \\ 0 & n + \delta n < \ell', \\ \theta_H(b_n < b_{n-1}) \delta_{b_{n-1}, j} + \delta_{b_n, j} & b_{n+1} < b_n, \ell' = n + \delta n, \end{cases} \quad (\text{A9}) $$

and, hence,

$$ \begin{aligned} \langle n | \sigma_{y\ell'-j}^{z} | n \rangle - \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle \\ = \langle 0 | \sigma_{y(\ell'-1)-j}^{z} | 0 \rangle \theta_H(n + \delta n + 1 < \ell') - \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle \theta_H[n + \delta n - \theta_H(b_{n+1} < b_n) < \ell'] \\ + \delta_{\ell', n+\delta n+1} [\theta_H(b_{n+1} < b_n) \delta_{b_{n+1}, j-1}] + \delta_{\ell', n+\delta n} [\theta_H(b_{n+1} < b_n < b_{n-1}) \delta_{b_{n-1}, j} - \theta_H(b_{n+1} < b_n)(1 - \delta_{b_n, j})]. \end{aligned} \quad (\text{A10}) $$

### 3. More impurities: Reduction to the single-impurity case

From Eqs. (3) and (A1), it follows that

$$ \begin{aligned} \langle n_1, n_2, \dots, n_\nu | \sigma_{y\ell'-j}^{z} | n_1, n_2, \dots, n_\nu \rangle &= \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle - 2 \sum_{\tilde{j} \in B} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j}, \ell'} + 2 \sum_{\tilde{j} \in B} \delta_{b_j, j} \delta_{\tilde{j} + \delta\tilde{j}} + \sum_{k=1}^{\nu} \theta_H(n_k+1-k<\tilde{j}), \ell' \\ &= \langle 0 | \sigma_{y\ell'-j}^{z} | 0 \rangle + 2 \sum_{k=1}^{\nu} \sum_{\tilde{j} > n_k+1-k} \delta_{b_j, j} (\delta_{\tilde{j}+\delta\tilde{j}, \ell'-k} - \delta_{\tilde{j}+\delta\tilde{j}, \ell'-k+1}). \end{aligned} \quad (\text{A11}) $$
We can then use

$$ \langle n | \sigma_{y'\ell' - j}^z | n \rangle - \langle 0 | \sigma_{y'\ell' - j}^z | 0 \rangle = 2 \sum_{j>n} \delta_{b_j, j} (\delta_{j+\delta_j, \ell'-1} - \delta_{j+\delta_j, \ell'}) $$

to recast the expectation value as follows:

$$ \langle n_1, n_2, \dots, n_\nu | \sigma_{\ell'}^z : | n_1, n_2, \dots, n_\nu \rangle = \sum_{k=1}^{\nu} \langle n_k + 1 - k | \sigma_{\ell'-y(k-1)}^z : | n_k + 1 - k \rangle, \quad (A12) $$

where we use the notation introduced in Eq. (5) in the text, i.e.,

$$ :O:= O - \langle 0; \underline{b} | O | 0; \underline{b} \rangle. \quad (A13) $$

Note that a similar decomposition holds in standard U(1) symmetric models in which $n_j$ simply represents the position of the jth spin up; the main difference is that here both the position of the operator and the pseudpositions of the impurities are shifted depending on the number of the impurities on their left.

## APPENDIX B: MAPPING OF DIAGONAL OPERATORS

In this section, we discuss the mapping of diagonal operators (in a basis diagonalizing $\sigma^z$) onto the pseudospace. It is convenient to split them in classes depending on the minimal number of consecutive domains of spins down that are required to characterize their local action.

### 1. Class I

The simplest operators are determined by their action on each fragment of a product state of the form

$$ |\dots \uparrow_j \downarrow \dots \downarrow (\uparrow_{j+n+1}) \dots \rangle. \quad (B1) $$

The corresponding contribution to their eigenvalue is summarized in Table I. The duality mapping (from spins to pseudospins) acts as follows:

<table><caption>TABLE I. Local contribution to the eigenvalues of some class I operators ordered by their range.</caption><thead><tr><th></th><th>|...↑<sub>j</sub>↓...↓(↑<sub>j+n+1</sub>)...⟩</th></tr></thead><tbody><tr><td>1</td><td>n + 1</td></tr><tr><td>σ<sup>z</sup> ⊗ σ<sup>z</sup></td><td>1 - n</td></tr><tr><td>σ<sup>z</sup> ⊗ 1/2 σ<sup>z</sup> ⊗ σ<sup>z</sup></td><td>n - 3 + 4δ<sub>n,0</sub></td></tr><tr><td>σ<sup>z</sup> ⊗ (1/2 σ<sup>z</sup>) ⊗ 2 σ<sup>z</sup></td><td>n - 4 + 4(δ<sub>n,0</sub> + δ<sub>n,1</sub>)</td></tr><tr><td>σ<sup>z</sup> ⊗ (1/2 σ<sup>z</sup>) ⊗ 3 σ<sup>z</sup></td><td>n - 5 + 4(δ<sub>n,1</sub> + δ<sub>n,2</sub>) + 5δ<sub>n,0</sub></td></tr><tr><td>σ<sup>z</sup> ⊗ (1/2 σ<sup>z</sup>) ⊗ 4 σ<sup>z</sup></td><td>n - 6 + 4(δ<sub>n,3</sub> + δ<sub>n,4</sub>) + 5δ<sub>n,1</sub> + 6δ<sub>n,0</sub></td></tr></tbody></table>

$$ |\dots \uparrow_j \downarrow \dots \downarrow (\uparrow_{j+n+1}) \dots \rangle \mapsto |\dots \downarrow^b \uparrow \dots \uparrow (\downarrow^{b-n-1} (\text{mod } y)) \dots \rangle, \quad (B2) $$

where, on the right-hand side, the superscripts of the down arrows label the species of particles, which we call $b$ for the first spin down. On the contrary, if the image of the state is $|\dots \downarrow^b \uparrow \dots \uparrow (\downarrow^{b-j+1}) \dots \rangle$, the number $n$ of spins down in the original fragment is

$$ \begin{aligned} n &= ym + (b_j - b_{j+1} - 1) (\text{mod } y) \\ &= ym + \frac{y-1}{2} - \frac{1}{2} \sum_{n=1}^{y-1} \frac{\sin\left(\frac{2\pi n(b_j - b_{j+1} - \frac{1}{2})}{y}\right)}{\sin\left(\frac{n\pi}{y}\right)}, \end{aligned} \quad (B3) $$

where we use

$$ x (\text{mod } y) = \frac{y-1}{2} - \frac{1}{2} \sum_{n=1}^{y-1} \frac{\sin\left(\frac{2\pi n(x+\frac{1}{2})}{y}\right)}{\sin\left(\frac{n\pi}{y}\right)}. \quad (B4) $$

The mapping of the operators in Table I can be worked out using Eq. (B3) and the action of the diagonal operators in the pseudospace reported in Table II. We recognize three kinds of contributions to the corresponding eigenvalues.

**Additive constant:** The contribution of an additive constant $c$ in the expressions appearing in the second column in Table I is given by the number of spins up multiplied by $c$. This is represented in the pseudospace by

$$ c \rightarrow c \sum_j \frac{1 - \tau_j^z}{2}. \quad (B5) $$

In fact, the total pseudospin in the $z$ direction is conserved and depends only on the sequence of species. In the following, we label by $B(\sim\mathbb{Z})$ the set of the relative positions and call “configuration space” the associated space. In order to reduce confusion, we instead call $I(\sim\mathbb{Z})$ the set of the pseudopositions. Using these notations, we see that

<table><caption>TABLE II. Local contribution to the eigenvalues of some diagonal operators in the pseudospace.</caption><thead><tr><th></th><th>|...↓<sup>b</sup>↑...↑(↓<sup>b-n-1</sup> (mod y))...⟩</th></tr></thead><tbody><tr><td>1/2 π̇<sup>z</sup></td><td>1</td></tr><tr><td>1/2 π̇<sup>z</sup></td><td>1/2 π̇<sup>z</sup></td></tr><tr><td>1/2 π̇<sup>z</sup> ⊗ (1/2 π̇<sup>z</sup>) ⊗ k</td><td>δ<sub>[y],k</sub></td></tr></tbody></table>
$$ \frac{1 - \tau_{j-k_1-1}^z}{2} \left( \prod_{j'=1}^{k_1} \frac{1 + \tau_{j-j'}^z}{2} \right) \frac{1 - \tau_j^z}{2} \left( \prod_{j''=1}^{k_2} \frac{1 + \tau_{j+j''}^z}{2} \right) \frac{1 - \tau_{j+k_2+1}^z}{2} |\cdots(\downarrow^{b+n_-+1 \pmod y})\underbrace{\uparrow\cdots\uparrow}_{[\frac{n_-}{y}]}\downarrow^b\uparrow\cdots\uparrow(\downarrow^{b-n_+-1 \pmod y})\cdots\rangle \\ = \delta_{[\frac{n_-}{y}], k_1} \delta_{[\frac{n_+}{y}], k_2} |\cdots(\downarrow^{b+n_+-1 \pmod y})\underbrace{\uparrow\cdots\uparrow}_{[\frac{n_-}{y}]}\downarrow^b\uparrow\cdots\uparrow(\downarrow^{b-n_+-1 \pmod y})\cdots\rangle. \tag{B15} $$

*Product of two Kronecker deltas* $\delta_{n_-,k_1}\delta_{n_+,k_1}$: This kind of contribution can be worked out by properly multiplying the representation of the delta in Eq. (B8):

$$ \delta_{n_-, k_1} \delta_{n_+, k_2} = \delta_{[\frac{n_-}{y}], [\frac{k_1}{y}]} \delta_{[\frac{n_+}{y}], [\frac{k_2}{y}]} \frac{1}{y^2} \sum_{m,m'=0}^{y-1} e^{2\pi i \frac{m_1(b_k-b_{k_1+1}-s_1)+m_2(b_{k_1+1}-b_{k_2+2-s_2})}{y}} e^{2\pi i \frac{b_1-b_{j+1}-1-k_2}{y}} e^{2\pi i \frac{b_1-b_{j-1}-1-k_1}{y}}. \tag{B16} $$

Thus, we find

$$ \delta_{n_-, k_1} \delta_{n_+, k_2} \rightarrow \sum_{j \in I} (Z_{1+k_1 \pmod y, 1+k_2 \pmod y}^{([\frac{k_1}{y}, \frac{k_2}{y}])})_j, \tag{B17} $$

where $(Z_{s_1 \pmod y, s_2 \pmod y}^{(n_1, n_2)})_j$ is the semilocal operator

$$ (Z_{s_1 \pmod y, s_2 \pmod y}^{(n_1, n_2)})_j = \int_{-\pi}^{\pi} \frac{d\omega}{2\pi} \sum_{k \in B} \sum_{m_1, m_2=0}^{y-1} \frac{e^{2\pi i \frac{m_1(b_k-b_{k_1+1}-s_1)+m_2(b_{k_1+1}-b_{k_2+2-s_2})}{y}}}{y^2} e^{i\omega(k+\sum_{j=1}^0 1-\sum_{j'=1}^{y-1} \frac{1-\tau_j^z}{y})} \\ \times \frac{1-\tau_j^z}{2} \left( \prod_{m=1}^{n_1-1} \frac{1+\tau_{j+m}^z}{2} \right) \frac{1-\tau_{j+n_1}^z}{2} \left( \prod_{m=1}^{n_2-1} \frac{1+\tau_{j+n_1+m}^z}{2} \right) \frac{1-\tau_{j+n_1+n_2}^z}{2}. \tag{B18} $$

We point out that identity (B11) is generalized in that the sum over all $s$ results in a (strictly) local operator.

### 3. Class III and beyond

The explicit calculations for the operators in classes I and II have revealed the general structure of the mapping. Specifically, an operator with a density consisting of the tensor product of projectors $P_\pm$ sandwiched between two $\sigma^z$ is class $M$ if there is a number $M-1$ of $P_\pm$. The eigenvalue will then be written as a linear combination of Kronecker deltas and constants, the highest product involving $M$ Kronecker deltas. The species associated with the deltas appearing in the corresponding representation (B8) should be chosen according to which species surround the associated domains of spins down. The product of deltas is then represented by generalizations of Eq. (B10) along the line of Eq. (B18).

## APPENDIX C: INTEGRABLE REPRESENTATION IN THE PRESENCE OF TWO IMPURITIES

In this section, we provide more details about the two-particle sector of the symmetric Hamiltonian shown in Eq. (13). We have

$$ \begin{aligned} :\tilde{H}: |n_1, n_2\rangle ={}& 2J[|n_1 - 1, n_2\rangle + (1 - \delta_{n_2, n_1+1})(|n_1 + 1, n_2\rangle \\ & + |n_1, n_2 - 1\rangle) + |n_1, n_2 + 1\rangle] \\ & + 4J\Delta(\delta_{n_2, n_1+1} - 2)|n_1, n_2\rangle \\ & + 2Jg\delta_{n_2, n_1+1}(|n_1 - 1, n_1\rangle \\ & + |n_1 + 1, n_1 + 2\rangle), \end{aligned} \tag{C1} $$

where the sequence $\underline{b}$ is understood. By inspecting the single-impurity problem, we identify the following single-particle eigenstates:

$$ |\Psi(p)\rangle = \sum_{n} e^{inp} |n\rangle \tag{C2} $$

normalized according to $\langle\Psi(p_1)|\Psi(p_2)\rangle = 2\pi\delta(p_1 - p_2)$; the corresponding excitation energy reads

$$ E(p) = 4J(\cos p - \Delta). \tag{C3} $$

The Bethe ansatz for the two-impurity problem can be written as follows:
$$|\Psi(p_1, p_2)\rangle \sim \sum_{n_1<n_2} [e^{in_1p_1+in_2p_2} + S(p_1, p_2)e^{in_1p_2+in_2p_1}] \\ \times |n_1, n_2\rangle, \qquad (C4)$$

where $S(p_1, p_2)$ is the scattering phase. Enforcing the secular equation gives

$$S(p_1, p_2) = - \frac{1 - 2[\Delta + g \cos(p_1 + p_2)]e^{ip_2} + e^{i(p_1+p_2)}}{1 - 2[\Delta + g \cos(p_1 + p_2)]e^{ip_1} + e^{i(p_1+p_2)}}. \qquad (C5)$$

We note that, for $g \neq 0$, it is not possible to define rapidities $\lambda_{1,2}$ so as to make the scattering phase of the difference form $S(\lambda_1 - \lambda_2)$. For the sake of transparency, we apply an artificial finite-volume regularization with periodic boundary conditions. The momenta satisfy, in turn, the Bethe equations

$$e^{iLp_1} = S(p_2, p_1), \quad e^{iLp_2} = S(p_1, p_2). \qquad (C6)$$

The first effect of interaction is that these equations exhibit both real and complex solutions, the latter ones describing bound states of two particles. Specifically, the real solutions

differ from the noninteracting ones in a negligible $O(1/L)$ correction that is nevertheless responsible for the incompleteness of the spanned space (it leaves out a one-dimensional space of solutions in the two-dimensional space of momenta).

In the following, we evaluate the late time behavior of the two cumulative distributions

$$F_t^{(1)}(n) = \sum_{n_1 \le n} p_{n_1,t}^{(1)} \\ = \sum_{n_1 \le n} \sum_{n_2 > n_1} |\langle n_1, n_2 | e^{-i\tilde{H}_{\text{eff}}t} \tilde{P} | \Psi_0 \rangle|^2, \qquad (C7)$$

$$F_t^{(2)}(n) = \sum_{n_2 \le n} p_{n_2,t}^{(2)} \\ = \sum_{n_1 < n_2 \le n} |\langle n_1, n_2 | e^{-i\tilde{H}_{\text{eff}}t} \tilde{P} | \Psi_0 \rangle|^2. \qquad (C8)$$

## 1. Real solutions

We start by considering the contribution from the real solutions. Let us define the non-normalized states in a symmetric way:

$$|\Psi(p_1, p_2)\rangle = \sum_{-\frac{L}{2} < n_1 < n_2 \le \frac{L}{2}} [A_{p_1, p_2} e^{in_1 p_1 + in_2 p_2} - A_{p_2, p_1} e^{in_1 p_2 + in_2 p_1}] |n_1, n_2\rangle, \qquad (C9)$$

where

$$A_{p_1, p_2} = 1 - 2[\Delta + g \cos(p_1 + p_2)]e^{ip_1} + e^{i(p_1+p_2)}, \qquad (C10)$$

$$A_{-p_1, -p_2} = A_{p_2, p_1} e^{-i(p_1+p_2)}. \qquad (C11)$$

In the thermodynamic limit such states are normalized according to

$$\langle \Psi(p_1, p_2) | \Psi(p'_1, p'_2) \rangle = (2\pi)^2 \delta(p'_1 - p_1) \delta(p'_2 - p_2) |A(p_1, p_2)|^2, \qquad (C12)$$

with finite-size corrections that scale as $L^{-1}$. In addition, their density is asymptotically flat; indeed, for a smooth test function $f(p_j)$ and $\epsilon$ approaching zero more slowly than $1/L$, we have $(\bar{j} = 3 - j)$

$$\sum_{\substack{p_j \\ \text{Bethe eqs}}} f(p_j) = \oint_{[-\pi, \pi]_e} \frac{dp}{2\pi} \frac{L - i\partial_p \log S(p, p_j)}{1 - e^{-iLp} S(p_j, p)} f(p) \to L \int_{-\pi}^{\pi} \frac{dp}{2\pi} \left[ 1 - \frac{i\partial_p \log S(p, p_j)}{L} \right] f(p), \qquad (C13)$$

where $[-\pi, \pi]_e$ is the rectangular contour surrounding the interval $[-\pi, \pi]$ with imaginary part $\pm ie$. We can conclude that the density is flat, observing that the correction associated with the scattering phase approaches zero as $O(L^{-1})$. We are now in a position to compute the contribution $F_{t,m}^{(j)}(n)$ to $F_{t,m}^{(j)}(n)$ of the excitations with real momenta. The antisymmetry of the state under the interchange of $p_1$ and $p_2$ allows us to express them as follows:
$$ F_{t,m}^{(1)}(n) = \frac{1}{4} \sum_{m_1=0}^{\infty} \sum_{m=1}^{\infty} \left| \int \frac{d^2 p}{(2\pi)^2} e^{-i[E(p_1)+E(p_2)]t} \frac{\langle -m_1+n, -m_1+n+m | \Psi(p_1, p_2) \rangle \langle \Psi(p_1, p_2) | n_1^{(0)}, n_2^{(0)} \rangle}{|A(p_1, p_2)|^2} \right|^2 \\ = \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p d^2 p'}{(2\pi)^4} \frac{e^{-i[E(p_1)+E(p_2)-E(p_1')-E(p_2')]t} e^{in(p_1+p_2-p_1'-p_2')} \frac{\langle \Psi(p_1, p_2) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi(p_1', p_2') \rangle}{|A(p_1, p_2)|^2 |A(p_1', p_2')|^2} \\ \times \left[ \frac{A(p_1, p_2)A(-p_1', -p_2')}{e^{-i(p_2-p_2'+i\epsilon)} - 1} - \frac{A(p_2, p_1)A(-p_1', -p_2')}{e^{-i(p_1-p_2'+i\epsilon)} - 1} \right]. \quad (C14) $$

It is convenient to shift the integration variables in such a way that only one momentum appears in the phases at the denominators. This, in turn, allows us to carry them out separately; specifically, we get

$$ F_{t,m}^{(1)}(n) = \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p d^2 p'}{(2\pi)^4} \frac{e^{-i[E(p_1)+E(p_2)-E(p_1'+p_1-p_2')-E(p_2'+p_2)]t} e^{-in p_1'}}{1 - e^{i(p_1'+i\epsilon)}} \\ \times \frac{\langle \Psi(p_1, p_2) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi(p_1' + p_1 - p_2', p_2' + p_2) \rangle}{|A(p_1, p_2)|^2 |A(p_1' + p_1 - p_2', p_2' + p_2)|^2} \frac{A(p_1, p_2)A(-p_1' - p_1 + p_2', -p_2' - p_2)}{e^{i(p_2'-i\epsilon)} - 1} \\ - \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p d^2 p'}{(2\pi)^4} \frac{e^{-i[E(p_1)+E(p_2)-E(p_1'+p_2-p_2')-E(p_2'+p_2)]t} e^{-in p_1'}}{1 - e^{i(p_1'+i\epsilon)}} \\ \times \frac{\langle \Psi(p_1, p_2) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi(p_1' + p_2 - p_2', p_2' + p_1) \rangle}{|A(p_1, p_2)|^2 |A(p_1' + p_2 - p_2', p_2' + p_1)|^2} \frac{A(p_2, p_1)A(-p_1' - p_2 + p_2', -p_2' - p_1)}{e^{i(p_2'-i\epsilon)} - 1} \\ \sim \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p}{(2\pi)^2} \theta_H(n - \min(v(p_1), v(p_2))t) \frac{\langle \Psi(p_1, p_2) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi(p_1, p_2) \rangle}{|A(p_1, p_2)|^2}. \quad (C15) $$

Thus, we have

$$ F_{t,m}^{(1)}(n) \sim \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p}{(2\pi)^2} \theta_H \{n - \min[v(p_1), v(p_2)]t\} |e^{i(n_1^{(0)} - n_2^{(0)})(p_1 - p_2)} + S(p_1, p_2)|^2 \quad (C16) $$

and analogously for $F_{t,m}^{(2)}(n)$

$$ F_{t,m}^{(2)}(n) \sim \frac{1}{2} \int_{-\pi}^{\pi} \frac{d^2 p}{(2\pi)^2} \theta_H \{n - \max[v(p_1), v(p_2)]t\} |e^{i(n_1^{(0)} - n_2^{(0)})(p_1 - p_2)} + S(p_1, p_2)|^2. \quad (C17) $$

The incompleteness of the real solutions is evident if we take the limit $n \to \infty$, which gives

$$ \lim_{n \to \infty} F_{t,m}^{(j)}(n) \sim 1 + \operatorname{Re} \int_{-\pi}^{\pi} \frac{d^2 p}{(2\pi)^2} e^{i(n_2^{(0)} - n_1^{(0)})(p_1 - p_2)} S(p_1, p_2). \quad (C18) $$

We see before long that the bound-state contribution cancels out the displacement from 1.

## 2. Bound states

Except for the region $g = \Delta \in (-\frac{1}{2}, \frac{1}{2})$, the Bethe equations exhibit always complex solutions. In the simple case we are considering, such solutions can be easily worked out

imposing that, in the Bethe equations, $|e^{iLp_j}|$ approaches either 0 or $\infty$. In particular, we find that the bound state can have any momentum within the region

$$ 1 + \cos p < 2(\Delta + g \cos p)^2, \quad (C19) $$

i.e.,

$$ \cos p > \frac{1 - 4g\Delta + \sqrt{1 + 8g^2 - 8g\Delta}}{4g^2} \lor \cos p < \frac{1 - 4g\Delta - \sqrt{1 + 8g^2 - 8g\Delta}}{4g^2}. \quad (C20) $$
Note that this condition is always satisfied when $g\Delta > \frac{1}{8} + g^2$ or $1 < (\Delta + g)^2 < 4g(\Delta + g) < 1$. The corresponding bound state reads

$$|\Psi_2(p)\rangle \sim \sqrt{\left(\frac{\Delta + g \cos p}{\cos \frac{p}{2}}\right)^2 - 1} \times \sum_{n_1 < n_2} e^{i[(n_1+n_2)/2]p} \left(\frac{\cos \frac{p}{2}}{\Delta + g \cos p}\right)^{n_2-n_1} |n_1, n_2\rangle, \quad (C21)$$

where we enforce the standard normalization $\langle\Psi(p')|\Psi(p)\rangle = 2\pi\delta(p - p')$. Its energy is given by

$$E_2(p) = 4g \cos p + 4J \left[ \frac{\cos^2 \frac{p}{2}}{\Delta + g \cos p} - \Delta \right]. \quad (C22)$$

We can now compute the bound-state contribution to the cumulative distributions. We have

$$\begin{aligned} F_{t,b}^{(1)}(n) &= \sum_{m_1=0}^{\infty} \sum_{m=1}^{\infty} \left| \int_{-\pi}^{\pi} \frac{dp}{2\pi} \theta_H[1 + \cos p < 2(\Delta + g \cos p)^2] e^{-iE_2(p)t} \langle -m_1 + n, -m_1 + n + m | \Psi_2(p) \rangle \langle \Psi_2(p) | n_1^{(0)}, n_2^{(0)} \rangle \right|^2 \\ &= \int_{-\pi}^{\pi} \frac{dp}{(2\pi)^2} \langle \Psi_2(p) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi_2(p' + p) \rangle \theta_H[\cos^2 \frac{p}{2} < (\Delta + g \cos p)^2] \\ &\quad \times \theta_H \left( \cos^2 \frac{p' + p}{2} < [\Delta + g \cos(p' + p)]^2 \right) \\ &\quad \times \sqrt{\left(\frac{\Delta + g \cos p}{\cos \frac{p}{2}}\right)^2 - 1} \sqrt{\left(\frac{\Delta + g \cos(p' + p)}{\cos \frac{p' + p}{2}}\right)^2 - 1} \frac{e^{-inp'}}{1 - e^{i(p'+ie)}} \frac{e^{-i[E_2(p) - E_2(p'+p)]t}}{\frac{[\Delta + g \cos p][\Delta + g \cos(p'+p)]}{e^{-ip'} \cos^2 \frac{p}{2} \cos \frac{p'+p}{2}}} - 1 \\ &\sim \int_{-\pi}^{\pi} \frac{dp}{2\pi} \theta_H[1 + \cos p < 2(\Delta + g \cos p)^2] \theta_H[n - v_2(p)] \langle \Psi_2(p) | n_1^{(0)}, n_2^{(0)} \rangle \langle n_1^{(0)}, n_2^{(0)} | \Psi_2(p) \rangle. \end{aligned} \quad (C23)$$

Thus, we have

$$F_{t,b}^{(1)}(n) \sim \int_{-\pi}^{\pi} \frac{dp}{2\pi} \theta_H[n - v_2(p)t] \\ \times \max \left[ \left( \frac{\Delta + g \cos p}{\cos \frac{p}{2}} \right)^2 - 1, 0 \right] \\ \times \left( \frac{\cos^2 \frac{p}{2}}{[\Delta + g \cos p]^2} \right)^{n_2^{(0)} - n_1^{(0)}}. \quad (C24)$$

Concerning $F_{t,b}^{(2)}(n)$, as one could expect from physical arguments based on the fact that the impurities are bounded, we find $F_{t,b}^{(2)}(t) \sim F_{t,b}^{(1)}(n)$, which we can therefore simply call $F_{t,b}(n)$.

Importantly, despite the different appearance of the integrals, the limit $n \to \infty$ exactly compensates the displacement from 1 that we get in Eq. (C18) considering only the real solutions, supporting the completeness of the Bethe states identified in this sector.

[1] T. Giamarchi, *Quantum Physics in One Dimension* (Oxford University, New York, 2003).

[2] F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, *The One-Dimensional Hubbard Model* (Cambridge University Press, Cambridge, England, 2005).

[3] B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen, *Quantum Information Meets Quantum Matter* (Oxford University, New York, 2019).

[4] I. Bloch, J. Dalibard, and W. Zwerger, *Many-body physics with ultracold gases*, Rev. Mod. Phys. **80**, 885 (2008).

[5] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, *One dimensional bosons: From condensed matter systems to ultracold gases*, Rev. Mod. Phys. **83**, 1405 (2011).

[6] I. Buluta and F. Nori, *Quantum simulators*, Science **326**, 108 (2009).

[7] J. I. Cirac and P. Zoller, *Goals and opportunities in quantum simulation*, Nat. Phys. **8**, 264 (2012).

[8] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, *Probing many-body dynamics on a 51-atom quantum simulator*, Nature (London) **551**, 579 (2017).

[9] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, *Colloquium: Nonequilibrium dynamics of closed interacting quantum systems*, Rev. Mod. Phys. **83**, 863 (2011).

[10] C. Gogolin and J. Eisert, *Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems*, Rep. Prog. Phys. **79**, 056001 (2016).

[11] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, *Relaxation and prethermalization in an isolated quantum system*, Science **337**, 1318 (2012).
[12] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, *Generalized hydrodynamics on an atom chip*, Phys. Rev. Lett. **122**, 090601 (2019).

[13] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, *Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1D lattice hard-core bosons*, Phys. Rev. Lett. **98**, 050405 (2007).

[14] F. H. L. Essler and M. Fagotti, *Quench dynamics and relaxation in isolated integrable quantum spin chains*, J. Stat. Mech. (2016) 064002.

[15] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, *Quasilocal charges in integrable lattice systems*, J. Stat. Mech. (2016) 064008.

[16] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, *Complete generalized Gibbs ensembles in an interacting theory*, Phys. Rev. Lett. **115**, 157201 (2015).

[17] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, *Emergent hydrodynamics in integrable quantum systems out of equilibrium*, Phys. Rev. X **6**, 041065 (2016).

[18] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, *Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents*, Phys. Rev. Lett. **117**, 207201 (2016).

[19] S.-J. Ran, E. Tirrito, C. Peng, X. Chen, L. Tagliacozzo, G. Su, and M. Lewenstein, *Tensor Network Contractions Methods and Applications to Quantum Many-Body Systems*, 1st ed., Lecture Notes in Physics Vol. 964 (Springer Nature, Cham, 2020).

[20] M. Moeckel and S. Kehrein, *Interaction quench in the Hubbard model*, Phys. Rev. Lett. **100**, 175702 (2008).

[21] M. Rigol, *Breakdown of thermalization in finite one-dimensional systems*, Phys. Rev. Lett. **103**, 100403 (2009).

[22] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson, *Prethermalization and thermalization in models with weak integrability breaking*, Phys. Rev. Lett. **115**, 180601 (2015).

[23] M. Babadi, E. Demler, and M. Knap, *Far-from-equilibrium field theory of many-body quantum spin systems: Prethermalization and relaxation of spin spiral states in three dimensions*, Phys. Rev. X **5**, 041005 (2015).

[24] T. Langen, T. Gasenzer, and J. Schmiedmayer, *Prethermalization and universal dynamics in near-integrable quantum systems*, J. Stat. Mech. (2016) 064009.

[25] D. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, *A rigorous theory of many-body prethermalization for periodically driven and closed quantum systems*, Commun. Math. Phys. **354**, 809 (2017).

[26] V. Alba and M. Fagotti, *Prethermalization at low temperature: The scent of long-range order*, Phys. Rev. Lett. **119**, 010601 (2017).

[27] P. Reimann and L. Dabelow, *Typicality of prethermalization*, Phys. Rev. Lett. **122**, 080603 (2019).

[28] J. Durnin, M. J. Bhaseen, and B. Doyon, *Nonequilibrium dynamics and weakly broken integrability*, Phys. Rev. Lett. **127**, 130601 (2021).

[29] M. Fagotti, *On conservation laws, relaxation and prerelaxation after a quantum quench*, J. Stat. Mech. (2014) P03016.

[30] B. Bertini and M. Fagotti, *Pre-relaxation in weakly interacting models*, J. Stat. Mech. (2015) P07012.

[31] L. Zadnik and J. P. Garrahan, *Slow heterogeneous relaxation due to constraints in dual XXZ models*, Phys. Rev. B **108**, L100304 (2023).

[32] A. H. MacDonald, S. M. Girvin, and D. Yoshioka, *t/U expansion for the Hubbard model*, Phys. Rev. B **37**, 9753 (1988).

[33] L. Zadnik and M. Fagotti, *The folded spin-1/2 XXZ model: I. Diagonalisation, jamming, and ground state properties*, SciPost Phys. Core **4**, 010 (2021).

[34] S. Moudgalya and O. I. Motrunich, *Hilbert space fragmentation and commutant algebras*, Phys. Rev. X **12**, 011050 (2022).

[35] G. Biroli, C. Chamon, and F. Zamponi, *Theory of the superglass phase*, Phys. Rev. B **78**, 224306 (2008).

[36] Z. Nussinov, P. Johnson, M. J. Graf, and A. V. Balatsky, *Mapping between finite temperature classical and zero temperature quantum systems: Quantum critical jamming and quantum dynamical heterogeneities*, Phys. Rev. B **87**, 184202 (2013).

[37] C. Artiaco, F. Balducci, G. Parisi, and A. Scardicchio, *Quantum jamming: Critical properties of a quantum mechanical perceptron*, Phys. Rev. A **103**, L040203 (2021).

[38] G. Biroli, *A new kind of phase transition?*, Nat. Phys. **3**, 222 (2007).

[39] C. Weitenberg, M. Endres, J. F. Sherson, M. Cheneau, P. Schauß, T. Fukuhara, I. Bloch, and S. Kuhr, *Single-spin addressing in an atomic Mott insulator*, Nature (London) **471**, 319 (2011).

[40] T. Došlić, M. Puljiz, S. Šebek, and J. Žubrinić, *Complexity function of jammed configurations of Rydberg atoms*, arXiv:2302.08791.

[41] L. Gotta, L. Mazza, P. Simon, and G. Roux, *Two-fluid coexistence in a spinless fermions chain with pair hopping*, Phys. Rev. Lett. **126**, 206805 (2021).

[42] Ž. Krajnik, J. Schmidt, V. Pasquier, T. Prosen, and E. Ilievski, *Universal anomalous fluctuations in charged single-file systems*, Phys. Rev. Res. **6**, 013260 (2024).

[43] M. Borsi, L. Pristyák, and B. Pozsgay, *Matrix product symmetries and breakdown of thermalization from hard rod deformations*, Phys. Rev. Lett. **131**, 037101 (2023).

[44] B. Pozsgay, T. Gombor, and A. Hutsalyuk, *Integrable hard-rod deformation of the Heisenberg spin chains*, Phys. Rev. E **104**, 064124 (2021).

[45] R. Z. Bariev, *Integrable spin chain with two- and three-particle interactions*, J. Phys. A **24**, L549 (1991).

[46] K. Bidzhiev, M. Fagotti, and L. Zadnik, *Macroscopic effects of localized measurements in jammed states of quantum spin chains*, Phys. Rev. Lett. **128**, 130603 (2022).

[47] L. Zadnik, S. Bocini, K. Bidzhiev, and M. Fagotti, *Measurement catastrophe and ballistic spread of charge density with vanishing current*, J. Phys. A **55**, 474001 (2022).

[48] S. Bocini and M. Fagotti, *Growing Schrödinger's cat states by local unitary time evolution of product states*, arXiv: 2210.15585.

[49] As it is written, this Hamiltonian is symmetric only in the weak sense that the dependency on the species is a trivial constant fixed by the sequence *b*, the exact symmetry
requiring one also to add an irrelevant constant and a coupling to $S^z$ (omitted here for the sake of compactness).

[50] H. Bethe, Z. Theor. Met. Z. Phys. **71**, 205 (1931).

[51] L. Bonnes, F. H. L. Essler, and A. M. Läuchli, "Light-cone" dynamics after quantum quenches in spin chains, Phys. Rev. Lett. **113**, 187203 (2014).

[52] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. **109**, 1492 (1958).

[53] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, Commun. Math. Phys. **28**, 251 (1972).

[54] F. Fröwis and W. Dür, Measures of macroscopicity for quantum spin systems, New J. Phys. **14**, 093039 (2012).

[55] V. Zauner, M. Ganahl, H. G. Evertz, and T. Nishino, Time evolution within a comoving window: Scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains, J. Phys. Condens. Matter **27**, 425602 (2015).

[56] V. Eisler and F. Maislinger, Front dynamics in the XY chain after local excitations, SciPost Phys. **8**, 037 (2020).

[57] M. Fagotti, Global quenches after localized perturbations, Phys. Rev. Lett. **128**, 110602 (2022).

[58] M. Fagotti, V. Marić, and L. Zadnik, Nonequilibrium symmetry-protected topological order: emergence of semi-local Gibbs ensembles, Phys. Rev. B **109**, 115117 (2024).

[59] S. Bocini (private communication).