\subsection{Analytic derivation of level statistics from the Brownian motion}
\subsubsection{Wigner-Dyson statistics}
Let us begin with reproducing the WD statistics from the Brownian motion model. We assume the Brownian motion in the space of random matrices, 
\begin{equation} \label{Eq:dH}
 H(\tau) = H_0+\int_0^\tau d\tau\, V(\tau),
\end{equation}
where matrix $V(\tau)$ satisfies the following properties:
\begin{eqnarray} \label{Eq:Random-diag}
\langle V_{nm}(\tau) \rangle &=& -H_{nm}(\tau),\\
\langle V_{nm}(\tau) V_{mn}(\tau') \rangle &=&  \delta(\tau-\tau'), \\ \label{Eq:Random-offdiag}
\langle V_{nn}(\tau) V_{mm}(\tau') \rangle &=&  \frac{2}{\beta_{RM}} \delta(\tau-\tau')\delta_{nm}.
\end{eqnarray}
Parameter $\beta_{RM}$ specifies the symmetry class, $\beta_{RM}=1$ for GOE, and $\beta_{RM}=2$ for GUE. 

For convenience, we fix the basis to coincide with the (instantaneous) eigenbasis of $H(\tau)$. We  apply the perturbation theory to calculate correction to eigenvalues of $H(\tau)$, $\{s_n\}$ induced by the change in the matrix
\begin{equation} \label{Eq:dV}
\Delta V  = \int_\tau^{\tau+\Delta \tau} d\tau \, V(\tau) .
\end{equation}
Resulting correction to $s_n$ reads:
\begin{equation} \label{Eq:sn-2opt}
  \Delta s_n = \Delta V_{nn}+\sum_{m\neq n} \frac{|\Delta V_{nm}|^2}{s_n-s_m}.
\end{equation}
Averaging this equation over $V$ using the fact that $H_{nn} = s_n$ in the eigenbasis, we get:
\begin{eqnarray} \label{Eq:ds}
  \langle \Delta s_n \rangle
  &=&
  \mu_n(\{s\})\Delta \tau,
  \\\label{Eq:drift}
    \mu_n[\{s\}]
    &=&
    - s_n + \sum_{m\neq n} \frac{1}{s_n-s_m}.
\end{eqnarray} 
Using the second moment 
\begin{equation} \label{Eq:Dl2}
  \langle \Delta s_n \Delta s_m \rangle
  =
 \frac{2}{\beta_{RM}}\delta_{nm} \Delta \tau,
\end{equation}
we find that $s_n(\tau)$ obeys the following Langevin equation:
\begin{equation} \label{Eq:Langevin}
  \frac{d s_n}{d \tau}
  =
  \mu_n[\{s\}]
  + \xi_n (\tau),
\end{equation}
where drift term is given by Eq.~(\ref{Eq:drift}) and white noise is specified by
\begin{eqnarray} \label{Eq:noise}
  \langle \xi_{n}(\tau )\xi_m (\tau')\rangle
  &=&
 \frac{2}{\beta_{RM}} \delta_{nm}\delta(\tau-\tau').
\end{eqnarray}

Translating the Langevin equation into the Fokker-Planck  equation governing the evolution of the joint probability distribution function $P(\{s_n\},\tau)$, we get:
\begin{multline} \label{Eq:FP}
\frac{\partial P(\{s\},\tau)}{\partial \tau}
=
\sum_n \left\{
-\frac{\partial }{\partial s_n} [ \mu_n(\{s\}) P(\{s\},\tau)] \right.
\\ \left.
+ \frac{1}{\beta_{RM}} \frac{\partial^2  }{\partial s_n^2}  P(\{s\},\tau)]
\right\}.
\end{multline}
Using explicit expression for the drift, Eq.~(\ref{Eq:drift}) we see that the joint probability distribution
\begin{equation} \label{Eq:PWD}
P(\{s\})
=
C \prod_{m<n} |s_n-s_m|^{\beta_{RM}} \exp\left(-\frac{\beta_{RM}}{2} \sum_n s_n^2\right)
\end{equation}
is a stationary solution of the Fokker-Planck equation~(\ref{Eq:FP}). This distribution function corresponds to Wigner-Dyson statistics. It is equivalent to the partition function of plasma with logarithmic repulsion and a parabolic confining potential. 
\begin{figure*}
\begin{center}
\includegraphics[width=1.99\columnwidth]{figS1.pdf}\\
\caption{ \label{Fig:dnm} Normalized diagonal correlator $\langle d_{nm}\rangle /\langle d_{nn}\rangle$ as a function of $|m-n|$ for different values of disorder. Although for intermediate disorders $\langle d_{nm}\rangle$ decays as a power-law with $|n-m|$, the value of $d_{nn}$ is significantly enhanced, compared to the power-law tail. }
\end{center}
\end{figure*}
\subsubsection{Plasma model}
After discussing the random matrix example, we move on to the many-body case. There the Brownian motion is induced by a random walk over different realizations of disorder. For a specific case of XXZ spin chain in a random magnetic field, we have
\begin{equation} \label{Eq:Vt}
  V(\tau) = \sum_{i=1}^L h_i(\tau)S^z_i,
  \quad
  \langle h_i(\tau) h_j(\tau)\rangle = v^2\delta(\tau-\tau')\delta_{ij}.
\end{equation}
Repeating the procedure outlined above for the case of random matrices, we get the following Langevin equation:
\begin{equation} \label{Eq:LangevinMB}
  \frac{d s_n}{d \tau}
  =
\mu_n[\{s\},\tau]
  +  \xi_n (\tau),
\end{equation}
with drift and noise terms reading:
\begin{eqnarray} \label{Eq:driftMB}
 \mu_n[\{s\},\tau]
    &=&
   \frac{v^2}{\delta}\sum_{i=1}^L \sum_{m\neq n} \frac{c^i_{nm}(\tau)}{s_n-s_m},
  \\   \label{Eq:noiseMB}
  \langle  \xi_{n}(\tau )  \xi_m (\tau')\rangle
  &=&
   \frac{v^2}{\delta} \sum_{i=1}^L d^i_{nm}(\tau)\delta(\tau-\tau'),
\end{eqnarray} 
where the level spacing $\delta$ appeared since we assume that $s$ represents unfolded energy spectrum with mean level spacing of one. 
The correlation functions $c_{nm}$ and $d_{nm}$ are defined by the matrix elements of the operator $S^z_i$ on a site $i$ (which couples to the random magnetic field) in the instantaneous eigenbasis:
\begin{eqnarray} \label{Eq:c}
 c^i_{nm}(\tau) &=& \frac{1}{\delta}\langle n |S^z_i | m \rangle \langle m |S^z_i | n \rangle,
  \\ \label{Eq:d}
 d^i_{nm}(\tau) &=& \frac{1}{\delta}\langle n |S^z_i | n \rangle \langle m |S^z_i | m \rangle.
\end{eqnarray}
The mean field approximation employed in the main text amounts to replacing $c^i_{nm}(\tau)$ and $d^i_{nm}(\tau)$ with their average over ensemble, which becomes only a function of $n-m$. Fixing $v^2=\delta/L$ to cancel the level spacing from resulting equations, we get the following Lagnevin dynamics:
\begin{equation} \label{Eq:LangevinMF}
  \frac{d s_n}{d \tau}
  =
 \mu_n[\{s\},\tau]
  + \xi_n (\tau),
\end{equation}
with  drift and noise terms
\begin{eqnarray} \label{Eq:driftMF}
    \mu_n[\{s\},\tau]
    &=&
    \sum_{m\neq n} \frac{c(s_n-s_m)}{s_n-s_m},
  \\   \label{Eq:noiseMF}
  \langle \xi_{n}(\tau ) \xi_m (\tau')\rangle
  &=&
   d(s_n-s_m)\delta(\tau-\tau'), 
\end{eqnarray} 
 expressed via rescaled correlators:

Assuming a delta-function form of $d(\omega)$, and a power-law ansatz for $c(\omega)$,
\begin{eqnarray} \label{Eq:cw}
 c(\omega) = \frac{A}{|\omega|^\gamma}, 
 \quad 
0<\gamma<1,
\end{eqnarray}
the stationary solution of corresponding Fokker-Planck equation reproduces the partition function of plasma with interaction potential given by:
\begin{eqnarray} \label{Eq:Up}
  U(s) = \frac{C}{|s|^\gamma}.
\end{eqnarray}

\subsubsection{Numerical results for $d(\omega)$}
In order to check the validity of approximating $d(\omega)$ by a delta-function, we calculate $d_{nm}$ numerically for the XXZ spin chain. Figure~\ref{Fig:dnm} reveals the evolution of average $\langle d_{nm}\rangle$ for different system sizes and disorders. At disorder $W=0.5$, when system is in the ergodic phase, we see that $d_{nm}$ when $|m-n|>0$ does not depend on $|m-n|$, which is a manifestation of the eigenstate thermalization  hypothesis(ETH).  Nevertheless, values of $\langle d_{nm}\rangle$ for $n\neq m$ are suppressed compared to $d_{nn}$. Even nearby eigenstates can have average spin of different sign, this does not contradict to the ETH which rigorously applies to coarse-grained observables. 

For larger disorder, $W>1.5$, the $d_{nm}$ decays approximately as power-law. Nevertheless, the value of $d_{nn}$ is still considerably enhanced, hence a proper approximation for the $d(\omega)$ is $d(\omega) = c_1 \delta(\omega)+ c_2/\omega^{\gamma'}$ [scaling of the relative weight of delta function in $d(\omega)$ in the thermodynamic limit is an interesting and open question]. Presence of the delta-function contribution in $d(\omega)$ is sufficient to make the  spectral statistics different compared to the  case of Anderson transition. In particular, repeating the mean-field treatment of Ref.~\cite{ChalkerRW}, we do not get constant spectral form-factor $K(t,\tau)$ in the limit $t\to 0$. This suggests that level compressibility is vanishing, consistent with plasma model~\cite{KravtsovPlasma}. 


\subsection{Spectral rigidity and different estimates for $\gamma$ across the MBL transition}
To probe the spectral rigidity, we study the behavior of the variance of number of levels in the box of size $N$. The variance  $\mathop{\mathrm{var}} N$ as a function of the box size is shown in Fig.~\ref{Fig:varN} for $L=16$ spins, along with the power-law fits of its behavior. The exponent extracted from the fits, $\gamma_{\text{var}}$ below is compared to exponents extracted by other means.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.9\columnwidth]{F3.pdf}
\caption{ \label{Fig:varN} Variance of the level number in the box of size $N$ becomes essentially linear in $N$ for disorder above $W\geq 2.5$. To minimize the influence of the unfolding procedure, we show and use the part where results agree for both local and global unfolding. All data is for $L=16$ spins.
}
\end{center}
\end{figure}
\begin{figure}[h]
\includegraphics[width=0.9\columnwidth]{Fexp.pdf}\\
\caption{ \label{Fig:gamma-app} Exponent $\gamma$ extracted from $c(\omega)$ agrees reasonably with the exponent $\gamma_\text{var}$ from the tails of the level statistic. The exponent $\gamma_P$ governing the tails of the $P(s)$ is consistently smaller, but agrees with $\gamma'= (d \ln c(\omega)/d\ln \omega)_{\omega=0}$.}
\end{figure}


Three exponents, $\gamma, \gamma_P, \gamma_\text{var}$, obtained from matrix elements, level statistics, and level rigidity respectively,  are expected to coincide in the region of applicability of plasma model. In Fig.~\ref{Fig:gamma-app} we plot these exponents as a function of disorder for spin chains with $L=16$ spins.  We observe the reasonable agreement between the values of $\gamma_\text{var}$ and $\gamma$ measured from the tails of level statistics. On the other hand, $\gamma_P$ is consistently smaller, and matches much better the $\gamma'$, logarithmic derivative of $c(\omega)$ at $\omega\to 0$, confirming that $P(s)$ is more strongly influenced by a non-universal part in the $c(\omega)$. Both $\gamma$ and $\gamma'$ become larger than one at $W\approx 2.5$, indicating that for larger disorders the level statistics enters a semi-Poisson regime. 




\end{document}

