# On intermediate statistics across many-body localization transition

Bitan De  Institute of Theoretical Physics, Jagiellonian University in Kraków,
Łojasiewicza 11, 30-348 Kraków, Poland  Piotr Sierant  The Abdus Salam
International Center for Theoretical Physics, Strada Costiera 11, 34151,
Trieste, Italy  Institute of Theoretical Physics, Jagiellonian University in
Kraków, Łojasiewicza 11, 30-348 Kraków, Poland  Jakub Zakrzewski
<EMAIL_ADDRESS>Institute of Theoretical Physics, Jagiellonian
University in Kraków, Łojasiewicza 11, 30-348 Kraków, Poland  Mark Kac Complex
Systems Research Center, Jagiellonian University in Krakow, Kraków, Poland.

(August 27, 2024)

###### Abstract

The level statistics in the transition between delocalized and localized
phases of many body interacting systems is considered. We recall the joint
probability distribution for eigenvalues resulting from the statistical
mechanics for energy level dynamics as introduced by Pechukas and Yukawa. The
resulting single parameter analytic distribution is probed numerically via
Monte Carlo method. The resulting higher order spacing ratios are compared
with data coming from different quantum many body systems. It is found that
this Pechukas-Yukawa distribution compares favorably with $\beta$–Gaussian
ensemble – a single parameter model of level statistics proposed recently in
the context of disordered many-body systems. Moreover, the Pechukas-Yukawa
distribution is also only slightly inferior to the two-parameter $\beta$-h
ansatz shown earlier to reproduce level statistics of physical systems
remarkably well.

##  I Introduction

Fritz Haake made fundamental contributions first to quantum optics (see e.g.
Bonifacio  _et al._ (1971); Haake  _et al._ (1979)) then to the fast
developing in the eighties of the last century area of Quantum Chaos not only
with his seminal monograph on the subject Haake (2010) but also with many
original works from introducing a celebrated kicked top model Haake  _et al._
(1987) to providing a link between Gutzwiller’s periodic orbit theory
Gutzwiller (1971) and random matrix statistics Müller  _et al._ (2004, 2005);
Heusler  _et al._ (2007).

Those late works provided a highlight of Haake’s fascination of the link
between spectral properties of physical models and random matrix theory. Of
particular interest for him has always been level dynamics (see Chapter.6 of
Haake (2010)). The dependence of the energy levels of a Hamiltonian
$H(\lambda)$ on some parameter $\lambda$ may be viewed as the motion of
interacting fictitious particles (levels) with “time” $\lambda$ following the
original formulation of Pechukas Pechukas (1983) and Yukawa Yukawa (1985). For
large matrices of size $N\times N$, in $N\rightarrow\infty$ limit, statistical
mechanics is applicable to those fictitious particles. Haake devoted
particular interest to periodically driven (Floquet) systems also in this
context, in particular since a transition between level clustering for
integrable systems to level repulsion (for quantally chaotic system) may be
viewed as a relaxation toward equilibrium Haake and Kuś (1988).

For this transition, following the original Wigner $2\times 2$ matrices
approach, Lenz and Haake found an interpolating spacing distribution Lenz and
Haake (1991) which compared well with random matrix simulations adding a
significant contribution to the topic which originated with the early work of
Rosenzweig and Porter Rosenzweig and Porter (1960). The Rosenzweig-Porter
distribution was an important twist on random matrices. As it is well known
Mehta (1990); Haake (2010); Stöckmann (1999) for generalized time reversal
invariant systems (the case we shall solely concentrate on) the Gaussian
orthogonal ensemble of random matrices (GOE) may be constructed by considering
matrices with independent Gaussian distributed entries with the variance of
the diagonal elements being twice the variance of the off-diagonal elements.
GOE ensemble faithfully represent statistical properties of quantally chaotic
system as conjectured by Bohigas, Giannoni, and Schmidt Bohigas  _et al._
(1984). The Rosenzweig-Porter ensemble interpolates between GOE and Poisson
cases using a single parameter $\chi=\sigma^{2}/N$, where $\sigma^{2}$ is the
variance of off-diagonal elements and $N$ the matrix size Rosenzweig and
Porter (1960). Lenz and Haake solution Lenz and Haake (1991) gives the nearest
spacing distribution for such a $2\times 2$ matrices ensemble.

This is by no means a unique solution for the transition between GOE and
Poisson level statistics. The different models that were proposed over the
years are reviewed in Section II. The renewed interest in the problem comes
from a realization that the transition between spectral properties of ergodic
many-body interacting systems and the so called many-body localized phase
Nandkishore and Huse (2015); Alet and Laflorencie (2018); Abanin  _et al._
(2019) may be described in the statistical sense by an ensemble interpolating
between GOE and Poissonian case as in the case of chaotic systems or single
particle Anderson localization. In Section III we review the statistical
approach of Yukawa Yukawa (1985) following the excellent textbook of H.-J.
Stockman Stöckmann (1999). Such a Pechukas-Yukawa statistical mechanics model
provides another, to the best of our knowledge not tested yet, proposition for
the interpolating ensemble. This approach has certain aesthetic advantages
over other propositions opening also possibilities of further studies going
beyond the eigenvalue statistics. A critical comparison of this ensemble
predictions with numerical data obtained for disordered Heisenberg chain and
further examples for other systems are the content of Section IV. We conclude
giving future perspectives of the presented approach in Section VI.

##  II The interpolating ensembles and formulae

The early approaches to level statistics often concentrated, as the work of
Lenz and Haake Lenz and Haake (1991), on the nearest neighbor level spacing
distribution. The distribution of normalized spacings $s_{i}=e_{i+1}-e_{i}$
where $e_{i}$ are the Hamiltonian matrix eigenvalues renormalized (in the
process called the level unfolding Haake (2010); Stöckmann (1999)) in such a
way that their mean density is unity. Some proposed an ad hoc expressions as a
popular Brody distribution Brody (1973) which surprisingly well fitted low
resolution experimental data – see e.g. Stöckmann (1999). In contrast, Berry
and Robnik Berry and Robnik (1984) proposed a distribution that was based on a
clear physical assumption of the separation between “chaotic” wavefunctions
faithful to GOE and those localized on the regular parts of the phase space.
This simple foundation resulted in the distribution which was shown to work
well for mixed phase space situations in the deep semiclassical limit Prosen
(1998). Importantly, long range correlations of energy levels can also be
addressed within this model Prosen and Robnik (1999).

Over the years, the most common approach to construct an ensemble
interpolating between GOE and uncorrelated Poisson level statistics has been
to modify the GOE distribution resigning from its beautiful property of being
invariant under orthogonal transformations. The original Rozenzweig-Porter
proposition Rosenzweig and Porter (1960) is just one of them. In the original
version a single parameter $\lambda=\sigma^{2}/N$ defined the ensemble of
matrices of rank $N$ where $\sigma^{2}$ was the variance of off-diagonal
matrix elements of the Hamiltonian matrix. For $\lambda\rightarrow 0$ one
recovers the Poisson ensemble while for $\lambda=1$ the GOE ensemble is
reproduced provided that the variance of the diagonal elements is equal to 2.
Recently, modifications of the Rosenzweig-Porter ensemble allowing for a
broad, log-normal distribution of the off-diagonal matrix elements, have been
shown to be valid for various types of many-body and Anderson localization
transitions Khaymovich  _et al._ (2020); Kravtsov  _et al._ (2020); Khaymovich
and Kravtsov (2021).

A second proposition goes back to Seligman and coworkers Seligman  _et al._
(1984) who postulated that the variance of off-diagonal elements $a_{ij}$
should scale as $\exp\left[-(i-j)^{2}/\sigma^{2}\right]$. For
$\sigma\rightarrow 0$ one recovers the Poisson case while for
$\sigma\rightarrow\infty$ the matrix belongs to GOE.

Another well-known approach is that of Guhr Guhr (1996) who used
supersymmetric techniques to express the two-level correlation function in the
Poisson-GOE ensemble in terms of a double integral. In the meantime Casati and
coworkers Casati  _et al._ (1990, 1991) as well as Fyodorov and Mirlin
Fyodorov and Mirlin (1991) considered banded Gaussian random matrices as a
useful tool in describing the transition, the corresponding parameter was
$y=b^{2}/N$ with $b$ being the matrix bandwidth and $N$ its rank.

For a long time, the nearest neighbor spacings were treated as essential
entities to capture the GOE-to-Poisson transition in the context of quantum
chaos. An another relevant parameter is the number-variance, defined as
$\Sigma^{2}(L)=\langle n(L)^{2}\rangle-\langle n(L)\rangle^{2}$, where $n(L)$
is the number of (unfolded) eigenvalues in an interval of length $L$. It was
demonstrated that for a large system size $L$, the number variance scales as
$\chi L$, where $\chi$ varies from $0$ to $1$ during the Possion-to-GOE
transitions Bogomolny, E.  _et al._ (2001); Bogomolny and Giraud (2011);
Bertrand and García-García (2016). However, the unfolding process, implemented
to compute $e_{i}$ and the number variance, is based on a separation of the
density of states $\rho(E)$ into smooth and fluctuating parts. The distinction
between these two parts of $\rho(E)$ is to some extent arbitrary and may lead
to dubious outcomes see, for instance, the differences between local unfolding
and Gaussian broadening Gómez  _et al._ (2002). Over the past years, there is
a constant effort to tackle such problems from the mathematical and
computational viewpoints Morales  _et al._ (2011); Torres-Vargas  _et al._
(2017).

To simplify the problem, Huse and Oganesyan Oganesyan and Huse (2007)
introduced an important measure known as gap ratio, defined as
$r_{n}=min[\delta_{n},\delta_{n-1}]/max[\delta_{n},\delta_{n-1}]$, where
$\delta_{n}=E_{n+1}-E_{n}$ is the energy gap between the consecutive energy
levels. The dimensionless gap-ratio is independent of the local density-of-
states, $\rho(E)$, and the problem of unfolding is avoided. Employing an exact
diagonalization, they also established that a spin system with time reversal
invariance undergoes a phase transition from ergodic (with GOE level
statistics) to the many-body localised (MBL) phase (signified by the
uncorrelated Possonian statistics of energy levels) as the magnitude of
disorder increases (see also Santos  _et al._ (2004)). In a further endeavor,
Atas et.al. Atas  _et al._ (2013) proposed a Wigner-like surmise for the gap
ratio statistics. While the level statistics across the ergodic-to-MBL
transition was still studied with standard level spacing distributions Santos
and Rigol (2010), the gap ratio became gradually a more popular tool Pal and
Huse (2010); Mondaini and Rigol (2015); Luitz  _et al._ (2015).

The quest for a proper level statistics model to describe ergodic-to-MBL
transition speeded up following the work of Serbyn and Moore Serbyn and Moore
(2016) who argued that the flow of level statistics between the GOE and
Poisson limits occurs in two stages: (1) A power-law type interaction between
the eigenvalues Kravtsov and Lerner (1995) of range that decreases with
increase of disorder strength and (2) an appearance of semi-poissonian level
statistics between the energy levels as an offshoot of local interactions
between the energy levels (see also Bogomolny  _et al._ (1999); Bogomolny, E.
_et al._ (2001)).

Among the latest developments a so called $\beta$-Gaussian($\beta-G$) model
was introduced Buijsman  _et al._ (2019) as a model to describe level
statistics across ergodic-to-MBL transition. The $\beta-G$ model is dependent
on a single parameter $\beta$ that characterizes the pairwise interaction of
the energy levels taking all the energy pairs in consideration. A real value
of $\beta\in[0,1]$ allows one to describe the whole crossover between GOE and
Poisson level statistics. For a comparison of the performance of different
models see Sierant and Zakrzewski (2019). In our recent work Sierant and
Zakrzewski (2020) , we have proposed a modified 2-parameter model, namely, the
$\beta-h$ model, where the interaction between the energy levels is limited to
a $h$ neighboring energy levels. This $2$-parameter model was shown to be
superior to its $\beta-G$ counterpart since it reproduced not only the nearest
neighbor spacing ratio but also the higher order spacing ratios defined by
Chavda  _et al._ (2014); Tekur  _et al._ (2018):

$r_{n}^{i}=min\bigg{[}\frac{E_{i+2n}-E_{i+n}}{E_{i+n}-E_{i}},\frac{E_{i+n}-E_{i}}{E_{i+2n}-E_{i+n}}\bigg{]},$
(1)

where $E_{i}$ stands for the $i$-th energy level. In the present work we
compare predictions of these recent models with the so called Pechukas-Yukawa
($P-Y$) level-dynamics approach and verify its underlying efficiency to
reproduce the level statistics in terms of the higher order spacing ratios
(1). In the next section we remind the derivation of the joint-probability-
distribution of the $P-Y$ model (which may be found, e.g. in Stöckmann (1999))
to make the paper self contained while later we compare the statistics
generated with different models in the transition between delocalized and MBL
cases.

##  III Statistical approach to level dynamics

Consider a Hamiltonian dependent on parameter $\lambda$:
$H(\lambda)=H_{0}+\lambda V$ for arbitrary $H_{0}$ and $V$ and the fate of its
eigenvalues as $\lambda$ is varied. Differentiating the eigenvalue equation

$\displaystyle H(\lambda)|a(\lambda)\rangle=x(\lambda)_{a}|a(\lambda)\rangle,$
(2)

with respect to $\lambda$ (above $x(\lambda)_{a}$ is the eigenvalue
corresponding to eigenvector $|a(\lambda)\rangle$) and taking appropriate
scalar products one immediately gets

$\displaystyle\frac{d}{d\lambda}x_{a}\equiv\dot{x}_{a}=\langle
a|V|a\rangle\equiv V_{aa}.$ (3)

Defining $p_{a}\equiv\dot{x}_{a}$ and looking at its $\lambda$ derivative
yields

$\displaystyle\dot{p}_{a}=2\sum_{b\neq
a}\frac{V_{ab}V_{ba}}{x_{a}-x_{b}}=2\sum_{b\neq
a}\frac{|f_{ab}|^{2}}{(x_{a}-x_{b})^{3}}$ (4)

where the second equality involves an additional definition
$f_{ab}=V_{ab}(x_{a}-x_{b})$. Most interestingly, differentiating $f_{ab}$
over $\lambda$, the set of differential equations closes:

$\dot{f}_{ab}=\sum_{r\neq
a,b}f_{ar}f_{rb}\bigg{[}\frac{1}{(x_{a}-x_{r})^{2}}-\frac{1}{(x_{b}-x_{r})^{2}}\bigg{]}.$
(5)

The set of equations (3)-(5) is known as the Pechukas-Yukawa equations.
Identifying $\lambda$ as a fictitious time, $x_{a}$’s and $p_{a}$’s can be
interpreted as positions and momenta of fictitious particles moving under the
force decaying as $1/x^{3}$ with the distance (assuming $f_{ab}$’s are
constant). This highly nonlinear set of of equations is neverthless integrable
(for a discussion see Haake (2010); Stöckmann (1999)).

Clearly, for $\lambda$ large the eigenvalues of $H=H_{0}+\lambda V$ will be
$V$ dominated and the variation of eigenvalues (and eigenvectors) trivializes.
For that reason Haake Haake (2010) introduces a slightly different $\lambda$
dependence (equivalent for small $\lambda$):
$H(\lambda)=\sqrt{f}(H_{0}+\lambda V)$ with $f=(1+\lambda^{2})^{-1}$ while we
shall use the form adopted from Zakrzewski and Delande (1993); Zakrzewski  _et
al._ (1993); Stöckmann (1999) $H=H_{0}\cos(\lambda)+V\sin(\lambda)$.

Then the equations above become slightly modified with

$\displaystyle\dot{x}_{a}=\langle a|\dot{H}|a\rangle=p_{a}$ (6)

and

$\displaystyle\dot{p}_{a}=-x_{a}+2\sum_{b\neq
a}\frac{|f_{ab}|^{2}}{(x_{a}-x_{b})^{3}}$ (7)

with $f_{ab}=\langle a|\dot{H}|b\rangle(x_{a}-x_{b})$. Thus the additional
harmonic force binds the fictitious particles preventing their escape to
infinity.

Pechukas-Yukawa equations corresponding to integrable system call for an
appropriate statistical mechanics description which should be realized within
the generalized Giggs ensemble Vidmar and Rigol (2016). This requires an
identification of the complete set of nontrivial constants of motion being in
convolution. Such a possible approach is a song of the future. Here we rather
follow Yukawa approach and consider the simplest integrals of motion. One of
them will be the total energy of the system of interacting classical particles

$E=\frac{1}{2}\sum_{n}(p_{n}^{2}+x_{n}^{2})+\frac{1}{2}\sum_{n,m}\frac{|f_{nm}|^{2}}{(x_{n}-x_{m})^{2}}.$
(8)

and the other $Q=\frac{1}{2}Tr(F^{2})$ called a total angular momentum, for a
justification see Stöckmann (1999). Let us note that those are the only two
second order constants of the motion.

With this two constants the phase-space density is given by the Gibbs ensemble
as

$\rho=\frac{1}{Z}\exp(-\beta E-\gamma Q),$ (9)

which in the presence of the harmonic binding potential reads explicitly

$\begin{split}\rho=\frac{1}{Z}\exp(-\beta\bigg{[}\frac{1}{2}\sum_{n}(p_{n}^{2}+x_{n}^{2})+\frac{1}{2}\sum_{n,m}\frac{|f_{nm}|^{2}}{(x_{n}-x_{m})^{2}}\bigg{]}\\\
-\gamma\frac{1}{2}\sum_{n,m}\bigg{|}f_{nm}\bigg{|}^{2})\end{split}$ (10)

Next, we can compute the joint probability distribution (JPD) of eigenvalues
Stöckmann (1999) by integrating out the variables $p_{n}$ and $f_{nm}$ from
the phase space density obtaining

$\begin{split}P(x_{1},x_{2},.....,x_{n})\sim\prod_{n<m}\bigg{|}\frac{(x_{n}-x_{m})^{2}}{1+\frac{\gamma}{\beta}(x_{m}-x_{n})^{2}}\bigg{|}^{\nu/2}\\\
exp\bigg{(}-\frac{\beta}{2}\sum_{n}x_{n}^{2}\bigg{)},\end{split}$ (11)

where $\nu=1,2,4$ corresponds to three possible ensembles of Dyson. All three
cases appear due to different properties of (integrated over) $f_{nm}$
variables. We shall consider the generalized time reversal invariant case of
$\nu=1$ only. From the expression above, one can easily reach the Possonian
distribution by putting $\gamma/\beta\gg 1$. The distribution becomes
proportional to

$P(x_{1},x_{2},.....,x_{n})\sim\exp\big{(}-\frac{\beta}{2}\sum_{n}x_{n}^{2}\big{)},$
(12)

while in the opposite limit, $0<\gamma/\beta<<1$, the distribution mimics the
Gaussian ensemble.

$P(x_{1},x_{2},.....,x_{n})\sim\prod_{n>m}\bigg{|}x_{n}-x_{m}\bigg{|}^{\nu}\exp\big{(}-\frac{\beta}{2}\sum_{n}x_{n}^{2}\big{)}.$
(13)

To obtain an ensemble that interpolates between GOE and Poisson level
statistics we fix $\beta=\nu=1$ and denote $\gamma/\beta=10^{p}$. The
distribution (11) becomes

$\begin{split}P(x_{1},x_{2},.....,x_{n})\sim\prod_{n<m}\bigg{|}\frac{(x_{n}-x_{m})^{2}}{1+10^{p}(x_{m}-x_{n})^{2}}\bigg{|}^{1/2}\\\
exp\bigg{(}-\frac{1}{2}\sum_{n}x_{n}^{2}\bigg{)}\end{split},$ (14)

where $p=\log_{10}\frac{\gamma}{\beta}$ is the single free parameter of the
distribution. The first term in the (14) signifies the pairwise interaction
between the particles and the exponential term provides the harmonic binding
of the eigenvalues. This single parameter distribution, resulting from
Pechukas-Yukawa statistical mechanics for time reversal invariant case, will
be tested against other existing propositions as well as numerical data from
interacting many-body models.

##  IV Comparison with numerical data

For the simulation purpose, the eigenvalues are drawn by sampling (14) with
Metropolis-Hastings algorithm Metropolis  _et al._ (1953); Hastings (1970) for
$p\in[-10,+15]$ with an initial sample size $N=500$ distributed over a linear
chain. This moderate system size is chosen for convenience and speed. We
consider about 100 eigenvalues in the middle of the spectrum. Since we
consider high order gap ratios too, this sample size may be insufficient for
$r_{n}$ with $n$ large. We consider the issue of the system size in the next
Section. The obtained eigenvalue distributions interpolate between a faithful
semi-circle histogram of the density-of-states in the GOE limit (when
$p\rightarrow-\infty$) and a Gaussian distribution for $p$ large and positive.
In the following we extract the higher order spacing ratios from (14) and fit
different physical models results with that distribution. For comparison we
use the single parameter $\beta$-Gaussian model Buijsman  _et al._ (2019) as
well as the two-parameter $\beta-h$ model shown to yield accurate
representation of many body data Sierant and Zakrzewski (2020).

###  IV.1 Random Heisenberg chain

Figure 1: Probability distributions of the higher order spacing ratios,
$r_{n}$, for $n=1,2,3,4,5,8$ for a disordered Heisenberg spin chain with
$L=18$ and the closest distributions given by $P-Y$ model. To facilitate the
comparison the numerical data for the Heisenberg chain are shown by bullets
and the corresponding fits with the $P-Y$ model by lines of the same color.
$W=1.8$ data (black bullets) are fitted with P-Y model with $p=2.3$; blue set
is for $W=2.0$ ($p=3.1$); orange set is for $W=2.4$ ($p=4.4$); light green set
for $W=2.8$ ($p=5.4$); dark green set for $W=3.0$ ($p=5.8$); and the red set
for $W=4.0$ ($p=7.4$).

As a first example we take a disordered Heisenberg spin-$\frac{1}{2}$ chain, a
paradigmatic model for MBL studies Luitz  _et al._ (2015, 2016); Serbyn and
Moore (2016); Alet and Laflorencie (2018); Macé  _et al._ (2019); Sierant and
Zakrzewski (2019). The Hamiltonian of the system reads

$\hat{H}_{H}=J\sum_{i=1}^{L}\vec{S}_{i}{\cdot}\vec{S}_{i+1}+\sum_{i=1}^{L}h_{i}{S}_{i}^{z},$
(15)

where the first term represents the exchange interaction between the
neighboring 1/2-spins with exchange coupling $J$ (normalized to unity). The
last term constitutes the on-site disorder potential at site $i$ with disorder
$h_{i}$ drawn from uniform random distribution, $h_{i}\in[-W,W]$. It is well
established that this system undergoes a transition from ergodic to MBL phase
Luitz  _et al._ (2015, 2016); Serbyn and Moore (2016); Alet and Laflorencie
(2018); Sierant and Zakrzewski (2019) for sufficiently large $W$.

The eigenvalues of the Heisenberg chain are obtained for the length $L=18$
with periodic boundary conditions (PBC) via the shift-and-invert method. In
the due course, we use 2000 disorder realizations for each value of $W$. The
gap ratios are calculated from approximately $500$ eigenvalues from the middle
of the spectrum. In Fig. 1, we present a comparison between the obtained gap-
ratio distributions $P(r_{n})$ and those resulting from the closest
predictions coming from the $P-Y$ model (14). Different $n$ values allow us to
explore correlations on different distances between eigenvalues. Let us stress
that we use a single $p$ value and the resulting $P-Y$ distribution to fit all
$P(r_{n})$ for $n\in[1,10]$ by minimizing the cumulative error between the
gap-ratios of Heisenberg chain and the gap-ratios predicted by $P-Y$ model. As
may be appreciated in Fig. 1, the $P-Y$ distribution fits the data of the spin
model remarkably well, in particular on the delocalized side of the transition
(smaller $W$ values). Only close to the fully localized case some small
deviations are seen.

Figure 2: Same as Fig. 1 but for the $\beta-G$ ensemble. The value of
${}^{\prime}beta$ is fitted for $P(r_{1})$ distribution as in Buijsman  _et
al._ (2019). Higher order gap ratios are clearly not reproduced
satisfactorily. The numerical data for the Heisenberg chain are shown by
bullets and the fits by $\beta-G$ model by lines of the same color. $W=1.8$
data (black bullets) are fitted with $\beta-G$ model with $\beta=1.0$; blue
set is for $W=2.0$ ($\beta=0.94$); orange set is for $W=2.4$ ($\beta=0.68$);
light green set for $W=2.8$ ($\beta=0.38$); dark green set for $W=3.0$
($\beta=0.26$); and the red set for $W=4.0$ ($\beta=0.06$).

To see how competitive $P-Y$ model is we compare its predictions with those
given by a single parameter $\beta-G$ approach, in which sampling of the
distribution is much easier Buijsman  _et al._ (2019) as well as a two
parameter family of $\beta-h$ distributions Sierant and Zakrzewski (2020). For
a most direct comparison with data in Fig. 1 we present similar figures for
$\beta-G$ in Fig. 2 and for $\beta-h$ ensemble in Fig. 3. Even a casual
inspection of Fig. 2 reveals that $\beta-G$ ensemble fitted for $r_{1}$ (where
it performs remarkably well) does not represent properly higher gap ratios.
Clearly this ensemble does not reproducer short-range spacings and long-range
correlations equally well. On the other hand $\beta-h$ ensemble performs at
least as well as the P-Y distribution considered.

Figure 3: Same as Fig. 1 and Fig. 2 but for $\beta-h$ distridution which fits
Heisenberg spin data remarkably well. $W=1.8$ data (black bullets) are fitted
with $\beta-h$ model with $\beta=1$, $h=10$; blue set is for $W=2.0$
($\beta=0.96$, $h=6$); orange set is for $W=2.4$ ($\beta=0.74$, $h=2.7$);
light green set for $W=2.8$ ($\beta=0.46$, $h=1.8$); dark green set for
$W=3.0$ ($\beta=0.3$, $h=1.4$); and the red set for $W=4.0$ ($\beta=0.08$,
$h=1.05$). Some of the data plotted were presented in Fig. 3 of Sierant and
Zakrzewski (2020) where $\beta-h$ model was introduced.  Figure 4: (a) The
comparison of the cumulative error of the fit $\epsilon$, (16) for all three
studied distributions. Red dashed line corresponds to $\beta-G$ model, blue
solid line for $\beta-h$ while green dashed-dot line for $P-Y$ model. (b) The
performance of different models compared for different mean gap ratios
$\langle r_{n}\rangle$ as compared with the data for XXZ spin chain for
$W=2.0$ (squares), $W=2.8$ (circles) and $W=4.0$ (diamonds) We plot the
difference between the average $\langle r_{n}\rangle$ and the corresponding
value for the Poisson distribution as it better represents the quality of
fits. The lines correspond to different models consistently with panel (a).

For a global, quantitative comparison, we construct a measure of the fit
quality, i.e., the average deviation

$\epsilon\equiv\frac{1}{10}\sum_{n=1}^{10}\bigg{[}\langle
r_{n}\rangle_{M}-\langle r_{n}\rangle_{H}\bigg{]},$ (16)

where $\langle r_{n}\rangle_{H}$ is the numerical average for the Heisenberg
chain while $\langle r_{n}\rangle_{M}$ is the corresponding value coming from
a fit to a given model, be it $P-Y$, $\beta-h$ or $\beta-G$ version. The
obtained values of $\epsilon$ are plotted in Fig. 4(a), for $\beta-G$ model
$\beta$ is fitted to reproduce as closely as possible $\langle
r_{1}\rangle_{H}$. Close to the delocalized GOE-like regime the resulting
$\epsilon$ for that ensemble is comparable to both $P-Y$ and $\beta-h$ values.
Closer to the transition and on the localized side one may clearly observe
that $\beta-G$ model performs poorly when compared with $P-Y$ and $\beta-h$
models. The latter is only slightly superior but it involves fitting of two
parameters instead of one.

Fig. 4(b) presents the average gap ratios $\langle r_{n}\rangle$ for
$n\in[1,40]$ or rather, for a better visualization, a difference
$\Delta\langle r_{n}\rangle=\langle r_{n}\rangle-\langle r^{\rm P}_{n}\rangle$
where the latter corresponds to an analytic prediction for the mean gap ratio
for a Poissonian ensemble Sierant and Zakrzewski (2020) (note that
normalization factor is missing there). On a first glance it is clear, that
the $\beta-G$ model is inaccurate beyond $r_{1}$ as it severely overestimates
$\Delta\langle r_{n}\rangle$, which means that the long-range correlations
between eigenvalues of $\beta-G$ model are much stronger than the correlations
of eigenvalues of the disordered Heisenberg spin chain . Both $\beta-h$ and
$P-Y$ reproduce higher order spacing ratios much more accurately. The single
parameter $P-Y$ model is only slightly inferior to the $\beta-h$ model.

Figure 5: The distributions of the higher order gap-ratios of disordered
quasi-periodic spin chain with $L=16$ and the best fit obtained in the $P-Y$
model. The numerical data for the quasi-periodic chain for a given disorder
$W$ are presented by bullets and the resulting fits by $P-Y$ model are shown
by lines. $W=2.0$ data (black bullets) are fitted with $P-Y$ model with
$p=3.4$; blue set is for $W=2.5$ ($p=5.0$); orange set is for $W=3.0$
($p=6.2$); green set for $W=3.5$ ($p=7.2$); and the red set for $W=4.0$
($p=7.8$).

###  IV.2 Quasi-periodic Heisenberg spin chain

Consider the same Heisenberg chain with the Hamiltonian (15) but now $h_{i}$
are not random but taken as $h_{i}=(W/2)\cos(2\pi\chi i+\alpha)$, where $i$ is
the site index. Different realizations of the disorder correspond to different
choices of $\alpha$ drawn from random uniform distribution on $[0,2\pi)$
interval. Contrary to the previous case, this quasi-periodic disorder is fully
correlated. Still, for sufficiently large $W$ (dependent on the value of
$\chi$ Guarrera  _et al._ (2007); Doggen and Mirlin (2019)) the system
undergoes a transition to MBL phase. Such a quasiperiodic disorder was
implemented in experiments of the Munich group Schreiber  _et al._ (2015);
Lüschen  _et al._ (2017) by placing an additional weak standing wave on top of
the primary one forming the optical lattice. Then $\chi$ is the ratio of
wavevectors of both laser wavelengths.

In our study we take $\chi=(\sqrt{5}-1)/2$ and diagonalize the Hamiltonian for
$L=16$ and open boundary conditions. Like in the preceding case, we implement
the shift-invert method to collect eigenvalues from $2000$ disorder
realizations and the higher order spacing ratios are computed form $500$
eigenvalues taken from the middle of the spectrum.

Figure 6: Same as Fig. 5 but a comparison of the quasi-periodic chain data
with their fits for $\beta-G$ ensemble is shown. As previously color
differentiates between different disorder amplitudes $W$ (bullets) and
different $beta-G$ distributions (lines) obtained by fitting $P(r_{1})$.
$W=2.0$ data (black bullets) are fitted with $\beta-G$ model with
$\beta=0.92$; blue set is for $W=2.5$ ($\beta=0.58$); orange set is for
$W=3.0$ ($\beta=0.18$); green set for $W=3.5$ ($\beta=0.06$); and the red set
for $W=4.0$ ($\beta=0.01$).

The whole analysis is performed similarly to the random disorder case. We
concentrate on the localized side of the crossover (for the system size
studied) where the differences between the considered distributions show more
clearly. Fig. 5 shows the resulting fits for $P-Y$ ensemble showing that
indeed it works very well also for this case. For comparison results obtained
for $\beta-G$ and $\beta-h$ ensembles are shown in Fig. 6 and Fig. 7,
respectively. The conclusions are very similar to the previous case (although
limited to smaller interval of $W$ taken). Firstly, the single parameter $P-Y$
model outperforms $\beta$-Gaussian proposition considerably, especially for
higher order gap ratios. Similarly to the uniform random disorder case,
$\beta$-Gaussian model works reasonably well only very close to the GOE limit
(data corresponding to $W=1.5$ in Fig. 5. Secondly, $P-Y$ distribution is
clearly less effective than the phenomenological, two parameter $\beta-h$
approach in the whole interval of $W$ values. But, let us stress, this
difference is quite small.

Figure 7: Same as Fig. 5 but a comparison of the quasi-periodic chain data
with their fits for $\beta-h$ model is presented. Different color corresponds
to disorder amplitudes $W$ (bullets) and different fits (lines) for the
$\beta-h$ model. $W=2.0$ data (black bullets) are fitted with $\beta-h$ model
with $\beta=0.92$; blue set is for $W=2.5$ ($\beta=0.58$); orange set is for
$W=3.0$ ($\beta=0.18$); green set for $W=3.5$ ($\beta=0.06$); and the red set
for $W=4.0$ ($\beta=0.01$).

The global error comparison presented in Fig. 8 additionally confirms the
above conclusions.

It is worthwhile to comment more on the comparison of random and quasiperiodic
cases. The crossover to MBL present in both models is quite well reproduced by
$P-Y$ or $\beta-h$ statistical models despite the fact that the systems with
random and quasiperiodic disorder behave quite differently on a microscopic
level. This was observed in [64] inspecting the entanglement entropy
statistics and further analysed, on the level of level statistics, in our
earlier work [49]. We have shown there that one may distinguish the MBL
transition in random and quasiperiodic disorder by examining the so called
intersample variances which are the sample-to-sample variances of gap ratios.
Those show a peak in the crossover region for a random disorder case. This
fact was attributed to the presence of rare Griffiths-type regions in systems
with random disorder. The corresponding peak was absent in the quasiperiodic
data supporting rare regions interpretation. One may pose a question: how
inter sample variances of gap ratio behave in the statistical models such as
$P-Y$ or $\beta-h$ models? We have verified that inter-sample variations of
gap ratio practically do not change in the GOE-Poisson transition for the
$P-Y$ model by a direct evaluation. Our explanation of this discrepancy
between behavior of a system with random disorder at the MBL transition and
$P-Y$ model is the following. The rare regions affect variances of the former,
for which the number of independent random variables determining the
Hamiltonian is the same as the system size, $L$, very small compared to the
Hilbert space dimension. For the latter, the number of random variables in
simulated matrices is much larger and scales as a square of the matrix size
making the sample-to-sample variations very small across the whole GOE-Poisson
crossover.

Figure 8: Same as Fig. 4 but for quasiperiodic disorder. The cumulative error
of $\beta-G$ distribution is much larger in the transition regime than for the
other two distributions. Note that $P-Y$ again performs comparably to
$\beta-h$ model.

###  IV.3 Bose-Hubbard model

Figure 9: $P(r_{n})$ for $n=1,2,3,4,5,8$ for a Bose-Hubbard model with $L=8$
as fitted by the $P-Y$ model. The numerical data for BH model are shown by
dots and the fits of $P-Y$ model by lines. $W=7$ data (black bullets) are
fitted with $P-Y$ model with $p=3.0$; blue set is for $W=10$ ($p=4.4$); orange
set is for $W=12$ ($p=5.1$); green set for $W=15$ ($p=5.8$); and the red set
for $W=25$ ($p=7.3$).

As the last example we test $P-Y$ model on a different disordered system
namely the Bose-Hubbard model with the Hamiltonian Sierant  _et al._ (2017);
Sierant and Zakrzewski (2018); Orell  _et al._ (2019); Hopjan and Heidrich-
Meisner (2020)

$\hat{H}_{BH}=-J\sum_{<i,j>}\hat{a}_{i}^{\dagger}\hat{a}_{j}+\frac{U}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i}-1)+\sum_{i}\mu_{i}\hat{n}_{i},$
(17)

where $a_{k}^{\dagger}(a_{k})$ are bosonic creation and annhilation operators
at site $k$, $J$ is the tunneling and $U$ is the onsite interaction strength.
We assume further on $J=U=1$ while the chemical potential $\mu_{i}$ is assumed
to be uniformly randomly distributed within an interval $[-W,W].$ As the
considered previously spin chain a transition from delocalized phase to MBL is
observed for such a disordered Bose-Hubbard model with increasing $W$ Sierant
and Zakrzewski (2018).

Figure 10: Same as Fig. 9 but a comparison of the Bose-Hubbard model data with
their fits for $\beta-G$ ensemble is presented. Colors differentiate between
different disorder amplitudes $W$ (bullets) and different $beta-G$
distributions (lines) obtained by fitting $P(r_{1})$. $W=7$ data (black
bullets) are fitted with $\beta-G$ model with $\beta=0.94$; blue set is for
$W=10$ ($\beta=0.62$); orange set is for $W=12$ ($\beta=0.42$); green set for
$W=15$ ($\beta=0.22$); and the red set for $W=25$ ($\beta=0.04$).  Figure 11:
Similar comparison to Fig. 10 but fitting Bose-Hubbard data with the $\beta-h$
ensemble. $W=7$ data (black bullets) are fitted with $\beta=0.98$ and $h=6.5$;
blue set is for $W=10$ ($\beta=0.70$, $h=2.6$); orange set is for $W=12$
($\beta=0.46$, $h=2.4$); green set for $W=15$ ($\beta=0.24$, $h=2.4$); and the
red set for $W=25$ ($\beta=0.06$, $h=2$). Some of the data plotted were
presented in Fig. 8 of Sierant and Zakrzewski (2020) where $\beta-h$ model was
introduced.

We consider a chain of $L=8$ sites at unit filling and eigenvalues are
computed via exact-diagonalization method Sierant and Zakrzewski (2019). As
previously, the gap ratios are determined from the eigenvalues collected from
the middle of the spectrum for no less than $500$ disorder realizations. The
distribution of gap ratio $P(r_{n})$ obtained are fitted with the $P-Y$ model
as shown in Fig. 9 for $n=1,2,3,4,5,8$. One may observe a rather outstanding
agreement between the data and the statistical $P-Y$ model except for slight
deviations close to full MBL (Poisson) cases.

Figure 12: Comparison of the error of all three distributions for Bose-Hubbard
model. For that system $\beta-h$ model clearly outperforms $P-Y$ result, the
latter is still acceptable as may be seen in the previous figure. $\beta-G$
model is performing much worse

The corresponding comparison of numerical data for the disordered Bose-Hubbard
model with $\beta-G$ and $\beta-h$ models are shown in Fig. 10 and Fig. 11,
respectively. Similarly to previously considered spin chains , $P-Y$ model is
clearly superior over the $\beta-G$ model in reproducing the data while being
slightly outperformed by the two parameters $\beta-h$ model, as further
summarised in Fig. 12.

##  V Size considerations and the universal distribution

The results presented for all three cases discussed show convincingly that the
single parameter Pechukas-Yukawa distribution (14) reproduces remarkably well
distributions of higher order spacing ratios in the transition between GOE and
Poisson regime. Let us recall the fact that we have consistently used data for
$N=500$ ensemble. It is clear, from the form of the distribution (14) that in
the Poisson limit the density of states becomes a Gaussian with a unit
variance independent of the system size. In the other, GOE limit the density
of states follows a semicircle law in $(-\sqrt{2N},\sqrt{2N})$ interval being
strongly $N$ dependent. In effect for any given set of data, the value of the
fitted parameter $p$ in (14) depends on $N$ chosen.

This estetic drawback may be cured by rewriting (14) as

$\begin{split}P(x_{1},x_{2},.....,x_{n})\sim\prod_{n<m}\bigg{|}\frac{(x_{n}-x_{m})^{2}}{1+10^{Y}(x_{m}-x_{n})^{2}/\Delta^{2}}\bigg{|}^{1/2}\\\
exp\bigg{(}-\frac{1}{2}\sum_{n}x_{n}^{2}\bigg{)},\end{split}$ (18)

where in the denominator the factor $10^{p}$ is replaced by
$10^{Y}/\Delta^{2}$, where $\Delta$ is the mean spacing of $\\{x_{n}\\}$. This
defines a new single parameter $Y$ characterizing the distribution, while
$\Delta=\Delta(N,Y)$ can be easily obtained from generated sequence of
$x_{n}$. We have verified that such an approach makes the obtained
distributions $N$ independent as visualised in Fig. 13 which presents the
average gap ratios $\langle r_{i}\rangle$ for $i=1,3$ for $N=500$ and
$N=2000$.

Figure 13: The mean gap ratios $\left\langle r_{1}\right\rangle$ and
$\left\langle r_{3}\right\rangle$ offset by their Poisson value
$r_{1}^{P}\approx 0.3862$ and $r_{3}^{P}\approx 0.5625$ for $N=500$ and
$N=2000$ obtained across the MBL transition using the universal distribution
(18).

##  VI Conclusions

The aim of this work was to investigate how the distributions obtained from
JPD of eigenvalues resulting from standard Pechukas-Yukawa statistical
approach compare with currently discussed models of level statistics
interpolating between GOE and integrable cases. The nice feature of $P-Y$
model, as given by (14), is that a single parameter, denoted by $p$ allows for
an interpolation in the whole interval between GOE and Poisson limits. The
distribution studied results in statistics closely following the numerical
data and in this respect is clearly superior to other single parameter model
such as $\beta$-Gaussian proposition Buijsman  _et al._ (2019). The latter has
a clear advantage that probing it does not require the Metropolis algorithm,
so probing it is relative simple. We also made a comparison with the so called
$\beta-h$ model Sierant and Zakrzewski (2020) which quite accurately
represents numerical data. The difficulty of applying the latter are
comparable. It is quite rewarding to see that a single parameter $P-Y$ model
yields results only slightly inferior to two-parameters $\beta-h$ approach.

It is worth stressing that while we concentrated on eigenvalue statistics, the
Pechukas-Yukawa JPD contains also information on distribution of matrix
elements $f_{nm}$ of a specific operator $F=[H,\dot{H}]$. The latter is
manifestly time independent - thus matrix elements $f_{nm}$ carry information
on eigenstates dependence on the parameter. Indeed Nakamura and Lakshmanan
Nakamura and Lakshmanan (1986) rephrased the dynamics in terms of eigenvalues
and components of eigenvectors. Thus JPD integrated over eigenvalues may yield
prediction on eigenvector statistics in the transition regime (or
alternatively on matrix elements of a generic operator).

Last but not least let us mention that one could aim at more “proper” analysis
of integrable model given by (3)-(5) with the appropriate generalized Gibbs
ensemble Vidmar and Rigol (2016) provided a proper identification of all
independent constants of the motion is made. A more modest approach will just
include some of them going beyond the presented here quadratic terms approach.

The presented study may be generalized to cases with broken generalized time-
reversal invariance as well as to an analogous approach possible for time-
periodic Floquet systems.

###### Acknowledgements.

J.Z. thanks Yan Fyodorov for suggestions on the literature of the subject and
Mariusz Gajda and Krzysztof Sacha for discussions on Metropolis algorithm
implementation. The support of PL-Grid Infrastructure is acknowledged. This
research has been supported by National Science Centre (Poland) under project
2019/35/B/ST2/00034 (B.D., J.Z.). P.S. acknowledges the support of Foundation
for Polish Science (FNP) through scholarship START.

## References

  * Bonifacio  _et al._ (1971) R. Bonifacio, P. Schwendimann, and F. Haake, Phys. Rev. A 4, 302 (1971).
  * Haake  _et al._ (1979) F. Haake, H. King, G. Schröder, J. Haus, and R. Glauber, Phys. Rev. A 20, 2047 (1979).
  * Haake (2010) F. Haake,  _Quantum Signatures of Chaos_ (Springer, Berlin, 2010).
  * Haake  _et al._ (1987) F. Haake, M. Kuś, and R. Scharf, Zeitschrift für Physik B Condensed Matter 65, 381 (1987).
  * Gutzwiller (1971) M. C. Gutzwiller, Journal of Mathematical Physics 12, 343 (1971).
  * Müller  _et al._ (2004) S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Phys. Rev. Lett. 93, 014103 (2004).
  * Müller  _et al._ (2005) S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Phys. Rev. E 72, 046207 (2005).
  * Heusler  _et al._ (2007) S. Heusler, S. Müller, A. Altland, P. Braun, and F. Haake, Phys. Rev. Lett. 98, 044103 (2007).
  * Pechukas (1983) P. Pechukas, Phys. Rev. Lett. 51, 943 (1983).
  * Yukawa (1985) T. Yukawa, Phys. Rev. Lett. 54, 1883 (1985).
  * Haake and Kuś (1988) F. Haake and M. Kuś, Europhysics Letters (EPL) 6, 579 (1988).
  * Lenz and Haake (1991) G. Lenz and F. Haake, Phys. Rev. Lett. 67, 1 (1991).
  * Rosenzweig and Porter (1960) N. Rosenzweig and C. E. Porter, Phys. Rev. 120, 1698 (1960).
  * Mehta (1990) M. L. Mehta,  _Random Matrices_ (Elsevier, Amsterdam, 1990).
  * Stöckmann (1999) H. Stöckmann,  _Quantum Chaos: An Introduction_ (Cambridge University Press, Cambridge, 1999).
  * Bohigas  _et al._ (1984) O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52, 1 (1984).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
  * Alet and Laflorencie (2018) F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018), quantum simulation / Simulation quantique.
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Brody (1973) T. A. Brody, Lettere al Nuovo Cimento (1971-1985) 7, 482 (1973).
  * Berry and Robnik (1984) M. V. Berry and M. Robnik, Journal of Physics A: Mathematical and General 17, 2413 (1984).
  * Prosen (1998) T. Prosen, Journal of Physics A: Mathematical and General 31, 7023 (1998).
  * Prosen and Robnik (1999) T. Prosen and M. Robnik, J. Phys. A: Math. Gen. 32, 1863 (1999).
  * Khaymovich  _et al._ (2020) I. M. Khaymovich, V. E. Kravtsov, B. L. Altshuler, and L. B. Ioffe, Phys. Rev. Research 2, 043346 (2020).
  * Kravtsov  _et al._ (2020) V. Kravtsov, I. Khaymovich, B. Altshuler, and L. Ioffe, arXiv preprint arXiv:2002.02979 (2020).
  * Khaymovich and Kravtsov (2021) I. Khaymovich and V. Kravtsov, arXiv preprint arXiv:2106.01965 (2021).
  * Seligman  _et al._ (1984) T. H. Seligman, J. J. M. Verbaarschot, and M. R. Zirnbauer, Phys. Rev. Lett. 53, 215 (1984).
  * Guhr (1996) T. Guhr, Annals of Physics 250, 145 (1996).
  * Casati  _et al._ (1990) G. Casati, L. Molinari, and F. Izrailev, Phys. Rev. Lett. 64, 1851 (1990).
  * Casati  _et al._ (1991) G. Casati, F. Izrailev, and L. Molinari, Journal of Physics A: Mathematical and General 24, 4755 (1991).
  * Fyodorov and Mirlin (1991) Y. V. Fyodorov and A. D. Mirlin, Phys. Rev. Lett. 67, 2405 (1991).
  * Bogomolny, E.  _et al._ (2001) Bogomolny, E., Gerland, U., and Schmit, C., Eur. Phys. J. B 19, 121 (2001).
  * Bogomolny and Giraud (2011) E. Bogomolny and O. Giraud, Phys. Rev. Lett. 106, 044101 (2011).
  * Bertrand and García-García (2016) C. L. Bertrand and A. M. García-García, Phys. Rev. B 94, 144201 (2016).
  * Gómez  _et al._ (2002) J. M. G. Gómez, R. A. Molina, A. Relaño, and J. Retamosa, Phys. Rev. E 66, 036209 (2002).
  * Morales  _et al._ (2011) I. O. Morales, E. Landa, P. Stránský, and A. Frank, Phys. Rev. E 84, 016203 (2011).
  * Torres-Vargas  _et al._ (2017) G. Torres-Vargas, R. Fossion, C. Tapia-Ignacio, and J. C. López-Vieyra, Phys. Rev. E 96, 012110 (2017).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Santos  _et al._ (2004) L. F. Santos, G. Rigolin, and C. O. Escobar, Phys. Rev. A 69, 042304 (2004).
  * Atas  _et al._ (2013) Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
  * Santos and Rigol (2010) L. F. Santos and M. Rigol, Phys. Rev. E 81, 036206 (2010).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Mondaini and Rigol (2015) R. Mondaini and M. Rigol, Phys. Rev. A 92, 041601 (2015).
  * Luitz  _et al._ (2015) D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
  * Serbyn and Moore (2016) M. Serbyn and J. E. Moore, Phys. Rev. B 93, 041424 (2016).
  * Kravtsov and Lerner (1995) V. E. Kravtsov and I. V. Lerner, J. Phys. A: Mat. Gen. 28, 3623 (1995).
  * Bogomolny  _et al._ (1999) E. B. Bogomolny, U. Gerland, and C. Schmit, Phys. Rev. E 59, R1315 (1999).
  * Buijsman  _et al._ (2019) W. Buijsman, V. Cheianov, and V. Gritsev, Phys. Rev. Lett. 122, 180601 (2019).
  * Sierant and Zakrzewski (2019) P. Sierant and J. Zakrzewski, Phys. Rev. B 99, 104205 (2019).
  * Sierant and Zakrzewski (2020) P. Sierant and J. Zakrzewski, Phys. Rev. B 101, 104201 (2020).
  * Chavda  _et al._ (2014) N. Chavda, H. N. Deota, and V. K. B. Kota, Phys. Lett. A 378, 3012 (2014).
  * Tekur  _et al._ (2018) S. H. Tekur, U. T. Bhosale, and M. S. Santhanam, Phys. Rev. B 98, 104305 (2018).
  * Zakrzewski and Delande (1993) J. Zakrzewski and D. Delande, Phys. Rev. E 47, 1650 (1993).
  * Zakrzewski  _et al._ (1993) J. Zakrzewski, D. Delande, and M. Kuś, Phys. Rev. E 47, 1665 (1993).
  * Vidmar and Rigol (2016) L. Vidmar and M. Rigol, Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
  * Metropolis  _et al._ (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).
  * Hastings (1970) W. K. Hastings, Biometrika 57, 97 (1970).
  * Luitz  _et al._ (2016) D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 93, 060201 (2016).
  * Macé  _et al._ (2019) N. Macé, N. Laflorencie, and F. Alet, SciPost Phys. 6, 50 (2019).
  * Guarrera  _et al._ (2007) V. Guarrera, L. Fallani, J. E. Lye, C. Fort, and M. Inguscio, New J. Phys. 9, 107 (2007).
  * Doggen and Mirlin (2019) E. V. H. Doggen and A. D. Mirlin, Phys. Rev. B 100, 104203 (2019).
  * 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).
  * 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).
  * Sierant  _et al._ (2017) P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. A 95, 021601 (2017).
  * Sierant and Zakrzewski (2018) P. Sierant and J. Zakrzewski, New Journal of Physics 20, 043032 (2018).
  * Orell  _et al._ (2019) T. Orell, A. A. Michailidis, M. Serbyn, and M. Silveri, Phys. Rev. B 100, 134504 (2019).
  * Hopjan and Heidrich-Meisner (2020) M. Hopjan and F. Heidrich-Meisner, Phys. Rev. A 101, 063617 (2020).
  * Nakamura and Lakshmanan (1986) K. Nakamura and M. Lakshmanan, Phys. Rev. Lett. 57, 1661 (1986).