# Ultraslow Growth of Number Entropy in an $\ell$-Bit Model of Many-Body Localization

David Aceituno Chávez, Claudia Artiaco, Thomas Klein Kvorning, Loïc Herviou, Jens H Bardarson

► To cite this version:

David Aceituno Chávez, Claudia Artiaco, Thomas Klein Kvorning, Loïc Herviou, Jens H Bardarson. Ultraslow Growth of Number Entropy in an $\ell$-Bit Model of Many-Body Localization. Physical Review Letters, 2024, 133 (12), pp.126502. 10.1103/physrevlett.133.126502. hal-04726996

HAL Id: hal-04726996

https://cnrs.hal.science/hal-04726996v1

Submitted on 9 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Ultraslow Growth of Number Entropy in an $\ell$-bit Model of Many-Body Localization

David Aceituno Chávez,¹ Claudia Artiaco,¹ Thomas Klein Kvorning,¹ Loïc Herviou,² and Jens H. Bardarson¹

¹Department of Physics, KTH Royal Institute of Technology, Stockholm, 106 91 Sweden

²Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

We demonstrate that slow growth of the number entropy following a quench from a local product state is consistent with many-body localization. To do this we construct a random circuit $\ell$-bit model with exponentially localized $\ell$-bits and exponentially decaying interactions between them. We observe an ultraslow growth of the number entropy starting from a Néel state, saturating at a value that grows with system size. This suggests that the observation of such growth in microscopic models is not sufficient to rule out many-body localization.

The many-body localized insulator [1, 2] is widely believed to be an emergent phase of interacting quasilocal integrals of motion, often referred to as $\ell$-bits [3–5]. The $\ell$-bits in spin-1/2 models are dressed physical onsite Pauli spin operators [3, 4], while in fermionic models they are Anderson orbitals dressed by local electron-hole excitations [6, 7]. Deep in the localized phase, the $\ell$-bits have large overlaps with the physical spins. Consequently, in a global quench from a local product state, such as the Néel state, the anti-ferromagnetic order remains nonzero at all times [8], as the overlaps of the $\ell$-bits with the initial state are constants of the motion [3]. In contrast, in thermalizing systems the anti-ferromagnetic order decays to zero due to the diffusion of magnetization (equivalently, transport of particles) [9]. Nonetheless, in localized systems the entanglement entropy grows slowly towards an extensive sub-thermal volume law [10, 11]. This behavior has been explained by quasilocal (exponentially decaying) interactions between the $\ell$-bits, which induce entanglement through dephasing between different energy states [12]. In this picture, there is no transport of particles, apart from an initial expansion of local densities within the single-particle localization length. This is consistent with early numerical work on many-body localization, which explored half-chain particle number fluctuations and found it to not increase with system size [13].

More recent work [14] investigating particle transport via the number entropy $S_N$—the Shannon entropy of the particle number distribution in half of the chain—has found evidence for it increasing extremely slowly with time and saturating in finite systems at a value that increases with system size, albeit in a subextensive way. The authors of this work interpreted their findings to imply the presence of transport and the absence of many-body localization. While these first results were likely obtained either on the ergodic side of the transition [15], or in the purported pre-thermal finite-size crossover region [16–19], later simulations deeper in the localized phase obtained consistent results [20]. These observations have been suggested to arise from rare resonances and their associated transport [21], though this was subject to debate [22, 23].

Other recent work based on finite-size numerics has ei-

ther cast doubt on the existence of many-body localization [24, 25], or reinterpreted the localized phase diagram to have a large nonlocalized crossover region, shifting the transition to the localized phase to a much larger disorder strength than hitherto believed [16–19]. It can, however, be difficult to interpret finite-size numerics without a definite scaling theory of the localization transition [26, 27]. It also seems clear that, even if there is no many-body localization in a large region of disorder strengths, both finite-size numerics and finite-time experiments are dominated by the physics of many-body localization, albeit in some cases with an evident creep towards more extensive states [28–32]. With more detailed observables and numerics, it seems pertinent to also refine our understanding of the many-body localization phenomenology.

Here we investigate to what extent the growth of number entropy can be explained directly within the phenomenology of many-body localization. Instead of working with microscopic models, we work directly with an effective $\ell$-bit model of localized interacting integrals of motion. The $\ell$-bit model may be sufficiently rich in features to host phenomena beyond those initially understood. In the $\ell$-bit model, particle number fluctuations arise both in the exponentially decaying support of the $\ell$-bits and from their exponentially decaying interactions. While it is expected that particle transport in the localized phase is minimal, in line with earlier findings, it is not a priori obvious that transport is entirely absent. We therefore construct a comprehensive $\ell$-bit model of many-body localization that takes into account both of these exponential tails, and measure the number entropy through numerical simulations following a quench from a local (in the physical spins) product state. Since the Hamiltonian in the $\ell$-bit basis is diagonal, we can simulate to large times and for relatively large systems—computational limitations arise from the basis change from the $\ell$-bits to the physical spins, required for calculating observables. Focusing on the quench dynamics, we show that the $\ell$-bit model exhibits ultraslow growth of the number entropy, which in a finite system saturates at a subextensive value that nevertheless increases with system size. We find a growth in time consistent with $\ln \ln t$, though other functional forms can not be ruled out. We
interpret this as resulting from the two exponential tails in the $\ell$-bit model, both of which, as already explained, contain particle fluctuations. Our observation of slow growth of number entropy implies that the $\ell$-bit model is richer than previously understood, and observation of such slow growth in microscopic models may therefore still be consistent with many-body localization.

**Number entropy.**— Our main observable of interest is the number entropy $S_N^A$, which is the contribution to the entanglement entropy $S_E^A$ arising from the particle number fluctuation between subsystems A and $\bar{A}$ (where $A \cup \bar{A}$ is the whole system) [33–35]. We consider pure states $|\psi\rangle$ describing L qubits labeled by $n_i \in \{0, 1\}$ on each site $i$ of a one-dimensional lattice. For $A = [0...i]$, the number entropy is defined as $S_N^A = -\sum_n p_i(n) \ln p_i(n)$, where

$$p_i(n) = \sum_{\substack{n_0...n_{L-1}, \\ \text{such that} \\ \sum_{k=0}^{i} n_k = n}} |\langle n_0...n_{L-1} | \psi \rangle|^2, \quad (1)$$

is the probability that n qubits in A are measured to be in the state $|1\rangle$. From $p_i(n)$, we also obtain the probability $q_i(n)$ that the particle number grows from $n-1$ to n at site i (counting from the left):

$$q_i(n) = \sum_{m=n}^{L} p_i(m) - p_{i-1}(m), \quad (n \ge 1), \quad (2)$$

where we set $p_{-1}(m) = 0$. This satisfies $\sum_{i=0}^{L-1} q_i(n) = 1$ and $\sum_{n=1}^{L} q_i(n) = \langle n_i \rangle$, and can be used to calculate the moments of position $x$ for the nth particle, that is, $\langle x_n \rangle = \sum_{i=1}^{L} iq_i(n)$ and $\sigma_{x_n}^2 = \sum_{i=1}^{L} (i - \langle x_n \rangle)^2 q_i(n)$ [36].

We denote by $S_E$ and $S_N$ the half-chain entropies for the subsystem $\bar{A} = [0...L/2)$. Furthermore, we will refer to each $|1\rangle$ as particle, charge, or magnetic excitation interchangeably, and focus exclusively on systems where their number is conserved on $A \cup \bar{A}$. Lastly, we interpret $\sigma_{x_n}^2 > 0$ as the spread of the nth particle $|1\rangle$ onto multiple sites by becoming dressed with particle-hole excitations, in which case we refer to it as a *quasiparticle*.

**Model.**—We consider Hamiltonians $H$ in the many-body-localized regime describing spin-1/2 chains of length $L$ with short-range interactions, in the sector of zero magnetization (half-filling). This includes the disordered Heisenberg XXZ spin chain, which maps to the $t-V$-model of interacting fermions via the Jordan-Wigner transformation [10, 11, 37]. In such Hamiltonians, physical spins are the set of Pauli spin operators $\{\vec{\sigma}_i\}$ acting on the Hilbert space of site $i$. A corresponding set of $\ell$-bits $\tau_i^z \equiv U\sigma_i^z U^\dagger$ emerges in the many-body localized phase, which are physical spins dressed by a quasilocal finite-depth unitary transformation $U$ that diagonalizes the Hamiltonian [3, 4, 38, 39].

To study the growth of number entropy, we need to time evolve large spin chains to large times with high

precision. We achieve this by time evolving matrix product states directly in the diagonal basis $H' = U^\dagger H U$, thereby avoiding Trotterization errors [36, 40],

$$|\psi(t)\rangle = U e^{-iH't} U^\dagger |\psi(0)\rangle. \quad (3)$$

We thus model an *effective* Hamiltonian $H$ by two random components: the diagonal $\ell$-bit Hamiltonian $H'$, which models the interactions between $\ell$-bits, and the transformation $U$ expressed as a finite-depth circuit of two-site gates. The $\ell$-bit Hamiltonian is defined as [4]

$$H' = \sum_{i} h_i \tau_i^z + \sum_{i<j} J_{ij} \tau_i^z \tau_j^z + \sum_{i<j<k} J_{ijk} \tau_i^z \tau_j^z \tau_k^z + \dots, \quad (4)$$

where the $h$'s are random fields and $J$'s are random quasilocal couplings. On average, the couplings decay exponentially with their index diameter, such that $J_{ij} \propto e^{-|i-j|/\xi_J^j}$, $J_{ijk} \propto e^{-|i-k|/\xi_J^{jk}}$, and so on, where $\xi_J^j, \xi_J^{jk} \ll L$ are distinct characteristic length scales. We implement $H'$ in the $\ell$-bit basis (where $\tau_i^z$ are diagonal Pauli matrices) keeping terms up to third order [41]:

$$h_i = R_i \qquad i = 0, 1, \dots, L-1, \quad (5)$$

$$J_{ij} = R_{ij} e^{-|i-j|/\xi_J}, \qquad \text{when } 1 \le |i-j| < L, \quad (6)$$

$$J_{ijk} = R_{ijk} e^{-|i-k|/\xi_J}, \qquad \text{when } |i-k| = 2, \quad (7)$$

with $\xi_J = 1$ and $R$'s random variables drawn from a normal distribution $N(0, \sigma_R^2 = 1)$. We expect $\xi_J$ and the variance $\sigma_R^2$ to be related in a microscopic model, and find that setting them equal strikes a good balance between two detrimental extremes: setting $\xi_J < \sigma_R^2$ results in a step-like entropic growth because sites at distance $d$ and $d+1$ entangle at distinct time scales, whereas $\xi_J > \sigma_R^2$ leads to excessively overlapping time scales, making the growth process less clear and too short-lived [42].

To be consistent with the $\ell$-bit picture, and to draw direct parallels with studies of Heisenberg-like models, we require that $U$ is quasilocal and that it preserves the $U(1)$ symmetry. The first condition means that the support of $\tau_i^z$ should decay exponentially in the physical basis. The second implies that $U$ conserves the total magnetization (particle number). These requirements are met if we express $U$ as a random unitary circuit with finite depth $d_u$, as depicted in Fig. 1a, where each gate $u_i$ preserves the $U(1)$ symmetry. We use $(L-1) \times d_u$ unique gates $u_i$ defined as unitary operators acting on sites $i$ and $i+1$

$$u_i = e^{-i f w_i M_i}, \quad (8)$$

$$w_i = e^{-2|h_i - h_{i+1}|}, \quad (9)$$

$$M_i = \sum_{p,q=\pm 1} \frac{\theta_{pq}}{4} (1 + p\sigma^z) \otimes (1 + q\sigma^z) \\ + c\sigma^+ \otimes \sigma^- + c^*\sigma^- \otimes \sigma^+, \quad (10)$$

where $f$ is a *mixing factor*, $w_i$ are weights discussed momentarily, and $M_i$ are random Hermitian matrices.
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) An $\ell$-bit $\tau_i^z = U\sigma_i^z U^\dagger$ on $L = 8$ sites with $U$ expressed as a “brickwork” circuit of two-site gates $u_i$ [white boxes, Eq. (8)]. This circuit has depth $d_u = 2$ because $U$ acts twice on each bond $i, i + 1$. The range of the support (shaded gray) grows as $4d_u$. (b) The support of $\tau_{i=10}^z$ on $L = 20$ sites for 25 circuit realizations with fixed $h_i$ (blue bars). This demonstrates the staggered decay caused by $w_i$ (black) where $|h_i - h_{i+1}|$ is large. Without $w_i$ (gray), the decay is smooth. (c) Supports for $\tau_{L/2}^z$ averaged over 10000 circuit realizations for each $L$, including the $w_i$ on and not fixing $h_i$. The parameters $\xi_\tau$ and $\beta$ have been extracted from stretched-exponential fits on all sites $i$. In both (b) and (c) $f = 0.2$ and $d_u = 16$.</figcaption></figure>

The mixing factor controls the exponential length scale $\xi_\tau$ with which $\tau_i^z$ decays (on average) away from site $i$. We set $f = 0.2$ and $d_u = 16$ and use these values unless stated otherwise. This circuit depth allows $\ell$-bits to spread edge-to-edge up to $L = 32$.

In Eq. (10), $\sigma^z$ and $\sigma^\pm = \sigma^x \pm i\sigma^y$ are $2 \times 2$ Pauli matrices and $\theta_{pq}$, Re($c$), Im($c \in \mathcal{N}(0, 1)$) are six normal random variables drawn independently for each gate. The first four $\theta$-terms along the diagonal mix the superposition of different configurations at fixed magnetization on sites $i, i + 1$, while the last two off-diagonal $c$-terms exchange the magnetization of those two sites. The weights $w_i$ in Eqs. (8) and (9) are introduced such that the spread of the $\ell$-bits is suppressed when the neighboring fields $h_i$ from Eq. (5) differ substantially. The effect of $w_i$ is readily seen in Fig. 1b, where we have used the operator overlap $O(i, j) = \text{Tr}(\tau_i^z \sigma_j^z)/2^L$ to measure the support of $\ell$-bits in the physical basis [36, 43, 44]. Upon disorder-averaging, Fig. 1c shows that weights reproduce the expected decay from the $\ell$-bit picture, that is, an exponential ($\beta \approx 1$) with a characteristic length scale $\xi_\tau \approx 0.5$, comparable to previous estimations deep in the many-body localized phase [45]. In contrast, with the weights disabled, we typically get a larger stretching exponent $\beta$ [36]. We attribute the small drift with $L$ of both $\xi_\tau$ and $\beta$ to finite-size effects in the circuit.

The weights can be heuristically motivated by variable-range hopping in the context of Anderson localization [46–48]. In this conduction mechanism, the transition amplitude $T$ between two states localized at distance $d$ with energies $E_a$ and $E_b$ is said to decay exponentially as $T \sim e^{-d/\xi} e^{-|E_a - E_b|/k_B T}$, where $\xi$ is the characteris-

tic length-scale of localization. In our case, the circuit structure models the first exponential. To model the second exponential, we use the fields $h_i$ as a proxy for energies, since these are dominant in $H'$ and correspond to the single-particle energies in the noninteracting limit. The use of weights results in the distributions of saturation values of $S_E$ and $S_N$ becoming qualitatively similar to that of the disordered Heisenberg XXZ model [15]. We verified that our main conclusions, in particular on the growth of number entropy, do not depend on these weights.

We note that the quasilocality of $\ell$-bits manifests in two characteristic length scales $\xi_J$ and $\xi_\tau$ that are determined independently in $H'$ and $U$, respectively. In microscopic models, these two length scales are typically correlated and much smaller than $L$ in the many-body localized regime.

*Numerical results.*—In Fig. 2, we present numerical results from simulated quenches using our model of interacting $\ell$-bits for system sizes up to $L = 28$ and for at least $10^5$ random realizations of $H'$ and $U$. Every simulation is initialized in the physical basis in one of the two antiferromagnetic Néel states ($|0101\cdots1\rangle$ or $|1010\cdots1\rangle$) chosen randomly. Observables are then measured during time evolution at logarithmically spaced time points spanning the entire process of entropic growth, until saturating due to finite system size. The infinite-time saturation values, denoted with a superscript $\infty$, are estimated by time-averaging the corresponding time series after it has fully saturated. Symbols $\triangleright \cdot \Box \cdot \Box$ indicate the time scales for entropic growth and subsequent saturation (see the caption), as discussed further below.

Fig. 2a shows that our approach recovers the slow growth of disorder-averaged entanglement entropy consistent with $\bar{S}_E \sim \xi_J \ln t$, expected in the many-body localized regime [10, 11]. The entropy saturates fully when $t \approx t_\circ$ at values $\bar{S}_E^\infty \sim L$ (inset), indicating a sub-thermal volume-law growth [11, 49].

Our main result, presented in Fig. 2b, demonstrates ultraslow growth of disorder-averaged number entropy consistent with $\bar{S}_N \sim \ln \ln t$. The inset shows that the saturation values $\bar{S}_N^\infty$ (full lines) grow subextensively with $L$, though the functional dependence on $L$ cannot be extracted reliably from our data. The Néel initial state was selected specifically to avoid spurious nonmonotonic growth of $\bar{S}_N^\infty$ with $L$, observed with random initial product states [20, 21]. With random initial states our model shows the same nonmonotonicity, as we elaborate on in the supplemental material [36]. The inset also shows the disorder average of $S[p_{L/2}^\infty(n)] = - \sum_n p_{L/2}^\infty(n) \ln p_{L/2}^\infty(n)$ (dotted line), which is the Shannon entropy of the saturation values of the particle number probabilities at the half-chain, $i = L/2$ [15]. Since the entropy is a concave function of $p_i(n)$, we know from Jensen’s inequality that $\bar{S}[p_{L/2}^\infty(n)]$ bounds $\bar{S}_N$ from above. Following an initial
<figure><img src="image_3.png" /></figure>

FIG. 2. We time evolve the Néel initial state (in the physical basis) using our $\ell$-bit model for system sizes $L = 12$ up to 28. The legend in (a) applies in all panels. The half-chain entropies $S_E, S_N$, and variance $\sigma_{x_n}^2$ are measured at $N_t = 200$ logarithmically spaced time points $t$ for $L \le 20$, and $N_t = 100$ otherwise, spanning $0.1 \le t \le e^{L/\xi_J}$ (rounded up to a power of ten). Disorder-averages $\bar{S}_E, \bar{S}_N$, and $\bar{\sigma}_{x_n}^2$ are from $2 \cdot 10^5$ random realizations of $H'$, U-pairs for $L \le 20$ and $10^5$ for $L > 20$. Saturation values are denoted $S_E^\infty, S_N^\infty$ and $\sigma_{x_n}^{2\infty}$. Symbols indicate when $\bar{S}_E \sim \ln t$ and $\bar{S}_N \sim \ln \ln t$ both begin ($\triangleright$), and when $\bar{S}_N (\Box)$ and $\bar{S}_E (\circ)$ saturate. (a,b,e) Insets show saturation values as functions of $L$ with bars denoting the standard error of the disorder-average. (c) The approach to saturation values are fitted to $y = a - \kappa \ln \ln t$ (dashed) within the marked time intervals. (d) The distributions $p(S_N^\infty)$ are compared to results from Ref. 15 (gray) at fixed $L = 14$ and various disorder strengths $W$. We use 100 bins in $0 \le S_N^\infty \le 1.2$ for the histogram. (e) The disorder-averaged variance of position $\bar{\sigma}_{x_n}^2$ averaged for two central quasiparticles $n = L/4 - 1$ and $L/4$. (f) The saturation times $t_\Box(\bar{S}_E)$ and $t_\circ(\bar{S}_N)$—defined as the time when the running average $S_k = (N_t - k)^{-1} \sum_{x_n}^{N_t} \bar{S}_{E|N}(t_k)$ ceases to grow monotonically—are estimated from 100 bootstrap samples of the disorder realizations.

rise when $L \le 16$, we find that $\bar{S}[p_{L/2}^\infty(n)]$ then decreases slowly until $L = 28$, suggesting that $\bar{S}_N$ could eventually saturate with system size. However, as for $\bar{S}_N$, we cannot extract a functional dependence on $L$ for $\bar{S}[p_{L/2}^\infty(n)]$. Since both $\bar{S}_N$ and $\bar{S}[p_{L/2}^\infty(n)]$ are small averages from long-tailed distributions, these may require more disorder realizations to resolve the $L$-dependence precisely.

In Fig. 2c, the approach of $\bar{S}_N$ to its saturation value $\bar{S}_N^\infty$ strongly suggests that this growth is consistent with $\bar{S}_N \sim \ln \ln t$. The linear fits to $y = a - \kappa \ln \ln t$ (dashed lines) agree convincingly within the time interval $t_\triangleright < t < t_\Box$ with slope $\kappa \approx 10^{-2}$. In contrast, a power law fit $\sim t^{-\alpha}$, argued for in Ref. 21, only agrees for shorter time intervals [36].

In Fig. 2d we compare the distribution of $S_N^\infty$ with results from Ref. 15 (in gray) that studied the disordered isotropic Heisenberg model. Our results are in qualitative agreement with their high-disorder ($W \ge 6$) data. The use of weights in Eq. (9) has been crucial for attaining this agreement [36]. We recover the finite $L$ drift of a peak towards $\ln 2$, though our peaks increase slowly while becoming narrower as $L$ increases. The sharp drop

before $S_N = \ln 3$ has been observed previously, and has been interpreted as a sign that $S_N$ stems mainly from the oscillations of a single particle [15, 20].

In Fig. 2e we plot the disorder-averaged width of quasiparticles at the center of the chain over time, where the width is estimated with the variance of position $\bar{\sigma}_{x_n}^2$—see Eq. (2) and the surrounding discussion. We take the average width of two central quasiparticles, numbered $n = L/4 - 1$ and $n = L/4$, to account for the asymmetry of the particle positions inherent in the Néel states. The growth is consistent with $\bar{\sigma}_{x_n}^2 \sim \ln \ln t$ on a time interval that is commensurate with the growth of $\bar{S}_N$. This time dependence is further supported by the approach to saturation, which admits a linear fit similarly to Fig. 2c, albeit with a smaller slope $\kappa \approx 0.0075$ [36]. The inset shows that saturation values $\bar{\sigma}_{x_n}^{2\infty}$ grow monotonically but subextensively with $L$, similarly to $\bar{S}_N^\infty$. The case is markedly different for the position expectation values $\langle x_n \rangle$, whose disorder averages (and fluctuation thereof) are roughly constant over time, centered at the position given by the initial state [36].

Finally, we address the time scales involved during the growth and saturation of entropies. In the time inter-
val where $\overline{S}_E \sim \ln t$ holds, increasingly distant physical
spins become entangled as they interact through the ex-
ponentially decaying couplings in $H'$, giving rise to the
logarithmic entanglement lightcone [50, 51]. The time
interval is bracketed at $t < t_\Delta \approx 1$ by ballistic spreading
of entanglement within the localization length (approxi-
mately the inverse nearest-neighbor couplings) such that
$\overline{S}_E \sim t$, and much later, by the onset of saturation due
to finite size. Although we cannot determine precisely
the time point when $\overline{S}_E$ begins to saturate, Fig. 2a sug-
gests that it coincides with $t \approx t_\square$, the time when $\overline{S}_N$ has
saturated completely. From Fig. 2f, the time $t_\square \sim e^{L/2}$
is comparable to $t \sim J_{ij}^{-1}$ at distance $|i-j| = L/2$ [see
Eq. (6)]. From the lightcone perspective, this is the time
scale for entanglement to spread over a distance $L/2$,
from the half-chain to the edge, after which finite-size ef-
fects become noticeable. Figs. 2a and 2f show that the
saturation time for $\overline{S}_E$ is consistent with $t_\circ \sim e^L$ which,
following the same argument, is the time scale for entan-
glement to spread edge-to-edge.

In Figs. 2b and 2c, we see that the ultraslow growth
$\overline{S}_N \sim \ln \ln t$ occurs in the time interval $t_\Delta < t < t_\square$.
Similarly to $\overline{S}_E$, we again find ballistic growth $\overline{S}_N \sim t$
in the initial time period $t < t_\Delta \approx 1$. This is not a coin-
cidence, since $S_N$ constitutes the largest contribution to
$S_E$ during this time. In contrast to the lengthy satura-
tion process of $\overline{S}_E$, the ultraslow growth of $\overline{S}_N$ saturates
abruptly when $t \approx t_\square$.

We see in Fig. 2e that the time points $t_\Delta$ and $t_\square$ ap-
ply to the quasiparticle width $\overline{\sigma}_{x_n}^2$ as well. The initial
growth of both entropies can therefore be understood
from the perspective of a quasiparticle next to the half-
chain boundary: starting entirely localized on one site,
the quasiparticle spreads ballistically within the localiza-
tion length such that a small fraction crosses the bound-
ary, thus generating a proportional amount of number
entropy [36]. We observe the simultaneous saturation
of $\overline{\sigma}_{x_n}^2$ and $\overline{S}_N$ at time $t_\square \sim e^{L/2}$, which we take as
yet another indication that the ultraslow growth of both
quantities are directly related.

Discussion.—We have demonstrated that a circuit
model of interacting $\ell$-bits exhibits ultraslow growth in
time of the half-chain number entropy. The observed
growth is consistent with $\overline{S}_N \sim \ln \ln t$ within a time inter-
val that is exponential in system size. Additionally, we
have observed the monotonic and subextensive growth
of saturation values $\overline{S}_N^\infty$ with system size up to $L = 28$.
However, due to the limited range of available system
sizes, the functional form of $\overline{S}_N^\infty(L)$ remains elusive. The
observation that the entropy $\overline{S}[p_{L/2}^\infty(n)]$, which bounds
the number entropy from above, decreases in the range
$L \in [16, 28]$ indicates the possibility that $\overline{S}_N^\infty$ may even-
tually saturate with $L$.

We have studied the quench dynamics from the per-
spective of quasiparticles—particles dressed by particle-
hole excitations—which are characterized by their posi-

tion $\langle x_n \rangle$ and variance $\overline{\sigma}_{x_n}^2$. We discovered a robust con-
nection between the behaviors of number entropy and
quasiparticle variance. Specifically, quasiparticles widen
ultraslowly, with $\overline{\sigma}_{x_n}^2 \sim \ln \ln t$ (for half-chain particles)
simultaneously with $\overline{S}_N \sim \ln \ln t$, while $\langle x_n \rangle$ remains
nearly constant over time [36]. This suggests that the
main contribution to $S_N$ comes from a single particle
oscillating across the half-chain boundary—as noted in
Refs. 15 and 20—and crucially, that the growth of $\overline{S}_N$ is
caused primarily by the ultraslow widening of the central
quasiparticle. Notably, the applicability of the quasipar-
ticle picture is exclusive to many-body localized systems,
as generic states of ergodic systems are extended. This
suggests the possibility for distinct behaviors of $\overline{S}_N$ vs.
$L$ in localized and ergodic systems, irrespective of the
observed ultraslow temporal growth in both cases.

As our model is many-body localized by construction,
the ultraslow growth of number entropy is inherent to
many-body localized systems, at least for the range of
system sizes available to us. An observation of ultra-
slow growth of number entropy in microscopic models
is therefore not, by itself, sufficient evidence that such
systems must delocalize in the thermodynamic limit. In-
stead, it may stem from the truly many-body nature of
the $\ell$-bits. Future studies should explore whether other
slow dynamics traditionally attributed to the absence of
many-body localization can be captured by more compre-
hensive models of many-body localized $\ell$-bit dynamics.

This work received funding from the European Re-
search Council (ERC) under the European Union’s Hori-
zon 2020 research and innovation program (Grant Agree-
ment No. 101001902). Thomas Klein Kvorning’s re-
search is funded by the Wenner-Gren Foundations. Loïc
Herviou has been supported by the Swiss National Sci-
ence Foundation (FM, ZJ) grant 212082. The compu-
tations were enabled by resources provided by the Na-
tional Academic Infrastructure for Supercomputing in
Sweden (NAISS) and the Swedish National Infrastruc-
ture for Computing (SNIC) at Tetralith partially funded
by the Swedish Research Council through grant agree-
ments no. 2022-06725 and no. 2018-05973.

[1] D. Basko, I. Aleiner, and B. Altshuler, Ann. Phys. (N. Y.) **321**, 1126 (2006).

[2] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. **95**, 206603 (2005).

[3] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. **111**, 127201 (2013).

[4] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B **90**, 174202 (2014).

[5] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. **91**, 021001 (2019).

[6] V. Ros, M. Müller, and A. Scardicchio, Nucl. Phys. B **891**, 420 (2015).
[7] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Phys. Rev. Lett. 115, 046603 (2015).

[8] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).

[9] L. D'Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).

[10] M. Žnidarič, T. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).

[11] J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).

[12] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 110, 260601 (2013).

[13] R. Singh, J. H. Bardarson, and F. Pollmann, New J. Phys. 18, 023046 (2016).

[14] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. Lett. 124, 243601 (2020).

[15] D. J. Luitz and Y. B. Lev, Phys. Rev. B 102, 100202 (2020).

[16] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Phys. Rev. B 105, 174205 (2022).

[17] P. J. D. Crowley and A. Chandran, SciPost Phys. 12, 201 (2022).

[18] D. Sels, Phys. Rev. B 106, L020202 (2022).

[19] D. M. Long, P. J. D. Crowley, V. Khemani, and A. Chandran, Phys. Rev. Lett. 131, 106301 (2023).

[20] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. B 103, 024203 (2021).

[21] R. Ghosh and M. Žnidarič, Phys. Rev. B 105, 144203 (2022).

[22] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Comment on “Resonance-induced growth of number entropy in strongly disordered systems” (2022), arXiv:2203.06689.

[23] R. Ghosh and M. Žnidarič, Response to ‘comment on the’ory of growth of number entropy in disordered systems’ (2022), arXiv:2205.05768.

[24] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Phys. Rev. E 102, 062144 (2020).

[25] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Phys. Rev. B 102, 064207 (2020).

[26] D. A. Abanin, J. H. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. A. Parameswaran, F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur, Ann. Phys. (N. Y.) 427, 168415 (2021).

[27] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. Lett. 124, 186601 (2020).

[28] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Phys. Rev. B 98, 174202 (2018).

[29] F. Weiner, F. Evers, and S. Bera, Phys. Rev. B 100, 104204 (2019).

[30] T. L. M. Lezama, S. Bera, and J. H. Bardarson, Phys. Rev. B 99, 161106 (2019).

[31] P. Sierant and J. Zakrzewski, Phys. Rev. B 105, 224203 (2022).

[32] F. Evers, I. Modak, and S. Bera, Phys. Rev. B 108, 134204 (2023).

[33] H. M. Wiseman and J. A. Vaccaro, Phys. Rev. Lett. 91, 097902 (2003).

[34] R. G. Melko, C. M. Herdman, D. Iouchtenko, P.-N. Roy, and A. Del Maestro, Phys. Rev. A 93, 042336 (2016).

[35] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Science 364, 256 (2019).

[36] See Supplemental Material (attached) for details on the numerical implementation, time evolution, choice of initial states, and correspondence between the number entropy and particle width. See also the video lbits.mp4.

[37] V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).

[38] B. Bauer and C. Nayak, J. Stat. Mech. Theory Exp. 2013, P09005 (2013).

[39] D. Pekker and B. K. Clark, Phys. Rev. B 95, 035116 (2017).

[40] U. Schollwöck, Ann. Phys. (N. Y.) 326, 96 (2011).

[41] We have kept $J_{ijk}$ in the model, but our results are reproduced qualitatively by $h_i$ and $J_{ij}$.

[42] M. Žnidarič, Phys. Rev. B 97, 214202 (2018).

[43] A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Phys. Rev. B 91, 085425 (2015).

[44] L. Rademaker, M. Ortuño, and A. M. Somoza, Ann. Phys. (Berl.) 529, 1600322 (2017).

[45] A. K. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Simon, Phys. Rev. B 98, 184201 (2018).

[46] N. Mott, J. Non-Cryst. Solids 1, 1 (1968).

[47] S. Banerjee and E. Altman, Phys. Rev. Lett. 116, 116601 (2016).

[48] S. M. Girvin and K. Yang, Anderson localization, in *Modern Condensed Matter Physics* (Cambridge University Press, 2019) p. 252–300.

[49] A. Nanduri, H. Kim, and D. A. Huse, Phys. Rev. B 90, 064201 (2014).

[50] I. H. Kim, A. Chandran, and D. A. Abanin, (2014), arXiv:1412.3073.

[51] D.-L. Deng, X. Li, J. H. Pixley, Y.-L. Wu, and S. Das Sarma, Phys. Rev. B 95, 024202 (2017).
Supplemental Material for
Ultraslow Growth of Number Entropy in an $\ell$-bit Model of Many-Body Localization

David Aceituno Chávez,¹ Claudia Artiaco,¹ Thomas Klein Kvorning,¹ Loïc Herviou,² and Jens H. Bardarson¹

¹Department of Physics, KTH Royal Institute of Technology, Stockholm, 106 91 Sweden

²Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

S1. QUASILOCALITY OF $\ell$-BITS

The quasilocality of the $\ell$-bit operators $\tau_i^z$ implies that their support on the physical basis should decay, on average, as $\sim e^{-|i-j|/\xi_\tau}$, where $\xi_\tau$ is the $\ell$-bit localization length. Our circuit, which is defined in Eqs. (8-10) in the main text, sets $\xi_\tau$ indirectly through the mixing factor $f$. We estimate $\xi_\tau$ by measuring the decay of $\tau_i^z$ in the physical basis using the operator overlap [1, 2]

$$O(i, j) = \frac{\mathrm{Tr}(\tau_i^z \sigma_j^z)}{2^L}; \qquad (S1)$$

the trace is computed using the circuit as depicted in Fig. S1. Next, we compute the arithmetic disorder average $\overline{O(i, j)}$ over many random realizations of the circuit as well as field configurations $h_i$. Since the $\tau_i^z$ are localized, the matrix $\overline{O(i, j)}$ is peaked along the diagonal $i = j$. By averaging the decay to the left and right of the diagonal entry we obtain $\overline{O(i, |i-j|)}$, and finally, by averaging over $i$, we compute $\overline{O(|i-j|)}$. The localization length $\xi_\tau$ is obtained by fitting to a stretched exponential $|\overline{O}(x)| = Ce^{-(x/\xi_\tau)^\beta}$, where $x = |i-j|$ and $\beta$ is the stretching exponent. Examples of $|\overline{O}(i, j)|$ and fits are shown in Fig. S2 for a range of mixing factors $f$.

The calculation outlined above is used to guide the construction of the unitary circuit. In particular, we use it to ensure that $\beta$ is close to unity, since it would otherwise indicate a departure from the $\ell$-bit picture, where the operators should decay exponentially on average. This is an important motivation for using the field-dependent weights $w_i$ introduced in Eq. (9) in the main text. Indeed, we see in Fig. S2a that $\beta$ remains close to unity for a broad range of $f$ when weights are enabled. On the other hand, in Fig. S2b we see that $\beta$ grows rapidly beyond unity with increasing $f$ when the weights are disabled (that is, $w_i = 1$). The decay of $|\overline{O}(x)|$ also guides the choice of $f$ because of its relation to $\xi_\tau$. If $\xi_\tau$ is too small, we suppress the $\ell$-bit interactions whose dynamics we wish to investigate. On the other hand, if $\xi_\tau$ is too large, entanglement may proliferate to a point where it is no longer commensurate with the physics of localized systems, on top of escalating computational demands. From a previous work, we also expect that $\xi_\tau \lesssim 1 \ll L$ deep in the many-body localized phase [3]. Considering all these criteria, we find that $f = 0.2$ (with weights enabled) gives reasonable values for $\xi_\tau$ and $\beta$.

Finally, we can use $\overline{O(|i-j|)}$ to check that the tails of the $\ell$-bit supports are not overly truncated. This is only a concern in the case of shallow circuits, that is, with $d_u < L/2$,

<figure><img src="image_4.png" /><figcaption>FIG. S1. The operator overlap $O(i, j)$ [Eq. (S1)] measures the support of the $\ell$-bit $\tau_i^z$ at site $j$. The dotted lines show the tracing over free indices that takes place after contracting the circuit with $\sigma_j^z$. Only blue gates give a relevant contribution in this example because white gates contract to identities.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. S2. The disorder-averaged support of $\tau_L^z$, calculated with Eq. (S1), for system size $L = 28$ and circuit depth $d_u = 16$, while varying the mixing factor $f$. In (a) and (b) we compare using and not using, respectively, the field-dependent weights $w_i$ defined in Eq. (9) in the main text. For each choice of parameters, we average 10000 realizations of the random unitary circuit and onsite fields $h_i$. The dashed lines show the fit to a stretched exponential used to calculate $\xi_\tau$ and $\beta$.</figcaption></figure>

because deeper circuits can always accommodate the full diameter of $\ell$-bits on any site. Even though a shallow circuit with carefully chosen parameters could improve the computational efficiency, we have opted to eliminate this source of error entirely by setting $d_u = 16$, which exceeds the longest half-chain considered here.
## S2. TIME EVOLUTION

We are interested in studying the entire dynamical process following the quench of a disordered quantum system modeled by interacting $ℓ$-bits. Our approach is to time-evolve in the diagonal $ℓ$-bit basis and perform measurements in the physical basis, transforming between these bases using a finite-depth random unitary circuit. We use the matrix product states (MPS) ansatz to express quantum states at all times [4–6].

### A. Time step

Because the $ℓ$-bit Hamiltonian $H'$ is diagonal (see Eq. (4) in the main text), the time evolution can be performed in a single step without Trotterization and with an arbitrary time step $t$:

$$|\psi'(t)\rangle = e^{-iH't} |\psi'(0)\rangle. \qquad (S2)$$

More precisely, since all terms in $H'$ commute with each other, the exponential in Eq. (S2) is trivially factorized into operators that act independently on the Hilbert space of the corresponding sites, that is,

$$|\psi'(t)\rangle = \prod_{i} e^{-iH'_{i}t} \prod_{i<j} e^{-iH'_{ij}t} \dots |\psi'(0)\rangle. \qquad (S3)$$

After applying the time-evolution operators, the state is brought back to an MPS form using singular value decomposition (SVD), as done, for example, in the time-evolving block decimation (TEBD) algorithm [7–9]. In fact, except for the lack of Trotterization and arbitrary time step, this part of our algorithm is identical to TEBD.

A single simulation calculates a series of states at increasing time points $|\psi'(t_1)\rangle, |\psi'(t_2)\rangle, \dots, |\psi'(t_{\text{max}})\rangle$. Each time-evolved state $|\psi'(t_n)\rangle$ is obtained independently by applying the time evolution operator to the initial state $|\psi'(0)\rangle$ once. Following each time step, we use $U$ to transform the state back to the physical basis, where observables are measured. Thus, a simulation can be summarized as

$$|\psi(t_n)\rangle = U e^{-iH't_n} U^{\dagger} |\psi_{\text{init}}\rangle, \qquad (S4)$$

where $t_n$ are logarithmically spaced time points from $t_1 = 10^{-1}$ up to $t_{\text{max}} \gg e^{L/\xi_J}$. This choice of $t_{\text{max}}$ exceeds the time scale $J_{ij}^{-1}$ of the smallest interactions in $H'$, at range $|i-j| = L-1$, after which the entropies should have fully saturated.

### B. Long-range terms

We use the method of swap gates [10] to apply time-evolution operators for long-range terms in $H'$. In this method, a swap operation $s_{i,i+1}$ exchanges the physical

<figure><img src="image_6.png" /><figcaption>FIG. S3. A sequence of swap-gates $s_{1,2}, s_{2,3}$, drawn as crossings, moves the state $A_1$ on site 1 to site 3. The time-evolution operator can then act on a two-site MPS tensor $A_1 \otimes A_4$ with spin dimension $d^2$, before swapping the spins back to their original sites. For long-range terms, this operation reduces the memory cost in exchange for swap operations.</figcaption></figure>

indices of two neighboring sites of the MPS, such that indices $i$ and $i+1$ point to sites $i+1$ and $i$, respectively. After exchanging indices, the combined matrix of sites $i$ and $i+1$ is split with an SVD, which dominates the runtime cost of this operation. The core idea, depicted in Fig. S3, is to sandwich each operator $e^{-iH'_{ij}t}$ from Eq. (S3) between two sequences of swap gates. The first sequence moves site $i$ towards $j$, so that $e^{-iH'_{ij}t}$ can act on nearest neighbors (that is, on a subspace of size $2^2$ instead of $2^{j-i+1}$). The second sequence moves site $i$ in reverse to its original position. Note that, since repeated swaps cancel ($s_{i,i+1}^2 = 1$), it is possible to cancel sequences between successive long-range operators. Thus, by re-ordering the product in Eq. (S3), one can cancel most swaps and reduce the computational overhead dramatically.

### C. Precision

Our method avoids two sources of cumulative error found in conventional TEBD: the error due to the Trotter decomposition and the error from repeated truncation of singular values between incremental time steps. Instead, states at different time points have the same upper bound on numerical error because they are fully independent of each other. This means that each time-evolved state is produced by a fixed number of truncations, namely after each site swap operator, after each time evolution operator [Eq. (S3)], and after each two-site gate in the unitary circuit (Eq. (8) in the main text).

When truncating singular values, our strategy is to fix the truncation error threshold $\epsilon$, and let the bond dimension $\chi$ adapt freely to obtain constant precision at all times. More precisely, we require that the truncation error, defined as the norm of the discarded singular values $\lambda_\alpha$, be smaller than a fixed threshold $\epsilon$:

$$\sqrt{\sum_{\alpha=\chi}^{r} \lambda_{\alpha}^{2}} < \epsilon, \qquad (S5)$$
<figure><img src="image_7.png" /><figcaption>FIG. S4. (a) We compute the probability directly with $Q_n = N_n/N$, where N is the total number of random product states of size L (with zero magnetization) and $N_n$ is the subset of states having jump distance n before contributing to $S_N$, as described in the text. (b) We extrapolate the tendency for $L > 10$ to obtain $Q_n^\infty = Q_n(L \to \infty)$.</figcaption></figure>

where r is the rank of the SVD and $\chi$ is chosen to satisfy the inequality on every individual truncation. The value of $\epsilon$ translates to an observed error of the same magnitude in the entropies that we measure at each time point. We have found in trials that $\epsilon = 10^{-5}$ is sufficiently small: smaller values increase computational demand without any discernible differences in our results. Although $\chi$ saturates after some time due to finite system size, the level at which it saturates varies greatly for each disorder realization. Therefore, the main disadvantage of fixing $\epsilon$ is that the memory usage for a given simulation is uncontrolled. We have found empirically that the average bond dimension grows slowly enough with system size that the tradeoff is worthwhile for system sizes below $L \approx 32$, with present computational resources. For the largest system sizes, a small fraction of realizations will outgrow their memory allocation, but these can simply be run again separately with more memory.

Lastly, we remark that since each term in $H'$ is diagonal, one can avoid costly matrix exponentials in Eq. (S3). Instead, one can take the exponential of each diagonal element $x$ directly, as $e^{-i(xt \mod 2\pi)}$, where the modulo prevents numerical overflow in the exponential. We caution that the modulo yields fewer correct digits as $2\pi/xt$ gets closer to floating-point precision. For example, $10^7 \mod 2\pi$ gives 1 or 9 correct digits, respectively, when using 32-bit or 64-bit IEEE 754 floating-point numbers which, in turn, fit approximately 7 or 16 decimal digits of precision, respectively. We use 128-bit “quadruple” precision for all the numbers involved in generating the time evolution operator. Since 128-bit precision fits 34 decimal digits, the error from the modulo is negligible in the largest time scale that we target here.

## S3. INITIAL STATE

We concentrate on initial states $|\psi_{init}\rangle$ with an extensive number of spins available for magnetic exchange,

<figure><img src="image_8.png" /><figcaption>FIG. S5. Perturbation of $q_i(n+1)$ on a Néel-like state.</figcaption></figure>

in the sector of zero total magnetization such that $\sum_{i=1}^L \langle \sigma_i^z \rangle = 0$. This includes random product states (with zero magnetization) and Néel states (when L is even), but not, for instance, states with a single domain-wall in the center. While we observe $\bar{S}_N \sim \ln \ln t$ with both random product states and Néel states, we find with the former that the saturation values $\bar{S}_N^\infty$ grow nonmonotonically with L. This nonmonotonicity has been noticed previously [11, 12]; in Ref. 12 it has been attributed to the fact that an increase in L is accompanied by a change in the probability $Q_n$ that, for a given state, the shortest distance between any particle and a hole on the opposite half-chain will be n. The key idea is that states with shorter distance n will have larger $S_N$, on average. We have depicted $Q_n$ as a function of L in Fig. S4. While in the limit $L \to \infty$ we expect $Q_1 = 0.5$, the probability is larger and varies significantly in the range of system sizes that is accessible to our numerical analysis. Given that $\bar{S}_N^\infty$ is determined largely by how the central quasiparticle crosses the half-chain boundary, we have chosen the Néel initial state (which has $Q_1 = 1$ for all L) to eliminate the spurious nonmonotonicity in saturation values.

## S4. RELATION BETWEEN $S_N$ AND $\sigma_{x_n}^2$

In the following, we establish a relationship between the number entropy $S_N$ and the quasiparticle width $\sigma_{x_n}^2$ that was introduced in the discussion following Eq. (2) in the main text. We show that a perturbation that only widens the central quasiparticle, such that its variance $\sigma^2(x)$ changes by $-\epsilon$, will change $S_N$ by an amount $-\epsilon k + O(\epsilon^2)$, where k is a positive proportionality factor that depends on the number probabilities $p_i(n)$ before the perturbation.

Consider the following perturbation of the central quasiparticle $n+1$

$$ \Delta q_{L/2-1}(n+1) = -\epsilon/2 \quad (S6) $$

$$ \Delta q_{L/2}(n+1) = \epsilon \quad (S7) $$

$$ \Delta q_{L/2+1}(n+1) = -\epsilon/2, \quad (S8) $$
<table><caption>TABLE S1.</caption><thead><tr><th>Δp<sub>i</sub>(m)</th><th>m = n</th><th>m = n + 1</th><th>m = n + 2</th></tr></thead><tbody><tr><td>i = L/2 - 1</td><td>-ε</td><td>ε/2</td><td>ε/2</td></tr><tr><td>i = L/2</td><td>-2ε</td><td>ε</td><td>ε</td></tr><tr><td>i = L/2 + 1</td><td>-3ε/2</td><td>ε/2</td><td>ε</td></tr></tbody></table>

where $\Delta q = \tilde{q} - q$ is the difference between the perturbed and unperturbed values. This perturbation reduces the width of the particle without changing the expectation value of its position $\langle x \rangle$, as shown in Fig. S5. Note that $\tilde{q}_i(n + 1)$ remains normalized, since the changes cancel each other. The change in the expectation value for the position of this quasiparticle is

$$ \begin{aligned} \langle \tilde{x} \rangle - \langle x \rangle &= \sum_{i=0}^{L-1} i \Delta q_i(n+1) \\ &= \left(\frac{L}{2} - 1\right)\left(-\frac{\epsilon}{2}\right) + \frac{L}{2}\epsilon + \left(\frac{L}{2} + 1\right)\left(-\frac{\epsilon}{2}\right) = 0, \end{aligned} \quad (S9) $$

as expected. Given that $\langle \tilde{x} \rangle = \langle x \rangle$, the corresponding change in the variance is

$$ \begin{aligned} \tilde{\sigma}^2(x) - \sigma^2(x) &= \sum_{i=0}^{L-1} (i - \langle x \rangle)^2 \Delta q_i(n+1) \\ &= \sum_{i=0}^{L-1} (i^2 - 2i\langle x \rangle + \langle x \rangle^2) \Delta q_i(n+1), \end{aligned} \quad (S10) $$

where the first term gives

$$ \sum_{i=0}^{L-1} i^2 \Delta q_i(n+1) \quad (S11) $$

$$ = \left(\frac{L}{2} - 1\right)^2 \left(-\frac{\epsilon}{2}\right) + \left(\frac{L}{2}\right)^2 \epsilon + \left(\frac{L}{2} + 1\right)^2 \left(-\frac{\epsilon}{2}\right) \quad (S12) $$

$$ = -\epsilon, \quad (S13) $$

while [using Eq. (S9)] the second term gives

$$ 2\langle x \rangle \sum_{i=0}^{L-1} i \Delta q_i(n+1) = 2\langle x \rangle (\langle \tilde{x} \rangle - \langle x \rangle) = 0. \quad (S14) $$

Finally, the remaining term in Eq. (S10) gives

$$ \langle x \rangle^2 \sum_{i=0}^{L-1} \Delta q_i(n+1) = 0, \quad (S15) $$

where the sum vanishes because the perturbations in Eqs. (S6)-(S8) cancel each other. Thus, the change in variance is

$$ \tilde{\sigma}^2(x) - \sigma^2(x) = -\epsilon, \quad (S16) $$

confirming that the width reduces by $-\epsilon$.

Next, we calculate the change in the half-chain number entropy $\Delta S_N = \tilde{S}_N - S_N$ from the corresponding change

in number probabilities $\Delta p_i(m)$ for $m = n - 1, n$ and $n + 1$. There are many possible configurations $\Delta p_i(m)$ that would satisfy the present perturbation for $q_i(n+1)$, while leaving $q_i(n)$ and $q_i(n + 2)$ unchanged, and maintaining normalization such that $\sum_m p_i(m) = 1 \forall i$ and $\sum_i q_i(m) = 1 \forall m$. Table S1 gives a concrete example, where the normalization is maintained by virtue of each row adding to zero. For the half-chain number entropy, we only need the row for $i = L/2$. For brevity, we introduce the shorthands $p^m = p_{L/2}(m)$, $\tilde{p}^m = \tilde{p}_{L/2}(m)$ and $\delta^m = \Delta p_{L/2}(m) = \tilde{p}^m - p^m$:

$$ \Delta S_N = \sum_m -\tilde{p}^m \ln(\tilde{p}^m) + p^m \ln(p^m) \quad (S17) $$

$$ = \sum_m -p^m \ln \left( 1 + \frac{\delta^m}{p^m} \right) - \delta^m \ln(p^m + \delta^m), \quad (S18) $$

where both logarithms can now be Taylor expanded as $\ln(1+x) = x + O(x^2)$. Keeping the linear terms gives us

$$ \Delta S_N = \sum_m -p^m \frac{\delta^m}{p^m} - \delta^m \left( \ln(p^m) + \frac{\delta^m}{p^m} \right) \quad (S19) $$

$$ = \sum_m -\delta^m - \delta^m \ln(p^m) - \frac{(\delta^m)^2}{p^m}, \quad (S20) $$

The first term above vanishes because the changes cancel by construction. Lastly, we carry out the sum and insert the values of $\delta^m$ from the table

$$ \Delta S_N = -(-2\epsilon) \ln(p^n) - \frac{(-2\epsilon)^2}{p^n} - \epsilon \ln(p^{n+1}) - \frac{\epsilon^2}{p^{n+1}} \quad (S21) $$

$$ = -\epsilon \ln(p^{n+2}) - \frac{\epsilon^2}{p^{n+2}} \quad (S21) $$

$$ = -\epsilon \ln \left( \frac{p^{n+1} p^{n+2}}{(p^n)^2} \right) - \epsilon^2 \left( \frac{4}{p^n} + \frac{1}{p^{n+1}} + \frac{1}{p^{n+2}} \right) \quad (S22) $$

$$ = -\epsilon k - O(\epsilon^2), \quad (S23) $$

where the proportionality factor $k$ is positive under the assumption that the number probability is peaked for particle $n + 1$ at site $L/2$, that is $p^n \approx p^{n+2} \ll p^{n+1}$, which holds in the configuration illustrated in Fig. S5.

## S5. ADDITIONAL DATA

Aided by Fig. S6, we elucidate some remarks made in the main text. In the results section, we have argued that the onset of saturation for $\bar{S}_E$ coincides roughly with the time $t_\square$ for the full saturation of $\bar{S}_N$. Some support for this claim is found in the approach to saturation of $\bar{S}_E$ shown in Fig. S6a, which departs from the linear fit $y \sim \gamma \ln t$ at approximately the same time point as $t_\square$.

In Fig. S6b, we address the suitability of a power-law fit $\sim t^\alpha$ (dotted) to our data. By plotting the approach to saturation of $\bar{S}_N$ in log-log axes, it is clear that $\ln \ln t$
<figure><img src="image_9.png" /></figure>

FIG. S6. All panels are generated from the same dataset as Fig. (2) in the main text, except for the orange line in panel (f). All fits are taken in the time interval $t_{\triangle} < t < t_{\Box}$. (a) The approach to saturation for $\bar{S}_E$ is fitted to $y = a - \gamma \ln t$ (dashed). (b) The approach to saturation for $\bar{S}_N$ is fitted to both $y = a - \kappa \ln \ln t$ (dashed) and $y' = ct^{-\alpha}$ (dotted). (c) The approach to saturation for $\bar{\sigma}_{x_n}^2$ is fitted to $y = a - \kappa \ln \ln t$ (dashed) (d) We use 100 bins in $0 \le \bar{\sigma}_{x_n}^2 \le 1.2$ for the histogram. (e) The offset in the quasiparticle position expectation value from its initial position. For Néel states with $|0\rangle$ at site $i=0$ we subtract the initial position $L/2-1$, otherwise, if $|1\rangle$ is at site $i=0$, we subtract $L/2$. (f) The distribution $p(S_N^\infty)$ produced with and without weights $w_i$, at constant $\xi_\tau \approx 0.5$. The case without weights (orange) includes 10000 disorder realizations. Both are compared to results from Ref. 13 (gray) as in the main text.

(dashed) is a better functional form to describe the tran-
sient regime.

The approach to saturation for $\bar{\sigma}_{x_n}^2$, shown in Fig. S6c,
demonstrates that the growth $\bar{\sigma}_{x_n}^2 \sim \ln \ln t$ occurs within
the same time interval as for $\bar{S}_N \sim \ln \ln t$, namely during
$t_{\triangle} < t < t_{\Box}$, albeit with a smaller slope $\kappa \approx 0.0075$.

In Fig. S6d, we show the distribution of saturation
values $p(\sigma_{x_n}^2 \infty)$ for increasing system size $L$. Both the
peak near $\sigma_{x_n}^2 \approx 0.2$ and the “knee” near $\sigma_{x_n}^2 \approx 0.6$ show
a subtle drift with $L$ towards larger values, which results
in a slowly increasing disorder average. This resembles
the drift seen for $p(S_N^\infty)$ in Fig. (2d) in the main text,
and similarly, it is not clear whether this drift saturates
with $L$.

We have argued in the main text that, while the quasi-
particle width $\bar{\sigma}_{x_n}^2$ grows ultraslowly in time, their posi-
tion expectation values remain fixed near the position

given by the initial state. In Fig. S6e, we show the off-
set $\Delta\langle x_n \rangle$ of the disorder-averaged position of the cen-
tral quasiparticle from its initial position. Notably, since
there is no drift in this quantity over time, it cannot be
responsible for the growth of $\bar{S}_N$.

Finally, we provide an example comparing the distribu-
tion of saturation values $p(S_N^\infty)$ with (green) and without
(orange) weights $w_i$. Note that $f = 0.06$ in the latter,
in order to obtain matching localization lengths $\xi_\tau \approx 0.5$
in both cases. At low values $S_N^\infty \approx 0$, the non-weighted
case resembles the thermal (low $W$) case from Ref. 13.
In addition, it shows no discernible peak near $\ln 2$, and
lacks the tail that follows $S_N > \ln 2$ in the other cases.
All these discrepancies motivate the use of weights in
the circuit, although we have verified that the growth
$S_N \sim \ln \ln t$ is present in both cases.

[1] A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Phys. Rev. B 91, 085425 (2015).

[2] L. Rademaker, M. Ortuño, and A. M. Somoza, Ann. Phys. (Berl.) **529**, 1600322 (2017).

[3] A. K. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Simon, Phys. Rev. B **98**, 184201 (2018).

[4] F. Verstraete and J. I. Cirac, Phys. Rev. B **73**, 094423 (2006).

[5] B. Bauer and C. Nayak, J. Stat. Mech. Theory Exp. **2013**, P09005 (2013).

[6] M. Friesdorf, A. H. Werner, W. Brown, V. B. Scholz, and J. Eisert, Phys. Rev. Lett. **114**, 170505 (2015).

[7] G. Vidal, Phys. Rev. Lett. **91**, 147902 (2003).

[8] G. Vidal, Phys. Rev. Lett. **93**, 040502 (2004).
[9] U. Schollwöck, Ann. Phys. (N. Y.) **326**, 96 (2011).

[10] E. M. Stoudenmire and S. R. White, New J. Phys. **12**, 055026 (2010).

[11] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer,

and J. Sirker, Phys. Rev. B **103**, 024203 (2021).

[12] R. Ghosh and M. Žnidarič, Phys. Rev. B **105**, 144203 (2022).

[13] D. J. Luitz and Y. B. Lev, Phys. Rev. B **102**, 100202 (2020).