# Spectroscopic performance of Low-Gain Avalanche Diodes for different types
of radiation

Gabriele Giacomini  Wei Chen  Gabriele D’Amen  Enrico Rossi  Alessandro
Tricoli  Instrumentation Division, Brookhaven National Laboratory, Upton,
11973 New York, United States  Physics Department, Brookhaven National
Laboratory, Upton, 11973 New York, United States

###### Abstract

Low-Gain Avalanche Diodes are a type of silicon Avalanche Photo-Diodes
originally developed for the fast detection of minimum ionizing particles in
high-energy physics experiments. Thanks to their fast timing performance, the
Low-Gain Avalanche Diode paradigm enables detectors to accurately measure
minimum ionizing particles with a timing resolution of a few tens of
picoseconds. Such a performance is due to a thin substrate and the presence of
a moderate signal gain. This internal gain of a few tens is enough to
compensate for the reduced charge deposition in the thinner substrate and the
noise of fast read-out systems. While Low-Gain Avalanche Diodes are optimized
for the detection of minimum ionizing particles for high-energy particle
detectors, it is critical to study their performance for the detection of
different types of particle, such as X-rays, gamma-rays, or alphas. In this
paper, we evaluate the gain of three types of Low-Gain Avalanche Diodes: two
devices with different geometries and doping profiles fabricated by Brookhaven
National Laboratory, and one fabricated by Hamamatsu Photonics with a
different process.

Since the gain in LGADs depends on the bias voltage applied to the sensor,
pulse-height spectra have been acquired for bias voltages spanning from the
depletion voltage up to breakdown voltage. The signal-to-noise ratio of the
generated signals and the shape of their spectra allow us to probe the
underlying physics of the multiplication process.

###### keywords:

Low-Gain Avalanche Diodes , LGAD , Avalanche Diodes , APD , avalanche noise ,
gain , charge-sensitive amplifier

##  1 Introduction

Low-Gain Avalanche Diodes (LGADs) [1] differ from regular Avalanche Photo-
Diodes (APDs) as they are optimized for the detection of minimum ionizing
particles (mips). Their development originated from the need for precisely
measuring the timing of mips in High-Energy Physics (HEP) experiments at the
Large Hadron Collider (LHC) to disentangle the effect of multiple interactions
per bunch crossing. Arrays of LGADs will be the sensors of choice for the
timing sub-detectors of both the CMS [2] and the ATLAS [3] experiments.

The intrinsic timing resolution of LGADs, in the order of a few tens of
picoseconds, is achieved thanks to a combination of 1) the sensors being
fabricated on thin substrates (20 - 50 $\mu$m), and 2) the signals being
internally amplified with gain in the order of 5 - 100, considerably smaller
than that featured in standard APDs. Such gain compensates for the small
signal in the thin substrate, makes the jitter contribution to noise small,
and allows for a large signal-to-noise ratio even when the sensor is coupled
to fast read-out electronics [4].

The gain has however been accurately measured only with mips, which create
about 80 electron-hole pairs per micron along a track in the traversed silicon
substrate. The mechanism by which other types of radiation produce charge
clouds differs from the linear track produced by mips, as the energy
deposition can be more or less deep and localized in silicon and produces a
distinct charge density. The gain experienced by different types of radiation
needs therefore to be assessed by dedicated measurements, for instance using a
low-noise read-out such as a charge sensitive pre-amplifier (CSA).

In this work, we study the gain of three different types of LGADs produced in
two different batches by Brookhaven National Laboratory (BNL) and Hamamatsu
Photonics (HPK) in a different batch. All LGADs were exposed to X/gamma-rays
of different energies, as well as beta (close to mips) and alpha particles.
Measurements were performed at different bias voltages to determine the gain
and the excess noise as a function of this parameter.

This paper is organized as follows: in Section 2, we describe the setup, the
LGAD devices under test, and the employed radioactive sources; in Section 3,
we present the spectroscopic performance of the LGAD extracted from the pulse
height spectra, discussing the gain of all sensors in the various scenarios
and the extrapolated shape of the distribution of charge deposited by
radiation.

##  2 Experimental setup

###  2.1 Sensors under test

A total of three LGAD devices were tested:

Figure 1: The three devices under test: top left) the BNL-old LGAD, center)
the BNL-new LGAD, and bottom left) the HPK LGAD. The schematic (right)
represents the readout system connected to the sensors as detailed in text
body.

  * $\bullet$

An LGAD sensor, produced as part of the very first BNL fabrication in 2017
[5], has a depletion voltage of 150 V (due to the relatively low resistivity
of the epitaxial layer), an active thickness of 50 $\mu$m, a breakdown voltage
of 300 V, and an active area of 2$\times$2 mm2. It will be referred to as
"BNL-old" in the following.

  * $\bullet$

A second LGAD comes from a more recent production by BNL carried out in 2021,
it has a total depletion voltage of 30 V, an active thickness of 50 $\mu$m, a
breakdown voltage of 200 V, and an active area of 1.3$\times$1.3 mm2 to match
the LGAD pixels designed for high-luminosity upgrades of the ATLAS and CMS
experiments. It will be referred to as "BNL-new" in the following.

  * $\bullet$

The third device under test is an LGAD fabricated by HPK, with a total
depletion voltage of 50 V, an active thickness of 35 $\mu$m, a breakdown
voltage of 260 V and an active area of 1.3$\times$1.3 mm2. It will be referred
to as "HPK" in the following.

The three devices are shown in Fig. 1. A major difference among the three
LGADs is the depth of the gain layer, obtained via boron implantation: in the
BNL-old sensor the gain layer is shallower than the one in BNL-new, with a
depth of $\sim$400 nm for the former and $\sim$650 nm for the latter, as shown
in Fig. 2; the gain layer of the HPK LGAD is deeper than the one implanted in
either BNL devices.

Figure 2: Doping profile of the boron gain layer as extracted from
measurements of capacitance as a function of voltage for the two BNL LGADs
under test. The position at $x=0$ represents the position of the junction
between the n+ layer and the gain layer. Only BNL-fabricated devices are
shown.

###  2.2 Readout system

The setup used to characterize the spectroscopic performance of the three
LGADs under test is shown in Fig. 1. The sensors are glued on a dual-in-line
ceramic board which plugs into a general-purpose breadboard. From here, short
wires connect the substrate and the guard ring to the the high voltage and
ground lines, respectively. A third wire connects the LGAD pad to the input of
the readout system.

The signal pad of the LGAD is AC-coupled to the input of an Amptek A250 low-
noise CSA
111https://www.amptek.com/-/media/ametekamptek/documents/resources/products/specs/a250-specs.pdf
through a large-value capacitance. The leakage current is discharged to ground
through a 1 MOhm resistor. The AC-coupling is rendered necessary due to the
sensor leakage current, whose magnitude becomes very high at higher bias
voltages (when the gain becomes large). Without such AC-coupling scheme, the
leakage current would flow through a 300 MOhm feedback resistor and saturate
the CSA.

The output of the A250 CSA is directly fed to an Amptek PX5 Digital Pulse
Processor
222https://www.amptek.com/-/media/ametekamptek/documents/resources/products/specs/px5-spec.pdf
which provides a digital trapezoidal shaping with adjustable peaking times and
flat-top widths. For all the measurements, the minimum allowed shaping time of
0.8 $\mu$s has been used, together with a 0.2 $\mu$s flat top width. The PX5
also provides the measurement of the pulse height of acquired signals by means
of an internal multi-channel analyzer (MCA). The measured pulse-height
spectrum is then plotted in real-time on a PC and stored for off-line
analysis. The calibration of the CSA for gain calculation is made by
performing the same measurements using a BNL-made 2 mm$\times$2 mm, 500 $\mu$m
thick silicon PIN diode, which does not feature any gain.

###  2.3 Radiation sources

Four spontaneous radioactive sources have been used to generate different
types of radiation: a) a 90Sr source, emitting beta particles with energies up
to 2 MeV that deposit energies close to those of mips; b) a 55Fe source of 5.9
keV X-rays; c) an 241Am source of gamma-rays up to 60 keV and d) a 241Am
source of alpha particles with energies up to 5.46 MeV.

These radioactive sources, placed about a centimeter away from the device
under test, illuminate the junction side of the devices. The entire setup was
enclosed in a metallic light-proof box for electromagnetic shielding.

##  3 Results

Signal spectra generated in LGADs by particles emitted by the four sources
were acquired for the three LGADs under study. Figure 3 shows the spectra
obtained for 90Sr beta particles and 241Am alpha particles acquired with the
BNL-old LGAD, while Figure 4 shows the spectra obtained obtained for 241Am
gamma-rays. Spectra are acquired for several bias voltages, from depletion
(150 V) to the onset of breakdown (290 V). Spectra from 55Fe are not shown for
the BNL-old device as the signals were too small compared to the high noise
and the limited amplification of the device. All the spectra but the alpha
particle one have been acquired by setting the PX5 digital gain to 5; for the
alpha particle spectra, as the generated charge is larger, a gain of 0.75 was
used.

Figure 3: Pulse height spectra acquired by the BNL-old LGAD when exposed to
top) betas from 90Sr and bottom) alphas from 241Am. Counts are normalized for
better comparison of shapes. The charge on the x-axis is represented in units
of ADC "channels". The distributions are cut at lower ADC values to remove the
contribution from the sensor noise. Figure 4: Pulse height spectra acquired by
the BNL-old LGAD when exposed to gamma-rays from 241Am. Counts are normalized
for better comparison of the shapes. The charge on the x-axis is represented
in units of ADC "channels". The distributions are cut at lower ADC values to
remove the contribution from the sensor noise.

In Fig. 5 and 6 the spectra from the same sources as in Fig. 3 and 4 are shown
for the BNL-new LGAD, fabricated with a deeper gain layer. Voltages range from
50 V (just greater than the depletion voltage) to 200 V (close to the
breakdown voltage). In the only case of acquisitions with alpha, an additional
spectrum has been acquired for $V_{bias}=30~{}V$ as noticeably the gain is
measured to be larger at this voltage. Better noise performance are clearly
noticeable by looking at the low energy peaks of the 241Am gamma-rays which,
differing from the case of the BNL-old LGAD, are well resolved for the lowest
voltages. Noticeably, the spectra from 55Fe are also visible in this device.

Figure 5: Pulse height spectra acquired by the BNL-new LGAD when exposed to,
from top to bottom, beta from 90Sr, alphas from 241Am, and X-rays from 55Fe.
Counts are normalized for better comparison. The charge on the x-axis is
represented in units of ADC "channels". The distributions are cut at lower ADC
values to remove the contribution from the sensor noise. Figure 6: Pulse
height spectra acquired by a BNL-new LGAD when exposed to gamma-rays from
241Am. Counts are normalized for better comparison. The charge on the x-axis
is represented in units of ADC "channels". The distributions are cut at lower
ADC values to remove the contribution from the sensor noise.

In Fig. 7 and 8, the spectra acquired with the HPK LGAD are reported for
voltages ranging from 50 V (about depletion voltage) to 260 V (close to
breakdown voltage). Spectra for the HPK device are very similar to the ones
obtained by the BNL-new sensor, except that larger gains are observed.

Figure 7: Pulse height spectra acquired by the HPK LGAD when exposed to (from
top to bottom) beta from 90Sr, alphas from 241Am, and X-rays from 55Fe. Counts
are normalized for better comparison. The charge on the x axis is represented
in units of ADC "channels". The distributions are cut at lower ADC values to
remove the contribution from the sensor noise. Figure 8: Pulse height spectra
acquired by the HPK LGAD when exposed to gamma-rays from 241Am. Counts are
normalized for better comparison. The charge on the x-axis is represented in
units of ADC "channels". The distributions are cut at lower ADC values to
remove the contribution from the sensor noise.

###  3.1 Gain

The summary of the average gain measured on the three LGADs under test is
plotted in Fig. 9; the gain as a function of the applied bias voltage is
plotted for the different radioactive sources used for the evaluation. Gain
values are obtained by comparing the signal produced by the LGAD sensors to
that produced by a PIN diode exposed to the same radioactive sources and
connected to the same readout system.

Figure 9: Summary of the average gain for the three devices under test as a
function of bias voltage, when the devices are exposed to different
radioactive sources.

For all devices under test, LGADs exhibit the largest gain when detecting
betas from 90Sr. The largest gain can be attributed to the lowest density in
the charge track of electron-hole pairs created by the impinging particles.
Charges are distributed along a track rather than in a localized spot into the
silicon bulk. The lowest gain (about 2 for all LGADs) is observed in the
interaction with alpha particles, as a large number of electron-hole pairs is
generated locally. For alphas and X- or gamma-rays, inner carriers in the
charge cloud are shielded from the external electric field by the outermost
charges, which furthermore lowers the local magnitude of the electric field.
As a result, the electrons experience a limited multiplication.

Overall, for mips (as for beta particles from 90Sr) the maximum measured gain
has a broad range and is about 120, 60, and 30 for the HPK, BNL-new, and BNL-
old devices, respectively. Such a large variation is attributed to the
different depth of the gain layer, with the HPK sensor having the deepest and
BNL-old the shallowest. When detecting X-rays, LGADs feature roughly half the
gain observed for mips while they are operated in the same conditions. Larger
gains are observed in the detection of low energy X-rays than in the detection
of higher energy gamma-rays, as a result of the lower carrier density of the
former.

Particularly interesting is the case of the gain for alphas. For HPK and BNL-
new LGADs, the gain is not a monotonic function of the bias voltage. At the
lowest voltages the LGADs show maximum gain, which then decreases when
increasing the applied voltage only to eventually slightly increase again at
higher voltage. This can be explained by the low electric field present in the
bulk at low voltage, which causes a slow electron drift towards the
multiplication region. This effect allows electrons to laterally diffuse and
therefore a lower charge density is presented in the gain region [6]. In this
case the self-screening effect is therefore slightly mitigated when compared
to faster drifts occurring at higher voltages. Finally, at the highest
voltages the higher electric fields are responsible for an increase in gain,
which compensates this effect.

###  3.2 Signal-to-Noise Ratio

The signal-to-noise ratio (S/N), defined as ratio between the distribution
mode and its full width at half maximum (FWHM) and calculated using the 55Fe
spectra acquired for the BNL-new and the HPK devices, is shown in Fig. 10. For
the extraction of this quantity, we evaluated only the spectra from the 55Fe
source as it creates in silicon only one narrow line from its $K_{\alpha}$
decay at 5.9 keV (neglecting the less active $K_{\beta}$ line at 6.4 keV): it
is therefore straightforward to measure the S/N. In the lowest bias voltage
range, the S/N slightly increases for increasing voltages until a maximum is
reached, then it rapidly decreases at higher voltages, as the gain increases.
This can be better understood with the help of Fig. 11 [7], showing the signal
and the noise components as a function of the gain for a generic avalanche
photodiode. It shows how a gain increase brings the signal over the noise
threshold in the detector. The series and 1/f electronic noise sources of the
system are independent of the gain, while the shot noise and the
multiplication [8] noises increase more than linearly with the gain. There
exists therefore an optimal gain that maximizes the S/N while higher gains are
in fact detrimental as far as S/N is concerned.

The spectra are not symmetric around the mode but rather they show a skewness
in the lower energy shoulder. This effect suggests that the multiplication
process is biased towards smaller gain values. The skewness for the same 55Fe
datasets is calculated as the Pearson’s first skewness coefficient:

$\rm{skewness=\frac{mean-mode}{sigma}}$ (1)

where sigma is the standard deviation of the pulse height spectra; the
skewness of the dataset as a function of the bias voltage is also shown in
Fig. 10. The skewness is observed to rapidly increase as the bias voltage
increases, albeit a minor divergence from this trend can be observed in the
BNL-new dataset for high bias voltages due to the high noise of the sensor. It
can also be noticed that even at the smaller voltages (close to depletion) the
skewness differs from zero.

Figure 10: Top, signal-to-noise ratio (S/N) for BNL-new and HPK LGADs and,
Bottom, the calculated skewness (in absolute value) as functions of the bias
voltage using the 55Fe source. Figure 11: Sketch of the signal (S) and the
several noise sources (N) in an APD [7] as a function of the gain.There is an
optimal gain that maximize the S/N.

##  4 Conclusions

By acquiring the pulse height spectra with a variety of LGADs exposed to
different radioactive sources and at different bias voltages, we were able to
characterize gain and noise properties. The gain of LGAD sensors depends on
the nature of the detected particle, being maximum for mips, which create a
low density of electron-hole pairs along a particle track, and minimum for
alpha particles that create a large number of electron-hole pairs in a
localized cloud. Moreover, the gain changes drastically as a function of the
depth of the gain layer, being higher in sensors fabricated with deeper gain
implants. The spectroscopic performance, quantified by the S/N ratio, rapidly
deteriorates by increasing the bias voltage and thus the gain, as the
multiplication and shot noises dominate over the other noise sources and
overcompensates for the signal gain. It appears therefore that there is an
optimal gain, usually small, that maximize spectroscopic performance. We also
observed how another spectroscopic parameter, the skewness, increases as a
function of bias voltage as a lower energy shoulder is featured in the
spectra, indicating a tendency of the multiplication process towards smaller
than average amplification values.

## References

  * [1] G. Pellegrini, P. Fernández-Martínez, M. Baselga, C. Fleta, D. Flores, V. Greco, S. Hidalgo, I. Mandić, G. Kramberger, D. Quirion, M. Ullan, Technology developments and first measurements of Low Gain Avalanche Detectors (LGAD) for high energy physics applications, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 765 (2014) 12–16, hSTD-9 2013 - Proceedings of the 9th International "Hiroshima" Symposium on Development and Application of Semiconductor Tracking Detectors.  doi:https://doi.org/10.1016/j.nima.2014.06.008.
URL https://www.sciencedirect.com/science/article/pii/S0168900214007128

  * [2] CMS Collaboration, A MIP Timing Detector for the CMS Phase-2 Upgrade, Tech. rep., CERN, Geneva (2019).
URL https://cds.cern.ch/record/2667167

  * [3] Technical Design Report: A High-Granularity Timing Detector for the ATLAS Phase-II Upgrade, Tech. rep., CERN, Geneva (2020).
URL https://cds.cern.ch/record/2719855

  * [4] R. Heller, A. Abreu, A. Apresyan, R. Arcidiacono, N. Cartiglia, K. DiPetrillo, M. Ferrero, M. Hussain, M. Lazarovitz, H. Lee, S. Los, C. Moon, C. Peña, F. Siviero, V. Sola, T. Wamorkar, S. Xie, Combined analysis of hpk 3.1 lgads using a proton beam, beta source, and probe station towards establishing high volume quality control, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1018 (2021) 165828.  doi:https://doi.org/10.1016/j.nima.2021.165828.
URL https://www.sciencedirect.com/science/article/pii/S0168900221008135

  * [5] G. Giacomini, W. Chen, F. Lanni, A. Tricoli, Development of a technology for the fabrication of low-gain avalanche diodes at bnl, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 934 (2019) 52–57.  doi:https://doi.org/10.1016/j.nima.2019.04.073.
URL https://www.sciencedirect.com/science/article/pii/S0168900219305479

  * [6] M. Jaksic, A. Crnjac, G. Kramberger, M. Manojlovic, G. Laštovička-Medin, M. R. Ramos, Ion microbeam studies of charge transport in semiconductor radiation detectors with three-dimensional structures: An example of lgad, Front. Phys., 13 May 2022, Sec. Radiation Detectors and Imaging Volume 10 (2022).  doi:https://doi.org/10.3389/fphy.2022.877577.
  * [7] N. Cartiglia, Why shot noise in lgad does not degrade time resolution?, report number: AIDA-2020-SLIDE-2019-015 (2019).
URL https://cds.cern.ch/record/2667247

  * [8] R. McIntyre, Multiplication Noise in Uniform Avalanche Diodes, IEEE Transaction on Electron Devices (1966) 164–168doi:10.1109/T-ED.1966.15651.


pnasresearcharticle

Few particles in a quantum many-body system do not generally reach the realm
of the macroscopic world. We show that quantum jamming is a mechanism to make
microscopic phenomena visible at a macroscopic scale: a single particle
produces a permanent change in the state, and few particles leave an imprint
of their scattering properties and bound states. One of the most noteworthy
phenomena is that a single observation of a local observable can turn a jammed
separable state into a state in which quantum mechanics is necessary to
describe macroscopic properties.

<EMAIL_ADDRESS>
# Quantum jamming brings quantum mechanics to macroscopic scales

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

###### Abstract

A quantum spin-$\frac{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, realising the quantum analogue 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,
eigenstate localisation 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.

7]3Low 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,
topological order (3) in 2D. Since nowadays such systems have become
experimentally accessible (4, 5)—turning the concept of quantum simulator (6,
7, 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 settings (9, 10, 11, 12). On the one hand, they
can exhibit exotic behaviours such as relaxation to non-thermal states (13,
14, 15); on the other hand, there are interacting 1D models that allow for
non-perturbative 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, 21, 22, 23, 24, 25, 26, 27, 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 celebrated glass
transition, which could be approached by varying a thermodynamic quantity such
as the density.

The key question we address here is which phenomenology 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 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-$\frac{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^{s,-s}_{\ell,\ell+y}=\sigma_{\ell}^{s}\prod_{j=1}^{y-1}\frac{1-\sigma_{\ell+j}^{z}}{2}\sigma_{\ell+y}^{-s}\,,\qquad\sigma^{z}_{\ell}\,,$
(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 (44, 43), “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. Ref. (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}^{x}+\sigma_{\ell}^{y}\sigma_{\ell+1}^{y}$).
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.

We reconsider this nonequilibrium problem from a more general perspective and
show that the phenomenology characterised by macroscopic effects from local
perturbations, which we will compactly call MELP, survives both integrable and
non-integrable 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, setting
quantum jamming apart from its classical counterpart. The non-symmetric case
requires a separate treatment and the phenomenology turns out to strongly
depend on the inhomogeneity of the jammed state. Specifically, MELP can be
observed even in the absence of that symmetry whenever the jammed state is
(almost) homogeneous; the more disordered the sequence of species, however,
the more hindered the spreading of the perturbation.

## 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 will put 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-$\frac{1}{2}$ Hamiltonian with an axial symmetry.

### 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 distinguishing $y$ species of particles,
characterised by their position modulo $y$. We can then separate the
information about the species from the position $\ell$ by defining the
_macroposition_ $\ell^{\prime}$ as

$\ell=y\ell^{\prime}-b\,,$ (2)

where $b\in\\{0,1,\dots,y-1\\}$ labels the particle species.

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,
$\frac{1}{y}\sum_{n=0}^{y-1}T_{\ell-n,\ell-n+y}$ and
$\frac{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.

### Multispecies particle-hole duality transformation

Since the sequence of species $\underline{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 $\underline{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
between the corresponding particles minus one— see Fig. 1.

$\ket{\underbrace{\cdots\downarrow}_{jammed}\underbrace{\downarrow\downarrow}_{\ell^{\prime}=-2}\underbrace{\uparrow\downarrow}_{-1}\underbrace{\downarrow\uparrow_{\ell=0}}_{0}\underbrace{\downarrow\uparrow}_{1}\underbrace{\downarrow\downarrow}_{2}\underbrace{\uparrow\uparrow}_{3}\underbrace{\downarrow\cdots}_{jammed}}\sim\ket{\Downarrow\uparrow\downarrow^{(1)}\uparrow\downarrow^{(0)}\downarrow_{j=0}^{(0)}\uparrow\downarrow^{(1)}\downarrow^{(0)}\Downarrow}\equiv\ket{-4,-2,1;\dots,1,0,0,1_{\tilde{j}=0},0,\dots}$
Figure 1: Example of the duality mapping from spins to pseudospins for $y=2$.
The product state is almost jammed with three impurities. Under the braces on
the left hand side we displayed the macropositions $\ell^{\prime}$ of the
spins (position $\ell=0$ is highlighted as well). At the centre we showed the
corresponding sequence of pseudospins ($\Downarrow$ stands for a semi-infinite
string of pseudospins down), displaying the particle species of the
pseudospins down as superscripts. We have also highlighted pseudoposition
$j=0$. On the right hand side we represented the state by indicating the
pseudopositions of the pseudospins up followed by the sequence of species, as
used in the paper.

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 described 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
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,\dots,\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^{\prime}(\tilde{j})=\tilde{j}+\delta\tilde{j}+\sum_{k=1}^{\nu}\theta(n_{k}+1-k<\tilde{j})$
(3)

where $n_{1},\dots,n_{\nu}$ are the pseudopositions of the impurities and
$\delta\tilde{j}$ parametrises 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}\geq
0\\\
\theta_{H}(b_{-\tilde{n}}-b_{-\tilde{n}-1}<0)&\text{otherwise.}\end{cases}$
(4)

In the following we will denote the state with the sequence of species
$\underline{b}$ and impurities at $n_{1},\dots,n_{\nu}$ by
$\ket{n_{1},\dots,n_{\nu};\underline{b}}$ (the absence of impurities will be
denoted by $\ket{\emptyset;\underline{b}}$)—see Fig. 1.

For given $\underline{b}$, the dynamics can be projected onto the space of
pseudospins by applying the duality transformation to the Hamiltonian. A
procedure is described in ‘Material and Methods’. 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. Non-symmetric operators,
instead, have a sequence-dependent semilocal representation111We 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. [17]..

## Results

### Local magnetisation

We introduce the notation

$:O:=O-\braket{\emptyset;\underline{b}}{O}{\emptyset;\underline{b}}\,.$ (5)

In the basis $\ket{n_{1},\dots,n_{\nu};\underline{b}}$, the local
magnetisation is a diagonal operator. Using Eq. [3], we find the matrix
elements shown in Eq. [6] and Eq. [7]—see SI for more details.

$\braket{n;\underline{b}}{:\sigma_{y\ell^{\prime}-j}^{z}:}{n;\underline{b}}=\braket{\emptyset;\underline{b}}{\sigma_{y(\ell^{\prime}-1)-j}^{z}}{\emptyset;\underline{b}}\theta_{H}(n+\delta
n+1<\ell^{\prime})-\braket{\emptyset;\underline{b}}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset;\underline{b}}\theta_{H}(n+\delta
n-\theta_{H}(b_{n+1}<b_{n})<\ell^{\prime})\\\ +\delta_{\ell^{\prime},n+\delta
n+1}[\theta_{H}(b_{n+1}<b_{n})\delta_{b_{n+1},j}-1]+\delta_{\ell^{\prime},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})]$
(6)

$\braket{n_{1},n_{2},\dots,n_{\nu};\underline{b}}{:\sigma^{z}_{\ell}:}{n_{1},n_{2},\dots,n_{\nu};\underline{b}}=\sum_{k=1}^{\nu}\braket{n_{k}+1-k;\underline{b}}{:\sigma_{\ell-y(k-1)}^{z}:}{n_{k}+1-k;\underline{b}}$
(7)

Roughly speaking, the local magnetisation in a state with a single impurity is
identified with either one local magnetisation 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 have a similar
interpretation and can be reduced to those with a single impurity, Eq. [7].
Time evolution is readily obtained

$\braket{:\sigma^{z}_{\ell}:}_{t}=\sum_{k=1}^{\nu}\sum_{n}p_{n+k-1,t}^{(k)}\braket{n;\underline{b}}{:\sigma_{\ell-y(k-1)}^{z}:}{n;\underline{b}}$
(8)

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

$p^{(j)}_{n,t}=\sum_{n^{\prime}}|\braket{n_{1}^{\prime},\dots,n^{\prime}_{j-1},n,n^{\prime}_{j+1},\dots,n^{\prime}_{\nu}}{e^{-i\tilde{H}t}}{\tilde{\Psi}_{0}}|^{2}.$
(9)

Here $\ket{\tilde{\Psi}_{0}}$ 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^{\prime}_{1},\dots,n^{\prime}_{j-1},n^{\prime}_{j+1},\dots,n^{\prime}_{\nu}$
satisfying $n_{1}^{\prime}<\dots<n^{\prime}_{j-1}<n$ and
$n<n^{\prime}_{j+1}<\dots<n^{\prime}_{\upsilon}$. For finite $\nu$, the
qualitative behaviour of the magnetisation can be readily understood if
$\tilde{H}$ is local (i.e., the original Hamiltonian is symmetric) and doesn’t
exhibit exceptional localisation properties in the eigenstates. It is then
reasonable to expect $p^{(j)}_{n,t}$ to approach zero in time. Consequently,
the terms in the second line of Eq. [6] give a subleading contribution to Eq.
[8], and the argument of the step functions on the first line of Eq. [6] can
be shifted without affecting the asymptotic behaviour of the local
magnetisation. We then end up with the following simplified expression

$\braket{:\sigma^{z}_{\ell}:}_{t}\sim\sum_{k=0}^{\nu-1}\braket{\emptyset;\underline{b}}{\sigma_{\ell-y(k+1)}^{z}-\sigma_{\ell-
yk}^{z}}{\emptyset;\underline{b}}F_{t}^{(1+k)}(x_{\ell})$ (10)

where $F_{t}^{(k)}(x)=\sum_{n\lesssim x}p_{n,t}^{(k)}$ is the cumulative
distribution and $x_{\ell}$ is the effective pseudoposition associated with
site $\ell$, i.e., $\frac{2}{\ell}(x_{\ell}+\delta x_{\ell})\sim
1+O(\ell^{-1})$. For its generality and transparency, Eq. [10] is the most
important formula of this paper. Firstly, by modifying the underlying jammed
state, the coefficients of $F_{t}^{(k)}(x_{\ell})$ in Eq. [10] (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. Secondly, 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 magnetisation as a function of $x/t$, which
is the unusual lightcone behaviour pointed out in Refs (46, 47) in a specific
model with a single impurity. Corrections should be expected if there is a
finite number of localised eigenstates, whose more explicit effect would be to
create effective discontinuities in the aforementioned relevant scaling limit
of $F_{t}^{(k)}(x)$.

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

### Example: symmetric Hamiltonian

Ref. (47) has proven MELP in the special case in which the Hamiltonian is
local and noninteracting in the pseudospace. In order to prove the robustness
of the phenomenology under interacting and integrability breaking terms, we
consider the following Hamiltonian with $y=2$

$H=J\sum\nolimits_{\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}$
(11)

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
integrability222As 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 $\underline{b}$, the exact symmetry requiring also to add an
irrelevant constant and a coupling to $S^{z}$ (omitted here for the sake of
compactness).. The Hamiltonian of Ref. (47) is recovered by setting
$\Delta=g=0$. The projection of $H$ onto the pseudospace reads

$\tilde{H}\sim
J\sum\nolimits_{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^{x}_{j-1}\tau^{x}_{j+1}+\tau^{y}_{j-1}\tau^{y}_{j+1})\frac{1+\tau^{z}_{j}}{2}\,,$
(12)

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

#### 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 doesn’t have any
impact on the local magnetisation.

#### Two impurities: integrable representation

Figure 2: The cumulative distributions for the (symmetric) example in Eq. [11]
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 trivialising
the scattering phase. The region delimited by each curve represents the
contribution from the bound states.

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
(48, 2) 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})}}.$
(13)

Importantly, the corresponding integrable model exhibits also a bound state in
the 2-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$ travelling at a typical pseudodistance
$1/\log|\frac{\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,\mathrm{m}}^{(k)}(x)+F_{t,\mathrm{b}}^{(k)}(x)$, which
are reported in Eq. [14] and Eq. [15].

$\displaystyle
F_{t,\mathrm{m}}^{(j)}(x)\sim\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p}{8\pi^{2}}\theta_{H}(\frac{x}{t}-v^{(j)}(p_{1},p_{2}))|e^{i(n_{1}^{(0)}-n_{2}^{(0)})(p_{1}-p_{2})}+S(p_{1},p_{2})|^{2}\qquad
v^{(j)}(p_{1},p_{2})=\begin{cases}\min(E^{\prime}(p_{1}),E^{\prime}(p_{2}))&j=1\\\
\max(E^{\prime}(p_{1}),E^{\prime}(p_{2}))&j=2\end{cases}$ (14) $\displaystyle
F_{t,\mathrm{b}}^{(1)}(x)\sim
F_{t,\mathrm{b}}^{(2)}(x)\sim\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\theta_{H}(\frac{x}{t}-E_{2}^{\prime}(p))\max\Bigl{[}\bigl{(}\tfrac{\Delta+g\cos
p}{\cos\frac{p}{2}}\big{)}^{2}-1,0\Bigr{]}\bigl{(}\tfrac{\cos^{2}\frac{p}{2}}{(\Delta+g\cos
p)^{2}}\bigr{)}^{n_{2}^{(0)}-n_{1}^{(0)}}$ (15)

There, $n_{1}^{(0)}<n_{2}^{(0)}$ denote the pseudopositions of the impurities
at the initial time (note that we took the limit $x,t\rightarrow\infty$ for
given $n_{1}^{(0)},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 by trivialising
the scattering phase, i.e., $S(p_{1},p_{2})\rightarrow-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 lightcone.

### Example: non-symmetric Hamiltonian

The symmetry of the Hamiltonian shown in Eq. [11] is broken by the interaction

$H_{I}=JV\sum\nolimits_{\ell}\sigma_{\ell}^{z}\sigma_{\ell+1}^{z}\,,$ (16)

which has the following semilocal representation in the pseudospace

$\tilde{H}_{I}=4JV\sum\nolimits_{j}\sum\nolimits_{\tilde{j}}\frac{1-(-1)^{b_{\tilde{j}}-b_{\tilde{j}+1}}}{2}\int_{-\pi}^{\pi}\frac{\mathrm{d}\omega}{2\pi}\\\
e^{i\omega(\tilde{j}-j-1+\sum_{j^{\prime}=-\infty}^{j}\frac{1+\tau^{z}_{j}}{2})}\frac{1-\tau_{j}^{z}}{2}\frac{1-\tau_{j+1}^{z}}{2}\,.$
(17)

Here the integral over $\omega$ represents a Kronecker delta. In the following
we will 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\nolimits_{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}\,,$
(18)

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

#### Quasi-homogeneous jammed state

The semilocality of the interaction is masked if the sequence of particles
consists of a repeating pattern. 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\\}$ 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
$\ket{\cdots\uparrow\uparrow\downarrow\uparrow\uparrow\downarrow_{0}\downarrow\uparrow\downarrow\uparrow\uparrow\downarrow\cdots}\sim\ket{0;\underline{b}}$,where
the left (right) hand side shows the (pseudo)spins and (pseudo)position $0$
($\Downarrow$ stands for a semi-infinite string of pseudospins down). The
sequence of species reads

$\underline{b}=\dots,(1,0,0,1),(1,0,0,1_{0}),(0,0,1,1),(0,0,1,1),\dots$ (19)

where we displayed relative position $0$ and grouped 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. [18], 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
localised excited states. This is confirmed by our numerical investigation,
which shows that the initial state has a significant overlap with a localised
excited state. Nevertheless, there is a delocalised contribution that triggers
MELP. This is shown in the first row of Figure 3, where we can also see that
the localised excited state creates a discontinuity in the cumulative
distribution function (the sharp change of colour 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 MELP. We surmise however that the effects could
become less pronounced due to the localised inhomogeneities of the
interactions.

#### Disordered jammed state: localisation

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. [18] is disordered, which
is known to trigger Anderson’s localisation (49). As shown in the second row
of Figure 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 localised perturbation on the opposite side,
provided that the width of the barrier is large enough (with respect to
$1/V^{2}$).

Figure 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
(non-symmetric) example described by Eq. [18] starting from
$\ket{0;\underline{b}}$ with the sequence in Eq. [19]. 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 raw is for a configuration with a small disordered region delimited by
the two orange vertical lines.

### Macroscopic quantum states

The problem under investigation shares some similarities with the MELP studied
in Ref. (50), which considered local perturbations in low-entangled excited
states that cannot be interpreted as ground states of local conservation laws.
Following Ref. (50), let us then imagine to prepare the system in a jammed
product state $\ket{\Psi(0^{-})}$ (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\vec{\sigma}_{\ell_{0}}$,
with $\hat{n}$ different from $\hat{z}$. The state after the measurement
$\ket{\Psi(0^{+})}$ will be in a superposition of the original jammed state
and an almost identical quasi-jammed state with an impurity,
$\ket{\Psi(0^{+})}=\sqrt{p}\ket{\emptyset;\underline{b}}+e^{i\phi}\sqrt{1-p}\ket{n_{0};\underline{b}^{(-)}}$
($\phi$ and $p$ depend on $\hat{n}$), with probability $p$. By the Lieb-
Robinson bound (51), the speed at which the information about the impurity
spreads is bounded from above, therefore the state remains separable outside a
lightcone 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_{\vec{s}}^{z}=\frac{1}{2}\sum_{\ell}s_{\ell}\sigma_{\ell}^{z}$, with
$s_{\ell}=\pm 1$, which has the property not to fluctuate outside the
lightcone. Because of that, its variance matches the variance of its
restriction to $\Omega_{t}$. Denoting its expectation value and its variance
at time $t$ by $\braket{S^{z}_{\vec{s}}}_{t}$ and
$\braket{S^{z}_{\vec{s}},S^{z}_{\vec{s}}}_{t}$, respectively, we have

$\braket{S^{z}_{\vec{s}},S^{z}_{\vec{s}}}_{t}\sim\frac{p}{1-p}(\braket{S^{z}_{\vec{s}}}_{t}-\braket{S^{z}_{\vec{s}}}_{0^{-}})^{2}\\\
+(1-p)\braket{n_{0};\underline{b}^{(-)}}{e^{iHt}S^{z}_{\vec{s}},S^{z}_{\vec{s}}e^{-iHt}}{n_{0};\underline{b}^{(-)}}\,.$
(20)

Assuming Eq. [10], 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 find that even the fluctuations on the second line grow as $t^{2}$,
e.g., if
$s_{\ell}=\braket{\emptyset;\underline{b}}{\sigma_{\ell}^{z}}{\emptyset;\underline{b}}$—see
SI—hence the same conclusion applies also to the marginal case with $p=0$,
which is the basic problem studied before. Following Ref. (52), the square
root of the variance provides a lower bound to the “minimal system size for
which one has to assume validity of quantum mechanics in order to explain a
non-classical phenomenon”. Since such a lower bound is proportional to
$|\Omega_{t}|$, 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.

## Discussion

Macroscopic effects from local perturbations have been recently observed in
several settings: ground states in ordered phases with spontaneous symmetry
breaking (53, 54), symmetric initial states in the presence of semilocal
conservation laws (55, 56), separable quantum scars (50), 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 MELP 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.

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 quasi-jammed 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 to appear in the sector with three
impurities.

In the integrable case it is instead still unclear the relationship between
the lightcone 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 generalised hydrodynamics.

### Explicit construction of the Hamiltonian in the example

The Hamiltonians in the examples have both a linear and a quadratic dependency
on the building blocks in Eq. [1]. Specifically, up to a constant and a term
proportional to $S^{z}$, which is conserved and irrelevant, the Hamiltonian in
Eq. [11] reads

$H\sim
2J\sum\nolimits_{\ell}\bar{T}^{+-}_{\ell-1,\ell+1}+\bar{T}^{-+}_{\ell-1,\ell+1}+2\Delta\\{\bar{T}^{+-}_{\ell-1,\ell+2},\bar{T}^{-+}_{\ell-1,\ell+2}\\}\\\
+g(T_{\ell,\ell+2}^{\mp}T_{\ell-2,\ell}^{\mp}+T_{\ell-2,\ell}^{\pm}T_{\ell,\ell+2}^{\pm})$

where we expressed 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.

Concerning the non-symmetric example, the interaction in Eq. [16] was already
expressed in terms of the building blocks $\sigma_{\ell}^{z}$.

### Duality mapping of shift invariant operators

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}\frac{1-\sigma_{\ell+j}}{2}\sigma_{\ell+m}^{z}$
for $m\geq 1$ can by inferred by applying them to a sequence of spins starting
with a spin up and ending before the next spin up,
$\ket{\dots\uparrow_{j}\underbrace{\downarrow\dots\downarrow}_{n}(\uparrow_{j+n+1})\dots}$.
Such a fragment is indeed mapped into the elementary fragment
$\ket{\dots\downarrow^{b}\underbrace{\uparrow\dots\uparrow}_{\lfloor
n/y\rfloor}(\downarrow^{b-n-1\\!\\!\\!\pmod{y}})\dots}$, where arrows
represent here pseudospins. It is then easy to identify the pseudospin
operator that acts on the dual fragment as the original spin operator does on
the original fragment (see SI for explicit expressions).

### Bethe Ansatz with two impurities

In order to solve the two-impurity sector we have used the Bethe Ansatz

$\ket{p_{1},p_{2}}\propto\\!\\!\\!\\!\sum_{n_{1}<n_{2}}[e^{in_{1}p_{1}+in_{2}p_{2}}+S(p_{1},p_{2})e^{in_{1}p_{2}+in_{2}p_{1}}]\ket{n_{1},n_{2};\underline{b}}\,.$

In the symmetric example we obtain the scattering phase in Eq. [13] and energy
$E(p)=4J(\cos p-\Delta)$. By applying an artificial finite-volume
regularisation 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 2-particle sector. The corresponding excited states can be
parametrised as follows

$\ket{p}_{\mathrm{b}}=\sum_{n_{1}<n_{2}}e^{i\frac{n_{1}+n_{2}}{2}p}\Bigl{(}\frac{\cos\frac{p}{2}}{\Delta+g\cos
p}\Bigr{)}^{n_{2}-n_{1}}\ket{n_{1},n_{2};\underline{b}}$

provided that $1+\cos p<2(\Delta+g\cos p)^{2}$. Their energy reads
$E_{2}(p)=4J(\cos^{2}\frac{p}{2}+g^{2}\cos^{2}p-\Delta^{2})/(\Delta+g\cos p)$.
The reader can find more details in SI, which also includes the derivation of
Eq. [14] and Eq. [15]

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

## References

  * (1) T Giamarchi, Quantum Physics in One Dimension.  (Oxford University Press), (2003).
  * (2) FHL Essler, H Frahm, F Göhmann, A Klümper, VE Korepin, The One-Dimensional Hubbard Model.  (Cambridge University Press), (2005).
  * (3) B Zeng, X Chen, DL Zhou, XG Wen, Quantum Information Meets Quantum Matter.  (Oxford University Press), (2019).
  * (4) I Bloch, J Dalibard, W Zwerger, Many-body physics with ultracold gases.  Rev. Mod. Phys. 80, 885–964 (2008).
  * (5) MA Cazalilla, R Citro, T Giamarchi, E Orignac, M Rigol, One dimensional bosons: From condensed matter systems to ultracold gases.  Rev. Mod. Phys. 83, 1405–1466 (2011).
  * (6) I Buluta, F Nori, Quantum Simulators.  Science 326, 108–111 (2009).
  * (7) JI Cirac, P Zoller, Goals and opportunities in quantum simulation.  Nature Physics 8, 264–266 (2012).
  * (8) H Bernien, et al., Probing many-body dynamics on a 51-atom quantum simulator.  Nature 551, 579–584 (2017).
  * (9) A Polkovnikov, K Sengupta, A Silva, M Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems.  Rev. Mod. Phys. 83, 863–883 (2011).
  * (10) C Gogolin, J Eisert, Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems.  Reports on Progress in Physics 79, 056001 (2016).
  * (11) M Gring, et al., Relaxation and Prethermalization in an Isolated Quantum System.  Science 337, 1318–1322 (2012).
  * (12) M Schemmer, I Bouchoule, B Doyon, J Dubail, Generalized Hydrodynamics on an Atom Chip.  Phys. Rev. Lett. 122, 090601 (2019).
  * (13) M Rigol, V Dunjko, V Yurovsky, 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) FHL Essler, M Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains.  Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
  * (15) E Ilievski, M Medenjak, T Prosen, L Zadnik, Quasilocal charges in integrable lattice systems.  Journal of Statistical Mechanics: Theory and Experiment 2016, 064008 (2016).
  * (16) E Ilievski, et al., Complete Generalized Gibbs Ensembles in an Interacting Theory.  Phys. Rev. Lett. 115, 157201 (2015).
  * (17) OA Castro-Alvaredo, B Doyon, 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, M Fagotti, Transport in Out-of-Equilibrium $XXZ$ Chains: Exact Profiles of Charges and Currents.  Phys. Rev. Lett. 117, 207201 (2016).
  * (19) SJ Ran, et al., Tensor Network Contractions : Methods and Applications to Quantum Many-Body Systems, Lecture Notes in Physics, 964.  (Springer Nature, Cham), 1st ed. 2020. edition, (2020).
  * (20) M Moeckel, 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, FHL Essler, S Groha, NJ Robinson, Prethermalization and Thermalization in Models with Weak Integrability Breaking.  Phys. Rev. Lett. 115, 180601 (2015).
  * (23) M Babadi, E Demler, 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, J Schmiedmayer, Prethermalization and universal dynamics in near-integrable quantum systems.  Journal of Statistical Mechanics: Theory and Experiment 2016, 064009 (2016).
  * (25) D Abanin, W De Roeck, WW Ho, F Huveneers, A Rigorous Theory of Many-Body Prethermalization for Periodically Driven and Closed Quantum Systems.  Communications in Mathematical Physics 354, 809–827 (2017).
  * (26) V Alba, M Fagotti, Prethermalization at Low Temperature: The Scent of Long-Range Order.  Phys. Rev. Lett. 119, 010601 (2017).
  * (27) P Reimann, L Dabelow, Typicality of Prethermalization.  Phys. Rev. Lett. 122, 080603 (2019).
  * (28) J Durnin, MJ Bhaseen, B Doyon, Nonequilibrium Dynamics and Weakly Broken Integrability.  Phys. Rev. Lett. 127, 130601 (2021).
  * (29) M Fagotti, On conservation laws, relaxation and pre-relaxation after a quantum quench.  Journal of Statistical Mechanics: Theory and Experiment 2014, P03016 (2014).
  * (30) B Bertini, M Fagotti, Pre-relaxation in weakly interacting models.  Journal of Statistical Mechanics: Theory and Experiment 2015, P07012 (2015).
  * (31) L Zadnik, JP Garrahan, Slow heterogeneous relaxation due to constraints in dual XXZ models.  arXiv:2304.10394 (2023).
  * (32) AH MacDonald, SM Girvin, D Yoshioka, $\frac{t}{U}$ expansion for the Hubbard model.  Phys. Rev. B 37, 9753–9756 (1988).
  * (33) L Zadnik, 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, OI Motrunich, Hilbert Space Fragmentation and Commutant Algebras.  Phys. Rev. X 12, 011050 (2022).
  * (35) G Biroli, C Chamon, F Zamponi, Theory of the superglass phase.  Phys. Rev. B 78, 224306 (2008).
  * (36) Z Nussinov, P Johnson, MJ Graf, AV 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, 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?  Nature Physics 3, 222–223 (2007).
  * (39) C Weitenberg, et al., Single-spin addressing in an atomic Mott insulator.  Nature 471, 319–324 (2011).
  * (40) T Došlić, M Puljiz, S Šebek, J Žubrinić, Complexity Function of Jammed Configurations of Rydberg Atoms.  arXiv:2302.08791 (2023).
  * (41) L Gotta, L Mazza, P Simon, 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, E Ilievski, Universal anomalous fluctuations in charged single-file systems.  arXiv:2208.01463 (2023).
  * (43) M Borsi, L Pristyák, B Pozsgay, Matrix product symmetries and breakdown of thermalization from hard rod deformations (2023) arXiv:2302.07219.
  * (44) B Pozsgay, T Gombor, A Hutsalyuk, Integrable hard-rod deformation of the Heisenberg spin chains.  Phys. Rev. E 104, 064124 (2021).
  * (45) RZ Bariev, Integrable spin chain with two- and three-particle interactions.  Journal of Physics A: Mathematical and General 24, L549 (1991).
  * (46) K Bidzhiev, M Fagotti, 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, M Fagotti, Measurement catastrophe and ballistic spread of charge density with vanishing current.  Journal of Physics A: Mathematical and Theoretical 55, 474001 (2022).
  * (48) H Bethe, Zur Theorie der Metalle.  Zeitschrift für Physik 71, 205–226 (1931).
  * (49) PW Anderson, Absence of Diffusion in Certain Random Lattices.  Phys. Rev. 109, 1492–1505 (1958).
  * (50) S Bocini, M Fagotti, Growing Schrödinger’s cat states by local unitary time evolution of product states (2023) arXiv.2210.15585.
  * (51) EH Lieb, DW Robinson, The finite group velocity of quantum spin systems.  Communications in Mathematical Physics 28, 251 – 257 (1972).
  * (52) F Fröwis, W Dür, Measures of macroscopicity for quantum spin systems.  New Journal of Physics 14, 093039 (2012).
  * (53) V Zauner, M Ganahl, HG Evertz, T Nishino, Time evolution within a comoving window: scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains.  Journal of Physics: Condensed Matter 27, 425602 (2015).
  * (54) V Eisler, F Maislinger, Front dynamics in the XY chain after local excitations.  SciPost Phys. 8, 037 (2020).
  * (55) M Fagotti, Global Quenches after Localized Perturbations.  Phys. Rev. Lett. 128, 110602 (2022).
  * (56) M Fagotti, V Marić, L Zadnik, Nonequilibrium symmetry-protected topological order: emergence of semilocal Gibbs ensembles (2022) arXiv:2205.02221.

## Supplemental Information

In the next four sections we have collected some details that might help the
reader grasp the results reported in the main text.

## Expectation values of the local magnetisation

Equation 3 allows us to infer a representation of local operators in the
pseudospace. In particular we have

$\sigma_{\ell}^{z}=\sigma_{2\lceil\frac{\ell}{y}\rceil-[(-\ell)\\!\\!\\!\pmod{y}]}^{z}\rightarrow-1+2\sum_{\tilde{j}}\delta_{b_{\tilde{j}},(-\ell)\\!\\!\\!\pmod{y}}\int\frac{\mathrm{d}\tilde{\omega}}{2\pi}e^{i\tilde{\omega}(\ell^{\prime}(\tilde{j})-\lceil\frac{\ell}{y}\rceil)}$
(21)

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

#### Zero impurities: jammed states

In a jammed state (which is the vacuum $\ket{\emptyset}$ in the pseudospace)
we have

$\ell^{\prime}(\tilde{j})=\tilde{j}+\delta j\,.$ (22)

We then have

$\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}=-1+2\sum_{\tilde{j}}\delta_{\ell^{\prime},\tilde{j}+\delta\tilde{j}}\delta_{b_{\tilde{j}},j}\,.$
(23)

#### One impurity

We denote by $\ket{n}$ the state with an impurity at position $n$ in the
pseuodspace. We have

$\ell^{\prime}(\tilde{j})=\tilde{j}+\delta\tilde{j}+\theta_{H}(\tilde{j}>n)$
(24)

and hence

$\braket{n}{\sigma_{y\ell^{\prime}-j}^{z}}{n}=-1+2\sum_{\tilde{j}\in
B}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j}+\theta_{H}(\tilde{j}>n),\ell^{\prime}}$
(25)

We are going to express these matrix elements in terms of those in the
underlying jammed state using Eq. [23]. We start from the identities

$\displaystyle\sum_{\tilde{j}>n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j}+1,\ell^{\prime}}=\frac{1+\braket{\emptyset}{\sigma_{y(\ell^{\prime}-1)-j}^{z}}{\emptyset}}{2}-\sum_{\tilde{j}\leq
n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j}+1,\ell^{\prime}}$
(26) $\displaystyle\sum_{\tilde{j}\leq
n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}}=\frac{1+\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}}{2}-\sum_{\tilde{j}>n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}}$
(27)

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

$\delta(n+1)-\delta n=-\theta_{H}(b_{n+1}<b_{n})$

we finally get

$\displaystyle\sum_{\tilde{j}>n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j}+1,\ell^{\prime}}=\begin{cases}\frac{1+\braket{\emptyset}{\sigma_{y(\ell^{\prime}-1)-j}^{z}}{\emptyset}}{2}&n+\delta
n+1<\ell^{\prime}\\\ 0&n+\delta
n+1-\theta_{H}(b_{n+1}<b_{n})\geq\ell^{\prime}\\\
\delta_{b_{n+1},j}&b_{n+1}<b_{n},\ell^{\prime}=n+\delta n+1\end{cases}$ (28)
$\displaystyle\sum_{\tilde{j}\leq
n}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}}=\begin{cases}\frac{1+\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}}{2}&n+\delta
n-\theta_{H}(b_{n+1}<b_{n})\geq\ell^{\prime}\\\ 0&n+\delta n<\ell^{\prime}\\\
\theta_{H}(b_{n}<b_{n-1})\delta_{b_{n-1},j}+\delta_{b_{n},j}&b_{n+1}<b_{n},\ell^{\prime}=n+\delta
n\end{cases}$ (29)

and hence

$\braket{n}{\sigma_{y\ell^{\prime}-j}^{z}}{n}-\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}=\\\
\braket{\emptyset}{\sigma_{y(\ell^{\prime}-1)-j}^{z}}{\emptyset}\theta_{H}(n+\delta
n+1<\ell^{\prime})-\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}\theta_{H}(n+\delta
n-\theta_{H}(b_{n+1}<b_{n})<\ell^{\prime})\\\ +\delta_{\ell^{\prime},n+\delta
n+1}[\theta_{H}(b_{n+1}<b_{n})\delta_{b_{n+1},j}-1]\\\
+\delta_{\ell^{\prime},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})]\,.$
(30)

#### More impurities: reduction to the single impurity case

From Eq. [3] and Eq. [21] it follows

$\braket{n_{1},n_{2},\dots,n_{\nu}}{\sigma^{z}_{y\ell^{\prime}-j}}{n_{1},n_{2},\dots,n_{\nu}}=\\\
\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}-2\sum_{\tilde{j}\in
B}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta
j,\ell^{\prime}}+2\sum_{\tilde{j}\in
B}\delta_{b_{\tilde{j}},j}\delta_{\tilde{j}+\delta\tilde{j}+\sum_{k=1}^{\nu}\theta(n_{k}+1-k<\tilde{j}),\ell^{\prime}}=\\\
\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}+2\sum_{k=1}^{\nu}\sum_{\tilde{j}>n_{k}+1-k}\delta_{b_{\tilde{j}},j}(\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}-k}-\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}-k+1})$
(31)

We can then use

$\braket{n}{\sigma_{y\ell^{\prime}-j}^{z}}{n}-\braket{\emptyset}{\sigma_{y\ell^{\prime}-j}^{z}}{\emptyset}=2\sum_{\tilde{j}>n}\delta_{b_{\tilde{j}},j}(\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}-1}-\delta_{\tilde{j}+\delta\tilde{j},\ell^{\prime}})$

to recast the expectation value as follows

$\braket{n_{1},n_{2},\dots,n_{\nu}}{:\sigma^{z}_{\ell}:}{n_{1},n_{2},\dots,n_{\nu}}=\sum_{k=1}^{\nu}\braket{n_{k}+1-k}{:\sigma_{\ell-y(k-1)}^{z}:}{n_{k}+1-k}$
(32)

where we used the notation introduced in Eq. [5]. Note that a similar
decomposition holds in standard $U(1)$ symmetric models in which $n_{j}$
simply represents the position of the $j$-th spin up; the main difference is
that here both the position of the operator and the pseudopositions of the
impurities are shifted depending on the number of the impurities on their
left.

## Mapping of diagonal operators

In this section we discuss the mapping of diagonal operators (in a basis
diagonalising $\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 characterise their local action.

### Class I

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

$\ket{\dots\uparrow_{j}\underbrace{\downarrow\dots\downarrow}_{n}(\uparrow_{j+n+1})\dots}\,.$
(33)

The corresponding contribution to their eigenvalue is summarised in Table 1.

|
$\ket{\dots\uparrow_{j}\underbrace{\downarrow\dots\downarrow}_{n}(\uparrow_{j+n+1})\dots}$
---|---
$1$ | $n+1$
$\sigma^{z}$ | $1-n$
$\sigma^{z}\otimes\sigma^{z}$ | $n-3+4\delta_{n,0}$
$\sigma^{z}\otimes\frac{1-\sigma^{z}}{2}\otimes\sigma^{z}$ | $n-4+4(\delta_{n,0}+\delta_{n,1})$
$\sigma^{z}\otimes(\frac{1-\sigma^{z}}{2})^{\otimes 2}\otimes\sigma^{z}$ | $n-5+4(\delta_{n,1}+\delta_{n,2})+5\delta_{n,0}$
$\sigma^{z}\otimes(\frac{1-\sigma^{z}}{2})^{\otimes 3}\otimes\sigma^{z}$ | $n-6+4(\delta_{n,3}+\delta_{n,2})+5\delta_{n,1}+6\delta_{n,0}$
Table 1: Local contribution to the eigenvalues of some class I operators
ordered by their range.

The duality mapping acts as follows

$\ket{\dots\uparrow_{j}\underbrace{\downarrow\dots\downarrow}_{n}(\uparrow_{j+n+1})\dots}\mapsto\ket{\dots\downarrow^{b}\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n}{y}\rfloor}(\downarrow^{b-n-1\\!\\!\\!\pmod{y}})\dots}\,,$
(34)

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

$n=ym+(b_{\tilde{j}}-b_{\tilde{j}+1}-1)\\!\\!\\!\pmod{y}=ym+\frac{y-1}{2}-\frac{1}{2}\sum_{n=1}^{y-1}\frac{\sin(\frac{2\pi
n(b_{\tilde{j}}-b_{\tilde{j}+1}-\frac{1}{2})}{y})}{\sin(\frac{n\pi}{y})}\,,$
(35)

where we used

$x\\!\\!\\!\pmod{y}=\frac{y-1}{2}-\frac{1}{2}\sum_{n=1}^{y-1}\frac{\sin(\frac{2\pi
n(x+\frac{1}{2})}{y})}{\sin(\frac{n\pi}{y})}\,.$ (36)

|
$\ket{\dots\downarrow^{b}\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n}{y}\rfloor}(\downarrow^{b-n-1\\!\\!\\!\pmod{y}})\dots}$
---|---
$\frac{1-\tau^{z}}{2}$ | $1$
$\frac{1+\tau^{z}}{2}$ | $\lfloor\frac{n}{y}\rfloor$
$\frac{1-\tau^{z}}{2}\otimes\bigl{(}\frac{1+\tau^{z}}{2}\bigr{)}^{\otimes k}\otimes\frac{1-\tau^{z}}{2}$ | $\delta_{\lfloor\frac{n}{y}\rfloor,k}$
Table 2: Local contribution to the eigenvalues of some diagonal operators in
the pseudospace.

The mapping of the operators in Table 1 can be worked out using Eq. [35] and
the action of the diagonal operators in the pseudospace reported in Table 2.
We recognise 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 of Table 1 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}\,.$ (37)

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 the additive constant can also be represented in the configuration
space as follows

$c\rightarrow c\sum_{\tilde{j}\in B}1\,,$ (38)

which we prefer as it shows more explicitly that constants are conserved. In
the following we will always replace $\sum_{j\in I}\tau_{j}^{z}$ by
$\sum_{j\in I}1-2\sum_{\tilde{j}\in B}1$.

Number ($n$) of spins down:

This contribution can be worked out using Eq. [35], which links the number $n$
of spins down in between spins up with the number $\lfloor\frac{n}{y}\rfloor$
of pseudospins up in the corresponding domain in between pseudospins down. In
particular, there is also a term written solely in terms of the species, which
can be treated as a constant by virtue of the conservation of the
configuration. Putting all together we get

$n\rightarrow\sum_{j\in I}y1-\frac{y+1}{2}\sum_{\tilde{j}\in
B}1-\frac{1}{2}\sum_{\tilde{j}\in B}\sum_{n=1}^{y-1}\frac{\sin(\frac{2\pi
n(b_{\tilde{j}}-b_{\tilde{j}+1}-\frac{1}{2})}{y})}{\sin(\frac{n\pi}{y})}$ (39)

Kronecker deltas $\delta_{n,k}$:

These terms characterise the eigenvalues of the operators with densities that
consist of the tensor product of $k$ projectors $\frac{1-\sigma^{z}}{2}$ onto
spin down, preceded and followed by the projector $\frac{1+\sigma^{z}}{2}$
onto spin up (cf. Table 2 for the analogue in the pseudospin space). The
representation in the pseudospin space can then be derived from the identity

$\delta_{n,k}=\delta_{\lfloor\frac{n}{y}\rfloor,\lfloor\frac{k}{y}\rfloor}\delta_{(b_{\tilde{j}}-b_{\tilde{j}+1}-1-k)\\!\\!\\!\pmod{y},0}=\delta_{\lfloor\frac{n}{y}\rfloor,\lfloor\frac{k}{y}\rfloor}\frac{1}{y}\sum_{m=0}^{y-1}e^{2\pi
im\frac{b_{\tilde{j}}-b_{\tilde{j}+1}-1-k}{y}}$ (40)

which results in the representation

$\delta_{n,k}\rightarrow\sum_{j\in
I}(Z^{(1+\lfloor\frac{k}{y}\rfloor)}_{1+k\\!\\!\\!\\!\\!\pmod{y}})_{j}\,,$
(41)

where $(Z^{(n)}_{s\\!\\!\pmod{y}})_{j}$ is the following semilocal operator

$(Z^{(n)}_{s\\!\\!\\!\\!\\!\pmod{y}})_{j}=\int_{-\pi}^{\pi}\frac{\mathrm{d}\omega}{2\pi}\sum_{k\in
B}\sum_{m=0}^{y-1}\frac{e^{2\pi
im\frac{b_{k}-b_{k+1}-s}{y}}}{y}e^{i\omega(k-1-j+\sum_{j^{\prime}}^{j}\frac{1+\tau^{z}_{j}}{2})}\frac{1-\tau_{j}^{z}}{2}\Bigl{(}\prod_{m=1}^{n-1}\frac{1+\tau_{j+m}^{z}}{2}\Bigr{)}\frac{1-\tau_{j+n}^{z}}{2}\,.$
(42)

Here the integral over $\omega$ represents a Kronecker delta and the sum over
$m$ a condition on the species of particles modulo $y$.

Before considering other classes of diagonal operators, we point out the
identity

$\sum_{s=1}^{y}(Z_{s\\!\\!\\!\\!\\!\pmod{y}}^{(n)})_{j}=\frac{1-\tau_{j}^{z}}{2}\Bigl{(}\prod_{m=1}^{n-1}\frac{1+\tau_{j+m}^{z}}{2}\Bigr{)}\frac{1-\tau_{j+n}^{z}}{2}$
(43)

which shows that the sum over $s$ is a (strictly) local operator. This allows
us to identify some diagonal local operators that remain local after the
duality transformation. They are obtained by summing Eq. [41] over $k$ for
$my\leq k<(m+1)y$ and given integer $m$; that is to say, operators

$V_{m}=\sum_{\ell}\sum_{n=my+1}^{m(y+1)}P^{(+)}_{\ell}\Bigl{(}\prod_{j=1}^{n-1}P^{(-)}_{\ell+j}\Bigr{)}P^{(+)}_{\ell+n}\,,$
(44)

where

$P^{(\pm)}=\frac{1\pm\sigma^{z}}{2}\,,$ (45)

are symmetric under transmutation of species. The telescoping sum can be
readily simplified and we can also write

$V_{m}=\sum_{\ell}P^{(+)}_{\ell}\Bigl{(}\prod_{j=1}^{my}P^{(-)}_{\ell+j}\Bigr{)}\Bigl{[}1-\Bigl{(}\prod_{j^{\prime}=1}^{y}P^{(-)}_{\ell+my+j^{\prime}}\Bigr{)}\Bigr{]}$
(46)

### Class II

We now consider diagonal operators that are still characterised by their local
action around a particle, this time however the number of spins down both
after and before the particle are important. Few examples of their local
action is reported in Table 3.

|
$\ket{\dots(\uparrow)\underbrace{\downarrow\dots\downarrow}_{n_{-}}\uparrow_{j}\underbrace{\downarrow\dots\downarrow}_{n_{+}}(\uparrow)\dots}$
---|---
$\sigma^{z}\otimes\frac{1+\sigma^{z}}{2}\otimes\sigma^{z}$ | $(1-2\delta_{n_{-},0})(1-2\delta_{n_{+},0})$
$\sigma^{z}\otimes\frac{1-\sigma^{z}}{2}\otimes\frac{1+\sigma^{z}}{2}\otimes\sigma^{z}$ | $(2\delta_{n_{-},1}+\delta_{n_{-},0}-1)(2\delta_{n_{+},0}-1)$
$\sigma^{z}\otimes\frac{1+\sigma^{z}}{2}\otimes\frac{1-\sigma^{z}}{2}\otimes\sigma^{z}$ | $(2\delta_{n_{-},0}-1)(2\delta_{n_{+},1}+\delta_{n_{+},0}-1)$
$\sigma^{z}\otimes\frac{1-\sigma^{z}}{2}\otimes\frac{1-\sigma^{z}}{2}\otimes\frac{1+\sigma^{z}}{2}\otimes\sigma^{z}$ | $(2\delta_{n_{-},2}+\delta_{n_{-},1}+\delta_{n_{-},0}-1)(2\delta_{n_{+},0}-1)$
$\sigma^{z}\otimes\frac{1-\sigma^{z}}{2}\otimes\frac{1+\sigma^{z}}{2}\otimes\frac{1-\sigma^{z}}{2}\otimes\sigma^{z}$ | $(2\delta_{n_{-},1}+\delta_{n_{-},0}-1)(2\delta_{n_{+},1}+\delta_{n_{+},0}-1)$
$\sigma^{z}\otimes\frac{1+\sigma^{z}}{2}\otimes\frac{1-\sigma^{z}}{2}\otimes\frac{1-\sigma^{z}}{2}\otimes\sigma^{z}$ | $(2\delta_{n_{-},0}-1)(2\delta_{n_{+},2}+\delta_{n_{+},1}+\delta_{n_{+},0}-1)$
Table 3: Local contribution to the eigenvalues of some class II operators,
ordered by their range.

The duality transformation can be worked out using that, in the pseudospace,
Eq. [47] holds.

$\frac{1-\tau_{j-k_{1}-1}^{z}}{2}\Bigl{(}\prod_{j^{\prime}=1}^{k_{1}}\frac{1+\tau_{j-j^{\prime}}^{z}}{2}\Bigr{)}\frac{1-\tau_{j}^{z}}{2}\Bigl{(}\prod_{j^{\prime\prime}=1}^{k_{2}}\frac{1+\tau_{j+j^{\prime\prime}}^{z}}{2}\Bigr{)}\frac{1-\tau_{j+k_{2}+1}^{z}}{2}\\\
\ket{\dots(\downarrow^{b+n_{-}+1\\!\\!\\!\pmod{y}})\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n_{-}}{y}\rfloor}\downarrow_{j}^{b}\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n_{+}}{y}\rfloor}(\downarrow^{b-n_{+}-1\\!\\!\\!\pmod{y}})\dots}=\\\
\delta_{\lfloor\frac{n_{-}}{y}\rfloor,k_{1}}\delta_{\lfloor\frac{n_{+}}{y}\rfloor,k_{2}}\ket{\dots(\downarrow^{b+n_{-}+1\\!\\!\\!\pmod{y}})\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n_{-}}{y}\rfloor}\downarrow_{j}^{b}\underbrace{\uparrow\dots\uparrow}_{\lfloor\frac{n_{+}}{y}\rfloor}(\downarrow^{b-n_{+}-1\\!\\!\\!\pmod{y}})\dots}\,.$
(47)

Indeed, besides contributions that we worked out also for class I operators,
class II operators are characterised by an additional one (cf. Table 3):

Product of two Kronecker deltas $\delta_{n_{-},k_{1}}\delta_{n_{+},k_{1}}$:

This kind of contributions can be worked out by properly multiplying the
representation of the delta in Eq. [40]:

$\delta_{n_{-},k_{1}}\delta_{n_{+},k_{2}}=\delta_{\lfloor\frac{n_{-}}{y}\rfloor,\lfloor\frac{k_{1}}{y}\rfloor}\delta_{\lfloor\frac{n_{+}}{y}\rfloor,\lfloor\frac{k_{2}}{y}\rfloor}\frac{1}{y^{2}}\sum_{m,m^{\prime}=0}^{y-1}e^{2\pi
im\frac{b_{\tilde{j}}-b_{\tilde{j}+1}-1-k_{2}}{y}}e^{2\pi
im^{\prime}\frac{b_{\tilde{j}-1}-b_{\tilde{j}}-1-k_{1}}{y}}$ (48)

Thus we find

$\delta_{n_{-},k_{1}}\delta_{n_{+},k_{2}}\rightarrow\sum_{j\in
I}(Z^{(\lfloor\frac{k_{1}}{y},\lfloor\frac{k_{2}}{y}\rfloor)}_{1+k_{1}\\!\\!\\!\\!\\!\pmod{y},1+k_{2}\\!\\!\\!\\!\\!\pmod{y}})_{j}\,,$
(49)

where
$(Z^{(n_{1},n_{2})}_{s_{1}\\!\\!\\!\\!\\!\pmod{y},s_{2}\\!\\!\\!\\!\\!\pmod{y}})_{j}$
is the semilocal operator

$(Z^{(n_{1},n_{2})}_{s_{1}\\!\\!\\!\\!\\!\pmod{y},s_{2}\\!\\!\\!\\!\\!\pmod{y}})_{j}=\int_{-\pi}^{\pi}\frac{\mathrm{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}-s_{1})+m_{2}(b_{k+1}-b_{k+2}-s_{2})}{y}}}{y^{2}}e^{i\omega(k+\sum_{\tilde{j}}^{0}1-\sum_{j^{\prime}}^{j}\frac{1-\tau^{z}_{j}}{2})}\\\
\frac{1-\tau_{j}^{z}}{2}\Bigl{(}\prod_{m=1}^{n_{1}-1}\frac{1+\tau_{j+m}^{z}}{2}\Bigr{)}\frac{1-\tau_{j+n_{1}}^{z}}{2}\Bigl{(}\prod_{m=1}^{n_{2}-1}\frac{1+\tau_{j+n_{1}+m}^{z}}{2}\Bigr{)}\frac{1-\tau_{j+n_{1}+n_{2}}^{z}}{2}\,.$
(50)

We point out that identity Eq. [43] is generalised in that the sum over all
$s$ results in a (strictly) local operator.

### 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_{+}$. 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 Eq.
[40] should be chosen according to which species surround the associated
domains of spins down. The product of deltas is then represented by
generalisations of Eq. [42] along the line of Eq. [50].

## 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. [11].

We have

$:\tilde{H}:\ket{n_{1},n_{2}}=2J[\ket{n_{1}-1,n_{2}}+(1-\delta_{n_{2},n_{1}+1})(\ket{n_{1}+1,n_{2}}\\\
+\ket{n_{1},n_{2}-1})+\ket{n_{1},n_{2}+1}]+4J\Delta(\delta_{n_{2},n_{1}+1}-2)\ket{n_{1},n_{2}}\\\
+2g\delta_{n_{2},n_{1}+1}(\ket{n_{1}-1,n_{1}}+\ket{n_{1}+1,n_{1}+2})\,,$ (51)

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

$\ket{\Psi(p)}=\sum_{n}e^{inp}\ket{n}$ (52)

normalised according to
$\braket{\Psi(p_{1})}{\Psi(p_{2})}=2\pi\delta(p_{1}-p_{2})$; the corresponding
excitation energy reads

$E(p)=4J(\cos p-\Delta)\,.$ (53)

The Bethe Ansatz for the two-impurity problem can be written as follows

$\ket{\Psi(p_{1},p_{2})}\sim\sum_{n_{1}<n_{2}}[e^{in_{1}p_{1}+in_{2}p_{2}}+S(p_{1},p_{2})e^{in_{1}p_{2}+in_{2}p_{1}}]\ket{n_{1},n_{2}}$
(54)

where $S(p_{1},p_{2})$ is the scattering phase. Enforcing the secular equation
gives

$S(p_{1},p_{2})=-\frac{1-2[\Delta+\frac{g}{J}\cos(p_{1}+p_{2})]e^{ip_{2}}+e^{i(p_{1}+p_{2})}}{1-2[\Delta+\frac{g}{J}\cos(p_{1}+p_{2})]e^{ip_{1}}+e^{i(p_{1}+p_{2})}}\,.$
(55)

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 regularisation with periodic boundary conditions. The
momenta satisfy in turn the Bethe equations

$e^{iLp_{1}}=S(p_{2},p_{1})\qquad e^{iLp_{2}}=S(p_{1},p_{2})$ (56)

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 behaviour of the two cumulative
distributions

$F_{n}^{(1,2)}(t)=\sum_{n_{1}\leq n}p_{n_{1}}^{(1)}(t)=\sum_{n_{1}\leq
n}\sum_{n_{2}>n_{1}}|\braket{n_{1},n_{2}}{e^{-i\tilde{H}_{eff}t}\tilde{P}}{\Psi_{0}}|^{2}\,.$
(57) $F_{n}^{(2,2)}(t)=\sum_{n_{2}\leq
n}p_{n_{2}}^{(2)}(t)=\sum_{n_{1}<n_{2}\leq
n}|\braket{n_{1},n_{2}}{e^{-i\tilde{H}_{eff}t}\tilde{P}}{\Psi_{0}}|^{2}\,.$
(58)

### Real solutions

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

$\ket{\Psi(p_{1},p_{2})}=\sum_{-\frac{L}{2}<n_{1}<n_{2}\leq\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}}]\ket{n_{1},n_{2}}$
(59)

where

$A_{p_{1},p_{2}}=1-2[\Delta+\frac{g}{J}\cos(p_{1}+p_{2})]e^{ip_{1}}+e^{i(p_{1}+p_{2})}\,.$
(60) $(A_{-p_{1},-p_{2}}=A_{p_{2},p_{1}}e^{-i(p_{1}+p_{2})})$ (61)

In the thermodynamic limit such states are normalised according to

$\braket{\Psi(p_{1},p_{2})}{\Psi(p_{1}^{\prime},p_{2}^{\prime})}=(2\pi)^{2}\delta(p_{1}^{\prime}-p_{1})\delta(p_{2}^{\prime}-p_{2})|A(p_{1},p_{2})|^{2}\,,$
(62)

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_{p_{j}\atop\mathrm{Bethe\
eqs}}f(p_{j})=\oint_{[-\pi,\pi]_{\epsilon}}\frac{\mathrm{d}p}{2\pi}\frac{L-i\partial_{p}\log
S(p,p_{\bar{j}})}{1-e^{-iLp}S(p_{\bar{j}},p)}f(p)\rightarrow
L\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}[1-\frac{i\partial_{p}\log
S(p,p_{\bar{j}})}{L}]f(p)\,,$ (63)

where $[-\pi,\pi]_{\epsilon}$ is the rectangular contour surrounding the
interval $[-\pi,\pi]$ with imaginary part $\pm i\epsilon$. 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_{n,1}^{(j,2)}(t)$ to $F_{n}^{(j,2)}(t)$ 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_{n,1}^{(1,2)}(t)=\frac{1}{4}\sum_{m_{1}=0}^{\infty}\sum_{m=1}^{\infty}\Bigl{|}\int\frac{\mathrm{d}^{2}p}{(2\pi)^{2}}e^{-i[E(p_{1})+E(p_{2})]t}\frac{\braket{-m_{1}+n,-m_{1}+n+m}{\Psi(p_{1},p_{2})}\braket{\Psi(p_{1},p_{2})}{n_{1}^{(0)},n_{2}^{(0)}}}{|A(p_{1},p_{2})|^{2}}\Bigr{|}^{2}=\\\
\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}p^{\prime}}{(2\pi)^{4}}\frac{e^{-i[E(p_{1})+E(p_{2})-E(p_{1}^{\prime})-E(p_{2}^{\prime})]t}e^{in(p_{1}+p_{2}-p_{1}^{\prime}-p_{2}^{\prime})}}{1-e^{-i(p_{1}+p_{2}-p_{1}^{\prime}-p_{2}^{\prime}-i\epsilon)}}\frac{\braket{\Psi(p_{1},p_{2})}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi(p_{1}^{\prime},p_{2}^{\prime})}}{|A(p_{1},p_{2})|^{2}|A(p_{1}^{\prime},p_{2}^{\prime})|^{2}}\\\
\times\Bigl{[}\frac{A(p_{1},p_{2})A(-p_{1}^{\prime},-p_{2}^{\prime})}{e^{-i(p_{2}-p_{2}^{\prime}+i\epsilon)}-1}-\frac{A(p_{2},p_{1})A(-p_{1}^{\prime},-p_{2}^{\prime})}{e^{-i(p_{1}-p_{2}^{\prime}+i\epsilon)}-1}\Bigr{]}$
(64)

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_{n,1}^{(1,2)}(t)=\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}p^{\prime}}{(2\pi)^{4}}\frac{e^{-i[E(p_{1})+E(p_{2})-E(p_{1}^{\prime}+p_{1}-p_{2}^{\prime})-E(p_{2}^{\prime}+p_{2})]t}e^{-inp_{1}^{\prime}}}{1-e^{i(p_{1}^{\prime}+i\epsilon)}}\\\
\frac{\braket{\Psi(p_{1},p_{2})}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi(p_{1}^{\prime}+p_{1}-p_{2}^{\prime},p_{2}^{\prime}+p_{2})}}{|A(p_{1},p_{2})|^{2}|A(p_{1}^{\prime}+p_{1}-p_{2}^{\prime},p_{2}^{\prime}+p_{2})|^{2}}\frac{A(p_{1},p_{2})A(-p_{1}^{\prime}-p_{1}+p_{2}^{\prime},-p_{2}^{\prime}-p_{2})}{e^{i(p_{2}^{\prime}-i\epsilon)}-1}\\\
-\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}p^{\prime}}{(2\pi)^{4}}\frac{e^{-i[E(p_{1})+E(p_{2})-E(p_{1}^{\prime}+p_{2}-p_{2}^{\prime})-E(p_{2}^{\prime}+p_{1})]t}e^{-inp_{1}^{\prime}}}{1-e^{i(p_{1}^{\prime}+i\epsilon)}}\\\
\frac{\braket{\Psi(p_{1},p_{2})}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi(p_{1}^{\prime}+p_{2}-p_{2}^{\prime},p_{2}^{\prime}+p_{1})}}{|A(p_{1},p_{2})|^{2}|A(p_{1}^{\prime}+p_{2}-p_{2}^{\prime},p_{2}^{\prime}+p_{1})|^{2}}\frac{A(p_{2},p_{1})A(-p_{1}^{\prime}-p_{2}+p_{2}^{\prime},-p_{2}^{\prime}-p_{1})}{e^{i(p_{2}^{\prime}-i\epsilon)}-1}\\\
\sim\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p}{(2\pi)^{2}}\theta_{H}(n-\min(v(p_{1}),v(p_{2}))t)\frac{\braket{\Psi(p_{1},p_{2})}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi(p_{1},p_{2})}}{|A(p_{1},p_{2})|^{2}}$
(65)

Thus we have

$F_{n,1}^{(1,2)}(t)\sim\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{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}$
(66)

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

$F_{n,1}^{(2,2)}(t)\sim\frac{1}{2}\int_{-\pi}^{\pi}\frac{\mathrm{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}$
(67)

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

$\lim_{n\rightarrow\infty}F_{n,1}^{(j,2)}(t)\sim
1+\mathrm{Re}\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p}{(2\pi)^{2}}e^{i(n_{2}^{(0)}-n_{1}^{(0)})(p_{1}-p_{2})}S(p_{1},p_{2})$
(68)

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

### Bound states

Except for the region $g=J\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+\frac{g}{J}\cos p)^{2}\,,$ (69)

i.e.

$\cos
p>\frac{1-4\frac{g}{J}\Delta+\sqrt{1+8\frac{g^{2}}{J^{2}}-8\frac{g}{J}\Delta}}{4\frac{g^{2}}{J^{2}}}\quad\vee\quad\cos
p<\frac{1-4\frac{g}{J}\Delta-\sqrt{1+8\frac{g^{2}}{J^{2}}-8\frac{g}{J}\Delta}}{4\frac{g^{2}}{J^{2}}}\,.$
(70)

Note that this condition is always satisfied when
$\frac{g}{J}\Delta>\frac{1}{8}+\frac{g^{2}}{J^{2}}$ or
$1<(\Delta+\frac{g}{J})^{2}\wedge 4g(\Delta+\frac{g}{J})<1$. The corresponding
bound state reads

$\ket{\Psi_{2}(p)}\sim\sqrt{\Bigl{(}\frac{\Delta+g\cos
p}{\cos\frac{p}{2}}\Big{)}^{2}-1}\sum_{n_{1}<n_{2}}e^{i\frac{n_{1}+n_{2}}{2}p}\Bigl{(}\frac{\cos\frac{p}{2}}{\Delta+g\cos
p}\Bigr{)}^{n_{2}-n_{1}}\ket{n_{1},n_{2}}\,,$ (71)

where we enforced the standard normalisation
$\braket{\Psi(p^{\prime})}{\Psi(p)}=2\pi\delta(p-p^{\prime})$. Its energy is
given by

$E_{2}(p)=4g\cos p+4J[\frac{\cos^{2}\frac{p}{2}}{\Delta+\frac{g}{J}\cos
p}-\Delta]$ (72)

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

$F_{n,2}^{(1,2)}(t)=\sum_{m_{1}=0}^{\infty}\sum_{m=1}^{\infty}|\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\theta_{H}(1+\cos
p<2(\Delta+\frac{g}{J}\cos
p)^{2})e^{-iE_{2}(p)t}\braket{-m_{1}+n,-m_{1}+n+m}{\Psi_{2}(p)}\braket{\Psi_{2}(p)}{n_{1}^{(0)},n_{2}^{(0)}}|^{2}=\\\
\int_{-\pi}^{\pi}\frac{\mathrm{d}p\mathrm{d}p^{\prime}}{(2\pi)^{2}}\theta_{H}(1+\cos
p<2(\Delta+\frac{g}{J}\cos
p)^{2})\theta_{H}(1+\cos(p^{\prime}+p)<2(\Delta+\frac{g}{J}\cos(p^{\prime}+p))^{2})\\\
\sqrt{\Bigl{(}\frac{\Delta+g\cos
p}{\cos\frac{p}{2}}\Big{)}^{2}-1}\sqrt{\Bigl{(}\frac{\Delta+g\cos(p^{\prime}+p)}{\cos\frac{p^{\prime}+p}{2}}\Big{)}^{2}-1}\frac{e^{-inp^{\prime}}}{1-e^{i(p^{\prime}+i\epsilon)}}\frac{e^{-i[E_{2}(p)-E_{2}(p^{\prime}+p)]t}}{\frac{[\Delta+g\cos
p][\Delta+g\cos(p^{\prime}+p)]}{e^{-i\frac{p^{\prime}}{2}}\cos\frac{p}{2}\cos\frac{p^{\prime}+p}{2}}-1}\braket{\Psi_{2}(p)}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi_{2}(p^{\prime}+p)}\sim\\\
\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\theta_{H}(1+\cos
p<2(\Delta+\frac{g}{J}\cos
p)^{2})\theta_{H}(n-v_{2}(p))\braket{\Psi_{2}(p)}{n_{1}^{(0)},n_{2}^{(0)}}\braket{n_{1}^{(0)},n_{2}^{(0)}}{\Psi_{2}(p)}$
(73)

Thus we have

$F_{n,2}^{(1,2)}(t)\sim\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\theta_{H}(n-v_{2}(p)t)\max\left[\Bigl{(}\frac{\Delta+g\cos
p}{\cos\frac{p}{2}}\Big{)}^{2}-1,0\right]\Bigl{(}\frac{\cos^{2}\frac{p}{2}}{[\Delta+g\cos
p]^{2}}\Bigr{)}^{n_{2}^{(0)}-n_{1}^{(0)}}\,.$ (74)

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

Importantly, despite the different appearance of the integrals, the limit
$n\rightarrow\infty$ exactly compensates the displacement from $1$ that we got
in Eq. [68] considering only the real solutions, supporting the completeness
of the Bethe states identified in this sector.

## Macroscopicity of the quantum state

In this section we work out the variance of $S^{z}_{\vec{s}}$ when the system
is prepared in a product state with a single impurity. We have

$\braket{S_{\vec{s}}^{z},S_{\vec{s}}^{z}}_{t}=\frac{1}{4}\sum_{\ell}(1-\braket{\sigma_{\ell}^{z}}^{2}_{t})+\frac{1}{4}\sum_{\ell_{1}\neq\ell_{2}}s_{\ell_{1}}s_{\ell_{2}}[\sum_{n}p_{n,t}(1-p_{n,t})\braket{n}{:\sigma_{\ell_{1}}^{z}::\sigma_{\ell_{2}}^{z}:}{n}-\sum_{n\neq
n^{\prime}}p_{n,t}p_{n^{\prime},t}\braket{n}{:\sigma_{\ell_{1}}^{z}:}{n}\braket{n^{\prime}}{:\sigma_{\ell_{2}}^{z}:}{n^{\prime}}]$
(75)

In the symmetric case the asymptotic behaviour can be obtained using the
simplified expression

$\braket{n;\underline{b}}{:\sigma_{\ell_{1}}^{z}:}{n;\underline{b}}\sim\braket{\emptyset;\underline{b}}{\sigma_{\ell_{1}-y}^{z}-\sigma_{\ell_{1}}^{z}}{\emptyset;\underline{b}}\theta_{H}(n<x_{\ell_{1}})\,,$
(76)

indeed the correction is suppressed by $O(p)$, which is assumed to approach
zero in time. We then find

$\braket{S_{\vec{s}}^{z},S_{\vec{s}}^{z}}_{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}})]s_{\ell_{1}}s_{\ell_{2}}\braket{\emptyset;\underline{b}}{\sigma_{\ell_{1}-y}^{z}-\sigma_{\ell_{1}}^{z}}{\emptyset;\underline{b}}\braket{\emptyset;\underline{b}}{\sigma_{\ell_{2}-y}^{z}-\sigma_{\ell_{2}}^{z}}{\emptyset;\underline{b}}\,.$
(77)

Choosing, for example,
$s_{\ell}=\braket{\emptyset;\underline{b}}{\sigma_{\ell}^{z}}{\emptyset;\underline{b}}$
gives

$\braket{S_{\vec{s}}^{z},S_{\vec{s}}^{z}}_{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}})]\braket{\emptyset;\underline{b}}{\frac{1-\sigma_{\ell_{1}-y}^{z}\sigma_{\ell_{1}}^{z}}{2}}{\emptyset;\underline{b}}\braket{\emptyset;\underline{b}}{\frac{1-\sigma_{\ell_{2}-y}^{z}\sigma_{\ell_{2}}^{z}}{2}}{\emptyset;\underline{b}}\,.$
(78)

Generally this grows as $t^{2}$, indeed $F_{t}(x)$ approaches a function of
$x/t$ and
$\braket{\emptyset;\underline{b}}{\frac{1-\sigma_{\ell-y}^{z}\sigma_{\ell}^{z}}{2}}{\emptyset;\underline{b}}$
has a strictly positive local average. This implies that also the second line
of Eq. [20] grows quadratically in time.