Characterization of many-body mobility edges with random matrices

A single-particle mobility edge characterizing the separation in energy between extended and localized states occurs in a variety of non-interacting quantum models, as in one-dimensional ones with long-range hoppings Biddle and Das Sarma (2010); Biddle et al. (2011), certain incommensurate modulations of the potential Soukoulis and Economou (1982); Das Sarma et al. (1988, 1990); Ganeshan et al. (2015); Li et al. (2017a), or in three-dimensional disordered systems at a finite disorder strength Bul$\cdot$ka et al. (1985); Plyushchay et al. (2003). The fate of mobility edges in the presence of interactions is of great interest to experimentalists and theorists, but as of today, a large debate questions the existence of its many-body analogue. A variety of numerical results show the manifestation of many-body mobility edges in finite sizes, and some studies employing scaling analysis also argue that they should occur when approaching the thermodynamic limit Laumann et al. (2014); Luitz et al. (2015); Li et al. (2015); Kjäll et al. (2014); Modak and Mukerjee (2015); Baygan et al. (2015); Bar Lev et al. (2015); Serbyn et al. (2015); Mondragon-Shem et al. (2015); Nag and Garg (2017); Wei et al. (2019). Roeck et al., for example, contend this view, asserting that localized states will be eventually thermalized by ergodic ones, ruling out their concomitant appearance and consequently the occurrence of many-body mobility edges De Roeck et al. (2016); De Roeck and Huveneers (2017). Thus, due to the size limitations used in all numerical calculations, consensus has not been achieved.

To provide a new angle in this controversy, we will step back, and rather than following the standard procedure of analyzing a typical physical model, we will head towards the theory of random matrices, which provides the basis of our understanding of quantum chaotic behavior and its eventual absence Haake (2006); L. D’Alessio et al. (2016). These are naturally characterized by the rigidity of the matrix spectrum, with the former displaying a characteristic level repulsion that is not present in regular systems, which in turn, exhibit completely uncorrelated eigenvalues. When defining a gap between adjacent eigenvalues in the spectrum $S_{j}\equiv E_{j+1}-E_{j}$, chaotic and non-chaotic behaviors can be characterized by probability densities of the gaps $P_{\rm GOE}(S)=(\pi S/2)\exp{(-S^{2}\pi/4)}$ and $P_{\rm P}(S)=e^{-S}$, respectively, for unity mean spacing, when dealing with symmetric matrices Porter (1965). Often what is done in studying the many-body localization (MBL) transition Nandkishore and Huse (2015); Altman (2018); Abanin et al. (2019); Alet and Laflorencie (2018) for physical systems is to combine both types of matrices, gradually changing their relative weight until the statistics of the spectrum suddenly changes ¹¹1There have been discussions of whether this is ever possible, with some arguing that ergodic behavior survives at any finite value of disorder, i.e., at any finite weight of the Poisson matrix Suntajs et al. (2019), while others claim that finite size effects are yet large to draw definitive conclusions Abanin et al. (2019); K. Panda et al. (2019); Sierant and Zakrzewski (2019).. If increasing the weight of the matrix displaying Poisson statistics, this signals the onset of non-ergodic behavior once the eigenvalues of the combined matrix become uncorrelated, and their eigenfunctions display support not scaling with the matrix dimension L. D’Alessio et al. (2016). When purely dealing with random matrices, this type of combined ensemble mixing matrices of different symmetry classes has been studied in the past Lenz and Haake (1990, 1991); Lenz et al. (1991); Kota and Sumedha (1999); Schierenberg et al. (2012); Chavda et al. (2014), with corresponding surmises for the level spacing distributions being specifically defined in terms of the weights, and interpolating between $P_{\rm GOE}(S)$ and $P_{\rm P}(S)$, if, e.g., adding random matrices from the Gaussian orthogonal ensemble (GOE) and diagonal random matrices.

In this Letter, we employ the mixed ensemble model to study the tuning process from GOE to Poisson statistics, unlike others with intermediate ensemble models Serbyn and Moore (2016); Buijsman et al. (2019); Sierant et al. (2019), providing strong evidence on the existence of the mobility edges and the stability of the coexistence of localized and delocalized states in approaching the thermodynamic limit.

We construct a combined random matrix ${\sf M}$, which reads

where ${\sf M}_{e}$ is a GOE random matrix constructed by means of a normal random distribution matrix ${\sf A}$ by $({\sf A}+{\sf A}^{T})/2$, ${\sf M}_{i}$ is a diagonal matrix with random elements, and $k\in[0,1]$ is a real coefficient. The limiting values $k=0$ (${\sf M}={\sf M}_{e}$) and $1$ (${\sf M}={\sf M}_{i}$) characterize chaotic and non-chaotic regimes, where the probability distributions of gaps in the spectrum follow $P_{\rm GOE}$ and $P_{\rm P}$, respectively. Our main interest is focused on the intermediate case $0<k<1$, where ${\sf M}$ can potentially describe the transition from chaos to localization. In this regime of mixed ensembles, surmises for the probability distribution have been derived Lenz and Haake (1991); Lenz et al. (1991); Schierenberg et al. (2012), in direct similarity with the Wigner surmises of pure ensembles Haake (2006). In what follows, we investigate two main indicators of chaotic matrices, as to classify both the eigenvalues spectrum as well as their eigenfunctions. For the first, we notice that the gap distributions albeit precise in characterizing the ergodic properties of the spectrum, suffer from technical difficulties related to the necessity of having mean level spacing $\langle S\rangle$ constrained to unity for the aforementioned distributions to be valid. Thus instead, we look at the distribution $P(r)$ of the ratio of two consecutive gaps in the spectrum $r_{\alpha}=\min\{S_{\alpha},S_{\alpha+1}\}/\max\{S_{\alpha},S_{\alpha+1}\}$, requiring no unfolding schemes  Atas et al. (2013).

In turn, a direct measure of the (de)localization of the eigenfunctions is obtained by the inverse participation ratio (IPR), ${\cal I}^{(\alpha)}=\sum_{n}|\psi_{n}^{(\alpha)}|^{4}$, where $\psi_{n}^{(\alpha)}$ is the $\alpha$-th eigenstate of the matrix and $n$ is the basis state index. The eigenstate properties are characterized according to its limiting values as

where $D$ is the matrix dimension. The IPR of a chaotic (extended) eigenstate is $D$-dependent, in contrast to a localized one; therefore, a finite-size scaling is necessary to differentiate the two types of states when approaching the limit $D\rightarrow\infty$. In that regime, it is important to extract a scaling parameter that can account for the modifications the mixed ensemble suffers. For example, the spectral variance of ${\sf M}_{e}$ is approximately proportional to $D$ for large enough matrix dimension Chavda et al. (2014). For that, we conjecture the following scaling parameter ²²2See Chavda et al. (2014) for a similar conjecture inspired by the preliminary derivation of probability distribution of eigenvalues of the mixed ensembles.,

where $\tilde{\sigma}$ and $\sigma$ are the standard deviation of the diagonal elements in ${\sf M}_{i}$ and ${\sf M}_{e}$, respectively. When $k$ changes from $0$ to $1$, $\xi$ changes from $0$ to infinity. Without loss of generality, we set $\tilde{\sigma}=\sigma=1$ and choose the standard deviation of non-diagonal elements $\sqrt{2}/2$ times of the diagonal elements in ${\sf M}_{e}$. We find that, by changing $k$ and $D$, but fixing $\xi$, the system keeps the same local and chaotic properties.

We test this scaling form by employing specific distributions for the diagonal elements in ${\sf M}_{i}$, Gaussian, uniform and cosine distributions in Fig. 1 [(a),(d)], 1 [(b),(e)] and 1 [(c),(f)], respectively. As shown in Fig. 1(d),  1(e),  1(f), the average inverse participation ratios over the whole spectrum, $\langle\cal I\rangle$, are completely coincident for $D=500$, $1000$ and $2000$ when using the scaling parameter $\xi$. This is indicative that for large matrices, the ${\cal I}^{(\alpha)}$ diagram for different distributions Fig. 1(a),  1(b),  1(c) (for $D=2000$) correctly reflect the phase diagrams for $D\rightarrow\infty$ (See Fig. 2 for details of IPR before the average for the different $D$s). Interestingly, the chaotic region shows a “D”-like shape in Fig. 1(a), and a “C”-like shape in Fig. 1(c). Since out of the closure of “C” and “D” it is a localized phase, both situations indicate a typical phase diagram of mobility edges. However, there is no obvious mobility edges for the uniformly distributed ${\sf M}_{i}$ [Fig. 1(b)]. Thus, those are model-dependent and do exist for ${\sf M}_{i}$ of Gaussian distribution or of cosine distribution for fixed $\xi$ even when the dimension of the matrix goes to infinity. Moreover, the transition from classical chaos to localization with increasing $\xi$ has characteristics of a crossover.

We further study the IPR and $\langle r\rangle$ with fixed $\xi$ but different $D$s in Fig. 2. The latter, apart from finite size fluctuations, is coincident for increasing dimensions $D$. For matrices ${\sf M}_{i}$ from the Gaussian distribution, the states at the edges of the spectrum ($\alpha\lesssim 0.05$ and $\alpha\gtrsim 0.95$, with $\xi=0.1$) do not follow GOE statistics, indicating that these are not chaotic. This is supported by the corresponding values of the IPR with increasing matrix sizes in Fig. 2(a): Only at the central part of the spectrum scaling with $D$ occurs. Conversely, the situation is opposite for ${\sf M}_{i}$ constructed from a cosine function distribution [See Figs. 2(b) and 2(d)]. At the edges of the spectrum, ${\cal I}^{(\alpha)}$ decreases with increasing dimensions and is accompanied by a mean ratio of adjacent gaps characteristic of chaotic matrics from the GOE ensemble for the value of $\xi=3$ used. Based on these results for mixed ensembles of random matrices, we apply similar formalism to Hamiltonian matrices of typical quantum systems.

In order to test our analysis and to determine whether a many-body mobility edge can exist in a one-dimensional quantum system on a lattice, we construct a many-body model, composed of $N$ spinless fermions in a chain with $L$ sites. The Hamiltonian is given by

where $\hat{T}$ is the (non-random) nearest-neighbor hopping matrix,

assuming $\hat{c}_{L+1}=\hat{c}_{1}$. The spectral variance of $\hat{T}$ is proportional to $N$ for large matrix sizes $D={L\choose N}$ and fixed filling $n=N/L$, as $\sigma_{E}^{2}=2N(1-n)$, which contrasts with a typical GOE matrix where $\sigma_{E}^{2}\propto D$. In turn, the second matrix is diagonal, and characterizes both the interaction and disorder, which is typically the case for Fock bases in real space. It assumes a form $\hat{V}=V{\sf G}$, with $V$ its strength, and ${\sf G}$ populated with random numbers selected from a particular distribution. Motivated by the cases where mobility edges were found when dealing with pure random matrices, we focus on two types: Gaussian and cosine distributions, both characterized by a standard deviation $\tilde{\sigma}$ of its elements. Other distributions, as power-law, can also result in mobility edges and are analyzed in the Supplemental Material SM . We notice that short-range homogeneous interactions result in diagonal matrices whose distribution of elements is given by a discrete Gaussian, and in the case the interactions are fully random, they approach a continuous distribution that can also lead to the onset of MBL Sierant et al. (2017); Li et al. (2017b). Furthermore, a typical Anderson-like disorder with uniform random local energies also result in a Gaussian distribution of entries in the diagonal elements for a many-body basis.

The scaling parameter of Eq. (3) in this model becomes,

which we use, as before, to scale the average values of IPR for increasing matrix sizes. Hereafter, the filling number is set to $n=N/L=0.5$, which results in the largest possible matrix sizes for a given $L$.

In Fig. 3, we calculate the phase diagrams of the IPR for $L=18$ and its average value, $\langle{\cal I}\rangle$, as a function of $\xi$ for different $L$s. The phase diagrams in Fig. 3(a),  3(b) are similar to those of Fig. 1(a),  1(c) for the pure random matrices, showing robust many-body mobility edges for Gaussian and cosine function distributions of diagonal matrix elements, respectively. In Fig. 3(c), $\langle{\cal I}\rangle$ shows converged values for different system (matrix) sizes for $\xi\gtrsim 7$ [details shown in Fig. 4(a)], indicating Eq. (6) serves as a scaling parameter, and that the system size $L$ here is no longer an important factor to trigger MBL. For $\xi<7$, $\langle{\cal I}\rangle$ decreases as $L$ increases, a signal that ergodic states start to populate the spectrum. Figure 3 (d) for the cosine distribution shows very similar behavior, but with a critical $\xi\simeq 16$ separating the regime with the presence of ergodicity and full MBL.

In Fig. 4, we investigate the system size effects on ${\cal I}^{(\alpha)}$ and $\langle r\rangle$ across the spectrum $\alpha$, for different typical values of $\xi$ in both distributions. For the many-body localized region [typical value $\xi=7$ (16) for Gaussian (cosine) distributions], Fig. 4 (a) and (b) show that completely coincident IPR for different sizes confirm localization across the whole spectrum. For smaller weights $\xi$ of the diagonal matrices, ergodic and delocalized states naturally occur in the central part (edges) of the spectrum for ${\sf G}$ of the Gaussian (cosine function) distribution. These regions are identified by an $L$-dependence in both ${\cal I}^{(\alpha)}$ and $\langle r\rangle$, with the former decreasing with the system size and the latter heading towards $\langle r\rangle_{\rm GOE}$ for increasing $L$s. When these quantities display system size independence, lack of ergodicity and localization ensues. We notice, however, that the onset of this behavior is not coincident in $\alpha$, i.e., it suggests that the many-body mobility edge occurs not as a sharp transition in the spectrum but rather as one characterized by an intermediate phase with non-ergodic but delocalized wavefunctions. This regime has been found in physical models and dubbed non-ergodic metal Li et al. (2015, 2016); Deng et al. (2017); Hsu et al. (2018).

Mixed ensembles provide further fundamental insights on the onset of ergodicity. Starting from a Poisson diagonal matrix, any non-zero weight of an added GOE matrix triggers a finite degree of level repulsion Lenz and Haake (1991); Schierenberg et al. (2012), a defining characteristic of ergodicity. However, the chaotic effects induced by the latter crucially depend on the density of states. That is, the effective weight of the GOE matrix on the resulting mixed matrix grows as denser the spectrum is Schierenberg et al. (2012), initially investigated via the analysis of the density-of-states-resolved statistics of the gaps $P(\rho,S)$, which requires unfolding of the spectrum. An $r$-statistic analysis avoids this inconvenience, but a surmise for $P(r)$ in the mixed GOE-Poisson ensemble is yet elusive Chavda et al. (2014). Nonetheless, we can infer similar information regarding the effect of the GOE random matrix by a working numerical definition. Starting from random $3\times 3$ matrices of the mixed ensemble form ${\sf M}={\sf M}_{i}+\lambda{\sf M}_{e}$, one can obtain the numerical surmise $P(r,\lambda)$ by using a set of 2,000 of such matrices for each $\lambda$. Subsequently, we extract the adjacent gap distribution for a much larger combined $D\times D$ matrix with a typical weight $\lambda=\Lambda$. By dividing the spectrum in a large number of homogeneous energy windows, we resolve the distribution by the density of states $\rho$, obtaining $P(r,\rho)$. Finally, the effect of the weight $\lambda$ of the chaotic matrix can be quantified as,

which is shown in Fig. 5, as a color plot. The minima of this work definition (white markers) assures that the effect of the $\lambda$ weight of the GOE matrix is linearly proportional to the density of states, a result also valid for the distributions of the gaps $P(\rho,S)$ Schierenberg et al. (2012). We now apply the same analysis in the case of a physical system [Eq. (4)]. Now, the different density of states results in a remarkable contrasting behavior for the effect of the thermalizing matrix. $\lambda$ approaches zero at finite density of states, which in practice results in the absence of ergodicity at the edges of the spectrum, that are typically less dense. This thus indicates that a sharp transition occurs in the ergodic properties of the eigenvalues in the spectrum, where a region with small density of states does not suffer from thermalizing effects and can be associated to the manifestation of a mobility edge.

We study the process from GOE to Poisson statistics by constructing a combined random matrix as to shed light on the long lasting debate on the existence of mobility edges. We find phase diagrams in the thermodynamic, determined by the invariant scaling parameter $\xi$, suggesting its existence for different distributions of the diagonal random matrix. Lastly, by employing a numerical analysis on the effect of the strength of the GOE matrix in different parts of the spectrum, we infer that less dense regions have an effective vanishing weight of the ergodic matrix, albeit denser ones do render ergodic behavior. Thus, we argue that mobility edges can exist in physical systems with short-ranged hoppings, even when approaching the thermodynamic limit. A closed form of the probability distribution of the ratio of adjacent gaps of mixed ensembles may result in an even more systematic way to define whether different quantum systems may or not display mobility edges and will be reserved for a future work.

In the main text, we have analyzed how the distribution of the random diagonal elements affect the overall shape of the localization-delocalization transition in the spectrum-disorder amplitude ($\alpha$ vs. $\xi$) space. In both cases where a mobility edge was observed, Gaussian and cosine function distribution, it is manifested mainly at the edges of the spectrum, where the density of states is small. Finite size effects are thus more dramatic at those regimes, as increasing $L$’s can significantly alter the sparsity of the levels at the ends of the spectrum. Since the ergodic-MBL transition shape significantly depends on the random distribution used to construct $\sf G$, we can proceed with the exercise of finding another distribution whose manifestation of a mobility edge avoids this problem. A natural choice, beyond the ones we have chosen, is a power-law function, which we select as $F(x)=(x^{4}-1/5)/(4/15)$. Here $x$ is a random number, picked in the interval $[-1,1]$, and $1/5$, and $4/15$ are the mean value and standard deviation of $F(x)=x^{4}$, respectively.

In Fig. S1(a), we report the phase diagram for the physical model (4), with $L=18$, and elements in $G$ selected according to the probability distribution $F(x)$ above. There is a clear many-body mobility edge with only one ergodic-to-MBL transition, and the critical point is around the middle of the spectrum for $\xi=9$. Fig. S1(b) indicates that the average IPR, $\langle{\cal I}\rangle$, possesses a small $L$-dependence in a large range of $\xi$, in similarity to the other distributions studied in Fig. 3. However, the spectrum resolved IPR, ${\cal I}^{(\alpha)}$, displays, on the other hand, a much more sharp transition in the spectrum approaching zero much faster than for other distributions studied in the main text. This sharp transition on the localization properties of the eigenfunctions is accompanied by a well marked transition on the level repulsion [Fig. S1(d)], displaying a GOE-to-Poisson transition for the statistics of the level spacings at similar values of $\alpha$.

In Ref. Chavda et al., 2014, Chavda et al. used a mixed ensemble of the form

where $\sf H_{0}$ is the Poisson matrix and $\sf V$ the corresponding GOE one. Similarly to what we did, they selected the variance of the off-diagonal elements of the GOE matrix being one-half of the variance of the diagonal ones, $\sigma^{2}$. After that, by arriving on a non-closed form of the joint probability of the eigenvalues of a 3$\times$3 ensemble, they notice that it depends on the transition parameter,

where $\overline{d_{0}}$ is the mean level spacing between nearest eigenvalues of the diagonal (unperturbed) Poisson matrix. To match our definition of the mixed ensemble, one can easily identify that $\lambda=\frac{1-k}{k}$ and that the mean level spacing of the Poisson matrix $\overline{d_{0}}=\frac{\tilde{\sigma}}{\sqrt{2}D}$. With this, the scaling parameter we used and theirs are related as,

As us, they identify the very small dependence of the average ratio of adjacent gaps over the whole spectrum on the matrix size $D$, when scaling with the transition parameter $\Lambda$, and we do so by looking at the average value of the IPR, scaling it with $\xi$ (see Fig. 1).

Figure 5 in the main text shows that for a mixed ensemble, the effective weight of the ergodic matrix on the eigenvalues of the mixed ensemble crucially depends on the density of the states of the latter. That is, for a denser region of the spectrum, the ergodicity sets in more strongly and is finite almost as any finite density is obtained for the eigenvalue spectrum. For the case of the physical system, where the dense GOE-matrix is substituted by a sparse non-random (but symmetric matrix), there is a critical density above which the ergodic aspects become relevant. We now analyze the matrix size dependence on these result in Fig. S2. For the case of the pure mixed ensemble, Fig. S2(a) shows that the linear relationship between the effective weight and the (normalized) density of states is valid for a wide range of matrix sizes $D$. Similarly, when using the physical model $\hat{H}=\hat{T}+\hat{V}$ [Fig. S2(b)], the finite size effects do not qualitatively alter the main aspect expressed in the main text, in which a critical level density is necessary for the effective ergodic weight to become relevant. It is worth mentioning that the critical density is similar in all system sizes analyzed.

Finally, to see that this threshold density of eigenlevels is not vanishingly small, we report in Fig. S3 the normalized density of states of both cases, marking with a horizontal dashed line the maximum density where the effective weight $\lambda$ goes to zero. For the physical model, it is clear that a substantial part of the spectrum is thus not affected by the ergodic properties, and one can consequently identify a clear separation in behavior in different parts of it, which we interpret as the manifestation of a mobility edge.