\section{Introduction}
Many-body localization (MBL) is understood as a distinct phase of matter that can not be described by conventional statistical physics \cite{Nandkishore15}. Driven by theoretical and experimental progress, there has been a surge of interest in the phenomenon over the last decade \cite{Abanin19}. MBL appears in sufficiently strongly disordered interacting quantum many-body systems, where the appearance of local integrals of motion \cite{Serbyn13-2, Ros15} leads to e.g. emergent integrability \cite{Huse14, Imbrie16}, the absence of thermalization \cite{Serbyn14}, and logarithmic growth of entanglement entropy in time after a quantum quench \cite{Znidaric08, Bardarson12, Serbyn13}.

Thermal and many-body localized phases are separated by an MBL transition \cite{Pal10, Luitz15, Khemani17}. At the localized side of the transition, eigenstates obey area-law scaling of entanglement entropy, while volume-law scaling is observed at the thermal side \cite{Dumitrescu17}. Entanglement entropies can be extracted from entanglement spectra \cite{Li08, Laflorencie16}. An entanglement spectrum encodes statistics beyond the entanglement entropy \cite{Vivo16}, of which several have been studied in the context of MBL \cite{Yang15, Geraedts16, Serbyn16, Yang17, Pietracaprina17, Geraedts17, Gray18}. 

The physical information encoded in the entanglement spectrum of a many-body localized eigenstate is almost fully carried by only a few elements, independent of system size \cite{Serbyn16}. This indicates the potential physical significance of the extreme value statistics \cite{Leadbetter83} of entanglement spectra in the context of MBL. In this work, we study the extreme value statistics of entanglement spectra of many-body localized eigenstates.

Extreme value statistics display universal characteristics over a wide range of physically relevant conditions \cite{Bramwell00, Antal01, Bertin05, Lakshminarayan08}. We report the surprising observation of Gumbel statistics \cite{Gumbel58, Leadbetter83}. These statistics, being observed in studies on various physical phenomena \cite{Bramwell01, Majumdar04, Katzgraber05, Hofferberth08, Lovas17}, apply to the extreme value of $n \to \infty$ independent samples drawn from a distribution with a faster than power-law asymptotic decay. Our result provides an analytical, parameter-free characterization of many-body localized eigenstates. 


\section{Gumbel statistics}
Following parts of Ref. \cite{Leadbetter83}, we here briefly discuss Gumbel statistics. Let $X_i$ ($i=1,2,\ldots, n$) be a sequence of independent and identically distributed random variables drawn from a distribution for which the distribution function (the probability that $X_i \le x$) is given by $F(x)$,
\begin{equation}
P \{ X_i \le x \} = F(x).
\end{equation}
Let $M_n$ denote the largest element of the sequence. It follows from the independence of the $X_i$ that the distribution function of $M_n$ is given by
\begin{equation}
P \{ M_n \le x \} = F^n(x).
\end{equation}
Two distribution functions $F_1$ and $F_2$ are said to be of the same type if, up to normalization,
\begin{equation}
F_2(x) = F_1(a x+b)
\label{eq: similar}
\end{equation}
for some $a > 0$ and $b$. From the extremal types theorem, it follows that if
\begin{equation}
\lim_{n \to \infty} F^n(a_n x + b_n) = G(x)
\end{equation}
for some $a_n$ and $b_n$, then $G$ is a distribution function of the same type as one of the three extreme value distributions. The distribution function of one of these extreme value distributions, relevant in the context of this work, is given by
\begin{equation}
G(x) = \exp \left( -e^{-x} \right).
\label{eq: G}
\end{equation}
This type emerges \emph{e.g.} for a density function $f = dF/dx$ asymptotically decaying faster than a power-law, i.e. as
\begin{equation}
f(x) \sim \exp(-x^\alpha)
\label{eq: tail}
\end{equation}
with $\alpha > 0$ a free parameter. Eq. \eqref{eq: tail} covers e.g. exponential ($\alpha = 1$) and Gaussian ($\alpha = 2$) decays.

The statistics of $G$ given in Eq. \eqref{eq: G} are referred to as Fisher-Tippett-Gumbel \cite{Antal01} or Gumbel \cite{Bertin05} statistics. For these statistics, the rate of convergence depends non-trivially on the structure of $F$ \cite{Gyorgyi10}. When comparing the statistics of a collection of extreme values with Gumbel statistics, it is convenient \cite{Bramwell00, Bramwell01, Antal01, Majumdar04, Katzgraber05, Bertin05, Lakshminarayan08, Hofferberth08, Lovas17} to take $a_n$ and $b_n$ such that the distribution has mean $0$ and standard deviation $1$. The distribution function of the same type as $G$ given in Eq. \eqref{eq: G} with these values for the mean and standard deviation is obtained through Eq. \eqref{eq: similar} with
\begin{equation}
a = \pi / \sqrt{6} , \qquad b = \gamma,
\end{equation}
where $\gamma \approx 0.577$ is Euler's constant. The corresponding density function is given explcitly in Eq. \eqref{eq: dens} below.


\section{Entanglement spectra}
Here, we review the concept of entanglement spectra. In the most general form, the setting is a quantum system divided into subsystems $A$ and $B$ with Hilbert space dimensions $m$ and $n$. A pure state $| \psi \rangle$ of the composite system can be expanded in basis states $| a_i \rangle$ and $| b_i \rangle$ of the respective subsystems as
\begin{equation}
| \psi \rangle = \sum_{i,j} X_{ij} \, | a_i \rangle \otimes | b_j \rangle,
\label{eq: X}
\end{equation}
where $X$ is an $m \times n$ matrix. Labeling the subsystems such that $m \ge n$, the Schmidt decomposition of $X$ uniquely expands $| \psi \rangle$ as a linear combination of product states over the subsystems,
\begin{equation}
| \psi \rangle = \sum_{i=1}^n \sqrt{\lambda_i} \; | \alpha_i \rangle \otimes | \beta_i \rangle,
\end{equation}
where $| \alpha_i \rangle$ and $| \beta_i \rangle$ are basis states for respectively subsystems $A$ and $B$, and the $\lambda_i$ ($\lambda_i \ge 0$) are the Schmidt coefficients. An element $\lambda_i$ can be interpreted as the physical weight of the product state $| \alpha_i \rangle \otimes | \beta_i \rangle$, providing a contribution of $-\lambda_i \ln(\lambda_i)$ to the entanglement entropy. The elements $e_i$ of the entanglement spectrum \cite{Li08, Laflorencie16} are given by
\begin{equation}
e_i = -\ln(\lambda_i).
\end{equation}
We remark that in some literature (e.g. Refs. \cite{Yang15, Serbyn16, Yang17}) the `quantum information definition' $e_i = \lambda_i$ of the entanglement spectrum is used. The smallest of the $e_i$ carry the largest physical weight. In this work, the focus is on the statistics of the smallest of the $e_i$.

The Schmidt coefficients for ergodic (`random') states \cite{Page93} obey the eigenvalue statistics of the fixed-trace Wishart-Laguerre random matrix ensemble \cite{Forrester10}. For this ensemble, the joint density function $ P \{\lambda_{1,2,\ldots, n} = x_{1,2,\ldots,n} \}$ of the eigenvalues is proportional to
\begin{equation}
\prod_{i=1}^n x_i^{\alpha \beta/2} \; \prod_{j<k} |x_j - x_k|^\beta \; \delta \left( \sum_{i=1}^n x_i - 1 \right),
\end{equation}
where $\alpha = (1+m-n)-2 / \beta$, and $\beta$ is the Dyson index given by $1$ or $2$ if the eigenstate is for a system with or without time-reversal symmetry, respectively. The elements are strongly correlated, due to which Gumbel statistics do not apply. For the values of $n$ relevant in the context of this work, the extreme value statistics of the smallest $e_i = -\ln(\lambda_i)$ are close to Gaussian up to $\sim 3$ standard deviations around the mean value (verified numerically).


\section{Physical setting}
We study the eigenstates of a spin chain with random onsite disorder. Let $\sigma^\alpha_i$ denote a Pauli matrix ($\alpha = x,y,z$) acting on site $i$, and let $S^\alpha_i = \sigma^\alpha_i / 2$ denote the corresponding spin-$1/2$ operator. The Hamiltonian of the model under consideration is given by
\begin{equation}
H = \sum_{i=1}^L \left( \vec{S}_i \cdot \vec{S}_{i+1} + h_i S_i^z + \Gamma S_i^x \right).
\label{eq: H}
\end{equation}
We impose periodic boundary conditions $S^\alpha_{L+1} \equiv S^\alpha_{1}$, set $\Gamma = 0.1$, and sample $h_i$ from the uniform distribution ranging over $[-W,W]$. We restrict the focus to eigenstates associated with eigenvalues close to the middle of the spectrum (quantified below) for system sizes $L=10, 12, 14$. The model is studied in Ref. \cite{Yang15}, where indications for an MBL transition at $W \approx 3.5$ are reported. 

The numerical analysis for systems of size $L=10$ or $L=12$ involves the $10$ eigenstates associated with energies closest to the middle $(\max(E_i) + \min(E_i))/2$ of the spectrum $E_i$ ($i = 1,2,\ldots, \dim(H)$), while for systems of size $L=14$ this number is set to $50$. Histograms are drawn from the data of at least $2.5 \times 10^5$ eigenstates, which corresponds to at least $25.000$ disorder realizations for $L=10$ and $L=12$, or at least $5.000$ disorder realizations for $L=14$.

For the calculation of entanglement spectra, we split the chain into subsystems $A$ and $B$ covering respectively the first and last $L/2$ sites, such that $n = 2^{L/2}$. Note that for $\Gamma = 0$ the Hamiltononian reduces to the `standard model of MBL' \cite{Oganesyan07, Pal10, Luitz15}. Then, the total spin projection
\begin{equation}
S^z = \sum_{i=1}^L S^z_i
\end{equation}
is a conserved quantity, due to which the entanglement spectrum is given by the union of independent subspectra labeled by either
\begin{equation}
S^z_A = \sum_{i \in A} S^z_i \qquad \text{or} \qquad S^z_B = \sum_{i \in B} S^z_i.
\end{equation}
This phenomenon is reflected in e.g. a block-diagonal structure of $X$ in Eq. \eqref{eq: X}.



\section{Results}
We here show the main result, namely the observation of Gumbel statistics for the entanglement spectra of many-body localized eigenstates. Let $e_\text{min} = \min_i(e_i)$ denote the smallest element of an entanglement spectrum. Because Gumbel statistics are formulated for the largest element of a sequence, we study the statistics of $-e_\text{min}$.

Let $\langle \cdot \rangle$ denote an expectation value, and let
\begin{equation}
\mu = \langle -e_\text{min} \rangle, \qquad \sigma^2 = \langle e_\text{min}^2 \rangle - \langle e_\text{min} \rangle^2
\end{equation}
denote respectively the mean and variance of the distribution of $-e_\text{min}$. We define $\tilde{e}_\text{min}$ as
\begin{equation}
\tilde{e}_\text{min} = \frac{-e_\text{min} - \mu}{\sqrt{\sigma^2}}.
\end{equation}
By construction, the distribution of $\tilde{e}_\text{min}$ has mean $0$ and standard deviation $1$. We compare the density $P \{ \tilde{e}_\text{min} = x \}$ of $\tilde{e}_\text{min}$ with the density function
\begin{equation}
g(x) = \frac{\pi}{\sqrt{6}} \exp \left[- \left( \frac{\pi}{\sqrt{6}} x + \gamma \right) - e^{-\left( \frac{\pi}{\sqrt{6}} x + \gamma \right)} \right]
\label{eq: dens}
\end{equation}
for Gumbel statistics having the same mean and standard deviation. For reference, we also compare it with the standard Gaussian, approximating the statistics of $\tilde{e}_\text{min}$ for ergodic states, for which the density function is given by
\begin{equation}
h(x) = \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{1}{2} x^2 \right).
\end{equation}

Fig. \ref{fig: Gumbel-W45} compares the density of $\tilde{e}_\text{min}$ for $W=4$ and $W=5$ at $L=10,12,14$ with $g(x)$ and $h(x)$. The density of $\tilde{e}_\text{min}$ is approximated by a histogram with bins of width $0.05$, which is normalized to unit area to allow for a direct comparison with the (normalized) probability densities. We observe good agreement with $g(x)$ for both disorder strengths at $L=12,14$. Deviations from Gumbel statistics can presumably be attributed to finite-size effects. These effects play a role in both the physics of the eigenstates (becoming stronger localized with increasing system size), as well as in the approach towards the limit $n \to \infty$ (required to observe Gumbel statistics).

\begin{figure}
\includegraphics[width= 0.9\columnwidth]{fig-Gumbel-W4}
\medskip
\includegraphics[width= 0.9\columnwidth]{fig-Gumbel-W5}
\caption{Density of $\tilde{e}_\text{min}$ for $W=4$ (top) and $W=5$ (bottom) at $L=10,12,14$, combined with the densities $g(x)$ and $h(x)$. Note the logarithmic scales on the vertical axes.}
\label{fig: Gumbel-W45}
\end{figure}

Fig. \ref{fig: Gumbel-L14} compares the density of $\tilde{e}_\text{min}$ with $g(x)$ and $h(x)$ for $W=2,3,4,5$ at $L=14$. Qualitative similarities between the density of $\tilde{e}_\text{min}$ and $g(x)$ can be observed at all disorder strengths. The eigenstate entanglement spectra of Hamiltonian \eqref{eq: H} are known to show statistics deviating from the expectation for ergodic states already at disorder strengths well below the MBL transition \cite{Yang15}. More generally, the thermal side of the MBL transition is not fully ergodic \cite{Luitz17}. We were not able to draw conclusions on the convergence of the distribution of $\tilde{e}_\text{min}$ towards Gumbel statistics with incrasing system size at the thermal side of the MBL transition, and remark that the statistics of $\tilde{e}_\text{min}$ can not be used as a probe for the MBL transition for the system sizes under consideration. Note that at the thermal side of the MBL transition the physical significance of the statistics of $\tilde{e}_\text{min}$ is presumably limited due to the vanishing physical weight in the thermodynamic limit $L \to \infty$.

\begin{figure}
\includegraphics[width= 0.9\columnwidth]{fig-Gumbel-L14} 
\caption{Density of $\tilde{e}_\text{min}$ for $W=2,3,4,5$ at $L=14$, combined with the densities $g(x)$ and $h(x)$. Note the logarithmic scales on the vertical axis.}
\label{fig: Gumbel-L14}
\end{figure}


\section{Discussion}
The observation of Gumbel statistics suggests that the largest elements of an entanglement spectrum $\tilde{e}_i$ are uncorrelated. To verify this, we study the presence of short-range correlations between the largest $\tilde{e}_i$. Short-range correlations can be probed through spacing statistics \cite{Atas13}. We order the $\tilde{e}_i$ in decreasing order (i.e. $\tilde{e}_i \ge \tilde{e}_{i+1}$), and focus on the ratio of consecutive spacings $r_1 \in [0,1]$ given by
\begin{equation}
r _1= \min \left(\frac{\tilde{e}_1 - \tilde{e}_2}{\tilde{e}_2 - \tilde{e}_3}, \frac{\tilde{e}_2 - \tilde{e}_3}{\tilde{e}_1 - \tilde{e}_2} \right).
\end{equation}
In the absence of short-range correlations, the distribution of $r_1$ obeys Poissonian statistics, for which
\begin{equation}
P \{ r_1 = x \} = \frac{2}{(1+x)^2}.
\end{equation}
Ergodic states obey Wigner-Dyson spacing statistics \cite{Mehta04}. For systems with time-reversal symmetry (Dyson index $\beta = 1$), the corresponding density of $r_1$ is well approximated \cite{Atas13} by
\begin{equation}
P \{ r_1 = x \} \approx \frac{27}{8} \frac{x+x^2}{(1+x+x^2)^{5/2}}.
\end{equation}
Fig. \ref{fig: spacing-L14} compares the density of $r_1$ for $W=2,3,4,5$ at $L=14$ with Poissonian and Wigner-Dyson spacing statistics. At all disorder strengths, the spacing statistics are close to Poissonian, indicating the (near) independence of the largest $\tilde{e}_i$.

\begin{figure}
\includegraphics[width= 0.9\columnwidth]{fig-spacing-L14} 
\caption{Densities of $r_1$ for $W=2,3,4,5$ at $L=14$, combined with the densities for Poissonian and Wigner-Dyson statistics ($\beta = 1$).}
\label{fig: spacing-L14}
\end{figure}


\section{Conclusion and outlook}
In summary, we have provided numerical evidence that the entanglement spectra of many-body localized eigenstates obey Gumbel statistics. Because the physical weight these spectra is almost fully carried by the few smallest elements, one might expect the extreme value statistics to have physical significance. Our result provides an analytical, parameter-free characterization of many-body localized eigenstates. We stress that no conclusions on the thermal side of the MBL transition can be drawn.

The main open question remaining is the physical mechanism responsible for the occurence of Gumbel statistics, and the way it can be explained in terms of phenomenological models \cite{Huse14}. A possible starting point for further investigations might be provided by the notion that an entanglement spectrum can be interpreted as the eigenvalue spectrum of the entanglement Hamiltonian
\begin{equation}
H_\text{ent} = -\ln(\rho_B),
\end{equation}
where $\rho_B = \Tr_A(| \psi \rangle \langle \psi |)$ is the partial trace of the density matrix $| \Psi \rangle \langle \Psi |$ over basis states of subsystem $A$ \cite{Li08}. One might hypothesize that the statistics of the eigenvector associated with $e_\text{min}$ carry relevant information. 


\begin{acknowledgments}
We thank Maksym Serbyn for useful discussions. This work is part of the Delta-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).
\end{acknowledgments}


