Model of level statistics for disordered interacting quantum many-body systems

TL;DR

This paper introduces a two-parameter plasma model that accurately captures the evolution of level statistics in disordered interacting quantum many-body systems across the ergodic to many-body localized transition, revealing universal spectral features.

Contribution

The paper develops and validates a novel two-parameter plasma model that describes spectral correlations and level statistics during the ergodic-MBL transition, including systems breaking time-reversal symmetry.

Findings

The $eta$-$h$ model reproduces level statistics across the transition.

Range $h$ correlates with the Thouless time and quantum chaos onset.

Model captures universal features in systems with broken time-reversal symmetry.

Abstract

We numerically study level statistics of disordered interacting quantum many-body systems. A two-parameter plasma model which controls level repulsion exponent $\beta$ and range $h$ of interactions between eigenvalues is shown to reproduce accurately features of level statistics across the transition from ergodic to many-body localized phase. Analysis of higher order spacing ratios indicates that the considered $\beta$-$h$ model accounts even for long range spectral correlations and allows to obtain a clear picture of the flow of level statistics across the transition. Comparing spectral form factors of $\beta$-$h$ model and of a system in the ergodic-MBL crossover, we show that the range of effective interactions between eigenvalues $h$ is related to the Thouless time which marks the onset of quantum chaotic behavior of the system. Analysis of level statistics of random quantum circuit which hosts chaotic and localized phases supports the claim that $\beta$-$h$ model grasps universal features of level statistics in transition between ergodic and many-body localized phases also for systems breaking time-reversal invariance.

Figures

Equations

$$r_i^{(n)}=\min \{\frac{E_{i+2n}-E_{i+n}}{E_{i+n}-E_i},\frac{E_{i+n}-E_i}{E_{i+2n}-E_{i+n}}\}$$

$${\mathcal{P}}_h^{\beta }(E_1,...,E_N)=Z_N^{-1}\prod _{i=0}^N|E_i-E_{i+1}{|}^{\beta }...|E_i-E_{i+h}{|}^{\beta }.$$

$${\Sigma }^2(L)=\chi L,$$

$$P(r^{(n)})={\mathcal{N}}_{n,\beta }\frac{(r^{(n)}{)}^{\beta +(n-1)(\beta +1)}}{(1+r^{(n)}{)}^{2(\beta +1)n}}$$

$$K(\tau )=\frac{1}{Z}⟨{\left|\sum _{j}g({ϵ}_j){\mathrm{e}}^{-i{E}_j\tau }\right|}^2⟩,$$

$$H=J\sum _{i=1}^L{\vec{S}}_i\cdot {\vec{S}}_{i+1}+\sum _{i=1}^L{h}_i{S}_i^z,$$

$$P(s)=C_1{s}^{\beta }{\mathrm{e}}^{-C_2{s}^{2-\gamma }}\mathrm{and}{\Sigma }_2(L)\propto L^{\gamma }$$

$$H_B=-J\sum _{⟨i,j⟩}{\hat{a}}_i^{†}{\hat{a}}_j+\frac{U}{2}\sum _i{\hat{n}}_i({\hat{n}}_i-1)+\sum _i{\mu }_i{\hat{n}}_i,$$

$$W_{a_1,\dots ,a_L;a_1^{\prime },\dots ,a_L^{\prime }}=U_{a_1,a_1^{\prime }}^{(1)}\dots U_{a_L,a_L^{\prime }}^{(L)}{\mathrm{e}}^{i{\sum }_n{\phi }_{a_n^{\prime },a_{n+1}^{\prime }}},$$

Full text

I Introduction

Many-body localization (MBL) Gornyi et al. (2005); Basko et al. (2006) manifesting ergodicity breaking in disordered interacting quantum many-body systems Huse et al. (2014); Nandkishore and Huse (2015); Alet and Laflorencie (2018); Abanin et al. (2019a) has attracted a vivid attention over the last decade. Important results include an emergent integrability of MBL phase due to the existence of local integrals of motion (LIOMs) Serbyn et al. (2013); Huse et al. (2014); Ros et al. (2015); Imbrie (2016); Mierzejewski et al. (2018) and the associated logarithmic growth of the bipartite entanglement entropy after a quench from a separable state Žnidarič et al. (2008); Bardarson et al. (2012). A wide regime of subdiffusive transport on the ergodic side of the transition was found Bar Lev et al. (2015); Luitz and Bar Lev (2016); Mierzejewski et al. (2016). Signatures of MBL have been observed experimentally in 1D Schreiber et al. (2015); Smith et al. (2016) and in 2D system Choi et al. (2016) see, however, De Roeck and Huveneers (2017). Recently, the very existence of MBL in the thermodynamic limit has been questioned Šuntajs et al. (2019) opening a new debate Abanin et al. (2019b); Sierant et al. (2019a); Panda et al. (2019). While the status of MBL in the thermodynamic limit is of utmost importance for the understanding of this phenomenon from purely theoretical viewpoint, the real systems studied in this respect are finite Kondov et al. (2015); Schreiber et al. (2015); Choi et al. (2016); Lüschen et al. (2017) often reaching very modest sizes that enable precise studies Roushan et al. (2017); Lukin et al. (2019); Rispoli et al. (2019). In this work we concentrate on systems of such a size.

Spectral statistics of ergodic systems with (without) time reversal invariance follow predictions of Gaussian orthogonal (unitary) ensemble (GOE, GUE, respectively) of random matrices Mehta (1990); Haake (2010) while eigenvalues of localized systems are uncorrelated resulting in Poisson statistics (PS). A ratio of consecutive spacings between energy levels

[TABLE]

was proposed as a simple probe of the level statistics in Oganesyan and Huse (2007) with $n=1$ and employed in investigation of ergodicity breaking in various settings Pal and Huse (2010); Santos and Rigol (2010); Mondaini and Rigol (2015); Luitz et al. (2015, 2016); Janarek et al. (2018); Wiater and Zakrzewski (2018); Macé et al. (2019); Sierant et al. (2019b). Higher order spacing ratios ($n>1$), studied in Chavda et al. (2014); Tekur et al. (2018); Kota and Chavda (2018); Bhosale (2019); Buijsman et al. (2019), are valuable tools to assess properties of level statistics. In contrast to standard measures such as level spacing distribution or number variance Mehta (1990); Haake (2010) they do not require the so called unfolding, i.e. the procedure of setting the density of energy levels $\rho(E)$ to unity which can lead to misleading results Gómez et al. (2002). Recently, an analytical understanding of an appearance of random matrix theory statistics in systems without a clear semiclassical limit have been developed in a periodically driven Ising models Kos et al. (2018); Bertini et al. (2018) or in random Floquet circuits Chan et al. (2018a). Variants of such systems have been argued to undergo ergodic-MBL transition Chan et al. (2018b); Braun et al. (2019).

In this work we introduce a two-parameter $\beta$-$h$ model which assumes a level repulsion determined by exponent $\beta$ between $h$ neighboring eigenvalues. Our model is a natural extension of the so called $\beta$-Gaussian model Buijsman et al. (2019) claimed to represent the level statistics in the transition to MBL. We show that the second parameter, the interaction distance $h$ is essential for understanding the transition and reproducing the numerical results obtained for various physical models. In particular, we demonstrate that distributions of higher order spacing ratios $r^{(n)}$ across the ergodic-MBL transition in disordered XXZ spin chain are faithfully captured by $\beta$-$h$ model and the obtained $\beta$ and $h$ parameters provide a simple perspective on short-range and long-range spectral correlations. The latter, captured effectively by the interaction range $h$, are further investigated by means of the SFF revealing a link between $h$ and the Thouless time. We demonstrate that $\beta$-Gaussian model fails to describe long-range spectral correlations. An analysis of a local Haar-random unitary nearest-neighbor quantum circuit system introduced in Chan et al. (2018b) indicates that also in such a generic system the spectral statistics can be grasped with the $\beta$-$h$ model demonstrating the robustness of observed features of level statistics.

The paper is organized as follows. In Section II we introduce the $\beta$-$h$ model and discuss properties of its level statistics. In Section III we show that the $\beta$-$h$ model accurately reproduces level statistics of disordered XXZ spin chain across the many-body localization transition. In Section IV we show that the $\beta$-$h$ model grasps also level statistics of disordered Bose-Hubbard model. In Section V we demonstrate that $\beta$-$h$ applies also to ergodic-MBL transition in systems with broken time-reversal symmetry and without local conservation laws by considering level statistics of a random quantum circuit. We conclude in Section VI.

II The $\beta$-$h$ model

The joint probablity density function (JPDF) of eigenvalues of matrix from GOE (GUE) with $\beta=1$ ($\beta=2$) can be written as a partition function of a fictitious 1D gas of particles $\mathcal{P}(E_{1},...,E_{N})=Z_{N}^{-1}\mathrm{e}^{-\beta\mathcal{E}(E_{1},...,E_{N})}$ where $Z_{N}$ is a normalization constant and the energy $\mathcal{E}$ includes a trapping potential $U(E)\propto E^{2}$ and pairwise logarithmic interactions $V(|E-E^{\prime}|)=-\log(|E-E^{\prime}|)$. Eigenvalues $E_{1}<...<E_{N}$ lie on a ring of length $N$ which confines them rendering the trapping potential $U(E)$ unnecessary. The JPDF can be written as

[TABLE]

The GOE (GUE) case is obtained when $h\rightarrow\infty$ with the appropriate value of $\beta$. The form of (2) suggests various models of intermediate level statistics between GOE (GUE) and PS. For instance, one can keep $h\rightarrow\infty$ and vary $\beta$, obtaining the so called $\beta$-Gaussian ensemble Buijsman et al. (2019). When $h$ is an integer number which sets the number of correlated neighboring eigenvalues one arrives at the so called short-range plasma model introduced in Bogomolny, E. et al. (2001) (see also Pandey et al. (2019); Kumar et al. (2019)).

In this work we extend this model by allowing $h$ to be a real number. Denoting by $\left\lfloor{.}\right\rfloor$ the floor function, the factor in (2) becomes $|E_{i}-E_{i+1}|^{\beta}...|E_{i}-E_{i+\left\lfloor{h}\right\rfloor}|^{\beta}|E_{i}-E_{i+\left\lfloor{h}\right\rfloor+1}|^{\beta(h-\left\lfloor{h}\right\rfloor)}$, hence defining the $\beta$-$h$ model where $h\in[1,\infty)$ and $\beta\in[0,1]$ ($\beta\in[0,2]$) for GOE(GUE)-PS transition. Varying continuously $h$ and $\beta$ allows us to capture spectral statistics of disordered quantum many-body systems across the ergodic-MBL transition, while a simple form of JDPF of $\beta$-$h$ model yields insight into correlations between eigenvalues. Semi-analytical results for $\beta-h$ model are available only for integer values of $h$ and $\beta$ Bogomolny, E. et al. (2001). In particular, the number variance, defined as the variance of the number of eigenvalues in an interval $(E,E+L)$ reads

[TABLE]

for $L\gg 1$ where $\chi=1/(\beta h+1)$. The spectral rigidity of GOE (GUE) which manifests itself in the logarithmic growth of the variance $\Sigma^{2}(L)$ is replaced by a finite spectral compressibility $\chi$. Thus, a profound change in long-range spectral correlations happens when $h<\infty$. Interestingly, we find that (3) is fulfilled with an excellent agreement for $\beta$-$h$ model as our Monte Carlo simulations (obtained sampling JPDF of $\beta$-$h$ model with the Metropolis-Hastings algorithm Hastings (1970)) shows for arbitrary real $\beta\in[0,2]$ and $h\in[1,40]$.

A straightforward application of the method of Bogomolny, E. et al. (2001) shows that distributions of higher order spacing ratios $P(r^{(n)})$ for $h=1$ are given by

[TABLE]

where $\mathcal{N}_{n,\beta}=[_{2}F_{1}(n(1+\beta),2n(1+\beta),1+n(\beta+1),-1)/(\beta+1)n]^{-1}$ is a normalization constant and ${}_{2}F_{1}$ is Gauss hypergeometric function. Such distributions of higher order spacing ratios $P(r^{(n)})$ at $h=1$ constitute a very good approximation for systems close to the MBL phase where $h\approx 1$ and provides analytical expressions for average higher order spacing ratios $\overline{r}^{(n)}$ (including PS for $\beta=0$). To obtain $P(r^{(n)})$ for arbitrary $h\in[1,\infty)$ and $\beta\in[0,1]$ we again sample JPDF of $\beta$-$h$ model with the Monte Carlo approach.

II.1 Spectral form factor of $\beta$-$h$ model

Consider the spectral form factor (SFF) Mehta (1990); Haake (2010):

[TABLE]

where $Z$ assures that $K(\tau)\stackrel{{\scriptstyle\tau\rightarrow\infty}}{{\rightarrow}}1$, the spectrum is unfolded (for remarks on unfolding see Appendix A) and $g(\epsilon)$ is a Gaussian function which vanishes at the edges of spectrum reducing their influence – see also Appendix. The SFF allows to identify two important time scales in disordered systems: the Heisenberg time $\tau_{H}$ equal to the inverse level spacing, beyond which the discrete nature of the energy spectrum manifests itself and the Thouless time $\tau_{Th}$, which is the time scale beyond which SFF admits universal GOE (GUE) form $K(\tau)\approx 2\tau$ Chan et al. (2018b); Šuntajs et al. (2019). The existence of two time scales is reflected in the JPDF of $\beta$-$h$ model, where the correlations between eigenvalues are of the GOE (GUE) form on energy scales smaller than $h$ level spacings so that $\tau_{Th}$, inversely proportional to $h$, (for $\beta=1,2$) provides a physical interpretation of the interaction range $h$ in $\beta$-$h$ model.

SFF of $\beta$-$h$ model is shown in Fig. 1. For $\beta=1$, SFF of $\beta$-$h$ model follows prediction for GOE down to Thouless time $\tau_{Th}$ which depends on the interaction range $h$ (roughly: $\tau_{Th}\approx 2/(h+1)$). SFF for $\beta<1$ shows that it is possible to have spectral statistics with $h>1$ and $\tau_{Th}=\tau_{H}=1$.

We note that (3) implies that $K(0)=1/(\beta h+1)$ – an analytical prediction for integer $\beta$ and $h$ which is very well confirmed by numerical data for arbitrary $\beta$ and $h$ as shown in Fig. 1.

Fig. 2 shows comparison of SFFs of $\beta-h$ model and of $\beta$-Gaussian ensemble. The two parameters of the $\beta-h$ model allow to reproduce the typical behavior of SFF of a disordered many-body system across ergodic-MBL transition Šuntajs et al. (2019); Sierant et al. (2019a) – the Thouless time $t_{Th}$ of a physical quantum many-body system increases with disorder strength $W$, which, in $\beta$-$h$ model is reflected by a decreasing range $h$ of eigenvalue interactions. The vanishing level repulsion in the vicinity of the localized phase is reflected by sufficiently small value of the exponent $\beta$. This is not the case for the $\beta$-Gaussian ensemble. As soon as $\beta<1$, the deviation from the SFF of GOE is observed for all $\tau<1$ as Fig. 2 illustrates. Thus, the $\beta$-Gaussian ensemble is unable of reproducing the typical behavior of SFF, $K(\tau)$, in disordered system in which $K(\tau)$ deviates from the universal GOE curve only for times smaller than the Thouless time $t_{Th}$. Moreover, for $\beta<1$ the predictions of $\beta$-Gaussian ensemble fail to reproduce a small $\tau$ behaviour showing a rapid decrease with decreasing $\tau$ instead of a saturation as expected closer to Poisson regime.

III XXZ spin chain

Le us go beyond a comparison of the statistical models among themselves and compare their predictions with different physical models. As a starting point for testing purposes we consider a standard disordered XXZ spin-1/2 chain with Hamiltonian given by

[TABLE]

where $\vec{S}_{i}$ are spin-1/2 matrices, $J=1$ is fixed as the energy unit, periodic boundary conditions are assumed $\vec{S}_{L+1}=\vec{S}_{1}$ and $h_{i}\in[-W,W]$ are independent, uniformly distributed random variables. The model (6) has been widely studied in the MBL context Žnidarič et al. (2008); Pal and Huse (2010); Berkelbach and Reichman (2010); Luitz et al. (2015); Agarwal et al. (2015); Bera et al. (2015); Enss et al. (2017); Bera et al. (2017); Herviou et al. (2019); Colmenarez et al. (2019), and its level statistics have been addressed in Serbyn and Moore (2016); Bertrand and García-García (2016); Kjäll (2018); Sierant and Zakrzewski (2019). Recently, the $\beta$-Gaussian ensemble was suggested to describe the ergodic-MBL transition Buijsman et al. (2019). As the analysis of Section II has shown this claim is questionable. Further, we shall show that the $\beta$-Gaussian ensemble reproduces level correlations only on a single level spacing scale while missing longer-range spectral correlations. Both aspects of level statistics are grasped by $\beta-h$ model.

Eigenvalues of the XXZ spin chain (6) are obtained by an exact diagonalization for small sizes or with shift-and-invert method Pietracaprina et al. (2018) for $L=18,20$. For each $W$ we accumulate eigenvalues from $2000$ ($400$) disorder realizations for $L\leq 18$ ($L=20$). The higher order spacing ratios (1) are calculated using $500$ eigenvalues from the middle of the spectrum. The resulting exemplary distributions $P(r^{(n)})$ of higher order spacing ratios for $n=1,3,5,8$ are shown in Fig. 3. Parameters for $\beta$-$h$ model are obtained by minimizing the deviation between $P(r^{(n)})$ distributions for XXZ spin chain and $\beta$-$h$ model. A very good agreement between the distributions obtained for the model (6) and predictions of $\beta$-$h$ model is observed in the whole transition region between the ergodic and MBL phases. Note that both parameters $\beta$ and $h$ are needed to reproduce $P(r^{(n)})$ distributions for $n\geq 1$. To demonstrate that the agreement between $\beta$-$h$ model and level statistics of XXZ spin chain in ergodic-MBL crossover persists to larger energy scales, we calculate $\Delta\overline{r}^{(n)}=\overline{r}^{(n)}-\overline{r}^{(n)}_{PS}$, where $\overline{r}^{(n)}$ is the average value of $n$’th order spacing ratio $\overline{r}^{(n)}$ and $\overline{r}^{(n)}_{PS}$ is the $n$’th order average gap ratio for PS. The resulting values of $\Delta\overline{r}^{(n)}$ as function of $n$ are shown in Fig. 4. Even though the parameters of $\beta$-$h$ model are determined by fit of $P(r^{(n)})$ for $n=1,3,5,8$ only (fits for all $n$ taken with the same weight), the good agreement between $\Delta\overline{r}^{(n)}$ for XXZ spin chain and for $\beta$-$h$ model persists up to $n=50$. Interestingly, on the ergodic side of the transition, for $W\leq 2.4$, the values of $\Delta\overline{r}^{(n)}$ for $n\geq 20$ predicted by $\beta$-$h$ model are consequently overestimating the values for XXZ spin chain. Since the larger value of $\Delta\overline{r}^{(n)}$ implies stronger level correlations at the scale determined by $n$, this means that energy levels of $\beta$-$h$ model, not coupled directly in JPDF (6), are still correlated more strongly than energy levels of the system across the ergodic-MBL transition. For comparison, we also show the predictions of $\beta$-Gaussian ensemble Buijsman et al. (2019) in Fig. 4. Only the values of $\Delta\overline{r}^{(1)}$ are well reproduced by this approach, for $n\geq 2$), the values of $\Delta\overline{r}^{(n)}$ are overestimated showing that finite $h$ is an essential feature of level statistics in the MBL transition.

III.1 Scaling of level repulsion exponent $\beta$ and range of interactions $h$ at the ergodic-MBL transition

The $\beta$ and $h$ parameters characterizing level statistics across the ergodic-MBL transition are shown in Fig. 5. In the ergodic phase, at small disorder strengths $W$, GOE describes level statistics well, hence $\beta=1$ and $h\rightarrow\infty$. Upon increase of $W$, the range of interactions $h$ and the level repulsion exponent $\beta$ decrease leading to PS for sufficiently strong disorder. Notably, the system size dependences of $h(W)$ and $\beta(W)$ are very different. The data for $\beta(W)$ collapse upon rescaling $W\rightarrow(W-W_{C})L^{1/\nu}$ with $W_{C}\approx 3.4$ and $\nu\approx 1$, similarly to the average gap ratio $\overline{r}^{(1)}$ Luitz et al. (2015) indicating that $\overline{r}^{(1)}(W)$ and $\beta(W)$ contain similar information. In particular, both measures lead to the exponent $\nu<2$ violating the Harris bound Harris (1974); Chayes et al. (1986); Chandran et al. (2015). On the other hand, data for $h(W)$ collapse upon rescaling $W\rightarrow W/L$. As inset in Fig. 5 d) demonstrates, the decrease of the level repulsion exponent $\beta$ in the transition region is accompanied by the interaction range $h$ increasing with $L$ for a given value of $\beta$. Therefore, our data indicate the presence of the transition to MBL phase with vanishing level repulsion exponent $\beta=0$ at disorder strength $W_{C}$ even though the interaction range admits a certain fixed value $h_{0}=h(W^{*})$ at disorder strength $W^{*}$ increasing linearly with the system size $L$. The linear dependence $W^{*}\sim L$ persists steadily up to the largest available system size $L=24$ but since the values of $h$ in the transition region do not exceed $10$ we cannot conclude whether $h$ diverges or stays finite at the transition in the thermodynamic limit. Nevertheless, in either case, the level repulsion vanishes at the transition in $L\rightarrow\infty$ limit, in accordance with recent phenomenological treatments Thiery et al. (2018); Morningstar and Huse (2019). We note that a linear with system size dependence for deviation of $\overline{r}^{(1)}$ from the value characteristic for GOE was found recently in Šuntajs et al. (2019). This observation is related to our finding that $W^{*}\sim L$ since disorder strength for which $h$ becomes of the order of e. g. $h_{0}=10$ is, at the same time, the moment for which $\overline{r}^{(1)}$ departs from the value characteristic for GOE. While $\beta-h$ model puts the long-range spectral statistics examined in Šuntajs et al. (2019) in another perspective, we emphasize that the presented data suggest transition to MBL phase at disorder strength $W_{C}$ in the thermodynamic limit.

III.2 The spectral form factor XXZ spin chain

Let us now consider the SFF of XXZ spin chain which is shown in Fig. 6 along with the predictions of $\beta$-$h$ model. Beyond the Heisenberg time $\tau_{H}$, $K(\tau)=1$. For smaller $\tau$, the SFF of XXZ spin chain follows the GOE prediction down to the Thouless time $\tau_{Th}$ which increases monotonically with disorder strength $W$. The behavior is captured by SFF of the $\beta$-$h$ model. For $\tau<\tau_{Th}$, an increase in the SFF of XXZ chain is observed for disorder strengths $W$ corresponding to the ergodic side of the transition, whereas SFF remains constant for the $\beta$-$h$ model. The latter behavior signals weak correlations between eigenvalues of $\beta$-$h$ model beyond energy scale determined by $h$, whereas the behavior of SFF of XXZ spin chain indicates even weaker correlations of its eigenvalues.

III.3 Level statistics and number variance across the ergodic-MBL transition

We revisit now the level spacing distribution and the number variance of the disordered XXZ spin chain in the ergodic-MBL crossover, shown in Fig. 7. Level spacing distributions are very faithfully reproduced by $\beta$-$h$ model in the whole crossover regime. There are, however slight deviations in the number variance $\Sigma^{2}(L)$ of XXZ spin chain and $\beta$-$h$ model. On the ergodic side of the crossover ($W<2.4$) the $\beta$-$h$ model underestimates number variance of XXZ spin chain indicating weaker long-range spectral correlations of the latter, in agreement with the analysis of $\Delta\overline{r}^{(n)}$ in this regime. For large $W$, the prediction of $\beta$-$h$ model overestimates the number variance of XXZ spin chain –that is probably related to effects of a finite number of eigenvalues $n_{e}$ from a single disorder realization which are known to contribute as $-L^{2}/n_{e}$ to the number variance $\Sigma^{2}(L)$ Sierant and Zakrzewski (2019).

We note that results for the number variance $\Sigma^{2}(L)$ are strongly dependent on the way the level unfolding is performed, in view of that we conclude that the higher order spacing ratios are more reliable in extracting informations about spectral correlations beyond single level spacing. Furthermore, when one inspects tails $s\gtrsim 4$ of the level spacing distribution on the logarithmic scale, deviations from predictions of $\beta$-$h$ model are found. As it was demonstrated in Sierant and Zakrzewski (2019), such a behavior at large $s$ is associated with the large inter-sample randomness associated with ergodic-MBL transition in random potentials.

The weighted short-range plasma model for level statistics in ergodic-MBL crossover was considered in Sierant and Zakrzewski (2019). This model takes into account inter-sample randomness, an important feature of MBL transition in random potentials Khemani et al. (2017). The JPDF of weighted short-range plasma model is a weighted superposition of JPDF of the form (2) and as such it is related to $\beta$-$h$ model. However, the necessity of reproducing the inter-sample randomness requires an introduction of many weight parameters. This makes use of the weighted short-range plasma model complicated. The simple picture of changes of interaction range between eigenvalues and its relation to Thouless time cannot be easily extracted due to the complexity of the model. It must be noted however, that taking into account the inter-sample randomness determined by a sample-averaged spacing ratio Sierant and Zakrzewski (2019) could diminish the (small) deviations in $P(r^{(n)})$ between $\beta$-$h$ model and XXZ spin chain for disorder strengths $W=2.4,2.6$ for which the inter-sample randomness is the largest at $L=18$.

A two-stage Serbyn and Moore (2016) picture of flow of level statistics between GOE and PS proposes that on the ergodic side of the crossover level statistics are described by a plasma model with power-law interactions between eigenvalues which yields the following expressions for the level spacing distribution and the number variance:

[TABLE]

with $C_{1,2}$ determined by normalization conditions $\langle 1\rangle=\langle s\rangle=1$. The exponent $\beta$ and $\gamma$ play a role similar to $\beta$ and $h$ parameters of the $\beta$-$h$ model. In the transition from extended to localized regime in the first stage $\gamma$ changes from 0 to 1 leading to Poissonian tail of $P(s)$ followed at the second stage by a change of level repulsion $\beta$. However, as demonstrated in Bertrand and García-García (2016), the predictions of (7) are not valid as the number variance $\Sigma^{2}(L)$ in ergodic-MBL transition grows linearly (or superlinearly) – see Fig. 7, contrary to prediction of (7) where $0<\gamma<1$ in the crossover regime. Moreover, (7) is obtained on the mean-field level Kravtsov et al. (1994), no other predictions for this model such as JPDF are available. The second stage of the flow Serbyn and Moore (2016) coincides with $\beta$-$h$ model with $h=1$. However, as shown above, the $h(\beta)$ dependence is such that the interaction range $h$ for fixed value of $\beta$ is increasing, hence $h$ becomes equal to unity only deep in the MBL regime.

IV Level statistics of disordered Bose-Hubbard model

To provide further evidence that $\beta$-$h$ model is able to reproduce level statistics of interacting disordered quantum many-body systems, we analyze higher order spacing ratios in ergodic-MBL transition in a disordered Bose-Hubbard model Sierant et al. (2017); Sierant and Zakrzewski (2018) with Hamiltonian:

[TABLE]

where $a^{{\dagger}}_{i},a_{i}$ are bosonic creation and annihilation operators respectively, the tunneling amplitude $J=1$ sets the energy scale, $U=1$ is interaction strength and the chemical potential $\mu_{i}$ is distributed uniformly in an interval $[-W;W]$. This model undergoes transition to MBL phase beyond critical disorder strength $W_{C}$ which depends on interaction strength $U$.

Distribution of higher order spacing ratios $(n=1,3,5,8)$ for the disordered Bose-Hubbard model are shown in Fig. 8. The $\beta$-$h$ model reproduces faithfully distributions $P(r^{(n)})$ in the whole crossover regime. We note that the dependence $h(\beta)$ is markedly different as compared to the XXZ spin chain – here we find $h=2$ even when the level repulsion exponent $\beta$ is close to [math]. Average higher order spacing ratios shown in Fig. 9 indicate that long-range spectral statistics are also well reproduced by the $\beta$-$h$ model. In particular, the tendency of $\beta$-$h$ model to overestimate long-range spectral correlations in XXZ spin chain is reversed in the case of Bose-Hubbard model, indicating that this is a model dependent feature.

V Random quantum circuit

Consider 1D chain of $q$-level systems of length $L$ with Floquet operator given by Friedman et al. (2019a)

[TABLE]

where $U^{(j)}$ are unitary matrices that generate rotations at each site, chosen independently from Haar distribution, $\varphi_{a_{n},a_{n+1}}$ are independent Gaussian random variables with zero mean and standard deviation $\epsilon$ that determine coupling between neighboring sites. The SFF is related to the Floquet operator via $K(t)=\left\langle\mathrm{Tr}[W^{t}]\mathrm{Tr}[(W^{{\dagger}})^{t}]\right\rangle$ where $t$ is an integer and (5) is recovered with $g(\epsilon)=1$ for $\tau\propto t$. Analytic calculation Friedman et al. (2019a) in the limit $q\rightarrow\infty$ shows that the system is chaotic in the thermodynamic limit and SFF follows the prediction for GUE: $K(\tau)=2\tau$. For $q=3$, numerical calculations indicate that the system undergoes a transition between ergodic phase at $\epsilon\gtrsim 0.25$ where the statistics of eigenphases $\theta_{j}$ are well described by GUE and MBL phase at $\epsilon\lesssim 0.25$ with PS statistics. We now turn to analysis of level statistics of (9) at finite $L$ and $q=3$.

Distributions of higher order spacing ratios (1) calculated for eigenphases $\theta_{j}$ are shown in Fig. 10. The $\beta$-$h$ model with level repulsion exponent $\beta\in[0,2]$ and appropriately chosen range of interactions $h$ reproduces the distributions of higher order spacing ratios $P(r^{(n)})$, despite the broken time-reversal symmetry in the system.

Average higher order gap ratios for the random quantum circuit are shown in Fig. 11. The $\beta$-$h$ model gives a good account for the spectral correlations reflected by $\Delta\overline{r}^{(n)}$. Notably, deviations at $n\gtrsim 20$ suggest also in this case that correlations between eigenphases of the Floquet operator $W$ in the crossover regime are weaker than correlations predicted by the $\beta$-$h$ model.

This suggests a similar behavior of level statistics at larger energy scales as in the case of XXZ spin chain. Fig. 12 shows the SFF of the considered Floquet operator (9) with predictions of $\beta$-$h$ model. The behavior of SFF is qualitatively very similar to the case of XXZ spin chain, $K(\tau)$ follows the prediction for GUE down to the Thouless time $\tau_{Th}$; for smaller $\tau$, $K(\tau)$ flattens matching SFF of $\beta$-$h$ model. On the ergodic side of transition SFF of the Floquet operator increases indicating weaker correlations between eigenvalues than in the $\beta$-$h$ model.

VI Discussion and outlook

We have analyzed level statistics across the ergodic-MBL transition. The proposed $\beta$-$h$ model provides a simple framework that allows one to reproduce universal features of level statistics of disordered interacting quantum many-body systems. The model captures the ergodic-MBL transition in the XXZ spin chain. Similarly, $\beta$-$h$ model is able to reproduce level statistics of disordered Bose-Hubbard models that undergo ergodic-MBL transition Sierant et al. (2017); Sierant and Zakrzewski (2018). The $\beta$-$h$ model grasps also level statistics of the random quantum circuit across the transition between ergodic and MBL phases in spite of broken time-reversal symmetry. Notably, the only feature encoded in the Floquet operator (9) is the locality of gates in the circuit and as such the random circuit can be regarded as a toy model of a generic disorder interacting quantum system. All this taken together allows us to conjecture that $\beta$-$h$ model grasps universal, robust features of level statistics of interacting disordered quantum many-body systems, independently, for instance, of local conservation laws Huang (2019); Friedman et al. (2019b).

The transition between chaotic and integrable regimes in systems with chaotic classical counterparts Bohigas et al. (1984); Haake (2010) is system specific as it is determined by the structure of underlying classical phase space Bohigas et al. (1993). Our analysis with $\beta$-$h$ model indicates that spectra of disordered interacting quantum many-body systems are effectively parametrized with good accuracy only by the level repulsion exponent $\beta$ and the range of interactions between eigenvalues $h$ suggesting an existence of a robust mechanism of delocalization of LIOMs that assure integrability and PS statistics in MBL phase. Detailed understanding of such a mechanism remains an open problem. Spectral properties of disordered interacting many-body systems resemble level statistics at the single particle Anderson localization transition Shklovskii et al. (1993); Evers and Mirlin (2008) which could be expected as MBL can be regarded as an Anderson localization in the Hilbert space Luca and Scardicchio (2013); Logan and Welsh (2019); Ghosh et al. (2019).

The $\beta$-$h$ model is capable of reproducing distributions of higher order spacing ratios, level spacing distributions, number variance of systems across ergodic-MBL transition. The range of interactions between eigenvalues $h$ sets the Thouless time $\tau_{Th}$ at which the SFF deviates from the universal RMT predictions. It would be interesting to compare this time scale to Thouless time extracted from matrix elements of local operators Luitz and Lev (2017); Serbyn et al. (2017) or from the return probability Schiulaz et al. (2019). The considered $\beta$-$h$ model can be also used to probe the entanglement spectrum in MBL systems Geraedts et al. (2016) or random fractonic circuits Pai et al. (2019) as it has been shown to hosts similar, local correlations between energy levels. It would be interesting to relate it to the associated multifractality observed deep in the MBL phase Macé et al. (2018) or to properties of level dynamics across the ergodic-MBL transition studied recently in Maksymov et al. (2019).

VII acknowledgments

We are most grateful to Fabien Alet for kindly sharing with us the eigenvalues for $L=22,24$ Heisenberg spin chain as well as discussions on subjects related to this work. We thank Dominique Delande for discussions and a careful reading of this manuscript. P. S. and J. Z. acknowledge support by PL-Grid Infrastructure. This research has been supported by National Science Centre (Poland) under projects 2015/19/B/ST2/01028 (P.S.), 2018/28/T/ST2/00401 (doctoral scholarship – P.S.) and 2016/21/B/ST2/01086 (J.Z.).

Appendix A Remarks on unfolding

One of the advantages of analysis of level statistics with higher order spacing ratios $r^{(n)}$ is that they do not require spectral unfolding, i.e. the level density $\rho(E)$ cancels out. This is of course valid only when $n$ is such that $\rho(E_{i})$ and $\rho(E_{i+2n})$ are not significantly different which seems to be a plausible assumption when dimension of Hilbert space is larger than few thousands.

The calculation of SFF of XXZ spin chain requires application of spectral unfolding. To this end we consider $40000$ of eigenvalues from the center of spectrum and fit the level staircase function Haake (2010) with a polynomial of degree $10$. To calculate $K(\tau)$ we use $g(E)\propto\exp\left(-(E-\bar{E})^{2}/(0.18\Delta E^{2})\right)$ (following Šuntajs et al. (2019)) where $\bar{E}$ is average of the ground state and highest excited state energies and $\Delta E$ is standard deviation of energy in given spectrum.

In order to obtain level spacing distribution and the number variance of XXZ spin chain we consider $500$ eigenvalues from the middle of the spectrum and we perform unfolding by fitting the level staircase function with a third order polynomial.

Eigenphases $\theta_{j}$ of the random quantum Haar-measured circuit are distributed uniformly in interval $[0,2\pi]$, hence no unfolding is required and SFF can be calculated directly from $K(t)=\left\langle\mathrm{Tr}[W^{t}]\mathrm{Tr}[(W^{{\dagger}})^{t}]\right\rangle$.