# Relaxation in dipolar spin ladders: from pair production to false-vacuum
decay

G. A. Domínguez-Castro  Institut für Theoretische Physik, Leibniz Universität
Hannover, Appelstr. 2, D-30167 Hanover, Germany  Thomas Bilitewski  Department
of Physics, Oklahoma State University, Stillwater, Oklahoma 74078, USA  David
Wellnitz  JILA, National Institute of Standards and Technology and Department
of Physics, University of Colorado, Boulder, CO, 80309, USA  Center for Theory
of Quantum Matter, University of Colorado, Boulder, CO, 80309, USA  Ana Maria
Rey  JILA, National Institute of Standards and Technology and Department of
Physics, University of Colorado, Boulder, CO, 80309, USA  Center for Theory of
Quantum Matter, University of Colorado, Boulder, CO, 80309, USA  Luis Santos
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstr. 2,
D-30167 Hanover, Germany

###### Abstract

Ultracold dipolar particles pinned in optical lattices or tweezers provide an
excellent platform for studying out-of-equilibrium quantum magnetism with
dipole-mediated couplings. Starting with an initial state with spins of
opposite orientation in each of the legs of a ladder lattice, we show that
spin relaxation displays an unexpected dependence on inter-leg distance and
dipole orientation. This intricate dependence, stemming from the interplay
between intra- and inter-leg interactions, results in three distinct dynamical
regimes: (i) ergodic, characterized by the fast relaxation towards equilibrium
of correlated pairs of excitations generated at exponentially fast rates from
the initial state; (ii) metastable, in which the state is quasi-localized in
the initial state and only decays at exceedingly long timescales, resembling
false vacuum decay; and, surprisingly, (iii) partially-relaxed, with
coexisting fast partial relaxation and very long-lived partial quasi-
localization. Realizing these intriguing dynamics is within reach of current
state-of-the-art experiments in dipolar gases.

Quantum simulators using ultracold gases in optical potentials have
dramatically improved our understanding of many-body quantum dynamics, as
highlighted by recent experiments on many-body localization [1, 2, 3], quantum
scars [4], or Hilbert-space fragmentation [5, 6]. Internal states of pinned
particles may be employed to simulate models of quantum magnetism. Up to
recently, nearest-neighbor spin interactions only resulted from relatively
weak super-exchange processes [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. In
contrast, the last years have witnessed rapid progress in the realization of
dipole-mediated spin interactions in various physical systems, including
magnetic atoms [18, 19, 20, 21, 22, 23, 24], Rydberg atoms [25, 26, 27, 28],
and polar molecules [29, 30, 31]. In parallel, the development of optical
tweezers [32, 33, 34, 35] and quantum gas microscopes for polar molecules [36]
opens fascinating perspectives for unveiling the intriguing physics of dipolar
spin models resulting from the long-range anisotropic dipole-dipole
interactions.

The full exploration of non-equilibrium dipolar spin models is further
facilitated by the capability to controllably prepare non-trivial initial spin
distributions. For polar molecules this was demonstrated in recent experiments
where each layer was initialized in opposite spin states (encoded in
rotational molecular levels) [37]. This case is particularly interesting
because dipolar exchange results in the correlated creation of spin
excitations in each layer [38, 39], resembling pair creation from vacuum
fluctuations, or parametric amplification and two-mode squeezing in quantum
optics [40].

Figure 1: (a) Ladder model considered. (b) Instability rate $\Gamma$ for a
ladder of $L=1001$ rungs, as a function of the orientation $\theta$ and the
separation $\Delta$. The colored region indicates where at least one momentum
mode is unstable. The white dashed curve is the effective instability
threshold for $L=11$. (c) Evolution of the imbalance for $\theta=0$ and
$\Delta=2$, $2.5$ and $6$; dashed, dash-dotted, and dotted red curves
correspond to the Bogoliubov predictions for each case, respectively.

In this Letter, we study the non-equilibrium dynamics of a ladder array of
spin-1/2 pinned dipoles (see Fig. 1 (a)). Ladder arrays provide a convenient
and highly relevant system to study non-equilibrium dynamics, as they are
numerically tractable, while still displaying interesting physics, as
illustrated by recent experiments on hydrodynamics in Hubbard models [41]. We
consider spins in each leg initialized in opposite spin orientations, which
may be then subsequently admixed by inter-chain dipolar exchanges. We focus on
whether and how the initial pattern relaxes towards an equilibrium state
characterized by an average zero magnetization in both legs. It might be
expected that relaxation is fastest when the inter-leg dipolar exchange is
maximal, i.e. for dipole moments oriented perpendicular to the ladder
direction ($\theta=0$ in Fig. 1 (a)), and that increasing the inter-chain
separation $\Delta$ only trivially enlarges the relaxation time scale.
Interestingly, none of these a priori reasonable expectations are correct,
illustrating the highly non-trivial nature of the dynamics of dipolar spins.

We show that spin relaxation is maximally favored at the dipole orientation
for which there is no spin exchange along the legs. Moreover, increasing
$\Delta$ does not only enlarge the relaxation time scale, but changes
qualitatively the long-time evolution, leading to three distinct regimes. The
ergodic regime, at low $\Delta$, is characterized by an exponential creation
of correlated pairs of spin excitations at short times, followed by a fast
relaxation towards equilibrium. At large $\Delta$, the system is in a highly
metastable, quasi-localized, regime, remaining in the initial state and only
relaxing in an exceedingly long timescale, resembling false vacuum decay [42].
Interestingly, at intermediate $\Delta$, the system is in a partially-relaxed
regime, with coexisting partial equilibration and very long-lived partial
quasi-localization. Using Bogoliubov theory we provide a microscopic
understanding of the dynamical regimes. Finally, we show that, also counter-
intuitively, positional disorder enhances relaxation.

#### Model.–

We consider two parallel spin-$1/2$ chains of $L$ sites, denoted as A and B.
The chains are characterized by a lattice spacing $a$, set to $a=1$ below, and
are separated by a distance $\Delta$ (Fig. 1 (a)). Due to resonant electric
dipole-dipole interactions, at zero electric field the spins undergo both
intra-leg and inter-leg exchange interactions. The system is then described by
the Hamiltonian:

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=A,B}\sum_{i>j}V_{ij}^{\alpha\alpha}(\hat{s}_{i\alpha}^{+}\hat{s}_{j\alpha}^{-}+\hat{s}_{i\alpha}^{-}\hat{s}_{j\alpha}^{+})$
(1) $\displaystyle+$
$\displaystyle\sum_{i,j}V_{ij}^{AB}(\hat{s}_{iA}^{+}\hat{s}_{jB}^{-}+\hat{s}_{iA}^{-}\hat{s}_{jB}^{+}),$

with $\hat{s}_{i}^{\alpha}=\hat{\sigma}_{i}^{\alpha}/2$ the spin operators in
terms of the Pauli matrices $\hat{\sigma}_{i}^{x,y,z}$ that act on the spin at
site $i$. We focus our attention below on exchange couplings of the form

$V_{ij}^{\alpha\alpha^{\prime}}=\frac{J}{|\mathbf{r}_{i}^{\alpha}-\mathbf{r}_{j}^{\alpha^{\prime}}|^{3}}\left(1-\frac{3[\hat{e}_{d}\cdot(\mathbf{r}_{i}^{\alpha}-\mathbf{r}_{j}^{\alpha^{\prime}})]^{2}}{|\mathbf{r}_{i}^{\alpha}-\mathbf{r}_{j}^{\alpha^{\prime}}|^{2}}\right),$
(2)

where $J$ is the spin-exchange rate (proportional to the dipole moment
squared), $\hat{e}_{d}=\hat{e}_{z}\cos\theta+\hat{e}_{x}\sin\theta$
characterizes the dipole orientation, which is fixed by an external magnetic
field, and $\mathbf{r}_{i}^{\alpha}$ is the position of site $i$ of chain
$\alpha$. Motivated by recent works on polar molecules in bilayers [37], we
consider all spins in A (B) initially in spin $\uparrow$ ($\downarrow$). We
are interested in the stability of this spin pattern, and whether and how the
system reaches equilibrium.

Figure 2: (a) Ratio $P(L/2)/P_{eq}(L/2)$ at a time $\bar{V}t=50$ as a function
of the dipole orientation $\theta$ and the inter-chain separation $\Delta$. In
the ergodic (E) regime (dark blue), $1-\delta<P(L/2)/P_{eq}(L/2)<1$, in the
the partially relaxed (P) regime (light blue),
$\delta<P(L/2)/P_{eq}(L/2)<1-\delta$, and in the quasi-localized (L) regime
(white), $0<P(L/2)/P_{eq}(L/2)<\delta$. In the figure we consider
$\delta=0.05$, but different reasonably small values of $\delta$ do not
significantly alter the results. The dashed red curve corresponds to the
Bogoliubov instability threshold. (b)-(e) Distribution of spin excitations at
four different points in panel (a), indicated by markers. (f)-(i) Evolution
towards thermalization of the distribution of spin excitations for $\Delta=2$
and $\theta=0$. Brown solid curves indicate the equilibrium distribution (5),
whereas dashed curves depict a re-scaled Gaussian characterizing the partially
relaxed regime. All plots are obtained using Chebyshev evolution for ladders
with $11$ rungs.

#### Stability analysis.–

Deviations from the initial condition are best analyzed using the Holstein-
Primakoff transformation
$\hat{s}_{A,i}^{z}=1/2-\hat{a}_{i}^{\dagger}\hat{a}_{i}$,
$\hat{s}_{A,i}^{+}=\hat{a}_{i}$,
$\hat{s}_{B,i}^{z}=-1/2+\hat{b}_{i}^{\dagger}\hat{b}_{i}$,
$\hat{s}_{B,i}^{-}=\hat{b}_{i}$, which maps the initial spin pattern to a
vacuum of spin excitations. The hard-core boson operator
$\hat{a}_{i}^{\dagger}$ ($\hat{b}_{i}^{\dagger}$) creates a spin excitation at
site $i$ in leg A (B). The Hamiltonian (1) transforms into:

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{i>j}\left[V_{ij}^{AA}\hat{a}_{i}^{\dagger}\hat{a}_{j}+V_{ij}^{BB}\hat{b}_{i}^{\dagger}\hat{b}_{j}+\mathrm{H.c.}\right]$
(3) $\displaystyle+$
$\displaystyle\sum_{i,j}V_{ij}^{AB}(\hat{a}_{i}^{\dagger}\hat{b}_{j}^{\dagger}+\hat{a}_{i}\hat{b}_{j}).$

For low density of spin excitations, the hard-core nature of the bosons can be
neglected, and the Hamiltonian may be re-expressed in quasi-momentum space as:

$\hat{H}=\sum_{k}\left[\varepsilon_{k}(\hat{a}_{k}^{\dagger}\hat{a}_{k}+\hat{b}_{k}^{\dagger}\hat{b}_{k})+\Omega_{k}\hat{a}_{k}^{\dagger}\hat{b}_{-k}^{\dagger}+\Omega_{k}^{*}\hat{b}_{-k}\hat{a}_{k}\right],$
(4)

where $\hat{a}_{k}=\frac{1}{\sqrt{L}}\sum_{j}e^{-ikj}\hat{a}_{j}$ and
$\hat{b}_{k}=\frac{1}{\sqrt{L}}\sum_{j}e^{-ikj}\hat{b}_{j}$. Neglecting the
hard-core constraint fails at longer times. Nevertheless, Eq. (4) provides
valuable insights. Inter-chain dipolar coupling results in the correlated
creation of pairs of spin excitations of opposite momentum in each leg, with a
momentum-dependent rate $\Omega_{k}=J\sum_{j}V_{0j}^{AB}e^{-ikj}$, whereas
intra-chain interaction leads to an effective band dispersion for the motion
of spin excitations along the legs, $\varepsilon_{k}=J\sum_{j\neq
0}V_{0j}^{\alpha\alpha}e^{-ikj}=J(1-3\sin^{2}\theta)\sum_{j\neq
0}e^{-ikj}/|j|^{3}$. The Hamiltonian (4) can be diagonalized by means of a
Bogoliubov transformation [43], yielding the eigenenergies
$\xi_{k}=\sqrt{\varepsilon_{k}^{2}-|\Omega_{k}|^{2}}$. Real $\xi_{k}$ lead to
oscillatory dynamics
$n_{k}^{A}(t)=n_{-k}^{B}(t)=(|\Omega_{k}|/\xi_{k})^{2}\sin^{2}(\xi_{k}t)$.
Crucially, $\xi_{k}$ may become imaginary for certain momenta $k_{c}$,
resulting in the dynamical instability of the spin pattern, triggering (at
short times) the exponential growth of correlated spin excitations in both
legs [39],
$n_{k_{c}}^{A}(t)=n_{-k_{c}}^{B}(t)\propto(|\Omega_{k_{c}}|/|\xi_{k_{c}}|)^{2}e^{2\Gamma_{k_{c}}t}$
with $\Gamma_{k_{c}}=\text{Im}[\xi_{k_{c}}]$ the growth rate of mode $k_{c}$.
We define the instability rate as $\Gamma=\max_{k_{c}}\Gamma_{k_{c}}$ (Fig. 1
(b)).

The momentum-dependent interplay between intra-leg dispersion and inter-leg
pair creation is crucial for the stability of the initial spin pattern and, as
shown below, also for the long-time evolution. Bogoliubov instability for a
given quasi-momentum $k$ demands $|\epsilon_{k}|<|\Omega_{k}|$. Note that,
whereas $|\Omega_{k}|$ decreases with the inter-chain separation $\Delta$, the
bandwidth associated with the intra-chain dispersion $\epsilon_{k}$ is
independent of $\Delta$. As a result, except in the vicinity of the magic
angle $\theta_{M}=\arcsin(1/\sqrt{3})$ discussed below, only $k\simeq
k_{0}\simeq 0.46\pi$, such that $\epsilon_{k_{0}}=0$, contribute to the
instability for large-enough $\Delta$. However, for a finite chain with $L$
sites, the quasi-momentum takes only discrete values $\\{k_{j}\\}$. Therefore
for system sizes where $k_{0}$ cannot be reached, an effective threshold
$\Delta_{c}(\theta,L)$ emerges (Fig. 1 (b)), such that for
$\Delta>\Delta_{c}$, $|\epsilon_{k_{j}}|>|\Omega_{k_{j}}|$ for all $k_{j}$,
and hence the initial spin pattern is Bogoliubov stable. The situation is
markedly different for $\theta\simeq\theta_{M}$. At the magic angle,
$\epsilon_{k}=0$ for all $k$, and hence $\xi_{k}$ is imaginary for all quasi-
momenta. Instability is hence strongly enhanced, being maximal for $k=0$, with
$\Gamma\simeq\frac{4J}{3\Delta^{2}}$.

#### Magnetization imbalance and spin distribution.–

We focus at this point on the long-time evolution, and in particular on
whether the system behaves ergodically, reaching equilibrium. To this aim, we
employ a Chebyshev expansion [44, 43], which provides a numerically exact
evolution for arbitrarily long times. Due to its numerical complexity, we
restrict to chains with $L=11$ rungs.

In order to characterize the relaxation dynamics of the initial spin
distribution, we monitor the imbalance $\eta=m_{A}-m_{B}$, where
$m_{\sigma=A,B}=\frac{1}{L}\sum_{i=1}^{L}\langle\hat{s}_{i\sigma}^{z}\rangle$
is the magnetization of chain A (B). Note that the imbalance is maximal,
$\eta=1$, for the initial condition, and reaches $\eta=0$ when the up and down
spins are evenly admixed between the two chains. In terms of the spin
excitations introduced above, $\eta=1-2\langle n_{A}\rangle$ with
$n_{A}=N_{A}/L$ the density of excitations in chain A
($N_{A}=\sum_{i}\langle\hat{a}_{i}^{\dagger}\hat{a}_{i}\rangle$).

A more detailed probe of the relaxation dynamics is provided by the
distribution of spin excitations in each chain. Due to energy conservation
($\langle\hat{H}\rangle=0$ at any time), and the fact that spin-exchange
preserves $m_{A}+m_{B}=0$, we expect that if ergodicity is reached, all states
with a given number of excitations $N_{A}$ along chain A, and $N_{B}=N_{A}$
excitations along chain B, are equally probable, irrespective of how the
excitations are spatially distributed. The probability $P(N_{A})$ that chain A
has magnetization $m_{A}=1/2-N_{A}/L$ would hence acquire the form [43]:

$P_{eq}(N_{A})=\frac{\binom{L}{N_{A}}^{2}}{\binom{2L}{L}}\simeq\frac{2}{\sqrt{\pi
L}}e^{-4(N_{A}-L/2)^{2}/L}.$ (5)

Figure 3: (a) Ratio $P(L/2)/P_{eq}(L/2)$ as a function of $\Delta$ for
$\theta=0$ and two different evolution times. The dashed red line corresponds
to the Bogoliubov stability threshold. Colored regions indicate the ergodic
(E), partially relaxed (P), and quasi-localized (L) regimes characterizing the
long-time evolution. (b)-(d) Imbalance $\eta$ after a time $Jt/\hbar=50$ as a
function of $\theta$ and $\Delta$ for a disorder-free fully filled lattice,
for sparse ($25\%$) random filling, and for a random on-site magnetic field
(with $W/J=0.15$, see text), respectively. For all plots, we employ Chebyshev
evolution for $11$ rungs, except in Fig. (c), where we average over $100$
realizations of $11$ randomly placed spins in each of the chains of a ladder
with $44$ rungs.

#### Relaxation dynamics.–

In order to properly compare the spin dynamics for different $\Delta$, we
gauge out the trivial stretching of the time-scale of the dynamics associated
to a growing inter-chain separation, by introducing the averaged inter-leg
interaction
$\bar{V}=\lvert(1/L)\sum_{i,j}V_{ij}^{AB}\rvert\simeq\frac{2J}{\Delta^{2}}\cos^{2}\theta$
and considering long evolution times, $\bar{V}t\gg 1$ [45]. Figure 1 (c) shows
the evolution of the imbalance for $\theta=0$ and different values of
$\Delta$. For short times, numerical results are in very good agreement with
the Bogoliubov prediction (red curves). At longer times, the hard-core nature
of the spin excitations, neglected in the Bogoliubov analysis, plays a
relevant role. For $\Delta=2$, for which the system is Bogoliubov unstable,
homogenization is quickly reached. For $\Delta=2.5$ and $6$, the initial
pattern is Bogoliubov stable, but beyond Bogoliubov physics eventually results
in relaxation, although with an exceedingly long-lived memory of the initial
condition.

The long-time evolution reveals three distinct relaxation regimes, as
illustrated in Figs. 2 (a–e) for $\bar{V}t=50$. In the ergodic regime, at low-
enough $\Delta$, the initial state relaxes into a fully balanced distribution
($\eta=0$), and $P(N_{A})$ reaches $P_{eq}(N_{A})$ (see Fig. 2 (b)). Figure 2
(a) shows $P(L/2)/P_{eq}(L/2)$ as a function of $\theta$ and $\Delta$;
$P(L/2)/P_{eq}(L/2)\simeq 1$ in the ergodic regime (dark blue region). As
illustrated in Figs. 2(f–i) for $\theta=0$ and $\Delta=2$, $P(N_{A})$ rapidly
approaches equilibrium, on a time scale of few tens of $\hbar/J$, well within
experimental reach. As expected from Bogoliubov analysis, spin relaxation at
short times is characterized by the exponential creation of correlated pairs
of spin excitations, and hence by an approximately exponentially decreasing
$P(N_{A})$ distribution (see Fig. 2(g)), and the formation of a thermo-field
double state [46, 47]. Eventually, further creation of excitations is arrested
by their hard-core nature, and $P(N_{A})$ transitions from an exponential to
the Gaussian $P_{eq}(N_{A})$ already at $Jt/\hbar\simeq 50$ (Fig. 2 (i)).

The long-time behavior presents similar qualitative features as the
instability rate of Fig. 1 (b). The enhanced Bogoliubov instability for
$\theta\simeq\theta_{M}$, translates into ergodic dynamics even for large
$\Delta$. Moreover, as expected from the stability analysis, for large-enough
$\Delta$ the system enters the quasi-localized regime, in which the spin
distribution remains fully imbalanced, with $\eta\simeq 1$ and
$P(N_{A}=0)\simeq 1$ (see Fig. 2 (e)). This regime is hence characterized by
$P(L/2)/P_{eq}(L/2)\simeq 0$ (white region in Fig. 2 (a)). We have checked
that this remains so up to an eventual relaxation at exponentially long times,
resembling the case of a metastable false vacuum [42].

Interestingly, at intermediate $\Delta$ values, the system enters a partially-
relaxed regime. In that regime $P(N_{A})$ acquires a peculiar bimodal
distribution (see Figs. 2(c,d)). Part of the probability reaches a Gaussian-
like dependence around $N_{A}=L/2$, as expected from ergodic evolution,
whereas the rest remains peaked at $N_{A}=0$, as in the quasi-localized
regime. As a result, this regime is characterized by $0<P(L/2)/P_{eq}(L/2)<1$
(light blue area in Fig. 2 (a)). This bimodal distribution, which is
maintained for very long evolution times [43], well over any present or
foreseeable experimental lifetimes, indicates that the excitation spectrum
presents eigenstates with very different degrees of ergodicity.

#### Experimental considerations.–

The abrupt ergodic to non-ergodic crossover as a function of $\Delta$
characteristic of short evolution times may be readily probed using polar
molecules or Rydberg atoms, and monitored using $\eta$ or $P(N_{A})$, as seen
for $\theta=0$ in Figs. 2 (f–i) and Fig. 3 (a). Probing the partially-relaxed
regime demands longer evolutions. For $\theta=0$ and $\Delta=3$, $\bar{V}t=50$
(depicted in Fig. 3 (a)) corresponds to $Jt/\hbar\simeq 200$. Recent
experiments on polar molecules [48], have demonstrated rotational coherence
times well over this time scale.

Experiments, especially those with polar molecules in optical lattices, are
characterized by a sparse lattice filling [29, 49], limited to less than
$25\%$. Sparse filling results in random positions of the spins, which amounts
to spatial disorder of the exchange couplings. Figure 3 (c) shows the
imbalance at $Jt/\hbar=50$ for a $25\%$ filling. Counter-intuitively,
positional disorder may enhance instability and homogenization compared to the
clean case depicted in Fig. 3 (b). Other sources of disorder may be provided
by the differential polarizability and imperfections of the lattice of the
tweezer array. These could result in a random effective local magnetic field,
$h_{j,\alpha}\hat{s}^{z}_{j\alpha}$. In Fig. 3 (d) we consider the effect of
uniformly distributed random $-W/2<h_{j,\alpha}<W/2$, with $W/J=0.15$. The
effect of the random magnetic field is almost negligible.

#### Conclusions.–

Dipolar spin ladders feature highly non-trivial dynamics, resulting from the
interplay between intra- and inter-leg dipolar interactions, which we have
illustrated for the particularly relevant case in which each leg is
initialized in the opposite spin state. The spin dynamics have a
counterintuitive dependence on the dipole orientation and the inter-chain
separation, being maximally ergodic at the magic angle in which the intra-
chain interactions vanish. Moreover, increasing the inter-leg spacing does not
only lead to the expected enlargement of the time scales, but to a qualitative
change in the stability of the initial spin pattern and in the long-term
evolution, which may fall into one of three possible regimes: ergodic, quasi-
localized and, interestingly, also partially relaxed. Surprisingly, the
Bogoliubov quasiparticle approach provides not only a very good description of
the short-time dynamics but also a qualitative understanding of the long-time
relaxation regimes. The predicted intriguing dynamical regimes shown may be
probed in current and near-future experiments, even in the presence of
positional or on-site disorder.

###### Acknowledgements.

L.S. and G.A.D.-C. acknowledge support of the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy –
EXC-2123 QuantumFrontiers – 390837967. A.M.R acknowledges support from the
AFOSR MURI, ARO W911NF-19-1-0210, NSF JILA-PFC PHY-2317149, NSF QLCI-2016244,
the DOE Quantum Systems Accelerator (QSA) and NIST.

## References

  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Kondov  _et al._ [2015] S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, Disorder-induced localization in a strongly correlated atomic hubbard gas, Phys. Rev. Lett. 114, 083002 (2015).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nature Physics 12, 907 (2016).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Scherg  _et al._ [2021] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted fermi-hubbard chains, Nature Communications 12, 4490 (2021).
  * Morong  _et al._ [2021] W. Morong, F. Liu, P. Becker, K. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. Gorshkov, and C. Monroe, Observation of stark many-body localization without disorder, Nature 599, 393 (2021).
  * Trotzky  _et al._ [2008] S. Trotzky, P. Cheinet, S. Folling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch, Time-resolved observation and control of superexchange interactions with ultracold atoms in optical lattices, Science 319, 295 (2008).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, A cold-atom fermi–hubbard antiferromagnet, Nature 545, 462 (2017).
  * Brown  _et al._ [2015] R. Brown, R. Wyllie, S. Koller, E. Goldschmidt, M. Foss-Feig, and J. Porto, Two-dimensional superexchange-mediated magnetization dynamics in an optical lattice, Science 348, 540 (2015).
  * Greif  _et al._ [2015] D. Greif, G. Jotzu, M. Messer, R. Desbuquois, and T. Esslinger, Formation and dynamics of antiferromagnetic correlations in tunable optical lattices, Phys. Rev. Lett. 115, 260401 (2015).
  * Hart  _et al._ [2015] R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse, and R. G. Hulet, Observation of antiferromagnetic correlations in the hubbard model with ultracold atoms, Nature 519, 211 (2015).
  * Drewes  _et al._ [2017] J. H. Drewes, L. A. Miller, E. Cocchi, C. F. Chan, N. Wurz, M. Gall, D. Pertot, F. Brennecke, and M. Köhl, Antiferromagnetic correlations in two-dimensional fermionic mott-insulating and metallic phases, Phys. Rev. Lett. 118, 170401 (2017).
  * Brown  _et al._ [2017] P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß, S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A. Huse, and W. S. Bakr, Spin-imbalance in a 2d fermi-hubbard system, Science 357, 1385 (2017).
  * Dimitrova  _et al._ [2020] I. Dimitrova, N. Jepsen, A. Buyskikh, A. Venegas-Gomez, J. Amato-Grill, A. Daley, and W. Ketterle, Enhanced superexchange in a tilted mott insulator, Phys. Rev. Lett. 124, 043204 (2020).
  * Nichols  _et al._ [2019] M. A. Nichols, L. W. Cheuk, M. Okan, T. R. Hartke, E. Mendez, T. Senthil, E. Khatami, H. Zhang, and M. W. Zwierlein, Spin transport in a mott insulator of ultracold fermions, Science 363, 383 (2019).
  * Lu  _et al._ [2021] R. Lu, H. Sun, S. Kumar, Y. Wang, M. Gu, M. Zeng, Y.-J. Hao, J. Li, J. Shao, X.-M. Ma, Z. Hao, K. Zhang, W. Mansuer, J. Mei, Y. Zhao, C. Liu, K. Deng, W. Huang, B. Shen, K. Shimada, E. F. Schwier, C. Liu, Q. Liu, and C. Chen, Half-magnetic topological insulator with magnetization-induced dirac gap at a selected surface, Phys. Rev. X 11, 011039 (2021).
  * Hirthe  _et al._ [2023] S. Hirthe, T. Chalopin, D. Bourgund, P. Bojović, A. Bohrdt, E. Demler, F. Grusdt, I. Bloch, and T. A. Hilker, Magnetically mediated hole pairing in fermionic ladders of ultracold atoms, Nature 613, 463 (2023).
  * de Paz  _et al._ [2013] A. de Paz, A. Sharma, A. Chotia, E. Maréchal, J. H. Huckans, P. Pedri, L. Santos, O. Gorceix, L. Vernac, and B. Laburthe-Tolra, Nonequilibrium quantum magnetism in a dipolar lattice gas, Phys. Rev. Lett. 111, 185305 (2013).
  * Chomaz  _et al._ [2023] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Dipolar physics: a review of experiments with magnetic quantum gases, Reports on Progress in Physics 86, 026401 (2023).
  * Alaoui  _et al._ [2022] Y. A. Alaoui, B. Zhu, S. R. Muleady, W. Dubosclard, T. Roscilde, A. M. Rey, B. Laburthe-Tolra, and L. Vernac, Measuring correlations from the collective spin fluctuations of a large ensemble of lattice-trapped dipolar spin-3 atoms, Phys. Rev. Lett. 129, 023401 (2022).
  * Gabardos  _et al._ [2020] L. Gabardos, B. Zhu, S. Lepoutre, A. M. Rey, B. Laburthe-Tolra, and L. Vernac, Relaxation of the collective magnetization of a dense 3d array of interacting dipolar $s=3$ atoms, Phys. Rev. Lett. 125, 143401 (2020).
  * Lepoutre  _et al._ [2019] S. Lepoutre, J. Schachenmayer, L. Gabardos, B. Zhu, B. Naylor, E. Maréchal, O. Gorceix, A. M. Rey, L. Vernac, and B. Laburthe-Tolra, Out-of-equilibrium quantum magnetism and thermalization in a spin-3 many-body dipolar lattice system, Nature communications 10, 1714 (2019).
  * Su  _et al._ [2023] L. Su, A. Douglas, M. Szurek, R. Groth, S. F. Ozturk, A. Krahn, A. H. Hébert, G. A. Phelps, S. Ebadi, S. Dickerson, F. Ferlaino, O. Marković, and M. Greiner, Dipolar quantum solids emerging in a hubbard quantum simulator, Nature 622, 724 (2023).
  * Patscheider  _et al._ [2020] A. Patscheider, B. Zhu, L. Chomaz, D. Petter, S. Baier, A.-M. Rey, F. Ferlaino, and M. J. Mark, Controlling dipolar exchange interactions in a dense three-dimensional array of large-spin fermions, Phys. Rev. Res. 2, 023050 (2020).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nat. Phys. 16, 132 (2020).
  * Scholl  _et al._ [2022] P. Scholl, H. J. Williams, G. Bornet, F. Wallner, D. Barredo, L. Henriet, A. Signoles, C. Hainaut, T. Franz, S. Geier, A. Tebben, A. Salzinger, G. Zürn, T. Lahaye, M. Weidemüller, and A. Browaeys, Microwave engineering of programmable $xxz$ hamiltonians in arrays of rydberg atoms, PRX Quantum 3, 020303 (2022).
  * Chew  _et al._ [2022] Y. Chew, T. Tomita, T. P. Mahesh, S. Sugawa, S. de Léséleuc, and K. Ohmori, Ultrafast energy exchange between two single Rydberg atoms on a nanosecond timescale, Nat. Photonics 16, 724 (2022).
  * Signoles  _et al._ [2021] A. Signoles, T. Franz, R. Ferracini Alves, M. Gärttner, S. Whitlock, G. Zürn, and M. Weidemüller, Glassy Dynamics in a Disordered Heisenberg Quantum Spin System, Phys. Rev. X 11, 011011 (2021).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Observation of dipolar spin-exchange interactions with lattice-confined polar molecules, Nature 501, 521 (2013).
  * Bohn  _et al._ [2017] J. L. Bohn, A. M. Rey, and J. Ye, Cold molecules: Progress in quantum engineering of chemistry and quantum matter, Science 357, 1002 (2017).
  * Moses  _et al._ [2017] S. A. Moses, J. P. Covey, M. T. Miecnikowski, D. S. Jin, and J. Ye, New frontiers for quantum gases of polar molecules, Nat. Phys. 13, 13 (2017).
  * Kaufman and Ni [2021] A. M. Kaufman and K.-K. Ni, Quantum science with optical tweezer arrays of ultracold atoms and molecules, Nature Physics 17, 1324 (2021).
  * Bao  _et al._ [2022] Y. Bao, S. S. Yu, L. Anderegg, E. Chae, W. Ketterle, K.-K. Ni, and J. M. Doyle, Dipolar spin-exchange and entanglement between molecules in an optical tweezer array (2022), arXiv:2211.09780 [physics.atom-ph] .
  * Holland  _et al._ [2022] C. M. Holland, Y. Lu, and L. W. Cheuk, On-demand entanglement of molecules in a reconfigurable optical tweezer array (2022), arXiv:2210.06309 [cond-mat.quant-gas] .
  * Zhang  _et al._ [2022] J. T. Zhang, L. R. B. Picard, W. B. Cairncross, K. Wang, Y. Yu, F. Fang, and K.-K. Ni, An optical tweezer array of ground-state polar molecules, Quantum Sci. Technol. 7, 035006 (2022).
  * Christakis  _et al._ [2023a] L. Christakis, J. S. Rosenberg, R. Raj, S. Chi, A. Morningstar, D. A. Huse, Z. Z. Yan, and W. S. Bakr, Probing site-resolved correlations in a spin system of ultracold molecules, Nature 614, 64 (2023a).
  * Tobias  _et al._ [2022] W. G. Tobias, K. Matsuda, J.-R. Li, C. Miller, A. N. Carroll, T. Bilitewski, A. M. Rey, and J. Ye, Reactions between layer-resolved molecules mediated by dipolar spin exchange, Science 375, 1299 (2022).
  * Bilitewski and Rey [2023] T. Bilitewski and A. M. Rey, Manipulating growth and propagation of correlations in dipolar multilayers: From pair production to bosonic kitaev models, Phys. Rev. Lett. 131, 053001 (2023).
  * Bilitewski  _et al._ [2023] T. Bilitewski, G. A. Domínguez-Castro, D. Wellnitz, A. M. Rey, and L. Santos, Tunable momentum pair creation of spin excitations in dipolar bilayers, Phys. Rev. A 108, 013313 (2023).
  * Agarwal [2013] G. Agarwal,  _Quantum Optics_, Quantum Optics (Cambridge University Press, 2013).
  * Wienand  _et al._ [2023] J. F. Wienand, S. Karch, A. Impertro, C. Schweizer, E. McCulloch, R. Vasseur, S. Gopalakrishnan, M. Aidelsburger, and I. Bloch, Emergence of fluctuating hydrodynamics in chaotic quantum systems (2023), arXiv:2306.11457 [cond-mat.quant-gas] .
  * Coleman [1977] S. Coleman, Fate of the false vacuum: Semiclassical theory, Phys. Rev. D 15, 2929 (1977).
  * [43] See the Supplemental Material (at the url provided by the publisher) for a detailed discussion of the Bogoliubov analysis, the Chebyshev method, the equilibrium distribution, and extended numerical results.
  * Fehske  _et al._ [2007] H. Fehske, R. Schneider, and A. Weisse,  _Computational many-particle physics_, Vol. 739 (Springer, 2007).
  * [45] In order to avoid very long time scales resulting from the vanishing angular dependence of $\bar{V}$, we restrict the long-time analysis to orientations with $0<\theta<0.9\pi$.
  * Takahashi and Umezawa [1996] Y. Takahashi and H. Umezawa, Thermo field dynamics, International Journal of Modern Physics B 10, 1755 (1996).
  * Chapman  _et al._ [2019] S. Chapman, J. Eisert, L. Hackl, M. P. Heller, R. Jefferson, H. Marrochio, and R. C. Myers, Complexity and entanglement for thermofield double states, SciPost Physics 6 (2019).
  * Gregory  _et al._ [2023] P. D. Gregory, L. M. Fernley, A. L. Tao, S. L. Bromley, J. Stepp, Z. Zhang, S. Kotochigova, K. R. A. Hazzard, and S. L. Cornish, Second-scale rotational coherence and dipolar interactions in a gas of ultracold polar molecules (2023), arXiv:2306.02991 [physics.atom-ph] .
  * Christakis  _et al._ [2023b] L. Christakis, J. S. Rosenberg, R. Raj, S. Chi, A. Morningstar, D. A. Huse, Z. Z. Yan, and W. S. Bakr, Probing site-resolved correlations in a spin system of ultracold molecules, Nature 614, 64 (2023b).

##  Appendix A Supplementary Information

This supplementary information contains additional details on the Bogoliubov
analysis, the Chebyshev method, the equilibrium distribution, and extended
numerical results.

###  A.1 Bogoliubov Analysis

In this section, we provide further details on the Bogoliubov analysis. We
write the Hamiltonian of the two-chain system in the quasi-momentum
representation

$\begin{split}\hat{H}&=\sum_{k}\varepsilon_{k}(\hat{a}_{k}^{\dagger}\hat{a}_{k}+\hat{b}_{k}^{\dagger}\hat{b}_{k})+\\\
&\sum_{k}[|\Omega_{k}|e^{-i\alpha_{k}}\hat{a}_{k}^{\dagger}\hat{b}_{-k}^{\dagger}+|\Omega_{k}|e^{i\alpha_{k}}\hat{b}_{-k}\hat{a}_{k}],\end{split}$
(S1)

where we have made explicit the complex nature of the inter-chain coupling
$\Omega_{k}=|\Omega_{k}|e^{-i\alpha_{k}}$. The Hamiltonian can be diagonalized
by means of a Bogoliubov transformation

$\begin{split}\hat{\gamma}_{k}&=u_{k}\hat{a}_{k}-v_{k}\hat{b}_{-k}^{\dagger},\\\
\hat{\beta}_{-k}^{\dagger}&=u_{k}^{*}\hat{b}_{-k}^{\dagger}-v_{k}^{*}\hat{a}_{k},\end{split}$
(S2)

where $\hat{\gamma}_{k}$ and $\hat{\beta}_{-k}^{\dagger}$ are the bosonic
operators associated with the Bogoliubov excitations, and $u_{k}$ and $v_{k}$
obey the Bogoliubov-de-Gennes equations

$\begin{split}\xi_{k}u_{k}&=\varepsilon_{k}u_{k}+|\Omega_{k}|e^{-i\alpha_{k}}v_{k}^{*},\\\
\xi_{k}v_{k}&=-|\Omega_{k}|e^{i\alpha_{k}}u_{k}-\varepsilon_{k}v_{k}^{*},\end{split}$
(S3)

where $\xi_{k}=\sqrt{\varepsilon_{k}^{2}-|\Omega_{k}|^{2}}$ is the Bogoliubov
spectrum. Real and imaginary Bogoliubov eigenenergies give rise to different
excitation dynamics. In the case of real eigenvalues, the time dependence of
the Bogoliubov operators is
$\hat{\gamma}_{k}(t)=e^{-i\xi_{k}t}\hat{\gamma}_{k}(0)$ and
$\hat{\beta}_{-k}^{\dagger}(t)=e^{i\xi_{k}t}\hat{\beta}_{-k}^{\dagger}(0)$. By
inverting the Bogoliubov transformation, one can obtain the time-dependence of
the bosonic operators

$\begin{split}\hat{a}_{k}(t)&=\hat{a}_{k}\left(\cos\xi_{k}t-\frac{i\varepsilon_{k}}{\xi_{k}}\sin\xi_{k}t\right)\\\
&-i\hat{b}_{-k}^{\dagger}\left(\frac{e^{i\alpha_{k}}|\Omega_{k}|}{\xi_{k}}\right)\sin\xi_{k}t,\\\
\hat{b}_{-k}^{\dagger}(t)&=\hat{b}_{-k}^{\dagger}\left(\cos\xi_{k}t+\frac{i\varepsilon_{k}}{\xi_{k}}\sin\xi_{k}t\right)\\\
&+i\hat{a}_{k}\left(\frac{e^{-i\alpha_{k}}|\Omega_{k}|}{\xi_{k}}\right)\sin\xi_{k}t.\end{split}$
(S4)

With the aid of the above equations, the vacuum expectation value of the
population of spin excitations in mode $k$ of the leg A is
$n_{k}^{A}(t)=\langle
0|(\hat{a}_{k}^{\dagger}\hat{a}_{k})(t)|0\rangle=[|\Omega_{k}|\sin(\xi_{k}t)/\xi_{k}]^{2}$.
Integration of $n_{k}^{A}(t)$ over the first Brillouin zone gives the total
density of excitations

$n^{A}(t)=\int_{BZ}\frac{dk}{2\pi a}\
|\Omega_{k}|^{2}\left[\frac{\sin|\xi_{k}|t}{|\xi_{k}|}\right]^{2}.$ (S5)

For Bogoliubov modes with imaginary eigenenergies, the time-dependence of the
bosonic operators is

$\begin{split}\hat{a}_{k}(t)&=\hat{a}_{k}\left(\cosh|\xi_{k}|t-\frac{i\varepsilon_{k}}{|\xi_{k}|}\sinh|\xi_{k}|t\right),\\\
&-i\hat{b}_{-k}^{\dagger}\left(\frac{e^{i\alpha_{k}}|\Omega_{k}|}{|\xi_{k}|}\right)\sinh|\xi_{k}|t\\\
\hat{b}_{-k}^{\dagger}(t)&=\hat{b}_{-k}^{\dagger}\left(\cosh|\xi_{k}|t+\frac{i\varepsilon_{k}}{|\xi_{k}|}\sin|\xi_{k}|t\right)\\\
&+i\hat{a}_{k}\left(\frac{e^{-i\alpha_{k}}|\Omega_{k}|}{|\xi_{k}|}\right)\sin|\xi_{k}|t.\end{split}$
(S6)

As mentioned in the main text, imaginary eigenvalues yield exponential growth
of spin excitations, that is
$n_{k}^{A}(t)=[|\Omega_{k}|\sinh(|\xi_{k}|t)/|\xi_{k}|]^{2}\propto(|\Omega_{k}|/|\xi_{k}|)^{2}e^{2|\xi_{k}|t}=(|\Omega_{k}|/|\xi_{k}|)^{2}e^{2\Gamma_{k}t}$,
where $\Gamma_{k}$ is the rate of the exponential growth. Since
$\sin(i|\xi_{k_{c}}|)/i|\xi_{k_{c}}|\rightarrow-\sinh(|\xi_{k_{c}}|)/|\xi_{k_{c}}|$,
one can safely use the expression in Eq. (S5) to obtain the time dependence of
the density of excitations in each leg.

##  Appendix B Chebyshev expansion

In the following, we briefly introduce the Chebyshev algorithm used to
calculate the time evolution of the spin ladder. A comprehensive description
of the method can be found in Ref. [44]. Considering that at time $t=0$ a
given quantum system is prepared in the state $|\psi_{0}\rangle$, its state at
a later time is given by $|\psi_{t}\rangle=e^{-i\hat{H}t}|\psi_{0}\rangle$
($\hbar=1$). Thus, by calculating an approximate form of the time evolution
operator $\hat{U}(t)=e^{-i\hat{H}t}$, one can estimate the evolution of the
quantum system. Such an approximation is given by the Chebyshev expansion

$U(t)=e^{-i(a\tilde{H}+b)t}\approx
e^{-ibt}\left(c_{0}+2\sum_{\nu=1}^{N}c_{\nu}T_{\nu}(\tilde{H})\right),$ (S7)

where $a=(E_{max}-E_{min})/2$, $b=(E_{max}+E_{min})/2$ with $E_{max}$
($E_{min}$) the highest (lowest) eigenstate energy of the Hamiltonian
$\hat{H}$, $\tilde{H}=(\hat{H}-b)/a$ is a rescaling of the original
Hamiltonian, and $T_{\nu}(x)$ is the Chebyshev polynomial of order $\nu$. The
coefficients $c_{\nu}$ of the Chebyshev expansion are given by
$c_{\nu}=(-i)^{\nu}J_{\nu}(at)$, where $J_{\nu}(x)$ denotes the Bessel
function of order $\nu$. The order $N$ of the Chebyshev expansion is fixed by
demanding a unitary time evolution.

Figure S1: Evolution of the distribution of spin excitations in chain A for
different dipole orientations. The brown curve depicts the expected thermal
distribution. Note the vastly longer time scales in the last two rows. Results
obtained from Chebyshev calculations for $11$ rungs and $\Delta=2$. Figure S2:
Evolution towards thermalization of the distribution of spin excitations in
chain A for $\theta=0$ and different separations of the chains. The brown
curve is associated with the expected thermal distribution. Note again the
much longer time scales in the last two rows. Results obtained from Chebyshev
calculations for $11$ rungs.

##  Appendix C Equilibrium distribution

In this section, we provide further details concerning the equilibrium
distribution that characterizes the ergodic regime. Following Ref. [22], the
long-time behavior of local observables, such as the leg magnetization
$m_{A}=\sum_{i=1}^{L}\hat{s}_{i,A}$, can be described by a thermal
distribution of the form

$\hat{\rho}_{th}=\frac{1}{Z}\text{exp}\left(-\beta\hat{H}-\mu_{1}\hat{S}_{z}-\mu_{2}(\hat{S}_{z}^{2}-\langle\hat{S}_{z}\rangle^{2})\right),$
(S8)

where

$Z=\text{Tr}\left[\text{exp}\left(-\beta\hat{H}-\mu_{1}\hat{S}_{z}-\mu_{2}(\hat{S}_{z}^{2}-\langle\hat{S}_{z}\rangle^{2})\right)\right]$
(S9)

is the partition function, $\beta$, $\mu_{1}$, and $\mu_{2}$ are the Lagrange
multipliers associated with the inverse temperature, total magnetization, and
the variance of the magnetization, respectively. Since in the two-chain system
spin flips are created in pairs, one per chain, the matrix elements of
$\hat{S}_{z}$ and $\hat{S}_{z}^{2}$ are zero. Furthermore, it is
straightforward to notice that the diagonal matrix elements of $\hat{H}$ are
also zero. Therefore, the density matrix is proportional to the identity
operator $\hat{I}$ and the system can be described using the microcanonical
ensemble, where each microstate has equal probability. In this scenario, the
partition function is simply given by the Hilbert space dimension

$Z=\binom{2L}{L}=\frac{(2L)!}{L!^{2}}.$ (S10)

The corresponding density matrix is $\hat{\rho}_{th}=\frac{\hat{I}}{Z}$.
Tracing over one chain we obtain the reduced density matrix. Since the matrix
elements of $\hat{\rho}^{th}$ are all equal, this is a straightforward task,
resulting in

$\rho^{th}_{A}(N_{A})=\frac{\binom{L}{N_{A}}^{2}}{\binom{2L}{L}}.$ (S11)

Using the Stirling approximation for large $L$ results in the Gaussian
distribution discussed in the main text.

##  Appendix D Extended results

In this section, we provide extended results for the time dependence of the
distribution of spin excitations. In Fig. S1, we show the evolution of the
distribution of spin excitations in chain A for different dipole orientations
and fixed separation of the chains $\Delta=2$. As one can readily observe,
regardless of the angle $\theta$, the distribution of excitations attains an
equilibrium state consistent with the thermal prediction. It is important to
note the role of the dipole orientation in the time scale to reach the thermal
state, while for $\theta=\pi/8$ and $\pi/4$ excitations are created rapidly,
bringing the system to a homogeneous state, for $\theta=3\pi/8$ and $\pi/2$
there is considerable overlap with the initial state up to very long time
scales.

To show the role of the separation of the chains on the dynamics of the
distribution of excitations, in Fig. S2, we plot $\rho_{A}(n)$ for three
different values of $\Delta$. As one can notice, thermalization is reached for
the three different values of the separations of the chains considered,
however, the time scale required to enter such steady-state increases
enormously. Note that cases well inside the partially-relaxed regime, e.g. the
lower panels of Fig. S1 and S2, keep the bimodal distribution for very long
times. The eventual relaxation towards the Gaussian distribution demands
evolution times typically well above the longest rotational coherent times
realized in experiments.