# Dark Matter Freeze-out via Catalyzed Annihilation

Chuan-Yang Xing<EMAIL_ADDRESS>Department of Physics and State Key
Laboratory of Nuclear Physics and Technology, Peking University, Beijing
100871, China  Shou-hua Zhu<EMAIL_ADDRESS>Department of Physics and State
Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing
100871, China  Collaborative Innovation Center of Quantum Matter, Beijing,
100871, China  Center for High Energy Physics, Peking University, Beijing
100871, China

###### Abstract

We present a new paradigm of dark matter freeze-out, where the annihilation of
dark matter particles is catalyzed. We discuss in detail the regime that the
depletion of dark matter proceeds via $2\chi\to 2A^{\prime}$ and
$3A^{\prime}\to 2\chi$ processes, in which $\chi$ and $A^{\prime}$ denote dark
matter and the catalyst respectively. In this regime, the dark matter number
density is depleted polynomially rather than exponentially (Boltzmann
suppression) as in classical WIMPs and SIMPs. The paradigm applies for a
secluded weakly interacting dark sector with a dark matter in the MeV-TeV mass
range. The catalyzed annihilation paradigm is compatible with CMB and BBN
constraints, with enhanced indirect detection signals.

##  I Introduction

The existence of dark matter (DM) is well established with ample evidence from
cosmological and astrophysical observations Bauer and Plehn (2019). Though,
the nature of dark matter is still unknown. To solve this puzzle, tremendous
efforts have been devoted to searching for dark matter candidates and studying
production mechanisms of dark matter in the early universe. Among all the
mechanisms that reproduce the observed abundance of dark matter, the
possibility of thermal dark matter where dark matter keeps in thermal
equilibrium with Standard Model (SM) particles in the early universe is
especially popular and compelling.

For massive thermal dark matter, DM particles remain in thermal and chemical
equilibrium while relativistic. As the universe cools down, DM particles are
depleted via some certain processes and the abundance of dark matter goes
down. These processes freeze-out when their interaction rate falls below the
expansion rate of the universe, and consequently DM abundance is settled.
There are essentially two kinds of process leading to depletion of DM
particles in the literature. The first one is that DM particles annihilate
into other particles, mostly SM particles. The other one is via number-
changing process of dark matter. For the former case, the most studied
scenario is self-annihilation process Lee and Weinberg (1977), e.g.
$2\text{DM}\to 2\text{SM}$. Especially, weakly interacting massive particles
(WIMPs) that naturally reproduce correct relic abundance attracted extensive
attentions Bertone  _et al._ (2005); Arcadi  _et al._ (2018); Roszkowski  _et
al._ (2018). Other variations on the self-annihilation case include co-
annihilation Griest and Seckel (1991); Ellis  _et al._ (1998); D’Agnolo  _et
al._ (2018), semi-annihilation D’Eramo and Thaler (2010) and so on Griest and
Seckel (1991); D’Agnolo and Ruderman (2015); Pospelov  _et al._ (2008);
Arkani-Hamed  _et al._ (2009); Feng and Kumar (2008); Belanger and Park
(2012); Dror  _et al._ (2016); D’Agnolo  _et al._ (2017); Garny  _et al._
(2019); Maity and Ray (2020); Kramer  _et al._ (2020); Dey  _et al._ (2017);
Cline  _et al._ (2017). Whereas, for the number-changing process, the most
studied process is $3\text{DM}\to 2\text{DM}$ annihilation, dubbed as strongly
interaction massive particles (SIMPs) Hochberg  _et al._ (2014, 2015).
Subsequently, other number-changing process are proposed and discussed,
including $Z_{2}$-symmetric SIMPs Bernal and Chu (2016); Bernal  _et al._
(2017), co-SIMPs Smirnov and Beacom (2020), etc.

In this Letter, we propose a new pattern of dark matter burning in the early
universe beside the two aforementioned kinds of process, where the abundance
of dark matter is determined by _catalyzed_ processes. In catalyzed processes,
there are some other particles beside DM that act as the catalyst. The
catalyst can _enhance the rate of dark matter burning_ , yet the catalyst
itself is _not consumed in the reaction_ Engel and Reid (2009); Rothenberg
(2017); Petrucci (2002). It provides with an alternative reaction pathway to
make the reaction happen without changing the reactants and products.
Specifically, in this work we study a simple regime of catalyzed annihilation
with two processes leading to depletion of DM particles: $2\chi\to
2A^{\prime}$ and $3A^{\prime}\to 2\chi$, where $\chi$ and $A^{\prime}$ denote
dark matter and the catalyst respectively. We shown in Figure 1 a depiction of
how these annihilation channels result in depopulation of DM particles, that
is, three $2\chi\to 2A^{\prime}$ processes together with two $3A^{\prime}\to
2\chi$ effectively deplete two DM particles. Note that the assisted
annihilation Dey  _et al._ (2017); Cline  _et al._ (2017); Fitzpatrick  _et
al._ (2020) are not catalyzed reactions since the assisters are consumed in
the reaction. The co-SIMP process $\text{SM}+\chi+\chi\to\text{SM}+\chi$
Smirnov and Beacom (2020) is not catalyzed reaction either, since
$\chi+\chi\to\chi$ is kinetically forbidden and it’s groundless to discuss
enhancement of this unphysical process. Same thing happens to Ref. Dolgov
(2017). We acknowledge that catalyzed processes are also considered in the Big
Bang Nucleosynthesis (BBN) Pospelov (2007).

Figure 1: Schematic illustration of the catalyzed annihilation of DM $\chi$
(red line) with a catalyst $A^{\prime}$ (blue line). Three $2\chi\to
2A^{\prime}$ processes plus two $3A^{\prime}\to 2\chi$ effectively deplete the
number of DM particles by two.

The observed DM abundance can be reproduced in the catalyzed annihilation
paradigm for a wide mass range of dark matter. We emphasis that the thermal
evolution in the catalyzed annihilation paradigm is unique. Different with
WIMPs and SIMPs, where DM number density $n_{\chi}$ tracks Boltzmann
distribution and shrinks exponentially before freeze-out, the catalyzed
annihilation could lead to a polynomial suppression of $n_{\chi}$ as the
universe cools down,

$n_{\chi}\propto s^{3/2}\propto T^{9/2},$ (1)

where $s$ and $T$ denote entropy density and temperature of the universe.
Thus, the catalyzed annihilation lasts longer and freezes-out at late times.
To reproduce correct relic abundance, the cross section of DM annihilation
$2\chi\to 2A^{\prime}$ should be enhanced since there is less time to redshift
to today Dror  _et al._ (2016), which corresponds to enhanced indirect
detection signals.

##  II Catalyzed Freeze-out

In order for the catalyzed annihilation paradigm to work, there are several
requirements listed as follows:

  * •

The dark sector is nearly secluded.

  * •

The catalyst is long-lived ($\gtrsim 10^{-9}\text{s}$).

  * •

The catalyst is slightly lighter than DM ($1<m_{\chi}/m_{A^{\prime}}\lesssim
2$).

  * •

Annihilation channels as in Figure 1.

The dark sector should be secluded so that the annihilation channels to SM
particles freeze-out before the catalyzed annihilation. The catalyst is
lighter than DM and long-lived so that its number density is large and the
$3A^{\prime}\to 2\chi$ process is not suppressed, which ensures that the
catalyzed annihilation happens. If $A^{\prime}$ decays fast, the paradigm
recovers to the secluded DM regime Pospelov  _et al._ (2008); Arkani-Hamed
_et al._ (2009). Besides, the catalyzed annihilation will heat up the dark
sector. For simplicity, we assume the dark sector could scatter with SM
particles intensely enough to maintain thermal equilibrium with SM particles.

Figure 2: Thermal evolution of DM $\chi$ (solid red) and the catalyst
$A^{\prime}$ (solid blue). The dashed colored curves denote equilibrium yields
$\bar{y}\equiv\bar{n}/s$, while the dashed black curve shows the approximate
DM yield $y^{app}_{\chi}\equiv n^{app}_{\chi}/s$ during the catalyzed
annihilation stage with
$n_{\chi}^{app}=\sqrt{n_{A^{\prime}}^{3}\left<\sigma_{3}v^{2}\right>/\left<\sigma_{2}v\right>}$
from Eq. 4. The parameters are taken for the model presented in Section IV.

We show in Figure 2 a typical thermal history of dark matter that freezes-out
via catalyzed annihilation. For now, we are focused on the regime that the
mass ratio of dark matter and the catalyst $r\equiv m_{\chi}/m_{A^{\prime}}$
is no larger than $1.5$. As is shown in the figure, there are four stages in
the thermal evolution:

  1. 1.

_Equilibrium stage._ Both $\chi$ and $A^{\prime}$ stay in chemical equilibrium
due to the number-changing processes in the dark sector. The dominate number-
changing process is $3A^{\prime}\to 2\chi$ for $r\gtrsim 1.1$. Other processes
with DM in the initial state, e.g. $\chi A^{\prime}A^{\prime}\to\chi
A^{\prime}$, are suppressed and negligible, since $n_{\chi}\ll
n_{A^{\prime}}$.

$n_{\chi}\simeq\bar{n}_{\chi},n_{A^{\prime}}\simeq\bar{n}_{A^{\prime}}.$ (2)

$n_{\chi,A^{\prime}}$ denote the number densities and
$\bar{n}_{\chi,A^{\prime}}$ are the equilibrium densities. In the non-
relativistic limit, we have
$\bar{n}_{\chi,A^{\prime}}=g_{\chi,A^{\prime}}\left(\frac{m_{\chi,A^{\prime}}T}{2\pi}\right)^{3/2}e^{-m_{\chi,A^{\prime}}/T}$,
where $g$ denotes number of internal degrees of freedom.

  2. 2.

_Chemical stage._ $\chi$ and $A^{\prime}$ are chemically decoupled from
equilibrium, but they can still maintain chemical equilibrium with each other
via the $2\chi\leftrightarrow 2{A^{\prime}}$ process.

$n_{\chi}/\bar{n}_{\chi}\simeq n_{A^{\prime}}/\bar{n}_{A^{\prime}}.$ (3)

  3. 3.

_Catalyzed annihilation._ As the rate of $2A^{\prime}\to 2\chi$ (inverse
process of $2\chi\to 2A^{\prime}$) descends exponentially at low temperature,
the $3A^{\prime}\to 2\chi$ process dominates over it. The evolution of DM
number density is now controlled by the catalyzed annihilation, i.e. $2\chi\to
2A^{\prime}$ and $3A^{\prime}\to 2\chi$. Before freeze-out, the rates of the
$2\chi\to 2A^{\prime}$ and $3A^{\prime}\to 2\chi$ reactions are much larger
than the rate of change of $n_{\chi}$, as well as the Hubble rate and rates of
other reactions (see in Eq. 10). Thus, neglecting the subdominant terms, we
get an approximate relation,

$\left<\sigma_{2}v\right>n_{\chi}^{2}\simeq\left<\sigma_{3}v^{2}\right>n_{A^{\prime}}^{3}.$
(4)

We used $\left<\sigma_{2}v\right>$ and $\left<\sigma_{3}v^{2}\right>$ to
denote the thermally averaged cross sections of $2\chi\to 2A^{\prime}$ and
$3A^{\prime}\to 2\chi$ respectively. In this stage, since
$y_{A^{\prime}}\equiv n_{A^{\prime}}/s$ is practically constant and
$\left<\sigma_{2}v\right>$ and $\left<\sigma_{3}v^{2}\right>$ are polynomial
functions of $T$, Eq. 4 indicates that $n_{\chi}$ is polynomially suppressed.
It is similar to the scaling of the number density of the assisting particle
after DM freeze-out in Ref. Cline  _et al._ (2017).

  4. 4.

_Freeze-out._ As the universe expands, the rate of the catalyzed annihilation
descends and dark matter freezes-out.

The equilibrium stage ends when the rate of $3A^{\prime}\to 2\chi$ falls below
Hubble constant $H$. The temperature of departure from equilibrium $T_{c}$ can
be determined approximately with,

$\left<\sigma_{3}v^{2}\right>\bar{n}_{A^{\prime}}^{3}\simeq
H(\bar{n}_{A^{\prime}}+\bar{n}_{\chi}).$ (5)

We note that the annihilation channels to SM particles or the $3A^{\prime}\to
2A^{\prime}$ process can also deplete dark sector particles and $T_{c}$ could
be altered if these channels freeze-out later. The ending of the chemical
stage is insignificant since the freeze-out temperature $T_{f}$ and relic
abundance can be estimated without it. Lastly, the catalyzed annihilation
freezes-out when the rate drops below $H$. $T_{f}$ is determined by,

$\left<\sigma_{2}v\right>n_{\chi}^{2}\simeq\left<\sigma_{3}v^{2}\right>n_{A^{\prime}}^{3}\simeq
Hn_{\chi}.$ (6)

The relic abundance of DM can be estimated approximately in the same spirit of
WIMPs Kolb and Turner (1990); Dror  _et al._ (2016),

$\Omega_{\chi}=\frac{m_{\chi}s_{0}H_{m}}{\rho_{c}s_{m}}\frac{\sqrt{g_{\star,m}}}{\sqrt{g_{\star,f}}}\frac{x_{f}}{\left<\sigma_{2}v\right>}.$
(7)

where $x_{f}\equiv m_{\chi}/T_{f}$. Since $x_{f}$ is dependent on $x_{c}\equiv
m_{\chi}/T_{c}$, we solve for $x_{c}$ with Eq. 5 first. Simplifying Eq. 5, we
find,

$x_{c}=\frac{r}{2}\log\left[0.0024\frac{g_{A^{\prime}}^{2}m_{\chi}^{4}M_{\text{Pl}}\left<\sigma_{3}v^{2}\right>_{c}}{g_{\star,c}^{1/2}r^{3}x_{c}}\right],$
(8)

where $M_{\text{Pl}}$ is Planck mass. With $x_{c}$ determined in Eq. 8, we can
solve for $n_{A^{\prime}}$. Substituting the result into Eq. 6, we get,

$x_{f}=\left(1.2\frac{g_{\star,f}^{2}}{g_{\star,c}^{9/4}}\frac{M_{\text{Pl}}^{1/2}\left<\sigma_{2}v\right>_{f}\left<\sigma_{3}v^{2}\right>_{f}}{m_{\chi}\left<\sigma_{3}v^{2}\right>_{c}^{3/2}}x_{c}^{6}\right)^{1/5}.$
(9)

The subscripts $m,c,f$ in the equations above mark the temperatures,
$T=m_{\chi},T_{c},T_{f}$, respectively, for the quantities, including entropy
density $s$, Hubble constant $H$, effective degrees of freedom
$g_{\star}$111We neglect the differences between effective entropy degrees of
freedom $g_{\star,s}$ and effective energy degrees of freedom $g_{\star}$ as
in Ref. Dror  _et al._ (2016) and the thermally averaged cross sections. Note
that if $x_{c}$ is delayed due to annihilation to SM particles or
$3A^{\prime}\to 2A^{\prime}$, Eq. 9 and Eq. 8 should be modified to include
these processes and DM will freeze-out earlier in this case with a smaller
relic abundance.

Based on the partial wave unitarity limit Griest and Kamionkowski (1990),
$\sigma_{2}v\leq\frac{4\pi}{m_{\chi}^{2}v}$, we can estimate the upper bound
of DM mass from Eq. 7 for the catalyzed annihilation paradigm. With
$x_{f}\gtrsim 100$, we deduce $m_{\chi}\lesssim 100\text{TeV}.$ Compared to
SIMP dark matter that lives in the MeV scale Hochberg  _et al._ (2014), it is
compelling to notice that $3\to 2$ process can apply to such a heavy dark
matter.

In order to study the thermal evolution and DM freeze-out in a quantitative
way, we turn to the Boltzmann equations. As is discussed above, we neglect the
subdominant $3\to 2$ annihilation channels, including $\chi
A^{\prime}A^{\prime}\to\chi A^{\prime}$, $\chi\chi A^{\prime}\to
A^{\prime}A^{\prime}$, $\chi\chi A^{\prime}\to\chi\chi$, $\chi\chi\chi\to\chi
A^{\prime}$ and assume $3A^{\prime}\to 2A^{\prime}$ is subdominant. If
$A^{\prime}$ decays to SM particles, the Boltzmann equations reads,

$\displaystyle\dot{n}_{\chi}$
$\displaystyle+3Hn_{\chi}=-\left<\sigma_{2}v\right>\left(n_{\chi}^{2}-\bar{n}_{\chi}^{2}\frac{n_{A^{\prime}}^{2}}{\bar{n}_{A^{\prime}}^{2}}\right)$
(10)
$\displaystyle+\left<\sigma_{3}v^{2}\right>\left(n_{A^{\prime}}^{3}-\bar{n}_{A^{\prime}}^{3}\frac{n_{\chi}^{2}}{\bar{n}_{\chi}^{2}}\right),$
$\displaystyle\dot{n}_{A^{\prime}}$
$\displaystyle+3Hn_{A^{\prime}}=+\left<\sigma_{2}v\right>\left(n_{\chi}^{2}-\bar{n}_{\chi}^{2}\frac{n_{A^{\prime}}^{2}}{\bar{n}_{A^{\prime}}^{2}}\right)$
$\displaystyle-\frac{3}{2}\left<\sigma_{3}v^{2}\right>\left(n_{A^{\prime}}^{3}-\bar{n}_{A^{\prime}}^{3}\frac{n_{\chi}^{2}}{\bar{n}_{\chi}^{2}}\right)-\left<\Gamma_{A^{\prime}}\right>(n_{A^{\prime}}-\bar{n}_{A^{\prime}}).$

The yield $y_{\chi,A^{\prime}}\equiv n_{\chi,A^{\prime}}/s$ can be solved
numerically and are shown in Figure 2.

##  III Mass Ratio

Figure 3: Curves of relic abundance $\Omega_{\chi}h^{2}$ for
$m_{\chi}=350\text{GeV}$ (solid red) and $m_{\chi}=1000\text{GeV}$ (solid
blue) with respect to different mass ratio $r$. For dashed colored curves, the
$4A^{\prime}\to 2\chi$ process is neglected. The dashed black curve denotes
the observed DM relic abundance Aghanim  _et al._ (2020a). The parameters are
taken for the model presented in Section IV.

In previous section, we concentrated on the mass ratio $r\leq 1.5$. In fact,
the catalyzed annihilation paradigm can go beyond this limit. Firstly, when
the mass ratio is slightly larger than $1.5$, i.e.
$3m_{A^{\prime}}<2m_{\chi}$, $\left<\sigma_{3}v^{2}\right>$ is exponentially
suppressed as the temperature goes down,

$\left<\sigma_{3}v^{2}\right>\propto e^{-(2r-3)x/r},$ (11)

where $x\equiv m_{\chi}/T$. During the catalyzed annihilation period, with
less DM particles produced via $3A^{\prime}\to 2\chi$ process since the cross
section is smaller, the DM number density shrinks more sharply. Consequently,
the catalyzed annihilation freezes-out much earlier.

As the mass ratio grows, when $r\lesssim 2$, it is intriguing to notice that
the $4A^{\prime}\to 2\chi$ process may play a part in the catalyzed
annihilation. To be specific, after a period of catalyzed annihilation
governed by $2\chi\to 2A^{\prime}$ and $3A^{\prime}\to 2\chi$ as usual, there
would be an extra stage of catalyzed annihilation predominated by $2\chi\to
2A^{\prime}$ and $4A^{\prime}\to 2\chi$, in which the non-suppressed
$4A^{\prime}\to 2\chi$ process takes over the role of converting $A^{\prime}$
to DM particles since the cross section of $3A^{\prime}\to 2\chi$ is
exponentially suppressed (Eq. 11), Similar to Eq. 4, we can deduce an
approximate relation that holds in this stage,

$\left<\sigma_{2}v\right>n_{\chi}^{2}\simeq\left<\sigma_{4}v^{3}\right>n_{A^{\prime}}^{4},$
(12)

where $\left<\sigma_{4}v^{3}\right>$ denotes the thermally averaged cross
section for $4A^{\prime}\to 2\chi$. The presence of $4A^{\prime}\to 2\chi$ is
essential. If it is neglected, as is discussed previously, $n_{\chi}$ shrinks
sharply and dark matter freezes-out early. Once $4A^{\prime}\to 2\chi$ takes
charge, the sharply falling of $n_{\chi}$ is bent and the polynomial
suppression recovers (compared to Eq. 1).

$n_{\chi}\propto s^{2}\propto T^{6}.$ (13)

Thus, the catalyzed annihilation freezes-out at later times, leading to
enhanced DM relic abundance.

For even larger mass ratio, we expect the processes with more catalysts
annihilating to two DM particles, e.g. $5A^{\prime}\to 2\chi$, to possibly
play a role in the catalyzed annihilation, especially when the dark sector is
strongly coupled.

We show in Figure 3 the variation of DM relic abundance $\Omega_{\chi}h^{2}$
with different mass ratio. When the mass ratio passes the critical value of
$1.5$, $\Omega_{\chi}h^{2}$ decreases rapidly. On the other hand, for
$r\lesssim 2$, relic abundance is uplifted if $4A^{\prime}\to 2\chi$ process
is included.

##  IV A Model

The requirements for realization of the catalyzed annihilation presented in
Section II can be easily met in many models. In this section, we simply
present a dark photon model Holdom (1986); Raggi and Kozhuharov (2015); Bauer
_et al._ (2020); Lao  _et al._ (2020); Fabbrichesi  _et al._ (2020) with a
Dirac fermion $\chi$ charged under a novel $U(1)^{\prime}$ gauge group and
$A^{\prime}$ being the gauge field. The Lagrangian for the dark sector is,

$\mathcal{L}_{\text{DS}}=-\frac{1}{4}F^{\prime}_{\mu\nu}F^{\prime\mu\nu}+\frac{1}{2}m_{A^{\prime}}^{2}A^{\prime}_{\mu}A^{\prime\mu}+\bar{\chi}(i\not{D}-m_{\chi})\chi,$
(14)

where $\not{D}=\not{\partial}-ig_{D}\not{A^{\prime}}$ and $g_{D}$ is the gauge
coupling constant. The mass of the dark photon can be generated via the Higgs
mechanism (or Stueckelberg mechanism Stueckelberg (1938); Ruegg and Ruiz-
Altaba (2004)). We assume the dark Higgs boson is heavy and can be neglected.
SM particles are neutral under the $U(1)^{\prime}$ gauge group. The dark
photon can be kinetically mixed with SM hypercharge field.

$\mathcal{L}_{\text{mix}}=-\frac{\epsilon}{2\cos{\theta_{W}}}F^{\prime}_{\mu\nu}B^{\mu\nu}.$
(15)

$\epsilon$ is the mixing constant and $\theta_{W}$ denotes the Weinberg angle.
$B^{\mu}$ is SM hypercharge field. Therefore, the dark sector can communicate
with SM particles via the mixing and the dark photon $A^{\prime}$ can decay to
SM particles. $\epsilon$ should be small so that the dark photon is long-lived
and acts as the catalyst. Additionally, the kinetic mixing could not keep the
dark sector in thermal equilibrium with SM since $\epsilon$ is small. In order
to thermalize the dark sector, we need another portal for the dark sector to
interact with SM particles, which might be the dark Higgs. Anyhow, we won’t
model this part and simply assume that the dark sector stays in thermal
equilibrium before freeze-out.

Figure 4: Phase diagram and constraints for the dark photon model in the
$(m_{\chi},\Gamma_{A^{\prime}})$ plane with $r=1.5$. The solid black curve
marks the boundary of the secluded phase and the catalyzed phase of the model.
Correct relic abundance can be reproduced for each point in the figure by
varying the value of $g_{D}$. Especially, we show in gray dashed curves for
five different values of $g_{D}$ that reproduce correct relic abundance. The
non-perturbative region is painted gray, while the color-shaded regions denote
the bounds from various experiments and observations correspondingly.

We show in Figure 4 different phases for the model in the calculation of DM
relic abundance. For short-lived dark photon, before dark matter freezes-out,
it simply stays in equilibrium with SM particles via the decay and inverse-
decay process. When DM particles annihilates into the dark photon, it
immediately decays. This is the secluded phase of the model. On the other
hand, when the dark photon width $\Gamma_{A^{\prime}}$ is small, the catalyzed
annihilation emerges. It is a continuous shift, since the decay of the dark
photon can occur during the catalyzed annihilation. When the dark photon
decays after DM freeze-out, $\Omega_{\chi}h^{2}$ is independent with
$\Gamma_{A^{\prime}}$. If $A^{\prime}$ is sufficiently long-lived
($\Gamma_{A^{\prime}}\lesssim 10^{-23}\text{GeV}$), it would come to dominate
the energy density of the universe. When it decays, considerable entropy is
produced and DM abundance is diluted. This effect can help to circumvent the
upper bound of DM mass Gelmini  _et al._ (2006); Bhupal Dev  _et al._ (2016);
Berlin  _et al._ (2016) (see in Figure 4). We show in dashed gray curves for
five different values of $g_{D}$ that reproduce the observed relic abundance
in Figure 4.

The catalyzed annihilation paradigm is constrained by numerous terrestrial and
celestial experiments and observations. Firstly, the residual annihilation of
$2\chi\to 2A^{\prime}$ after freeze-out will distort the anisotropy of the
Cosmic Microwave Background (CMB) since the decay products of $A^{\prime}$ are
electrically charged particles Padmanabhan and Finkbeiner (2005); Galli  _et
al._ (2009); Kawasaki  _et al._ (2016); Slatyer (2016); Liu  _et al._ (2016);
Aghanim  _et al._ (2020b); Cang  _et al._ (2020). Similarly, the signal of DM
annihilation at present is detectable in indirect detection experiments
Adriani  _et al._ (2013); Ackermann  _et al._ (2015); Albert  _et al._ (2017);
Aguilar  _et al._ (2016); Profumo  _et al._ (2018); Abdallah  _et al._ (2018).
The signal is enhanced compared to WIMPs with
$\left<\sigma_{2}v\right>_{0}\simeq 5\cdot 10^{-25}\text{cm}^{3}\text{s}^{-1}$
for $s$-wave annihilation. We used bounds from Fermi-LAT experiment to
constrain our model in Figure 4. The late time decay of $A^{\prime}$, on the
other hand, is also stringently constrained by CMB Chen and Kamionkowski
(2004); Acharya and Khatri (2019a, b) as well as BBN Kawasaki  _et al._
(2005); Jedamzik (2006); Kawasaki  _et al._ (2018). We note that these bounds
can be evaded in models that the catalyst decays into neutrinos or dark
radiations Ichiki  _et al._ (2004); Poulin  _et al._ (2016); Nygaard  _et al._
(2020). Since the dark sector is highly secluded, the constraints on
scattering with nucleons Cui  _et al._ (2017); Aprile  _et al._ (2018); Leane
and Smirnov (2021); Leane  _et al._ (2021) are evaded.

For light dark photon, beam dump and fixed target experiments provide great
sensitivity on the mixing coupling constant $\epsilon$ Bergsma  _et al._
(1986); Bjorken  _et al._ (1988); Davier and Nguyen Ngoc (1989); Blümlein and
Brunner (2014). There are also lots of new experiments Anelli  _et al._
(2015); Feng  _et al._ (2018); Berlin  _et al._ (2018); Adrian  _et al._
(2018) proposed in recent years that are focused on long-lived particles.
Besides, the long-lived dark photon can enhance the cooling of supernova and
the constraints from SN 1987A Chang  _et al._ (2017, 2018); Sung  _et al._
(2019) is widely discussed. These bounds on the dark photon model are
considered and presented in Figure 4.

##  V Conclusion and discussion

We proposed a novel paradigm for thermal relic dark matter, yielding the
observed relic abundance. The distinctive wisdom of the paradigm is that the
dark matter freeze-out proceeds via catalyzed annihilation. We discussed in
detail the scenario that the catalyzed annihilation includes $2\chi\to
2A^{\prime}$ and $3A^{\prime}\to 2\chi$, where $\chi$ and $A^{\prime}$ are
dark matter and the catalyst respectively. The paradigm applies for a wide
mass range of dark matter, from 1MeV to 100TeV, with a unique thermal history
compared with WIMPs and SIMPs. Besides, the paradigm offers rich phenomenology
including indirect DM search and long-lived particles.

We note that thermal decoupling effects can significantly modify dark matter
relic abundance Kuflik  _et al._ (2016); Pappadopulo  _et al._ (2016); Farina
_et al._ (2016); D’Agnolo  _et al._ (2017); Kuflik  _et al._ (2017);
Fitzpatrick  _et al._ (2020). We leave this to future works Xing and Zhu .
Additionally, catalyzed annihilation dominating DM abundance can go far beyond
the reactions considered here and should be investigated further.

## Acknowledgments

C.Y.X. would like to thank Yan-Fang Bai for encouragement. This work is
supported by the National Science Foundation of China under Grants No.
11635001, 11875072.

## References

  * Bauer and Plehn (2019) Martin Bauer and Tilman Plehn,  _Yet Another Introduction to Dark Matter: The Particle Physics Approach_, Lecture Notes in Physics, Vol. 959 (Springer, 2019) arXiv:1705.01987 [hep-ph] .
  * Lee and Weinberg (1977) Benjamin W. Lee and Steven Weinberg, “Cosmological Lower Bound on Heavy Neutrino Masses,” Phys. Rev. Lett. 39, 165–168 (1977).
  * Bertone  _et al._ (2005) Gianfranco Bertone, Dan Hooper, and Joseph Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept. 405, 279–390 (2005), arXiv:hep-ph/0404175 .
  * Arcadi  _et al._ (2018) Giorgio Arcadi, Maíra Dutra, Pradipta Ghosh, Manfred Lindner, Yann Mambrini, Mathias Pierre, Stefano Profumo, and Farinaldo S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78, 203 (2018), arXiv:1703.07364 [hep-ph] .
  * Roszkowski  _et al._ (2018) Leszek Roszkowski, Enrico Maria Sessolo, and Sebastian Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81, 066201 (2018), arXiv:1707.06277 [hep-ph] .
  * Griest and Seckel (1991) Kim Griest and David Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D 43, 3191–3203 (1991).
  * Ellis  _et al._ (1998) John R. Ellis, Toby Falk, and Keith A. Olive, “Neutralino - Stau coannihilation and the cosmological upper limit on the mass of the lightest supersymmetric particle,” Phys. Lett. B 444, 367–372 (1998), arXiv:hep-ph/9810360 .
  * D’Agnolo  _et al._ (2018) Raffaele Tito D’Agnolo, Cristina Mondino, Joshua T. Ruderman, and Po-Jen Wang, “Exponentially Light Dark Matter from Coannihilation,” JHEP 08, 079 (2018), arXiv:1803.02901 [hep-ph] .
  * D’Eramo and Thaler (2010) Francesco D’Eramo and Jesse Thaler, “Semi-annihilation of Dark Matter,” JHEP 06, 109 (2010), arXiv:1003.5912 [hep-ph] .
  * D’Agnolo and Ruderman (2015) Raffaele Tito D’Agnolo and Joshua T. Ruderman, “Light Dark Matter from Forbidden Channels,” Phys. Rev. Lett. 115, 061301 (2015), arXiv:1505.07107 [hep-ph] .
  * Pospelov  _et al._ (2008) Maxim Pospelov, Adam Ritz, and Mikhail B. Voloshin, “Secluded WIMP Dark Matter,” Phys. Lett. B 662, 53–61 (2008), arXiv:0711.4866 [hep-ph] .
  * Arkani-Hamed  _et al._ (2009) Nima Arkani-Hamed, Douglas P. Finkbeiner, Tracy R. Slatyer, and Neal Weiner, “A Theory of Dark Matter,” Phys. Rev. D 79, 015014 (2009), arXiv:0810.0713 [hep-ph] .
  * Feng and Kumar (2008) Jonathan L. Feng and Jason Kumar, “The WIMPless Miracle: Dark-Matter Particles without Weak-Scale Masses or Weak Interactions,” Phys. Rev. Lett. 101, 231301 (2008), arXiv:0803.4196 [hep-ph] .
  * Belanger and Park (2012) Genevieve Belanger and Jong-Chul Park, “Assisted freeze-out,” JCAP 03, 038 (2012), arXiv:1112.4491 [hep-ph] .
  * Dror  _et al._ (2016) Jeff Asaf Dror, Eric Kuflik, and Wee Hao Ng, “Codecaying Dark Matter,” Phys. Rev. Lett. 117, 211801 (2016), arXiv:1607.03110 [hep-ph] .
  * D’Agnolo  _et al._ (2017) Raffaele Tito D’Agnolo, Duccio Pappadopulo, and Joshua T. Ruderman, “Fourth Exception in the Calculation of Relic Abundances,” Phys. Rev. Lett. 119, 061102 (2017), arXiv:1705.08450 [hep-ph] .
  * Garny  _et al._ (2019) Mathias Garny, Jan Heisig, Marco Hufnagel, Benedikt Lülf, and Stefan Vogl, “Conversion-driven freeze-out: Dark matter genesis beyond the WIMP paradigm,” PoS CORFU2018, 092 (2019), arXiv:1904.00238 [hep-ph] .
  * Maity and Ray (2020) Tarak Nath Maity and Tirtha Sankar Ray, “Exchange driven freeze out of dark matter,” Phys. Rev. D 101, 103013 (2020), arXiv:1908.10343 [hep-ph] .
  * Kramer  _et al._ (2020) Eric David Kramer, Eric Kuflik, Noam Levi, Nadav Joseph Outmezguine, and Joshua T. Ruderman, “Heavy Thermal Relics from Zombie Collisions,” (2020), arXiv:2003.04900 [hep-ph] .
  * Dey  _et al._ (2017) Ujjal Kumar Dey, Tarak Nath Maity, and Tirtha Sankar Ray, “Light Dark Matter through Assisted Annihilation,” JCAP 03, 045 (2017), arXiv:1612.09074 [hep-ph] .
  * Cline  _et al._ (2017) James M. Cline, Hongwan Liu, Tracy Slatyer, and Wei Xue, “Enabling Forbidden Dark Matter,” Phys. Rev. D 96, 083521 (2017), arXiv:1702.07716 [hep-ph] .
  * Hochberg  _et al._ (2014) Yonit Hochberg, Eric Kuflik, Tomer Volansky, and Jay G. Wacker, “Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 113, 171301 (2014), arXiv:1402.5143 [hep-ph] .
  * Hochberg  _et al._ (2015) Yonit Hochberg, Eric Kuflik, Hitoshi Murayama, Tomer Volansky, and Jay G. Wacker, “Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 115, 021301 (2015), arXiv:1411.3727 [hep-ph] .
  * Bernal and Chu (2016) Nicolas Bernal and Xiaoyong Chu, “$\mathbb{Z}_{2}$ SIMP Dark Matter,” JCAP 01, 006 (2016), arXiv:1510.08527 [hep-ph] .
  * Bernal  _et al._ (2017) Nicolás Bernal, Xiaoyong Chu, and Josef Pradler, “Simply split strongly interacting massive particles,” Phys. Rev. D 95, 115023 (2017), arXiv:1702.04906 [hep-ph] .
  * Smirnov and Beacom (2020) Juri Smirnov and John F. Beacom, “New Freezeout Mechanism for Strongly Interacting Dark Matter,” Phys. Rev. Lett. 125, 131301 (2020), arXiv:2002.04038 [hep-ph] .
  * Engel and Reid (2009) Thomas Engel and Philip Reid,  _Physical chemistry_ , 2nd ed. (Pearson Education, Upper Saddle River, N.J. ; Harlow, 2009).
  * Rothenberg (2017) Gadi Rothenberg,  _Catalysis_ , 2nd ed. (Wiley, Newark, 2017).
  * Petrucci (2002) Ralph H Petrucci,  _General chemistry : principles and modern applications_ , 8th ed. (Prentice Hall, Upper Saddle River, N.J., 2002).
  * Fitzpatrick  _et al._ (2020) Patrick J. Fitzpatrick, Hongwan Liu, Tracy R. Slatyer, and Yu-Dai Tsai, “New Pathways to the Relic Abundance of Vector-Portal Dark Matter,” (2020), arXiv:2011.01240 [hep-ph] .
  * Dolgov (2017) A. D. Dolgov, “New Old Mechanism of Dark Matter Burning,” (2017), arXiv:1705.03689 [hep-ph] .
  * Pospelov (2007) Maxim Pospelov, “Particle physics catalysis of thermal Big Bang Nucleosynthesis,” Phys. Rev. Lett. 98, 231301 (2007), arXiv:hep-ph/0605215 .
  * Kolb and Turner (1990) Edward W. Kolb and Michael S. Turner,  _The Early Universe_ , Vol. 69 (1990).
  * Griest and Kamionkowski (1990) Kim Griest and Marc Kamionkowski, “Unitarity Limits on the Mass and Radius of Dark Matter Particles,” Phys. Rev. Lett. 64, 615 (1990).
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 [astro-ph.CO] .
  * Holdom (1986) Bob Holdom, “Two U(1)’s and Epsilon Charge Shifts,” Phys. Lett. B 166, 196–198 (1986).
  * Raggi and Kozhuharov (2015) Mauro Raggi and Venelin Kozhuharov, “Results and perspectives in dark photon physics,” Riv. Nuovo Cim. 38, 449–505 (2015).
  * Bauer  _et al._ (2020) Martin Bauer, Patrick Foldenauer, and Joerg Jaeckel, “Hunting All the Hidden Photons,” JHEP 18, 094 (2020), arXiv:1803.05466 [hep-ph] .
  * Lao  _et al._ (2020) Juebin Lao, Chengfeng Cai, Zhao-Huan Yu, Yu-Pan Zeng, and Hong-Hao Zhang, “Fermionic and scalar dark matter with hidden $\mathrm{U}(1)$ gauge interaction and kinetic mixing,” Phys. Rev. D 101, 095031 (2020), arXiv:2003.02516 [hep-ph] .
  * Fabbrichesi  _et al._ (2020) Marco Fabbrichesi, Emidio Gabrielli, and Gaia Lanfranchi, “The Dark Photon,” (2020), arXiv:2005.01515 [hep-ph] .
  * Stueckelberg (1938) E. C. G. Stueckelberg, “Interaction energy in electrodynamics and in the field theory of nuclear forces,” Helv. Phys. Acta 11, 225–244 (1938).
  * Ruegg and Ruiz-Altaba (2004) Henri Ruegg and Marti Ruiz-Altaba, “The Stueckelberg field,” Int. J. Mod. Phys. A 19, 3265–3348 (2004), arXiv:hep-th/0304245 .
  * Gelmini  _et al._ (2006) Graciela Gelmini, Paolo Gondolo, Adrian Soldatenko, and Carlos E. Yaguna, “The Effect of a late decaying scalar on the neutralino relic density,” Phys. Rev. D 74, 083514 (2006), arXiv:hep-ph/0605016 .
  * Bhupal Dev  _et al._ (2016) P. S. Bhupal Dev, Rabindra N. Mohapatra, and Yongchao Zhang, “Naturally stable right-handed neutrino dark matter,” JHEP 11, 077 (2016), arXiv:1608.06266 [hep-ph] .
  * Berlin  _et al._ (2016) Asher Berlin, Dan Hooper, and Gordan Krnjaic, “PeV-Scale Dark Matter as a Thermal Relic of a Decoupled Sector,” Phys. Lett. B 760, 106–111 (2016), arXiv:1602.08490 [hep-ph] .
  * Padmanabhan and Finkbeiner (2005) Nikhil Padmanabhan and Douglas P. Finkbeiner, “Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects,” Phys. Rev. D 72, 023508 (2005), arXiv:astro-ph/0503486 .
  * Galli  _et al._ (2009) Silvia Galli, Fabio Iocco, Gianfranco Bertone, and Alessandro Melchiorri, “CMB constraints on Dark Matter models with large annihilation cross-section,” Phys. Rev. D 80, 023505 (2009), arXiv:0905.0003 [astro-ph.CO] .
  * Kawasaki  _et al._ (2016) Masahiro Kawasaki, Kazunori Nakayama, and Toyokazu Sekiguchi, “CMB Constraint on Dark Matter Annihilation after Planck 2015,” Phys. Lett. B 756, 212–215 (2016), arXiv:1512.08015 [astro-ph.CO] .
  * Slatyer (2016) Tracy R. Slatyer, “Indirect dark matter signatures in the cosmic dark ages. I. Generalizing the bound on s-wave dark matter annihilation from Planck results,” Phys. Rev. D 93, 023527 (2016), arXiv:1506.03811 [hep-ph] .
  * Liu  _et al._ (2016) Hongwan Liu, Tracy R. Slatyer, and Jesús Zavala, “Contributions to cosmic reionization from dark matter annihilation and decay,” Phys. Rev. D 94, 063507 (2016), arXiv:1604.02457 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), “Planck 2018 results. V. CMB power spectra and likelihoods,” Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 [astro-ph.CO] .
  * Cang  _et al._ (2020) Junsong Cang, Yu Gao, and Yin-Zhe Ma, “Probing dark matter with future CMB measurements,” Phys. Rev. D 102, 103005 (2020), arXiv:2002.03380 [astro-ph.CO] .
  * Adriani  _et al._ (2013) O. Adriani _et al._ (PAMELA), “Cosmic-Ray Positron Energy Spectrum Measured by PAMELA,” Phys. Rev. Lett. 111, 081102 (2013), arXiv:1308.0133 [astro-ph.HE] .
  * Ackermann  _et al._ (2015) M. Ackermann _et al._ (Fermi-LAT), “Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data,” Phys. Rev. Lett. 115, 231301 (2015), arXiv:1503.02641 [astro-ph.HE] .
  * Albert  _et al._ (2017) A. Albert _et al._ (Fermi-LAT, DES), “Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT,” Astrophys. J. 834, 110 (2017), arXiv:1611.03184 [astro-ph.HE] .
  * Aguilar  _et al._ (2016) M. Aguilar _et al._ (AMS), “Antiproton Flux, Antiproton-to-Proton Flux Ratio, and Properties of Elementary Particle Fluxes in Primary Cosmic Rays Measured with the Alpha Magnetic Spectrometer on the International Space Station,” Phys. Rev. Lett. 117, 091103 (2016).
  * Profumo  _et al._ (2018) Stefano Profumo, Farinaldo S. Queiroz, Joseph Silk, and Clarissa Siqueira, “Searching for Secluded Dark Matter with H.E.S.S., Fermi-LAT, and Planck,” JCAP 03, 010 (2018), arXiv:1711.03133 [hep-ph] .
  * Abdallah  _et al._ (2018) H. Abdallah _et al._ (HESS), “Search for $\gamma$-Ray Line Signals from Dark Matter Annihilations in the Inner Galactic Halo from 10 Years of Observations with H.E.S.S.” Phys. Rev. Lett. 120, 201101 (2018), arXiv:1805.05741 [astro-ph.HE] .
  * Chen and Kamionkowski (2004) Xue-Lei Chen and Marc Kamionkowski, “Particle decays during the cosmic dark ages,” Phys. Rev. D 70, 043502 (2004), arXiv:astro-ph/0310473 .
  * Acharya and Khatri (2019a) Sandeep Kumar Acharya and Rishi Khatri, “New CMB spectral distortion constraints on decaying dark matter with full evolution of electromagnetic cascades before recombination,” Phys. Rev. D 99, 123510 (2019a), arXiv:1903.04503 [astro-ph.CO] .
  * Acharya and Khatri (2019b) Sandeep Kumar Acharya and Rishi Khatri, “CMB anisotropy and BBN constraints on pre-recombination decay of dark matter to visible particles,” JCAP 12, 046 (2019b), arXiv:1910.06272 [astro-ph.CO] .
  * Kawasaki  _et al._ (2005) Masahiro Kawasaki, Kazunori Kohri, and Takeo Moroi, “Big-Bang nucleosynthesis and hadronic decay of long-lived massive particles,” Phys. Rev. D 71, 083502 (2005), arXiv:astro-ph/0408426 .
  * Jedamzik (2006) Karsten Jedamzik, “Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles,” Phys. Rev. D 74, 103509 (2006), arXiv:hep-ph/0604251 .
  * Kawasaki  _et al._ (2018) Masahiro Kawasaki, Kazunori Kohri, Takeo Moroi, and Yoshitaro Takaesu, “Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles,” Phys. Rev. D 97, 023502 (2018), arXiv:1709.01211 [hep-ph] .
  * Ichiki  _et al._ (2004) Kiyotomo Ichiki, Masamune Oguri, and Keitaro Takahashi, “WMAP constraints on decaying cold dark matter,” Phys. Rev. Lett. 93, 071302 (2004), arXiv:astro-ph/0403164 .
  * Poulin  _et al._ (2016) Vivian Poulin, Pasquale D. Serpico, and Julien Lesgourgues, “A fresh look at linear cosmological constraints on a decaying dark matter component,” JCAP 08, 036 (2016), arXiv:1606.02073 [astro-ph.CO] .
  * Nygaard  _et al._ (2020) Andreas Nygaard, Thomas Tram, and Steen Hannestad, “Updated constraints on decaying cold dark matter,” (2020), arXiv:2011.01632 [astro-ph.CO] .
  * Cui  _et al._ (2017) Xiangyi Cui _et al._ (PandaX-II), “Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment,” Phys. Rev. Lett. 119, 181302 (2017), arXiv:1708.06917 [astro-ph.CO] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Leane and Smirnov (2021) Rebecca K. Leane and Juri Smirnov, “Exoplanets as Sub-GeV Dark Matter Detectors,” Phys. Rev. Lett. 126, 161101 (2021), arXiv:2010.00015 [hep-ph] .
  * Leane  _et al._ (2021) Rebecca K. Leane, Tim Linden, Payel Mukhopadhyay, and Natalia Toro, “Celestial-Body Focused Dark Matter Annihilation Throughout the Galaxy,” Phys. Rev. D 103, 075030 (2021), arXiv:2101.12213 [astro-ph.HE] .
  * Bergsma  _et al._ (1986) F. Bergsma _et al._ (CHARM), “A Search for Decays of Heavy Neutrinos in the Mass Range 0.5-{GeV} to 2.8-{GeV},” Phys. Lett. B 166, 473–478 (1986).
  * Bjorken  _et al._ (1988) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker, and P. Rassmann, “Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump,” Phys. Rev. D 38, 3375 (1988).
  * Davier and Nguyen Ngoc (1989) M. Davier and H. Nguyen Ngoc, “An Unambiguous Search for a Light Higgs Boson,” Phys. Lett. B 229, 150–155 (1989).
  * Blümlein and Brunner (2014) Johannes Blümlein and Jürgen Brunner, “New Exclusion Limits on Dark Gauge Forces from Proton Bremsstrahlung in Beam-Dump Data,” Phys. Lett. B 731, 320–326 (2014), arXiv:1311.3870 [hep-ph] .
  * Anelli  _et al._ (2015) M. Anelli _et al._ (SHiP), “A facility to Search for Hidden Particles (SHiP) at the CERN SPS,” (2015), arXiv:1504.04956 [physics.ins-det] .
  * Feng  _et al._ (2018) Jonathan L. Feng, Iftah Galon, Felix Kling, and Sebastian Trojanowski, “ForwArd Search ExpeRiment at the LHC,” Phys. Rev. D 97, 035001 (2018), arXiv:1708.09389 [hep-ph] .
  * Berlin  _et al._ (2018) Asher Berlin, Stefania Gori, Philip Schuster, and Natalia Toro, “Dark Sectors at the Fermilab SeaQuest Experiment,” Phys. Rev. D 98, 035011 (2018), arXiv:1804.00661 [hep-ph] .
  * Adrian  _et al._ (2018) P. H. Adrian _et al._ (HPS), “Search for a dark photon in electroproduced $e^{+}e^{-}$ pairs with the Heavy Photon Search experiment at JLab,” Phys. Rev. D 98, 091101 (2018), arXiv:1807.11530 [hep-ex] .
  * Chang  _et al._ (2017) Jae Hyeok Chang, Rouven Essig, and Samuel D. McDermott, “Revisiting Supernova 1987A Constraints on Dark Photons,” JHEP 01, 107 (2017), arXiv:1611.03864 [hep-ph] .
  * Chang  _et al._ (2018) Jae Hyeok Chang, Rouven Essig, and Samuel D. McDermott, “Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle,” JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Sung  _et al._ (2019) Allan Sung, Huitzu Tu, and Meng-Ru Wu, “New constraint from supernova explosions on light particles beyond the Standard Model,” Phys. Rev. D 99, 121305 (2019), arXiv:1903.07923 [hep-ph] .
  * Kuflik  _et al._ (2016) Eric Kuflik, Maxim Perelstein, Nicolas Rey-Le Lorier, and Yu-Dai Tsai, “Elastically Decoupling Dark Matter,” Phys. Rev. Lett. 116, 221302 (2016), arXiv:1512.04545 [hep-ph] .
  * Pappadopulo  _et al._ (2016) Duccio Pappadopulo, Joshua T. Ruderman, and Gabriele Trevisan, “Dark matter freeze-out in a nonrelativistic sector,” Phys. Rev. D 94, 035005 (2016), arXiv:1602.04219 [hep-ph] .
  * Farina  _et al._ (2016) Marco Farina, Duccio Pappadopulo, Joshua T. Ruderman, and Gabriele Trevisan, “Phases of Cannibal Dark Matter,” JHEP 12, 039 (2016), arXiv:1607.03108 [hep-ph] .
  * Kuflik  _et al._ (2017) Eric Kuflik, Maxim Perelstein, Nicolas Rey-Le Lorier, and Yu-Dai Tsai, “Phenomenology of ELDER Dark Matter,” JHEP 08, 078 (2017), arXiv:1706.05381 [hep-ph] .
  * (87) Chuan-Yang Xing and Shou-hua Zhu, In preparation.


# From ergodicity to Stark many-body localization in spin chains with single-
ion anisotropy

M. G. Sousa  Instituto de Física, Universidade Federal de Uberlândia,
Uberlândia, Minas Gerais 38400-902, Brazil.  Rafael F. P. Costa  Instituto de
Física, Universidade Federal de Uberlândia, Uberlândia, Minas Gerais
38400-902, Brazil.  G. D. de Moraes Neto<EMAIL_ADDRESS>Department of
Physics, Zhejiang Normal University, Jinhua 321004, Zhejiang, P. R. China  E.
Vernek<EMAIL_ADDRESS>Instituto de Física, Universidade Federal de Uberlândia,
Uberlândia, Minas Gerais 38400-902, Brazil.  Department of Physics, Zhejiang
Normal University, Jinhua 321004, Zhejiang, P. R. China

###### Abstract

The principles of ergodicity and thermalization constitute the foundation of
statistical mechanics, positing that a many-body system progressively loses
its local information as it evolves. Nevertheless, these principles can be
disrupted when thermalization dynamics lead to the conservation of local
information, as observed in the phenomenon known as many-body localization.
Quantum spin chains provide a fundamental platform for exploring the dynamics
of closed interacting quantum many-body systems. This study explores the
dynamics of a spin chain with $S\geq 1/2$ within the Majumdar-Ghosh model,
incorporating a non-uniform magnetic field and single-ion anisotropy. Through
the use of exact numerical diagonalization, we unveil that a nearly constant-
gradient magnetic field suppress thermalization, a phenomenon termed Stark
many-body localization (SMBL), previously observed in $S=1/2$ chains.
Furthermore, our findings reveal that the sole presence of single-ion
anisotropy is sufficient to prevent thermalization in the system.
Interestingly, when the magnitudes of the magnetic field and anisotropy are
comparable, they compete, favoring delocalization. Despite the potential
hindrance of SMBL by single-ion anisotropy in this scenario, it introduces an
alternative mechanism for localization. Our interpretation, considering local
energetic constraints and resonances between degenerate eigenstates, not only
provides insights into SMBL but also opens avenues for future experimental
investigations into the enriched phenomenology of disordered free localized
$S\geq 1/2$ systems.

##  I Introduction

In quantum mechanics, a given initial state of a closed interacting quantum
system will evolve unitarily according its Hamiltonian $H$. A fundamental
assumption in statistical physics is that a generic closed quantum many-body
systems thermalize under its own dynamics [1, 2]. This implies that, under
unitary time evolution, the reduced density matrix of any generic initial
state ${|{\Psi(0)}\rangle}$ tends to evolve towards the equilibrium Gibbs
state within that subsystem. Since any initial state can be spanned by
eigenstates of the Hamiltonian, then the reduced matrix constructed with the
eigenstates should also evolve toward equilibrium. This is the notion
underlying the known eigenstate thermalization hypotesis (ETH) [3, 4, 5]. At a
first glance, it seems that any out-of-equilibrium state evolves towards
equilibrium. However, this is not always the case. It was shown that disorder
can induce localization in a variety of interacting systems [6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], a phenomenon dubbed many-body-
localization (MBL), as a generalization of the well known Anderson
localization proposed by P. W. Anderson about fifty years ago [22, 23, 24].

The intriguing phenomenon of ergodicity break in many-body localized systems
has motivated great effort in the scientific community to understand the
mechanisms that lead to MBL [25, 26, 27, 28]. Indeed, the phenomena of
localization in interacting quantum systems is a puzzling problem still under
study [29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. From a theoretical point of
view, it is well established that the breakdown of the ETH in disorder induced
many-body localized systems can be captured by the level-spacing statistics
[2, 5, 39, 30] or entanglement properties of their eigenstates [15, 40, 10].
The ETH usually understood in two classes: strong sense, in which all states
thermalizes and weak sense, where almost all states thermalizes. The main
features of disorder-induced MBL phases have been observed in various
experimental platforms [41, 42, 43, 44, 45]. In spite of the complexities in
the MBL phenomena, theoretical and experimental results suggest that, at least
in one dimension, strong disorder induce the emergence of nearly conserved
local quantities, leading to integrability [46, 47, 48, 49, 50].

Recently, Pollman et al. showed that localization can occur in interacting
systems even in the absence of disorder [51]. The key ingredient here is the
presence of a nearly uniform gradient potential in an interacting system and
this is closely tied to the single-particle localization process known as
Wannier-Stark localization [52].This phenomenon was called Stark Many-Body
Localization (SMBL), and it shares similarities with the traditional MBL, such
as level statistics, but differs in aspects such as a strong dependence on the
initial conditions, as shown in [53] for a spin $1/2$ Heisenberg chain. The
mechanism leading to ETH violation in SMBL was believed to be Hilbert space
fragmentation [53, 54, 55, 56, 57], in which the Hilbert space fragments into
many disconnected subspaces, preventing themalization. However, this argument
have been questioned [58]. Striking experiments in quantum simulators were
performed to demonstrate Stark MBL [59, 60] and confirm that localization can
arise in disorder-free systems.

In a recent study, one of us have extended the results of Ref. [53] by
including exchange interaction up to second nearest neighbors [61]. It was
showed that the SMBL phenomenon is robust for arbitrary ration $J_{2}/J_{1}$,
a feature also found in spinless fermions with long-range interactions [62].
SMBL has been also investigated in the context of bosonic lattices [63, 64]
and it’s worth noting that other mechanisms can give rise to MBL, such as
quasi-periodic potentials[65, 66] and periodic driving[67, 68].

Nonetheless, it is remarkable that all existing studies on spin lattices have
exclusively focused on SMBL in spin $1/2$ systems, which does not allow
interactions between the systems’ elementary excitations, kinks or spinors.
Interactions arise naturally, in systems with spin $S\geq 1$, for instance,
magnon-magnon interactions are present in spin$-1$ FM chains, giving rise to a
nonlinear term $D\sum_{j=1}^{L}\left(S_{j}^{z}\right)^{2}$, where D represents
the uniaxial anisotropy strength. Interestingly, weak ergodicity breaking in
the form of many-body quantum scars appears in the spin$-1$ XY model [69]. We
should mention that investigations on magnon-bound states has been conducted
examining the impact of a weak anisotropy ($D$ much smaller than the exchange
interaction $J$) on thermalization time scales of spin$-1$ systems [70, 71].
Notably, tunable single-ion anisotropy have been realized experimentally in
trapped ions [72], ultra-cold atoms [73], and compounds such as NiNb2O6 [74]
and [Ni(HF2)($3-$Clpyradine)4]BF4 [75].

In this paper, we investigate localization of quantum states in a Majumdar-
Ghosh model [76] in the presence of a non-uniform magnetic field, and
incorporating an additional term to address uniaxial single-ion anisotropy,
which naturally arises in magnetic models with $S\geq 1/2$ [77, 78, 79]. By
employing numerically exact diagonalization, we calculate the temporal
evolution of complementary quantifiers, namely imbalance, entanglement
entropy, and participation entropy. These measures are commonly used to
characterize the ergodic-to-MBL phase diagram in quantum many-body systems.
Our results reveal that for some initial product states: (i) Similar to the
case of $S=1/2$ chain, MBL induced by finite $h_{0}$ is also capable of
localizing system of spin $S>1/2$; (ii) For $S\geq 1$, uniaxial anisotropy
alone localizes the system; (iii) If the system is localized for a particular
value of $h_{0}$ the presence both $D\sim h_{0}$ suppresses localization and
(iv) for $D,h_{0}\gg J$ the system localizes. These results suggests that the
localization phenomena for both $D$ and $h_{0}$ is of the same nature. Indeed,
they seem to result from a suppression of the spin dynamics by local energy
constraint. While our numerically exact approach is limited by the size of the
system, our results suggest a different way of inducing localization in a
quantum many-body system. Moreover, our results shed light on the nature of
SMBL currently under intense investigations and paves the way for future
experimental investigations of these phenomena.

The rest of this paper is organized as follow: In Sec. II we introduce the
model and numerical methods, in Sec. III we present our numerical results and
discussions. Finally, our work is summarized in Sec. IV

##  II Model and methods

For concreteness, we consider a spin chain modeled by an extended version of
Majundar-Ghosh Hamiltonian [80], which describes a physical system composed of
a one-dimensional chain of interacting spins. Explicitly, the Hamiltonian for
$L$ spins can be written as

$\displaystyle H=J_{1}\sum_{j=1}^{L-1}{\bf S}_{j}\cdot{\bf
S}_{j+1}+J_{2}\sum_{j=1}^{L-2}{\bf S}_{j}\cdot{\bf
S}_{j+2}+\sum_{j=1}^{L}\left(h_{j}+DS_{j}^{z}\right)S_{j}^{z},$ (1)

where the first two terms of our model Hamiltonian represent the exchange
interaction between neighboring spins, capturing the mutual influence of these
adjacent entities. $J_{1}$ and $J_{2}$ define the coupling strength between
nearest and next-nearest neighbors, respectively. These are the two terms that
corresponds to what is known as Majumdar-Ghosh Hamiltonian [76]. In the last
summation of Eq. (1), $h_{j}S_{j}^{z}$ (with $h_{j}=jh_{0}+\gamma
j^{2}/L^{2}$) introduces the effect of a non-uniform magnetic field along the
$z$-direction, allowing the study of Stark-many-body localization. Lastly, $D$
accounts for a single-ion anisotropy, which takes into account the directional
preference of spins along the z-axis. This particular form of the Hamiltonian
(1), provides a theoretical framework for investigating fundamental properties
of the system for spin $S=1/2$ and $S\geq 1$ as well. We should mention that
the effect of $D$ can only be observed in the dynamics of system of spins
$S>1/2$. In the special case of $S=1/2$, the effect of $D$ is just to shift
the spectrum of the Hamiltonian by an amount $DL/4$.

For convenience, following Burssil et al [81] and our previous work [61], we
parameterize the exchange couplings $J_{1}$ and $J_{2}$ as
$(J_{1},J_{2})=(J_{0}\cos\theta,J_{0}\sin\theta)$, with $0\leq\theta<2\pi$.
The ground state properties for $h_{0}=\gamma=0$ in Eq. (1) have received
great attention over the past decades [82, 83] However, our main focus here
are thermalization processes, which potentially involves all eigenstates of
the system consistent with the particular symmetry of the initial state. The
absence of translation-invariance due to a finite field gradient naturally
lead us to adopt open boundary conditions. This choice allows one to explore
static and dynamical properties of the system in its entirety, with the price
of finite sizes effects. Here we seek to understand how the system
thermalizes, i.e., how it evolves from a given initial non-equilibrium state
towards a thermal equilibrium state. More specifically, we will cover initial
states of the form ${|{\Psi_{0}}\rangle}=~{}\mid\cdots
m,m,m,\bar{m},\bar{m},\bar{m},\cdots\rangle$ in which $m=2S^{z}_{j}$ and
$\bar{m}=-m$. This particular type of states correspond to a state with a
single domains wall, but states with multiple domain walls will also be
considered. This type of initial state has already been analyzed in other
studies where it was possible to unveil several interesting characteristics of
Stark many-body localization phenomena [53, 61]. Here we will explore mainly
the case of $S>1/2$ by considering the effects of the single-ion anisotropy
$D$ and magnetic fields gradients $h_{0}$ in order to gain a better
understanding of both local and global effects over the collective behavior of
the spin dynamics.

###  II.1 Time-evolution analysis

To study the phenomena of localization and thermalization in our closed
quantum system, we adopt the Schrödinger representation and perform time
evolution of a given initial state
${|{\Psi_{0}}\rangle}\equiv{|{\Psi(t=0)}\rangle}$ at $t=0$ as
${|{\Psi(t)}\rangle}=U(t){|{\Psi_{0}}\rangle}$ where
$U(t)=\exp{\left(-iHt/\hbar\right)}$ is the time evolution operator. Having
the time-evolved quantum state ${|{\Psi(t)}\rangle}$, we are able to calculate
the relevant physical quantity that witnesses localization/thermalization
phenomenon, such as imbalance, participation and entanglement entropies. Since
they are in general individually inconclusive they are used complimentaryly.

#### Imbalance.–

Imbalance, denoted as ${\cal I}(t)$, is a key measure in investigating the
dynamics of quantum systems, particularly in the context of Heisenberg models
[53, 59, 84, 51]. It provides insights on whether the system retains or loses
information about local magnetization as it evolves in time. We define the
imbalance for a chain of $L$ spins $S$ as

$\displaystyle{\cal I}(t)=\frac{1}{LS^{2}}\sum_{j=1}^{L}\langle
S_{j}^{z}(t)\rangle\langle S_{j}^{z}(0)\rangle,$ (2)

where, ${\langle{S^{z}_{j}(t)}\rangle}=\langle\Psi(t)\mid
S_{j}^{z}\mid\Psi(t)\rangle$ represents the expectation value of $S^{z}$ of
the $j$-th spin of the chain within the evolving quantum state
$\mid\Psi(t)\rangle$. As such, we clearly have ${\cal I}(0)=1$ if local spins
are fully polarized, along either the positive or negative $S^{z}$
projections. Moreover, in the regime of strong localization, ${\cal I}(t)$
remains close to its value at t=0, indicating that information about local
magnetization is preserved during the evolution. On the other hand, in the
thermal regime, ${\cal I}(t)$ tends to some lower value at long time scales.
The precise asymptotic value of the imbalance in the thermal regime depends on
the initial state. For instance, let us assume an initial product state
${|{\Psi_{0}}\rangle}$ as defined above for a chain of $L$ spins, where
$N_{\uparrow}$ and $N_{\downarrow}$ represent the number of spins with maximum
and minimum $S^{z}$, respectively. Thus, the total spin projection along
$z$-axis is given by $S^{z}_{\rm tot}=(N_{\uparrow}-N_{\downarrow})S$. Since
in the thermal state, this quantity will be uniformly distributed along the
chain, each site will exhibit a magnetization $S^{z}_{\rm tot}/L$, if total
$S^{z}$ is conserved. With this, we can easily show that the thermal value for
the imbalance is given by

$\displaystyle{\cal I}_{\rm thermal}=\left(\frac{S^{z}_{\rm
tot}}{LS}\right)^{2}.$ (3)

In particular, for an initial state in the sector $S^{z}_{\rm tot}=0$ we
obtain ${\cal I}_{\rm thermal}=0$. We should remark that, as defined in Eq.
(2), the imbalance is not capable of detecting dynamics of spins that occurs
only on the $xy$-plane. In this case ${\cal I}(t)=0$, regardless the dynamics
of the initial state.

#### Entanglement entropy.–

Another important quantity to monitor the localization of quantum states is
the well known entanglement entropy $S_{\ell}(t)$ defined as

$\displaystyle S_{\ell}(t)=-\frac{1}{\ln 2}{\rm
Tr}\left[\rho^{\ell}_{A}(t)\ln\rho^{\ell}_{A}(t)\right],$ (4)

where $\rho^{\ell}_{A}(t)Tr_{B}[\rho^{\ell}_{B}]$ is the reduced density
matrix of a subsystem $A$. Here, given the 1D nature of a chain, $\ell$ simply
represents the site that defines the interface between two portions of the
chain which we call subsystem $A$ and $B$.

While the imbalance gives qualitative local data that is good for describing
the dynamics of individual states, the use of entanglement entropy for many-
body systems allows for the investigation of their collective and emergent
properties, providing insights into the nature of quantum correlations and how
information distributes throughout the entire system [85].

#### Participation entropy.–

Participation entropy (PE) has gained attention in studies of dynamics of
many-body systems. Here we follow Ref. [52] and defined this quantity as

$\displaystyle S_{2}=-\ln\left(\sum_{mn}\lvert c_{n}\rvert^{2}\lvert
c_{m}\rvert^{2}\delta_{E_{n},Em}\right),$ (5)

where $c_{n}={\langle{n}|}{{\Psi_{0}}\rangle}$ is the projection of the
initial state ${|{\Psi_{0}}\rangle}$ onto the eigenbasis
$\\{{|{n}\rangle}\\}$. For practical purpose, to calculate $S_{2}$ we first
define the quantity $f(t)={\langle{n}|}\exp(-iHt/\hbar){|{\Psi_{0}}\rangle}$,
which is projection of time evolved initial state onto the eigenstate
${|{n}\rangle}$. Upon Fourier transforming $f(t)$ to the energy domain, we
obtain we obtain $F(\omega)$. The inverse participation ratio (IPR) can be
obtained as

$\displaystyle IPR=\int F^{2}(\omega)d\omega=\sum_{mn}\lvert
c_{n}\rvert^{2}\lvert c_{m}\rvert^{2}\delta_{E_{n},Em}.$ (6)

With this, the participation entropy (5) is just $S_{2}=-\ln(IPR)$. For non
degenerate case, this quantity coincides with the quantity $S_{q}$ for $q=2$
defined in Ref. [86]. Defined as such, $S_{2}$ provides a measures of how the
projections of the initial state are distributed among the eigenstates of the
Hamiltonian. In the absence of symmetries, for a delocalized state,
$S_{2}\sim\ln({\cal N})$, where ${\cal N}$ the dimension of the Hilbert space.
This increases logarithmicaly with the system size. On the other hand, in a
strongly localized state, only a limited subset of coefficients contributes
significantly to the sum, resulting in a constant PE $S_{2}$. In other words,
the state is fully spanned on a restricted region of the Hilbert space,
indicating the presence of localization.

##  III Numerical results

To show our numerical results for the time evolution of the initial many-body
state, we first set $J_{0}=1$ as energy unit and assume $\hbar=1$, so that
time has units of inverse energy. Since the system is closed, exact numerical
unitary evolution can be performed for small chains. To this end, here we use
the python package Quspin [87, 88] that allows us to perform time evolution
and calculation all the physical quantities we need. Since the Hilbert space
increase as $(2S+1)^{L}$, where $L$ is the length of the chain, for a
spin-1/2, for instance we can handle chains up to $L=22$ spins or so,
depending on the symmetry of the relevant sector of the Hilbert space.

Figure 1: Localization in as spin $1/2$ chain, for different values of $h_{0}$
and $\theta$ with $D=0$, using the average of the imbalance $\langle{\cal
I}\rangle$ vs $h_{0}$ and the average entanglement entropy $\langle
S_{\ell}\rangle$ vs $h_{0}$. The system is configured in the space $S_{z}=0$
with $L=12$ and $L=20$ with a magnetization island of $1/2$ in the middle of
the chain and $-1/2$ at the ends.

###  III.1 Spin-$1/2$ system: effect of magnetic field

For the sake of completeness, we reproduce and elaborate on the main results
obtained by one of us in Ref. [61] for spin-1/2 SMBL. We validate that the
fundamental characteristics are qualitatively well captured at a scale
involving $12$ and $20$ spins. In Fig. 1 we use imbalance and entanglement
entropy to show how the system evolves from a thermalized to localized regime
as $h_{0}$ increases. We use an initial state within the Hilbert subspace with
$S_{\rm tot}^{z}=0$ consisting of an island of spin _up_ at the center of the
chain. For a given value of $h_{0}$ we perform time evolution up to
$J_{0}t=500$ just as in Ref. [61]. Since both ${\cal I}(t)$ and $S_{\ell}(t)$
oscillates around a fixed value for large $t$, we take their respective time
average, $\bar{\cal I}$ and $\bar{S}_{\ell}$, for $J_{0}t\in[400,500]$.

In Figs. 1(a) and 1(c), $\bar{\cal I}$ as a function of $h_{)}$ is depicted
for system sizes $L=12$ and $L=20$, respectively, for various value of
$\theta$. We observe that for $h_{0}=0$ $\bar{\cal I}\approx 0$ but increases
with increasing $h_{0}$ reaching almost the unity for $h_{0}=2$ for all value
of $\theta$. The vanishing $\bar{\cal I}$ for small $h_{0}$ reveals that the
local information of the initial state is lost at long time though the unitary
dynamics, while large value of $\bar{\cal I}$ indicates that the information
of the initial state is kept local at long times. There is an interesting
qualitative agreement between the results for $L=12$ and $L=20$, which shows
that localization can be observed for relatively small chains. Similarly the
entanglement entropy shown in Figs. 1(c) and 1(d) (for $L=12$ and $L=20$,
respectively) confirms localization in the system as $h_{0}$ increases. Since
the initial state is a product state, with entanglement entropy $S_{L}(t)=0$,
large values of $\bar{S}_{L}$ for $h_{0}\rightarrow 0$ indicate that the
system becomes entangled, signifying the spread of information throughout the
entire system in the thermalized regime. Conversely, vanishing $\bar{S}_{L}$
in the localized regime for large $h_{0}$ confirms that information is kept
local at long times. It is noteworthy to mention a symmetry pertaining to the
imbalances $\theta$ and $\theta+\pi$, as previously highlighted in Ref. [61].

In what follows we will turn our attention to systems with $S>1/2$ on which
$D$ exerts non-trivial effects. The natural progression leads us to the case
of $S=1$. However, in this instance, flipping a given spin from $S_{z}=-1$ to
$S_{z}=+1$ requires passing through $S_{z}=0$, which remains inert to the
magnetic field along the $z$-direction. To explore the interplay between
$h_{0}$ and $D$, the case of $S=1$ does not encompass the most general
scenario. Consequently, we shift our focus to the case of $S=3/2$, which
presents a richer set of dynamics. Further details regarding the $S=1$ case
can be found in Appendix A.

###  III.2 Spin-$3/2$ system: effect of magnetic field and single-ion
anisotropy

Figure 2: Imbalance as a function of time for the stark many-body localization
(for increasing $h_{0}$ and $D=0$) [panels (a)-(d)] and for $h_{0}=0$ and
increasind $D$ [panels (f)-(h)]. Different curves corresponds to different
values of $\theta$ (see legend). The system consists of $L=12$ spins and the
initial state corresponds to an island of of spins with $S^{z}_{j}=3/2$ in the
middle of the chain, while all other has $S^{z}_{j}=-3/2$.

The $S=3/2$ case offers a much richer Hilbert space as compared to $S=1/2$. In
this case each spins exhibit $S^{z}=(3/2,1/2,-1/2,-3/2)$ where the two _local_
subspaces defined by $|S_{j}^{z}|=1/2$ and $|S_{j}^{z}|=3/2$ are both active
to external magnetic field gradient $h_{0}$ and anisotropy $D$. This leads to
more complex dynamics of the system towards thermalization. We will first
study the localization in the presence of $h_{0}$ and contrast with the case
of $S=1/2$ and then investigate the effect of $D$ alone. In a third stage we
will study the interplay between $h_{0}$ and $D$ in the dynamics of the
system.

Following the same strategy used before, we will analyze the time-evolution
for an initial state where there is an island of $L/2$ spins with
$S_{j}^{z}=-3/2$ centered at $L/2$ while the rest of the system has
$S_{j}^{z}=3/2$. Like the previous initial state we used for $S=1/2$, this is
a non-equilibrium product state, hence non-entangled. The effect of $h_{0}$ in
the time evolution of this state is shown in Fig. 2(a)-2(d) for $h_{0}=0$,
$h_{0}=2.0$, $h_{0}=3.0$ and $h_{0}=4.0$, while keeping $D=0$. In Fig. 2(a),
corresponding to the case $h_{0}=0$, the Hamiltonian of the system is
exclusively governed by the Majumdar-Ghosh exchange terms. We observe rapid
thermalization in the system, leading to the decay of ${\cal I}(t)$. When the
magnetic field gradient is increased to $h_{0}=2$, as depicted in Fig. 2(b),
the system continues to undergo thermalization for all angles, albeit at a
slower pace for every $\theta$. Notably, ${\cal I}(t,\theta)={\cal
I}(t,\theta+\pi)$, a symmetry previously noted for spins $S=1/2$. Fig. 2(c)
and 2(d) exhibits the results $h_{0}=3$ and $h_{0}=4$, respectively. In both
situations, where the imbalances remain close to unity, indicating localized
states, the observed persistence of these imbalances in systems with spin
$3/2$ strongly suggests the occurrence of SMBL in lattices with higher-spin
configurations

We now turn our attention to the impact of single-ion anisotropy on the system
dynamics, utilizing the same initial states as in the previous analysis. In
Fig. 2(e)-2(h) we set $h_{0}=0$ and present the dynamics of ${\cal I}(t)$ for
$D>0$. Specifically, Fig. 2(e) displays the results for $D=0.5$ across various
values of $\theta$. All the curves exhibit rapid thermalization, except in the
specific case of $\theta=0$, where the imbalance remains finite for an
extended duration, up to $J_{0}t=500$. This observation strongly indicates
that anisotropy has the potential to induce localization, with a notable
dependence on the angle $\theta$. Indeed, for $D=4$, as shown in Fig. 2(h),
the largest $D$ considered here, the system remains localized, as all
imbalances remain close to unity. The localization observed for large $D$ can
be easily understood by observing that for $D\gg J$, the initial state becomes
_almost_ an eigenstate of the system, apart from some very small perturbation
introduced by $J_{1}$ and $J_{2}$. In this scenario, localization is trivially
achieved. The intriguing situation arises in the competing regime where $D\sim
J_{0}$. Let’s take a closer look at Figs. 2(f) (where $D=2$) in more detail.
Observe that the curves for $\theta=0$, $\theta=\pi/3$ and $\theta=5\pi/3$
remains close to unity, indicating strong localization. Conversely, for
$\theta=2\pi/3$, $\theta=\pi$ and $\theta=4\pi/3$ the system undergoes
thermalization. Note that $J_{2}$ is the same for both $\theta=0$ and
$\theta=\pi$, yet the dynamics is entirely distinct. This difference can be
attributed solely to the sign of $J_{1}$. In other words, thermalization is
primarily governed by nearest-neighbor spin-flip processes facilitated by
$J_{1}$. To comprehend the mechanisms leading to localization (thermalization)
$J>0$ ($J<0$) at large $t$, let us analyze the early stages of the demolishing
of the domain wall at a given interface between spins $S^{z}=+3/2$ and
$S^{z}=-3/2$. The very first process that destroys the interface is given by
the term $JS_{i}^{-}S_{i+1}^{+}$ acting on the interface ${|{\cdots
3,3,\bar{3},\bar{3}\cdots}\rangle}$, resulting in an smoother interface
${|{\cdots 3,1,\bar{1},\bar{3}\cdots}\rangle}$. Due to finite $D$, there is a
cost associated with smoothing the domain wall, as the local energy in sites
with $S_{j}^{z}=\pm 1/2$ significantly differs from that in sites with
$S_{j}^{z}=\pm 3/2$. This energy difference acts as a barrier that can only be
surpassed through virtual processes during the system’s time evolution.

Figure 3: Effect of a linearly varying magnetic field (a) and a constant
anisotropy (b) on the local energy of a spin-$S=3/2$ chain. (c) Representation
of the Hilbert space and typical initial state within the Hilbert space for a
spin-$3/2$ chain. White regions encloses only states in which all spins have
their spin projection on the $z$-axis is maximum or minimum, i. e., their
projections are $S_{j}^{z}\pm 3/2$. Darker region contains states in which
some of their spins has projections $S^{z}_{j}=\pm 1/2$. The red and violet
curves represent an initial states with single and double domains wall,
respectively. The dots marks the active regions upon time $t=0^{+}$ in the
time evolution operator $\exp(-iH0^{+}/\hbar)$.

To gain a deeper understanding of the physical origin of this localization,
one can explore the limit of very large $|D|$ or, equivalently, $J\rightarrow
0$. For simplicity, let us assume $D>0$. When $J=0$ our initial state is an
eigenstate of the Hamiltonian of highest energy, since all spins have the
maximum $|S^{z}_{j}|$. In fact, any product state constructed from the $S^{z}$
basis is an eigenstate of the Hamiltonian for $J=0$. The spectrum of the
Hamiltonian consists of a series of $delta$-peaks among which the initial
state corresponds to one of the states of highest energy (for $D>0$). For the
sake of clarity, in Figs. 3(a) and 3(b) we illustrate the energy cost
associated with placing a spin with different $S_{j}^{z}$ along the chain,
considering a constant-gradient magnetic field $h_{0}$ and a anisotropy $D$,
respectively. Figure 3(c) represents the entire Hilbert space of the system,
with the white region denoting the sector of the Hilbert space where the
states have spins fully aligned along the $z$-direction. Consider the simple
case of the initial state
${|{\cdots\bar{3},\bar{3},\bar{3},3,3,3\cdots}\rangle}.$ This state exhibits a
single domain wall, represented by the red line within the sector of the
Hilbert space indicated by the white region. The point where the curve touches
the border of the white region represents the domains wall. At $t=0^{+}$, only
the region represented by the red line is active. The purple line indicates a
state with a double wall, featuring two active regions. In the vicinity of the
white sector, there are states where the two spins around the wall have
projections $\pm 1/2$, such as the state
${|{\cdots\bar{3},\bar{3},\bar{1},1,3,3\cdots}\rangle}$. One can readily
confirm that for $J=0$, this state is trivially localized, but for finite
$J=0$, its temporal evolution becomes highly complex. Consider that at $t=0$,
the system is prepared in the initial state with a single domain wall. At
$t=0^{+}$ the only active region of the state is around the domain wall, since
$H$ acting on this state modifies only the spins at the vicinity of the domain
wall. Now, upon time evolution, in order to reach a thermal situation, the
mechanism involves states with $S_{j}^{z}\pm 1/2$ in the vicinity of the wall
(represented by the dark gray region). It turns out that, for $D>0$ theses
intermediate states has of much different energies for $D\gg J$. Therefore,
for a given $J$ fixed, by increasing $D$ these sector becomes isolated from
the rest of the Hilbert space.

In Figure 4(a)-4(d), we present the density of states for a system of $L=8$
spins $S=3/2$ for various values of $J$ and fixed $D=0.5$. In the absence of
$J$($J=0$) in 4(a), the initial state precisely corresponds to one of the
eigenstates associated with the rightmost degenerate energy peak. As $J$
increases, these peaks hybridize with other states, resulting in peak
splitting. For instance, at $J=0.1$, the peaks are split due the the
hybridization with the rest of the band. Therefore, the effect of $D>0$ is
evident in pushing the initial state to higher energy, creating a ”gap”
between this state and the rest. Consequently, the entire white region of the
Hilbert space, as illustrated in 3(c), becomes isolated from the rest. In
contrast to QMBS, where the Hilbert space fragmentation is driven by emergent
conserved quantities, the observed suppression of thermalization in this
system arises from the inhibition of spin flip processes due to energy
constraints.

Figure 4: Density of states as a function of energy for $D=0.5$ and $\theta=0$
for $J=0$ (a), $J=0.1$ (b), $J=1$ (c) and $J=0.4$ (d). For all of them, we use
the subspace of $S_{z}=0$, with the spin-up magnetization island in the middle
for L = 8.

The picture drown above is useful to understand why for a moderate value of
$D$, such as $D=2$ and $J=1$ shown in Fig. 2(f), by changing $\theta=0$ to
$\theta=\pi$ (or equivalently $J$ to $-J$), the system goes from localized to
thermalized. The reason is that, spectrum is not symmetric when
$J\rightarrow-J$. In fact, for $J>0$ the spectrum extends more towards
negative side of the energy axis. The opposite occurs $J<0$. Therefore, for a
$D$ fixed, the energy “gap” between the sector of the initial state is larger
for $J>0$, thus localization is obtained more easily. The evidence that
increasing $D$ suppresses dynamics of spins around the domain walls, places
the localization mechanism here in the same perspective of the SMBL. As
discussed in Ref. [58], in SMBL, only for infinite $h_{0}$ the spins a
completely frozen by “local energy” constraint. Likewise, here, the spin
dynamics will be fully frozen only at infinite $D$. Nevertheless, our results
suggest localization already for $D\sim 2J$ as in Fig. 2(b).

Since the demolishing of domain walls of the initial states occurs because of
the active regions around the interfaces, let us analyse how the number of
walls modifies the localization/thermalization. In Fig. 5 shows the time
evolution of the imbalance starting at an initial states containing one, two,
three and five domains walls, all states within the zero total magnetization.
Again, we use a chain of $L=12$ spins $S=3/2$ and set $h_{0}=0$. Panels 5(a)
and Fig. 5(b) show the case of $D=0$ and $\theta=0$ and $\theta=2\pi/3$. For
this we observe thermalization for all initial states shown. Now, for $D=4$
(lower panels), while for $\theta=0$ the system is still localized, as shown
in Fig. 5(c), for $\theta=2\pi/3$ localization becomes poorer as the number of
walls increases. For instance, for $\theta=2\pi/3$, note that for five
interfaces the system thermalizes, while remains localized for one and two
walls. This is because, the number of active regions increases with the number
of domain walls as illustrated in Fig. 3(a), thus favoring thermalization.
This shows that the larger is the number of interfaces, the larger is the
value of $D$ and $h_{0}$ necessary to produce the same localization, similar
to what was previously reported in Refs [53] for SMBL.

Figure 5: Time evolution of imbalance for different initial conditions
(different curves). Upper and lower panels correspond to $D=0$ and $D=4$,
while left and right panels correspond to $\theta=0$ and $\theta=2\pi/3$. For
all panels use a chain of spin $3/2$ and $L=12$ and $h_{0}=0$. All initial
states here lies on the $S_{z}=0$ subspace.

###  III.3 Interplay between anisotropy and magnetic field

Figure 6: Time evolution of imbalance (a)-(d) and entropy (e)-(h) for
different values of $h_{0}$ and $D$. The system and initial condition here is
the same as the one used in Fig. 2. The curves corresponds to different values
of $\theta$ (see legend).

Figure 7: (a) Participation entropy vs $\ln{\cal N}$ for $(D,h_{0})=(0,0)$
(red circles), $(D,h_{0})=(8,0)$ (black squares) and $(D,h_{0})=(0,8)$ (cyan
triangles). Dashed lines corresponds to a linear fitting of the curces and
serves to guide the eyes. To change the size of the Hilbert space, $L$
includes all even numbers form $4$ to $14$. (b) and (c) shows, respectively,
latter time average of imbalance and entropy as function of $D$ (for
$h_{0}=0$) (red squares) and $h_{0}$ (for $D_{0}=0$) (blue triangles) for
$L=12$. For all panels, we set $\theta=0$ ans use an initial state of type
${|{\cdots,\bar{3},3,3,\bar{3},\cdots}\rangle}$, which corresponds to an
island of spins $S_{j}^{z}=3/2$ in the middle of a chain of spins
$S_{j}^{z}=-3/2$.

So far we have analysed separately the SMBL, induce by the magnetic field
gradient, and the localization induced by anisotropy. Let us now address the
interplay between these two terms of the Hamiltonian. In Fig. 6(a)-6(d) we
repeat the same calculation of Fig. 2, but now for both $h_{0}$ and $D$
finite. Figure 6(a) shows the imbalance for $h_{0}=D=2$. We observe that while
localization were observed for some values of $\theta$ for $D=2$ and $h_{0}=0$
shown in Fig.2(f), now localization is no longer observed for any value of
$\theta$. Interestingly, if we make $h_{0}=2D=4$, as shown in Fig. 6(b),
localization is again recovered. However, for $D=2h_{0}=4$, localization is
observed only for some values of $\theta$, in fact for those closer to
$\theta=0$ of $\theta=2\pi$. Now, for $\Delta=h_{0}=4$, all the curves
approach to unity at large $t$. Despite the complexity of the behavior of the
various curves, results indicate that for moderate values of $D\approx h_{0}$,
they compete against each other and fail to localize the system. However, if
both are much larger than $J$ localization can be recovered again. The
competition between $D$ and $h_{0}$ observed here can be understood as
follows: the energy price that strongly suppresses spin dynamics around the
domain wall for $D$ finite (and $h_{0}=0$) is now, to some extent, compensated
by finite $h_{0}$. Likewise, the suppression of the dynamics for finite
$h_{0}$ (and $D=0$) is now energetically relaxed if $D$ is comparable to
$h_{0}$.

In Fig. 7(a) we show the participation entropy as function of ${\cal N}$ (the
dimension of the Hilbert space for $(D,h_{0}=0)=(0,0)$ (red circles),
$(D,h_{0})=(8,0)$ (black squares), and $(D,h_{0})=(0,8)$ (cyan triangles).
Here we use an initial state consisting of chains of even number of sites from
$L=4$ to $L=14$. The initial states are of the type
${|{\cdots\bar{3},3,3,\bar{3}\cdots}\rangle}$. For all system size, the
initial state contains an island of only two spin $S_{j}^{z}=3/2$. For
$(D,h_{0}=0)=(0,0)$ we note the $S_{2}$ increases linearly with $\ln{\cal N}$,
as expected for an ergodic regime. In contrast, for both cases
$(D,h_{0})=(0,8)$ and $(D,h_{0})=(8,0)$, $S_{2}$ remains nearly constant. This
indicates that for the ergodic regime, the number of eigenstates that
participates on the dynamics increases with the dimension of the Hilbert
space. In the localized regime, on the other hand, only a limited portion of
the eigenstates contributes to the dynamics.

Figure 7(b) shows $\bar{\cal I}$ vs $D$ for $h_{0}=0$ (green squares) and
$\bar{\cal I}$ vs $h_{0}$ for $D=0$. We use a chain of $L=12$ and the same
initial state as in Fig. 7(a) and set $\theta=0$. Again, $\bar{\cal I}$ is the
average of ${\cal I}(t)$ for $t\in[400,500]$. We note that both curves evolves
quite nicely from a thermalized ($\bar{\cal I}\approx 0.5$) to a localized
($\bar{\cal I}\approx 1$) regime as $D$ and $h_{0}$ increases from zero to
$8$. Interestingly, the crossover occur for similar values of $D$ and $h_{0}$.
This is corroborated by the entropy $\bar{S}_{L}$ depicted in Fig 7(c), where
it is evident that ${\bar{S}}_{L}$ vanishes if either $h_{0}$ and $D$
increase, indicating the localization of the states.

To further elucidate the numerical results, we employ the Holstein-Primakoff
(HP) transformation; mapping spin operators in a lattice system with spin$-S$
to bosonic operators: $S_{j}^{z}=(S-n_{j})$, $S_{j}^{+}=\sqrt{2S-n_{j}}a_{j}$,
and $S_{j}^{-}=a_{j}^{\dagger}\sqrt{2S-n_{j}}a_{j}$, where
$a_{j}^{\dagger}(a_{j})$ is a bosonic creation (annihilation) and the number
$n_{j}=a_{j}^{\dagger}a_{j}$ operator at site $j$. Subsequently, a $1/S$
expansion enables us to derive an effective description of the quantum
dynamics, expressed as $H_{HP}=H_{G}+H_{1}+H_{2}+\mathcal{O(}1/S)$. The term
$H_{G}$ characterizes the classical ground-state energy, while $H_{1}$
represents the magnon dispersion, both contingent on the parameters
$\left\\{J_{1},J\,_{2},D,h_{0}\right\\}$, yet not influencing localization.
The leading contribution for comprehending the interplay between $h_{0}$ and
$D$ emerges from $H_{2},$ since this term and higher-order ones, describes
magnon-mangnon interactions. To simplify, we consider only nearest-neighbor
interaction, i.e., $J_{1}=J>0,J_{2}=0$ and $D>0$. The HP representation leads
to the effective Hamiltonian

$\displaystyle H_{\rm HP}$ $\displaystyle\approx\sum_{j=1}^{L}\left(2DS+2JS-
h_{j}\right)n_{j}-JS\sum_{j=1}^{L-1}(a_{j}^{\dagger}a_{j+1}+H.c.)$ (7)
$\displaystyle+\frac{J}{8}\sum_{j}\left[\left(a_{j}^{\dagger}\right)^{2}a_{j}a_{j+1}+\left(a_{j+1}^{\dagger}\right)^{2}a_{j+1}a_{j}+H.c.\right]$
$\displaystyle+D\sum_{j}^{L}n_{j}(n_{j}-1)+\frac{J}{2}\sum_{j=1}n_{j}n_{j+1}+H_{G}.$

The form of $H_{\rm HP}$ is similar to the generalized Bose-Hubbard model in a
tilted optical lattice, with on-site chemical potential $\mu_{j}=2DS+2JS-
h_{j}$ (see Ref. [63]). The single-ion anisotropy $D$ is mapped into a
nonlinear repulsive on-site interaction $U$, symbolizing the additional energy
needed to accommodate more than one boson at each location. In the context of
a tilted lattice, where $\mu_{j}=jA$, in which $A$ represents the strength of
the tilt, the standard Bose-Hubbard model exhibits SMBL [63, 64], closely
resembling the behavior observed in our numerical results. Specifically, at
large values of $\ U$, high-energy states exhibit localized features at all
values of $A$. This behavior is inherently connected to the segmentation of
the Hilbert space into subbands, a characteristic that becomes more pronounced
in regions with strong interactions. It leads to the freezing of dynamics on
long time scales for initial states above a certain energetic threshold [89].
Another parallel emerges in the interplay between $U$ and $\mu_{j}$, where an
increase in $U$ corresponds to an increase in the critical Stark field $A_{c}$
necessary for the manifestation of the localization phenomenon. Additionally,
the symmetry noted in our results, encompassing both positive and negative
single-ion anisotropy $D$, is reflected in the bosonic scenario as an exact
symmetry between two Hamiltonians with the sign of $U$ changed and the energy
ordering of eigenvalues and eigenvectors reversed.

##  IV Conclusions

In summary, we have investigated the dynamics of a closed spin chain modeled
as a Majumdar-Ghosh Hamiltonian in the presence of finite-gradient magnetic
field $h_{0}$ inducing SMBL. An additional term accounting for single-ion
uniaxial anisotropy $D$ was included. Employing exact numerical calculations,
we compute various quantities, such as imbalance and entanglement entropy,
which are commonly used to probe localization phenomena. Our results show that
the known SMBL observed for $S=1/2$ chains is observed for $S>1/2$. Similarly,
in the absence of magnetic field, the presence of anisotropy alone is able to
induce localization in the system for $S>1/2$ for certain initial product
states where the spins have maximum $S_{j}^{z}$ projection. Our results
suggest that the localization induced by $D$ results from the suppression of
spin dynamics across the well defined domain walls of the initial states. This
can be understood as a energy barrier imposed by the anisotropy that has to be
overcome in order to scramble the spins from an product initial state where
all spins have their maximum or minimum $S^{z}_{j}$ projections. Even though
the sector of the Hilbert space containing these family of product states are
not completely decoupled from the rest under time evolution, they become
progressively isolated as $D$ increases. This is somewhat similar to what was
pointed out by Zisling et al, in Ref. [58] in the context of SMBL. More
interesting, we found that for $D\approx h_{0}\gtrsim J$ the magnetic field
and anisotropy compete to localized the system, favoring delocalization. We
interpret this as an interplay between energy costs payed to melt the domain
walls. In the competing regime, the energetic price charged by $D$ is somehow
compensated by the cost imposed by the $h_{0}$. While our exact numerical
calculations are currently applicable to few-site systems, we contend that
these initial findings hold relevance in the broader context of the dynamics
of closed many-body quantum systems. Particularly, they contribute to our
understanding of the connections between Hilbert space fragmentation and
disordered free localized systems. It is important to highlight that our
results do not rule out the possibility, for example, of the gradual
disappearance of localized regions in the Hilbert space in large systems or in
the long-time limit. Addressing these questions demands further exploration,
both experimentally using probes commonly employed in conventional MBL studies
and theoretically through investigations that could, for instance, extend the
theory of dynamical l-bits [90] to spin lattices with $S\geq 1$. Finally, it
is worth noting that the models examined in this study can be experimentally
realized using current quantum simulators, including cold atoms, trapped ions,
and superconducting qubits.

###### Acknowledgements.

The authors acknowledge financial support from CAPES, FAPEMIG and CNPq. EV
thanks FAPEMIG (Process PPM-00631-17) and CNPq (Process 311366/2021-0). EV
also acknowledge support and hospitality from the Zhejiang Normal University,
where part of this work was conducted. This work used resources of the “Centro
Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP).”

##  Appendix A Spin-$1$ chain

Figure 8: Imbalance (upper) and entanglement entropy (lower) as function of
time for a spin-$1$ chain for $D=0$ (left) and $D=6$ (right). The curves
correspond to different values of $\theta$. We used a chain of $L=12$ spins
with an initial state consisting of a island of six spins _up_ in the middle
of the chain, while all other spins are _down_ (see text).

Here, we briefly aim to demonstrate that the anisotropy $D$ also prevents
thermalization of the $S=1$ spin chain. We consider a chain of $L=12$ and set
$h_{0}=0$. Figure 8, the upper and lower panels display the imbalance and
entropy, respectively. We use an initial state
${|{\cdots\bar{2},2,2,2,2,2,2,\bar{2}\cdots}\rangle}$ and present the results
for various angles $\theta$. In Figs. 8(a) and 8(b) we show ${\cal I}(t)$ and
${\cal S}_{L}$ for $D=0$. It is noteworthy that the imbalance rapidly drops to
zero, while ${\cal S}_{L}(t)$ indicates the system’s tendency to thermalize
during time evolution. In Figs. 8(c) and 8(d), the corresponding results for
$D=6$ are presented. In this scenario, we observe that the imbalance remains
close to unity, while the entropy consistently stays close to zero. This
observation indicates that, in the presence of $D=6$, the system is localized.

## References

  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239–362 (2016).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Tasaki [1998] H. Tasaki, Phys. Rev. Lett. 80, 1373 (1998).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854–858 (2008).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498–525 (2018).
  * Chandran  _et al._ [2014] A. Chandran, V. Khemani, C. R. Laumann, and S. L. Sondhi, Phys. Rev. B 89, 144201 (2014).
  * Basko  _et al._ [2006] D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126–1205 (2006).
  * Zangara  _et al._ [2013] P. R. Zangara, A. D. Dente, A. Iucci, P. R. Levstein, and H. M. Pastawski, Phys. Rev. B 88, 195106 (2013).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Vasseur and Moore [2016] R. Vasseur and J. E. Moore, Journal of Statistical Mechanics: Theory and Experiment 2016, 064010 (2016).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. c. v. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
  * Oganesyan  _et al._ [2009] V. Oganesyan, A. Pal, and D. A. Huse, Phys. Rev. B 80, 115104 (2009).
  * Vasquez  _et al._ [2009] L. Vasquez, K. Slevin, A. Rodriguez, and R. Roemer, Annalen der Physik 521, 901–904 (2009).
  * Ioffe and Mézard [2010] L. B. Ioffe and M. Mézard, Phys. Rev. Lett. 105, 037001 (2010).
  * Aleiner  _et al._ [2010] I. L. Aleiner, B. L. Altshuler, and G. V. Shlyapnikov, Nature Physics 6, 900–904 (2010).
  * Monthus and Garel [2010] C. Monthus and T. Garel, Phys. Rev. B 81, 134202 (2010).
  * Berkelbach and Reichman [2010] T. C. Berkelbach and D. R. Reichman, Phys. Rev. B 81, 224429 (2010).
  * Anderson [1958] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Lagendijk  _et al._ [2009] A. Lagendijk, B. van Tiggelen, and D. S. Wiersma, Physics Today 62, 24 (2009).
  * Kramer and MacKinnon [1993] B. Kramer and A. MacKinnon, Reports on Progress in Physics 56, 1469 (1993).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Altman and Vosk [2015] E. Altman and R. Vosk, Annual Review of Condensed Matter Physics 6, 383 (2015).
  * Abanin and Papić [2017] D. A. Abanin and Z. Papić, Annalen der Physik 529, 1700169 (2017).
  * Smith  _et al._ [2016a] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nature Physics 12, 907–911 (2016a).
  * Šuntajs  _et al._ [2020] J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. E 102, 062144 (2020).
  * Huang [2022] Y. Huang, Annals of Physics 438, 168761 (2022).
  * Luitz [2016] D. J. Luitz, Phys. Rev. B 93, 134201 (2016).
  * Yu  _et al._ [2016] X. Yu, D. J. Luitz, and B. K. Clark, Phys. Rev. B 94, 184202 (2016).
  * Abanin  _et al._ [2021] D. Abanin, J. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. Parameswaran, F. Pollmann, A. Potter, M. Serbyn, and R. Vasseur, Annals of Physics 427, 168415 (2021).
  * Panda  _et al._ [2020] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, EPL (Europhysics Letters) 128, 67003 (2020).
  * Sierant  _et al._ [2020a] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. Lett. 124, 186601 (2020a).
  * Sierant  _et al._ [2020b] P. Sierant, M. Lewenstein, and J. Zakrzewski, Phys. Rev. Lett. 125, 156601 (2020b).
  * Kiefer-Emmanouilidis  _et al._ [2020] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. Lett. 124, 243601 (2020).
  * Kiefer-Emmanouilidis  _et al._ [2021] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. B 103, 024203 (2021).
  * Mondaini  _et al._ [2018] R. Mondaini, K. Mallayya, L. F. Santos, and M. Rigol, Phys. Rev. Lett. 121, 038901 (2018).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Smith  _et al._ [2016b] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nature Physics 12, 907 (2016b).
  * Xu  _et al._ [2018] K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song, W. Liu, Q. Guo, P. Zhang, D. Xu, H. Deng, _et al._ , Phys. Rev. Lett. 120, 050507 (2018).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016).
  * Lüschen  _et al._ [2017] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Phys. Rev. Lett 119, 260401 (2017).
  * Gong  _et al._ [2021] M. Gong, G. D. de Moraes Neto, C. Zha, Y. Wu, H. Rong, Y. Ye, S. Li, Q. Zhu, S. Wang, Y. Zhao, _et al._ , Phys. Rev. Res 3, 033043 (2021).
  * Serbyn  _et al._ [2013] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * Rademaker and Ortuño [2016] L. Rademaker and M. Ortuño, Phys. Rev. Lett. 116, 010404 (2016).
  * Imbrie [2016] J. Z. Imbrie, Journal of Statistical Physics 163, 998–1048 (2016).
  * O’Brien  _et al._ [2016] T. E. O’Brien, D. A. Abanin, G. Vidal, and Z. Papić, Phys. Rev. B 94, 144208 (2016).
  * Schulz  _et al._ [2019] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Phys. Rev. Lett. 122, 040606 (2019).
  * van Nieuwenburg  _et al._ [2019] E. van Nieuwenburg, Y. Baum, and G. Refael, Proceedings of the National Academy of Sciences 116, 9269 (2019).
  * Doggen  _et al._ [2021] E. V. H. Doggen, I. V. Gornyi, and D. G. Polyakov, Phys. Rev. B 103, L100202 (2021).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Phys. Rev. B 101, 174204 (2020).
  * Sala  _et al._ [2020] P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Phys. Rev. X 10, 011047 (2020).
  * Taylor  _et al._ [2020a] S. R. Taylor, M. Schulz, F. Pollmann, and R. Moessner, Phys. Rev. B 102, 054206 (2020a).
  * Herviou  _et al._ [2021] L. Herviou, J. H. Bardarson, and N. Regnault, Phys. Rev. B 103, 134207 (2021).
  * Zisling  _et al._ [2022] G. Zisling, D. M. Kennes, and Y. Bar Lev, Phys. Rev. B 105, L140201 (2022).
  * Morong  _et al._ [2021] W. Morong, F. Liu, P. Becker, K. S. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, and C. Monroe, Nature 599, 393 (2021).
  * Guo  _et al._ [2021] Q. Guo, C. Cheng, H. Li, S. Xu, P. Zhang, Z. Wang, C. Song, W. Liu, W. Ren, H. Dong, R. Mondaini, and H. Wang, Phys. Rev. Lett. 127, 240502 (2021).
  * Vernek [2022] E. Vernek, Phys. Rev. B 105, 075124 (2022).
  * Jiang  _et al._ [2023] X.-P. Jiang, R. Qi, S. Yang, Y. Hu, and G. Yang, arXiv preprint arXiv:2307.12376 (2023).
  * Yao and Zakrzewski [2020] R. Yao and J. Zakrzewski, Phys. Rev. B 102, 104203 (2020).
  * Taylor  _et al._ [2020b] S. R. Taylor, M. Schulz, F. Pollmann, and R. Moessner, Phys. Rev. B 102, 054206 (2020b).
  * Zhang and Yao [2018] S.-X. Zhang and H. Yao, Phys. Rev. Lett. 121, 206601 (2018).
  * Singh  _et al._ [2021] H. Singh, B. Ware, R. Vasseur, and S. Gopalakrishnan, Phys. Rev. B 103, L220201 (2021).
  * Bairey  _et al._ [2017] E. Bairey, G. Refael, and N. H. Lindner, Phys. Rev. B 96, 020201 (2017).
  * Yousefjani  _et al._ [2023] R. Yousefjani, S. Bose, and A. Bayat, Phys. Rev. Res. 5, 013094 (2023).
  * Schecter and Iadecola [2019] M. Schecter and T. Iadecola, Phys. Rev. Lett. 123, 147201 (2019).
  * Wu  _et al._ [2022] N. Wu, H. Katsura, S.-W. Li, X. Cai, and X.-W. Guan, Phys. Rev. B 105, 064419 (2022).
  * Sharma  _et al._ [2022] P. Sharma, K. Lee, and H. J. Changlani, Phys. Rev. B 105, 054413 (2022).
  * Senko  _et al._ [2015] C. Senko, P. Richerme, J. Smith, A. Lee, I. Cohen, A. Retzker, and C. Monroe, Phys. Rev. X 5, 021026 (2015).
  * Chung  _et al._ [2021] W. C. Chung, J. de Hond, J. Xiang, E. Cruz-Colón, and W. Ketterle, Phys. Rev. Lett. 126, 163203 (2021).
  * Chauhan  _et al._ [2020] P. Chauhan, F. Mahmood, H. J. Changlani, S. M. Koohpayeh, and N. P. Armitage, Phys. Rev. Lett. 124, 037203 (2020).
  * Pajerowski  _et al._ [2022] D. M. Pajerowski, A. P. Podlesnyak, J. Herbrych, and J. Manson, Phys. Rev. B 105, 134420 (2022).
  * Majumdar and Ghosh [2003] C. K. Majumdar and D. K. Ghosh, Journal of Mathematical Physics 10, 1388 (2003).
  * Li  _et al._ [2019] Y. Li, Z. Jiang, J. Li, S. Xu, and W. Duan, Phys. Rev. B 100, 134438 (2019).
  * Liu  _et al._ [2014] J. Liu, H.-J. Koo, H. Xiang, R. K. Kremer, and M.-H. Whangbo, The Journal of Chemical Physics 141, 124113 (2014).
  * Nardelli  _et al._ [2020] F. Nardelli, S. Borsacchi, L. Calucci, E. Carignani, F. Martini, and M. Geppi, Rendiconti Lincei. Scienze Fisiche e Naturali 31, 999–1010 (2020).
  * Bonechi  _et al._ [1992] F. Bonechi, E. Celeghini, R. Giachetti, E. Sorace, and M. Tarlini, Journal of Physics A: Mathematical and General 25, L939 (1992).
  * Bursill  _et al._ [1995] R. Bursill, G. A. Gehring, D. J. J. Farnell, J. B. Parkinson, T. Xiang, and C. Zeng, Journal of Physics: Condensed Matter 7, 8605 (1995).
  * Kumar [2002] B. Kumar, Phys. Rev. B 66, 024406 (2002).
  * Vernek  _et al._ [2020] E. Vernek, O. Ávalos-Ovando, and S. E. Ulloa, Phys. Rev. B 102, 174427 (2020).
  * 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, Science 349, 842 (2015).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Reviews of Modern Physics 82, 277 (2010).
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Phys. Rev. Lett. 123, 180601 (2019).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, SciPost Phys. 2, 003 (2017).
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, SciPost Phys. 7, 020 (2019).
  * Carleo  _et al._ [2012] G. Carleo, F. Becca, M. Schiró, and M. Fabrizio, Scientific Reports 2, 243 (2012).
  * Gunawardana and Buča [2022] T. M. Gunawardana and B. Buča, Phys. Rev. B 106, L161111 (2022).