Eigenstate properties of the disordered Bose−Hubbard chain

Abstract

Many-body localization (MBL) of a disordered interacting boson system in one dimension is studied numerically at the filling faction one-half. The von Neumann entanglement entropy $S_{\mathrm{v}\mathrm{N}}$ is commonly used to detect the MBL phase transition but remains challenging to be directly measured. Based on the $U(1)$ symmetry from the particle number conservation, $S_{\mathrm{v}\mathrm{N}}$ can be decomposed into the particle number entropy $S_N$ and the configuration entropy $S_C$. In light of the tendency that the eigenstate’s $S_C$ nears zero in the localized phase, we introduce a quantity describing the deviation of $S_N$ from the ideal thermalization distribution; finite-size scaling analysis illustrates that it shares the same phase transition point with $S_{\mathrm{v}\mathrm{N}}$ but displays the better critical exponents. This observation hints that the phase transition to MBL might largely be determined by $S_N$ and its fluctuations. Notably, the recent experiments [A. Lukin, et al., Science 364, 256 (2019); J. Léonard, et al., Nat. Phys. 19, 481 (2023)] demonstrated that this deviation can potentially be measured through the $S_N$ measurement. Furthermore, our investigations reveal that the thermalized states primarily occupy the low-energy section of the spectrum, as indicated by measures of localization length, gap ratio, and energy density distribution. This low-energy spectrum of the Bose model closely resembles the entire spectrum of the Fermi (or spin $XXZ$) model, accommodating a transition from the thermalized to the localized states. While, owing to the bosonic statistics, the high-energy spectrum of the model allows the formation of distinct clusters of bosons in the random potential background. We analyze the resulting eigenstate properties and briefly summarize the associated dynamics. To distinguish between the phase regions at the low and high energies, a probing quantity based on the structure of $S_{\mathrm{v}\mathrm{N}}$ is also devised. Our work highlights the importance of symmetry combined with entanglement in the study of MBL. In this regard, for the disordered Heisenberg $XXZ$ chain, the recent pure eigenvalue analyses in [J. Šuntajs, et al., Phys. Rev. E 102, 062144 (2020)] would appear inadequate, while methods used in [A. Morningstar, et al., Phys. Rev. B 105, 174205 (2022)] that spoil the $U(1)$ symmetry could be misleading.

Graphical abstract

Keywords

entanglement entropy decomposition / U(1) symmetry / thermalization-to-localization transition

1 Introduction

In the absence of many-body interactions, it is known that Anderson localization may occur in a disordered system depending upon dimensionality and disorder strength [1–3]. This localization phenomenon, essentially caused by interference of elastic scatterings, has been intensively investigated in various disordered circumstances for more than five decades [4–6]. In recent years, with the advances of numerical techniques, computational capacities, and experimental facilities, the interaction effect within many-body disordered systems has become of great interest in light of the novel competition between many-body interactions and quenched disorder. One would like to thoroughly explore the effects of disorder on the fundamental properties of many-body interacting systems associated with ergodicity breaking, eigenstate thermalization, and quantum criticality. Particularly, many-body localization (MBL) is expected to be intrinsically different from Anderson localization [7–10].

In the absence of disorder, a nonintegrable interacting system is usually considered to satisfy the so-called eigenstate thermalization hypothesis (ETH) [11, 12], which assumes that the ergodicity of many-body states exists, implying that eigenstate fluctuations in the expectation of a quantum mechanical observable is about the same order of magnitude as the statistical fluctuations in a microcanonical ensemble controlled by an appropriate energy density [7]. When sufficiently strong disorder is introduced into such many-body interacting systems, an ergodic many-body state may turn into an MBL state, in which local integrals of motion (LIOMs) emerge [13–15] such that the integrability is considerably restored, resulting in the Poisson distribution of the energy levels, rather than the Wigner−Dyson surmise for the typical ETH phases as is derived from the random matrix theory [7, 16].

Phenomenology of MBL can be relevant for many different but related settings, such as a Floquet disordered system [17–20], for which the driven system is considered to implement the so-called discrete time crystal [21, 22], the quantum frustrated Heisenberg spin chains [23], and the one-dimensional (1D) QED lattice version of the Schwinger model with a disordered gauge [24], as well as in several open quantum setups [25–30]. Physical properties of MBL have been widely studied numerically for 1D Heisenberg or similarly interacting hard-core boson systems with randomness [23, 31–36]. One finds that by exploiting entanglement entropy, a dynamical transition beyond the quantum phase transition occurs between the ETH and MBL states that is driven by disorder strength when varying the energy density [23, 31–34]. While an ETH state has an extended thermal nature, an MBL state is featured by its insulating nature, being close to a product of localized states. Correspondingly, entanglement entropy possesses a volume law in the ETH phase [23] but an area law in the MBL phase. Moreover, a slowly dephasing process is unveiled in MBL phase, engendering a logarithmical entropy increase upon the time evolution of a many-body state [37, 38].

Many experimental efforts are cast into the measurements of disorder effects for interacting many-body systems [21, 22, 36, 39–42]. Specifically, Lukin et al. [39] succeeded in realizing a 1D interacting Aubry−André model for bosons. By probing particle number entropy and configuration entropy, they detected that the signature of quantum thermalization, the particle’s finite localization length, and an area-law scaling of the particle number entropy can help to distinguish the MBL phase from the ergodic ETH phase. Moreover, they also found a slowly growing behavior of the configuration entropy that ultimately results in a volume-law scaling, demonstrating that the MBL state is qualitatively different from a noninteracting Anderson insulating state. Curiously, the key difference between the Bose model and the Fermi (or spin) model is that Bose statistics can result in the clustering of multiple particles, potentially giving rise to the novel phenomena of Bose type MBL. This distinctive character is missing in Refs. [39, 43].

In this paper, we extensively explore the nature of MBL phase in 1D Bose−Hubbard model with disorder by exact diagonalization (ED) [44] and several complementary techniques [45, 46]. For the present system, numerical calculations are limited to small sizes since the occupancy number of bosons can be any integer, depending on the filling factor under consideration. However, since Hamiltonian (1) below commutes with total number of particles, the ED calculation can be performed in the subspace with a given number of particles. For a chain with $L$ sites and $N$ bosons, the dimension for the matrix diagonalization is $(L+N-1)!/[N!(L-1)!]$ [44]. By implementing the separation of von Neumann entanglement entropy, one can further compute the particle number entropy and configuration entropy to access the so-called symmetry-resolved entropy. In our numerical calculations, we exclusively focus on the filling factor $1/2$. Considering the rapid growth of the Hilbert-space dimension, it is practically feasible that $10000$ samples were taken for the averages with $L=8,5000$ with $L=10,3000$ with $L=12$, and for $L=14,200$ samples are used for weak disorder, $400$ for medium disorder, and $600$ for strong disorder. We impose open boundary conditions in the calculations to compute the two-body density-density correlation function and periodic boundary conditions for other quantities. To access larger system sizes, we also employ time-evolving block decimation (TEBD) [45, 46] method to evaluate the time evolution of the two-body correlation function. Errors resulting from sample averages are given explicitly for most cases and are invisible otherwise.

In the following, we first introduce the disordered Bose−Hubbard (dBH) model and briefly discuss current status regarding the ground-state studies on this model in Section 2. In Section 3, we focus on the phase transition from ETH to MBL; at the beginning of this section, a detailed outline of its content is provided. In Section 4, we analyze different MBL phases within dBH model through the properties of its eigenstates; this section also starts with a detailed content overview. A summary and pertinent discussions are finally given in Section 5. For completeness, several corresponding supplementary results are included in Appendixes A−D.

2 Model

The 1D dBH model can be written as

(1)$\begin{array}{rl}\hat{H}=& -J\sum _i({\hat{a}}_i^{†}{\hat{a}}_{i+1}+\mathrm{H}.\mathrm{c}.)\\ & +\frac{U}{2}\sum _i{\hat{n}}_i({\hat{n}}_i-1)+\mu \sum _i{\delta }_i{\hat{n}}_i,\end{array}$

where ${\hat{a}}_i^{†}$ and ${\hat{a}}_i$ are particle creation and annihilation operators, respectively, and ${\hat{n}}_i={\hat{a}}_i^{†}{\hat{a}}_i$ is the particle occupation number operator at site $i$. Disorder is introduced through a random number ${\delta }_i$, which is uniformly distributed within $[-1,1]$ for all sites, thus chemical potential $\mu$ becomes a tuning parameter for the disorder strength. For convenience, we set $J=1$ as the energy-scale unit hereafter.

For Hamiltonian (1), we use the ED technique to find all its eigenvalues $E_i$ with an order of $E_1⩽E_2⩽\cdots ⩽E_D$ and $D$ is the dimension of the Hamiltonian matrix. The energy density $\epsilon$ used throughout this work is then defined by

(2)$\epsilon ={\displaystyle \frac{E_i-E_1}{E_D-E_1}}.$

At zero temperature, the BH model, in the absence of disorder, undergoes a quantum phase transition from a superfluid phase to the Mott-insulator state with an integer filling factor as the interaction strength $U$ increases [47]. When disorder is introduced, the nature of the ground state is anticipated to alter [48–51]. For a given diagonal disorder sample $\{{\delta }_i\}$ with $i\in [1,L]$, one finds that the system changes from the superfluid state first to the Bose-glass state and then to the Mott-insulator state under the rise of $U$ [49, 50]. In the case of an off-diagonal disorder configuration assigned to $J$, a phase diagram has also been established in terms of the interaction strength $U$ and the randomness strength in the hopping amplitude $J$ [49–51].

In this work, we fix $U=3J$ and always maintain the filling factor to be one-half per site, i.e., $N=\sum _i{\hat{n}}_i=\frac{L}{2}$, for which one has a superfluid ground state in the case of a clean chain [52].

3 Eigenstate phase transition

A plethora of physical quantities have been used to numerically detect the MBL phase transition, but most are experimentally inaccessible. Among these, the von Neumann entropy $S_{\mathrm{v}\mathrm{N}}$ is a common choice, yet remains challenging to be directly measured. Conceiving the experimentally accessible probes to assess the MBL transition therefore comprises one of the main goals of the present work.

In Section 3.1, based on the $U(1)$ symmetry, we decompose $S_{\mathrm{v}\mathrm{N}}$ into the particle number entropy $S_N$ (readily measurable) and the configuration entropy $S_C$ (not currently measurable). We first observe that the $S_C$ of the eigenstate nears zero in the localized phase, implying that $S_C$ may largely be phase transition irrelevant (or at least less relevant). Motivated by this, in Section 3.2, we introduce a new quantity that quantifies the deviation of $S_N$ from the ideal thermal distribution. Subsequent finite-size analysis shows that it shares the same phase transition point with $S_{\mathrm{v}\mathrm{N}}$ but bears a greater critical exponent, indicating that it is $S_N$ that primarily drives the phase transition. Most importantly, this deviation of $S_N$ can be preliminarily measured in the laboratory, as has already been demonstrated by recent Harvard experiments [39, 43] (especially in boson systems). Admittedly, here we have also asserted the assumption that the energy-resolved eigenstate measurements we consider can be qualitatively executed through the careful combinations of disorder realizations and initial state preparations for using the dynamics, as was successfully implemented in the recent programmable-superconducting-processor-based experiment on MBL [41]. Recent debates about whether the MBL transition is a true phase transition or a crossover [53, 54] prompt us to seriously inspect such experimental measurability, and these considerations can contribute to the settling of these ongoing controversies (see comments at the end of this subsection). Finally, in Section 3.3, to gain an overview of the eigenstate-transition scenario, we sketch the finite-system dynamical phase diagram obtained from the large-scale numerical simulations.

3.1 Entanglement entropy and particle number fluctuation

In Appendix A, we detail the decomposition of the total entanglement entropy $S_{\mathrm{v}\mathrm{N}}$ into the particle number entropy $S_N$ and the configuration entropy $S_C$ when the total number of particles $N$ in the system is conserved and discuss how to calculate the configuration entropy in an efficient way.

Now we mainly focus on the obtained numerical results for these two quantities, the particle number and configuration entropies. Fig.1(a) shows the size dependence of $S_N$. For each length $L$, $S_N$ decreases monotonically with the increasing disorder strength $\mu$. Likewise, for $S_C$, we can see from Fig.1(b) that it lessens as disorder enhances, and approaches zero once the disorder strength $\mu$ becomes significant. This is because if disorder is strong enough, the system is expected to be localized in several particular configurations. The configuration entropy due to the superposition of these configurations thereby tends to near zero. Under the circumstance of larger $\mu$, the particle number entropy $S_N$ appears to dominate the evolution of the total entropy. Besides, both entropies exhibit a volume law at weak disorder and an area law at strong disorder, consistent with the general trend of the total entropy found in previous MBL studies [23, 33].

Fig.1 (a) Particle number entropy $S_N$ versus disorder intensity $\mu$ for different sizes $L$ at an energy density $\epsilon =0.25$. (b) Results for configuration entropy $S_C$ at the same energy density. Companion results for $p_n$ and $S_{\mathrm{v}\mathrm{N}}^n$ with respect to differing sizes are shown in the insets of (a) and (b) for $\mu =2$.

Full size|PPT slide

As pointed out in Appendix A, $S_N$ and $S_C$ can be expressed in terms of $p_n$ and $S_{\mathrm{v}\mathrm{N}}^n$. Here, $n$ denotes the number of particles in the left half chain, marking the chunks of the reduced density matrix which we call channels. The value of $n$ ranges from $0$ to $N$ (there are a total of $N+1$ options for $n$). In the inset of Fig.1(a), we give the distribution of $p_n$ as a function of $n-N/2$ with respect to different chain lengths for a small disorder strength $\mu =2$. It can be seen that since $p_n$ needs to be normalized and there is no difference between tracing out the left or the right portion of the system at ETH, it is a symmetric distribution of the Gaussian lineshape. In the inset of Fig.1(b), the results of $S_{\mathrm{v}\mathrm{N}}^n$ are given, again symmetric. Overall, as the chain length increases, the dimension of the channel grows, so does the corresponding entropy.

$S_N$ arises from the superposition of states with varying particle numbers in half-chain subsystems, generated by particle motion across the boundary of the bipartite ($A$ and $B$). This fluctuation and its measurability can both be understood in terms of the particle number variances for part $A$, which is typically defined as follows:

(3)$\begin{array}{rl}{F}_N=& {\left\{⟨{\hat{N}}_A^2⟩-⟨{\hat{N}}_A{⟩}^2\right\}}_{\mathrm{a}\mathrm{v}}\\ =& -{\left\{\sum _{i<\frac{L}{2}}\sum _{j>\frac{L}{2}}\left[⟨{\hat{n}}_i{\hat{n}}_j⟩-⟨{\hat{n}}_i⟩⟨{\hat{n}}_j⟩\right]\right\}}_{\mathrm{a}\mathrm{v}},\end{array}$

where ${\hat{N}}_A$ is the total particle number operator for part $A$. $F_N$ accounts for the correlation between the two-particle densities on the two sites located inside portions $A$ and $B$, respectively, reflecting the particle-number-fluctuation induced entanglement between the bipartite.

Fig.2 shows the results of $F_N$ at $\epsilon =0.25$. One can see that it indeed behaves analogously to $S_N$. This implies that it also has a volume law for small $\mu$ in the ETH phase and tends toward the area law at large $\mu$ for the MBL phase. Analyzing the last expression of (3) when subject to the sufficiently large $\mu$, one may expect that the correlations would become exponentially small with increasing $|i-j|$. It turns out that in the summation of (3), those terms with small $|i-j|$ do dominate, i.e., only the contributions of those correlations at the interface between parts $A$ and $B$ remain finite for the MBL phase, thereby leading to an area law for $F_N$. This also reinforces the experimental observation on the area-law scaling of $S_N$ in the MBL phase as was reported in the Harvard experiment [39].

Fig.2 Half-chain particle number fluctuation $F_N$ versus $\mu$ at an energy density $\epsilon =0.25$ for various lengths $L$.

Full size|PPT slide

3.2 Number entropy dominates the MBL transition

As seen in the previous subsection, in the localized phase, the $S_C$ of the eigenstate nears zero and consequently $S_N$ dominates the total entropy. This implies that $S_C$ may be much less relevant to the phase transition and the transition is chiefly governed by $S_N$. Next, we will show numerical evidence to substantiate this claim.

Firstly, we calculate the total von Neumann entropy $S_{\mathrm{v}\mathrm{N}}$ for a given $\epsilon =0.25$ using the various sizes $L$, which is expected to capture the finer signal of the transition point according to the general scaling theory of the critical phenomena [55, 56]. As per the finite-size scaling theory, the phase transition point and the critical exponent thereof can be obtained when all data curves are collapsed into a single line describable by the following scaling form ansatz [23, 57],

(4)$Q(L,\mu )=g(L)f\mathbf{(}(\mu -{\mu }_c)L^{1/\nu }\mathbf{)},$

which gives rise to a collapsed scaling relation of $y_L=Q(L,\mu )/g(L)$ versus $x_L=(\mu -{\mu }_c)L^{1/\nu }$ and the critical value ${\mu }_c$ as well as the critical exponent $\nu$ for the transition. $g(L)$ is generally a polynomial of $L$. Here, we set $g(L)=L$. In Fig.3(a), we show a size dependence of $S_{\mathrm{v}\mathrm{N}}$ at $\epsilon =0.25$ for $L=6,8,10,12,14$. Using the scaling relation mentioned above, we plot $y_L=S_{\mathrm{v}\mathrm{N}}/L$ in the inset of Fig.3(a), where we can collapse the data into almost a single curve except for those points at the very low $\mu$. We are thus able to obtain a critical value ${\mu }_c\approx 7.8$ and the associated critical exponent $\nu \approx 1.2$, namely a phase transition point at $\epsilon =0.25$.

Fig.3 (a) Variation of the total entropy $S_{\mathrm{v}\mathrm{N}}$ with the disorder intensity for different sizes; the inset gives the result of its data collapse, yielding the estimate of the phase transition point at ${\mu }_c\approx 7.8$. (b) is the corresponding result of $|d(S_N)|$ for different chain lengths when disorder increases; its data collapse is given also in the inset, showing the same phase transition point at around ${\mu }_c\approx 7.8$. The energy density here is always fixed to be $\epsilon =0.25$.

Full size|PPT slide

Secondly, for an ideally thermalized system where the system has an equal probability to access various configurations, a classical probability can be obtained by dividing the dimension of each symmetry-channel (or symmetry-sector) block by the total dimension of the Hilbert space. For the present model, the total Hilbert space of a bosonic system with the chain length $L$ and the number of particles $N$ is $D_N^L=\frac{(L+N-1)!}{N!(L-1)!}$, and the dimension of the Hilbert subspace with particle number $n$ on the left half chain is $D_n=D_n^{L/2}{D}_{N-n}^{L/2}$, so $p_n^{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}=D_n/D_N^L$, and we can get $S_N^{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}=-\sum _{n=0}^N{p}_n^{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}{\log }_2{p}_n^{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}$. To characterize the extent to which the actual probability distribution deviates from this ideal thermalization situation, we next define the following quantity so as to potentially locate the phase transition more accurately than the bare $S_N$,

(5)$|d(S_N)|=|\{S_N{\}}_{\mathrm{a}\mathrm{v}}-S_N^{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}|.$

As shown in Fig.3(b), we take $g(L)=L^a$ for the case at hand, and one can see that the phase transition point thus obtained from $|d(S_N)|$ agrees with the phase transition point obtained from the total entropy $S_{\mathrm{v}\mathrm{N}}$ shown in Fig.3(a).

More encouragingly, the critical exponent obtained from $|d(S_N)|$, $\nu \approx 1.5$, is closer to satisfy the Harris bound $(\nu ⩾2)$ [58] compared to the critical exponent $\nu \approx 1.2$ obtained from $S_{\mathrm{v}\mathrm{N}}$. Previous numeric evaluations of various physical quantities have yielded critical exponents all around $\nu \approx 1.0$ [59], similar to that obtained from $S_{\mathrm{v}\mathrm{N}}$. This suggests that the inclusion of $S_C$ may weaken the critical phenomena in general and enhance the finite-size effects in particular. This could partially lead to the pronounced drifts in the finite-size calculations, making it challenging to distinguish between a phase transition and a crossover. Furthermore, it is noteworthy that we define the deviation $|d(S_N)|$ using $S_N$, and $S_N$ itself is computed by decomposing the total entropy as per the $U(1)$ symmetry of the system. It is this symmetry-based decomposition that allows us to separate out the phase-transition-less-relevant quantity $S_C$. This indicates that symmetry plays a crucial role in the MBL phase transition. Finally, the deviation $|d(S_N)|$ we defined can be accessed directly through the $S_N$ measurements, and recent experiments (especially in boson systems [39, 43]) have shown that $S_N$ might be available.

Currently, there has been substantial literature debating the existence of the MBL phase and whether the MBL transition is really a phase transition or merely a crossover [53, 54, 60]. We do not assess their conclusions. However, based on the above analyses, the devised physical quantities in these works may have the following limitations.

(i) The models they studied generally have the $U(1)$ symmetry, but the computed physical quantities do not incorporate this symmetry thoroughly. Even worse is the case that the computational scheme they used may explicitly break those important internal symmetries of the system, such as the open-system approach used in [54] which spoils the $U(1)$ symmetry of the Heisenberg chain. This lack of the symmetry embedding, particularly in a way parallel to the entropy decomposition already implemented here in Appendix A, may result in the dangerous inclusion of the phase-transition-irrelevant components or contributions, which probably distorts the revealing of the underlying critical phenomena.

(ii) Most quantities designed might not be experimentally measurable.

(iii) Various transport calculations require to involve the information of the whole energy spectrum. As will be illustrated shortly, the higher- and lower-energy behaviors of the generic bosonic systems can be notably distinct. This means that these regimes need to be considered and examined carefully and separately. For instance, when calculating the transport properties of a bosonic chain at infinite temperature, blindly mixing all eigenstates equally may pose a problem.

Following the recent theoretical framework proposed by Ref. [54], the observed MBL phase is likely akin to what they call the finite-size MBL regime, with the phase transition point corresponding to the landmark ${\mathcal{L}}_r$ within the finite-size accuracy. In essence, while debates over the MBL phase and transition persist, we provide an experimentally practical quantity to aid the investigation of these critical issues.

3.3 Finite-system dynamical phase diagram

Based on the scaling analyses discussed in the previous subsection, we analogously extracted the phase transition points for a spectrum of other energy densities, ultimately producing the eigenstate phase transition diagram for the dBH model, as shown by Fig.4(b) and Fig.5. It can be seen that $S_{\mathrm{v}\mathrm{N}}$ and $|d(S_N)|$ yield nearly identical estimates for the phase boundaries.

Fig.4 (a) Numerical results on the two-body density-density correlation function $G_2$ for $L=12$ from the ED calculations. $G_2$ versus $d$ shows an exponential decay for $\mu =3.5$ and $7.0$ at $\epsilon =0.3$, from which the localization lengths might be extracted. (b) Phase diagram on the plane of energy density $\epsilon$ and disorder strength $\mu$ is constructed based on the color contour of the localization length $\xi$ and consists of the ETH and MBL regimes, separated by a line derived from the scaling analyses of the total entropy $S_{\mathrm{v}\mathrm{N}}$ and $|d(S_N)|$.

Full size|PPT slide

Fig.5 Landscape of the phase diagram from the adjacent gap ratio $r$, using a dBH chain of $L=12$. The right color bar labels $r$ approximately between the GOE limit $r_{\mathrm{G}\mathrm{O}\mathrm{E}}\approx 0.536$ and the Poisson limit $r_{\mathrm{P}\mathrm{o}\mathrm{i}.}\approx 0.386$. A blue square in the bottom-left corner is an artifact that may occur from an insufficient number of eigenstates, particularly in such an area where both $\epsilon$ and $\mu$ are small, i.e., essentially originating from the finite-size effect. In addition, schematics of the different MBL regions are indicated: scatter MBL and cluster MBL. The red (from $S_{\mathrm{v}\mathrm{N}}$) and cyan (from $|d(S_N)|$) squared symbols with the horizontal error bars for $\mu$ mark the position estimates of the phase boundaries at $\epsilon =0.05,0.1,0.15,0.2,$$0.25,0.3,0.35,0.4$.

Full size|PPT slide

Inspired by the experimental pertinence, we may use the measurable localization length to illustrate that the system’s thermalization states predominantly reside in the low-energy-density region. Firstly, the two-body density-density correlation function $G_2(d)$ is introduced to measure the correlation of the particles located at the two sites separated by a distance $d$ in an eigenstate with the eigenenergy closest to a given energy density $\epsilon$. Specifically, we borrow the definition from the Harvard experiment [39]

(6)$G_2(d)=-{\left\{{\displaystyle \frac{{{\displaystyle \sum }}_{i+d⩽L}[⟨{\hat{n}}_i{\hat{n}}_{i+d}⟩-⟨{\hat{n}}_i⟩⟨{\hat{n}}_{i+d}⟩]}{N_{i+d⩽L}}}\right\}}_{\mathrm{a}\mathrm{v}},$

where $N_{i+d⩽L}$ indicates the number of pairs in the summation over $i$ and the given $d$, and ${\left\{...\right\}}_{\mathrm{a}\mathrm{v}}$ denotes an average over the disorder samples with a given set of $\mu$ and $\epsilon$. Since the total particle number is conserved for the present dBH model, one generally has $G_2(d)⩾0$. Further, one expects that $G_2(d)$ decreases exponentially with $d$ when the disorder strength is large enough. The localization length can thereby be extracted from this exponentially decaying behavior of $G_2(d)$, alternatively providing a distinguishable signal between the ETH and MBL phases.

In Fig.4(a), the two-body density−density correlation function at $\epsilon =0.3$ is shown for the four disorder strengths. In the region of small disorder strength, the correlation function decays very slowly with the distance $d$. It is nearly a constant, indicating that all the system’s parts are highly entangled and the many-body eigenstates are well thermalized. As the disorder strength increases, an exponential decay behavior becomes sharper and sharper so that it can be fitted by a formula of the form $G_2(d)=A{\mathrm{e}}^{-d/\xi }$ and the localization length $\xi$ can be extracted thereof. (Note, in comparison, the localization length cannot be gained via fitting to an exponential form when the disorder is weak and the energy density is small.) Accordingly, the emergence and the change of the localization length tunable by $\mu$ reveal an eigenstate phase transition from the ETH state to the MBL state.

As illustrated in Fig.4(b), by exploiting the extracted localization lengths, for the 1D dBH model with $U/J=3$ and the half-filling, we find an ETH phase in the bottom-left corner (yellow color) and the MBL phase (blue color) in the rest of the $\mu$−$\epsilon$ parameter plane. As already stressed, the thermalized states predominantly distribute within the low-energy, weakly disordered region of the phase diagram.

Moreover, as an aside, when comparing the phase transition boundary to the color map of the localization length, one notices that the localization lengths are generally much greater than several lattice spacings in the ETH phase. Nonetheless, we would like to remind the reader the procedure that for a thin area of the very weak disorder strengths in the ETH phase, the values of these large localization lengths that cannot be determined by using the exponential fit of $G_2$ have all been truncated to the $L-1$ lattice spacings, the maximal possible value of $d$, so that one can still represent the ETH phase in Fig.4(b) in a coherent fashion.

To further confirm the construct of the phase diagram, it is instructive to study the complementary energy-resolved adjacent gap ratio defined by [61, 62, 63]

(7)$r_n={\displaystyle \frac{\mathrm{m}\mathrm{i}\mathrm{n}[\delta E_n,\delta E_{n+1}]}{\mathrm{m}\mathrm{a}\mathrm{x}[\delta E_n,\delta E_{n+1}]}},$

where $\delta E_n=E_n-E_{n-1}$ is the difference between the two successive eigenenergies arranged in an ascending order. For each disorder sample, we calculate the $16$ adjacent $r_n$ to obtain the averaged value $\overline{r}$ subject to a given energy density $\epsilon$, and $r=\{\overline{r}{\}}_{\mathrm{a}\mathrm{v}}$ denotes the average over the disorder samples. The ratio $r$ reflects how the disorder strength $\mu$ determines the nature of the level statistics at a specific energy density $\epsilon$. Fig.5 shows the landscape of $r$ versus the energy density $\epsilon$ and the disorder strength $\mu$ with $L=12$. One can see that the ETH state is located in the yellow region with $r$ approaching the Gaussian orthogonal ensemble (GOE) limit $r_{\mathrm{G}\mathrm{O}\mathrm{E}}\approx 0.536$. By contrast, the MBL state, represented in blue, tends to reach the Poisson limit $r_{\mathrm{P}\mathrm{o}\mathrm{i}.}\approx 0.386$ [64]. The results of the phase boundaries, obtained from finite-size analyses of $S_{\mathrm{v}\mathrm{N}}$ and $|d(S_N)|$, are also overlaid in this phase diagram for comparison. From Fig.5, it is apparent that there appears to be an intermediate region between the yellow-region boundary and the phase-transition boundary. In Appendix B, we explore this low-energy and medium-$\mu$ region further from the perspective of the local compressibility distribution.

In the region where the energy density is less than $0.5$, the overall phase transition resembles that in the full spin or fermion models [32–34]. It is anticipated that when the energy density is relatively low, it is not easy to form the high-energy situation where multiple particles occupy the same position. While, in the high-energy-density region, particles can and do form clusters due to interaction and bosonic statistics, yet possibly with very different behaviors between the limits of small and large disorder strengths.

As shown in Fig.4(b) and Fig.5, there is basically no significant change in the color of the localized part of the diagram. But we are now equipped to elucidate the properties and distinctions of these phase regions in some detail by scrutinizing the eigenstate properties and the quench dynamics in the high- and low-energy portions of the phase diagram. These peculiarities of the dBH model, as will be described next in Section 4, exemplify the intricate interplay of interaction, quantum statistics, and disorder, and give the schematic delineation of the phase regions in Fig.5 as per the broad classification of scatter and cluster MBLs.

4 Many-body localized eigenstates

This section focuses on the properties of MBL eigenstates in the dBH chain. Our investigations reveal the following findings. i) In Section 4.1, we point out that the low-energy-density portion of the eigenspectrum of the dBH model closely resembles the entire eigenspectrum of the fermionic (or spin) models [32–34], featured by a transition from the thermalized to localized states. In view of the fact that the low-energy mobile identities are scatter particles, we name the resulting localized state the scatter MBL phase. In Section 4.1, we analyze the eigenstate particle number distribution for scatter MBL and summarize the associated quench dynamics characteristics. ii) Owing to the bosonic statistics, the high-energy-density portion of the eigenspectrum allows the formation of the distinct boson clusters in a random potential background. This situation is unique to the dBH-type model, and we name the resulting localized state the cluster MBL phase. After examining the eigenstate particle number distribution and the quench dynamics features, we show that the cluster MBL may differ in significant aspects from the scatter MBL. iii) Regarding the localized eigenstate properties, the traditional measures might not be directly sensitive to the differences between the high- and low-energy-density sections of the eigenspectrum. In Section 4.3, by exploring the structure of the symmetry-resolved entanglement entropy $S_{\mathrm{v}\mathrm{N}}^n$, we introduce another measure to help discriminate between the eigenstate phases emerging at the respective low and high energies.

4.1 Scatter MBL in the low-energy-density spectrum

First, we define the following two states that will be used below. Note, the subscripts in Eqs. (8) and (9) label the site indices.

(i) A line-shape initial product state, called the $l$-state, consists of one boson on each site of the left half chain. Pictorially, it is given by

(8)$|l⟩=|1_1,1_2,\dots ,1_{\frac{L}{2}},0_{\frac{L}{2}+1},\dots ,0_L⟩.$

(ii) A point-shape initial product state, called the $p$-state, is formed by putting all $N$ bosons onto the leftmost site. It is given by

(9)$|p⟩=|N_1,0_2,\dots ,0_L⟩.$

Their corresponding energy densities are shown in Appendix C. From there, it is clear that as disorder increases, the $l$-state always stays within the low-energy-density region, while the $p$-state, in comparison, always stays within the high-energy-density region.

Selecting a random sample, we first compute the energy of the $l$-state as per $E=⟨l|H|l{⟩}_{\mu =2}$. Next, we identify the 11 nearest eigenstates (5 above and 5 below the target eigenstate whose eigenenergy is closest to $E$ given above). Subsequently, we calculate the particle number distributions for these eigenstates at each lattice site. In Fig.6(a), the particle numbers at weak disorder $(\mu =2)$ are predominantly distributed around $0.5$. This is reasonable because the system is in a thermalized state whose particle numbers are distributed uniformly across the lattice sites. Consider that we are dealing with a half-filled system, then the particle number should indeed be at about $0.5$. In Ref. [65], we observe that for the weak disorder at $\mu =2$, starting from the $l$-state, $S_N$ grows with time as a function of $\ln t$, and its steady-state value grows with size as per $\ln L$. On the other hand, the companion $S_C$ grows with time linearly, and exhibits a volume-law scaling behavior at the long-time steady-state limit.

Fig.6 Particle number distributions of the 11 eigenstates, 5 above and 5 below the eigenstate closest in energy to the given initial state. (a) corresponds to the target energy $E=⟨l|H|l{⟩}_{\mu =2}$ realized in the region of low energy density and small disorder strength. (b) corresponds to $E=⟨l|H|l{⟩}_{\mu =20}$ within the low-energy and strong-disorder region. (c) gives the result of $E=⟨p|H|p{⟩}_{\mu =2}$ for the high energy density and weak disorder. (d) targets $E=⟨p|H|p{⟩}_{\mu =20}$ upon the high energy density and large disorder strength. The chain size here is fixed to be $L=14$.

Full size|PPT slide

We repeat the same calculation for the $l$-state but under the condition of strong disorder, i.e., $\mu =20$. As shown by Fig.6(b), in this scenario, some lattice sites exhibit particle clustering, while others remain nearly empty, indicating the particle localization within the chain. It is worth noticing that the maximum particle occupancy per site in Fig.6(b) remains low, and the situation where the majority of the particles occupy a single lattice site does not occur. We call this low-energy MBL phenomenon the scatter MBL. In Ref. [65], we show that for this strong-disorder condition at $\mu =20$, if starting from the $l$-state, the growth of $S_N$ slows down to $\ln \ln t$, and the scaling of its steady-state value obeys the area law. Correspondingly, the companion $S_C$ in this situation grows with time as per a function of $\ln t$, and its steady-state value fulfills a volume law.

Incidentally, with the channel-resolved analysis, we show that this slow $\ln \ln t$ growth of $S_N$ is mainly owing to the localization phenomenon and its final steady-state scaling remains to be an area law [66]. These observations are in contradiction with the prediction of the eventual thermalization for the similar situation in the disordered spin-chain model as was claimed in recent literature [67].

In the recent seminal experimental work on the dBH model [39], it was observed that within the MBL phase, $S_N$ follows an area-law scaling behavior, while concurrently it was implicit that $S_C$ follows a volume-law scaling behavior ($S_C$ cannot be measured directly, so a closely related correlator was measured experimentally as an alternative). This suggests that the low-energy scatter MBL behavior has been experimentally detectable, but to our knowledge, the subsequent discussion regarding the high-energy cluster MBL behaviors has never been reported before.

4.2 Cluster MBL in the high-energy-density spectrum

We do the same calculation by switching the initial state from the low-energy $l$-state to the high-energy $p$-state. For weak disorder $(\mu =2)$, as depicted in Fig.6(c), we observe that the 11 nearest eigenstates at the target energy $E=⟨p|H|p{⟩}_{\mu =2}$ have a common feature in their particle number distributions. Concretely, almost all the particles are concentrated to occupy a single lattice site, while the particle numbers at other sites decay exponentially. This signifies a strong degree of localization within the system; we call it the hard-type cluster MBL. This is because, in the high-energy-density regime, Bose statistics allows for the formation of the multi-particle clusters. Without disorder, these clusters can superpose to form the cat like states, but disorder breaks the translational symmetry, thus obstructing their coherent superposition and consequently giving rise to the localized, isolated clusters, ultimately resulting in the robust cluster MBL.

In Ref. [65], we find that, for the situation of low disorder $(\mu =2)$, starting from the $p$-state, both $S_N$ and $S_C$ exhibit almost no growth over time, with the steady-state values following the area law. Loosely speaking, these behaviors resemble the Anderson localization in some aspects.

Similar calculations upon the $p$-state are performed under the strong-disorder condition $(\mu =20)$ as well. It can be seen from Fig.6(d) that, at this circumstance, the system also exhibits the distinct clustering behavior. But it is no longer as extreme as in Fig.6(c), where nearly all particles are localized onto a single point, nor is it close to Fig.6(b), where particles are localized to only reach the relatively low local concentrations. Therefore, we refer to this situation as the soft-type cluster MBL.

In Ref. [65], we demonstrate that for this strong disorder at $\mu =20$, starting from the $p$-state, $S_N$ exhibits the diminished time-dependent growth, with the steady-state values following an area law. By contrast, the companion $S_C$ keeps increasing slowly over time in a tentative $\ln \ln t$ format, but its steady-state values seem to still obey the area law.

4.3 Fine structure of $S_{\mathrm{v}\mathrm{N}}$

Curiously, is there a single quantity that can both reveal the thermalized-state distribution and differentiate between the behaviors within the low- and high-energy-density sections of the phase diagram? Through exploiting the channel reflection symmetry of the set $\{S_{\mathrm{v}\mathrm{N}}^n\}$, we introduce a new physical quantity that effectively addresses this question. Complementarily, this might provide a different means to probe the phase transition zone from the pure quantum entanglement $\{S_{\mathrm{v}\mathrm{N}}^n\}$ via its intrinsic structure.

As claimed in Appendix A, for the symmetry-resolved entropy $\{S_{\mathrm{v}\mathrm{N}}^n\}$, the relation $S_{\mathrm{v}\mathrm{N}}^{n_A=n}=S_{\mathrm{v}\mathrm{N}}^{n_B=N-n}$ always holds. In the absence of disorder, the preserved spatial reflection symmetry guarantees the general satisfaction of $S_{\mathrm{v}\mathrm{N}}^{n_A=n}=S_{\mathrm{v}\mathrm{N}}^{n_A=N-n}$, which generates a channel reflection symmetry of $\{S_{\mathrm{v}\mathrm{N}}^n\}$, as already illustrated in the inset of Fig.1(b) above. It turns out that we can determine the degree of deviation from this symmetry by the summed difference of the corresponding pairs of channels, which is equivalent to determining the degree of deviation from the ideal thermalization. Thus, we try to define the pertinent physical quantity as follows:

(10)$|d(S_{\mathrm{v}\mathrm{N}}^n)|={\left\{{\displaystyle \sum _{n=1}^M|S_{\mathrm{v}\mathrm{N}}^n-S_{\mathrm{v}\mathrm{N}}^{N-n}|}\right\}}_{\mathrm{a}\mathrm{v}},$

where

(11)$M=\{\begin{array}{l}(N-2)/2\mathrm{f}\mathrm{o}\mathrm{r}N=\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\\ (N-1)/2\mathrm{f}\mathrm{o}\mathrm{r}N=\mathrm{o}\mathrm{d}\mathrm{d}\end{array}.$

Since $S_{\mathrm{v}\mathrm{N}}^{n_A=0}=S_{\mathrm{v}\mathrm{N}}^{n_A=N}=0$, the summation starts from $n=1$. In the ETH region, the left and right parts are approximately symmetric, giving rise to a small value of $|d(S_{\mathrm{v}\mathrm{N}}^n)|$. In the scatter MBL region, although the relative difference is nonzero, the value of $|d(S_{\mathrm{v}\mathrm{N}}^n)|$ remains small because of the absolute smallness of the values of $\{S_{\mathrm{v}\mathrm{N}}^n\}$ themselves. Taken together, these trends produce a peak in the middle, corresponding to the critical regime of the underlying phase transition.

Further, from Fig.7, we can see that the ETH region corresponds to the valley in the 3D figure of $|d(S_{\mathrm{v}\mathrm{N}}^n)|$. The closer to the bottom of the valley, the lower the degree of the channel-reflection-symmetry breaking is and the more thorough the thermalization of the system assumes. There are two peaks at the valley’s boundaries, corresponding to the location of the energy-resolved phase transition. When the two peaks merge and disappear, $|d(S_{\mathrm{v}\mathrm{N}}^n)|$ becomes flattened, signalling that the system begins to enter the scatter MBL region. Between the high- and low-energy-density sections, a distinctive dip effectively separates the cluster MBL at high energies and the scatter MBL at low energies. In this regard, $|d(S_{\mathrm{v}\mathrm{N}}^n)|$ may possess the capacity to distinguish between these MBL phases and aptly identify the thermalized region at the low energies.

Fig.7 3D plot of $|d(S_{\mathrm{v}\mathrm{N}}^n)|$ at $L=12$, where the phase transition boundary and the three different phases can be perceived.

Full size|PPT slide

The above content indicates that the internal structure of the entanglement entropy effectively distinguishes between different phase regions, upon resort to the concept of symmetry decomposition. Once again, one can appreciate the importance of “symmetry combined with entanglement” in the study and elucidation of the MBL phases and transitions.

5 Summary and discussion

This paper provides a comprehensive study on the eigenstate properties of the dBH model in one dimension.

Firstly, we explore the transition from the thermalized state to the MBL state at low energies. The von Neumann entropy, $S_{\mathrm{v}\mathrm{N}}$, is often used to detect the MBL transition, even though it cannot be directly measured. Utilizing $U(1)$ symmetry, the von Neumann entropy is decomposed into the particle number entropy, $S_N$, and the configurational entropy, $S_C$. Notably, $S_C$ approaches zero within the localized phase. We take advantage of this trend and define $|d(S_N)|$ to quantify the deviation of $S_N$ from the ideal thermal distribution. Finite-size scaling analysis shows that $|d(S_N)|$ shares the same critical point with $S_{\mathrm{v}\mathrm{N}}$ but exhibits a larger critical exponent, suggesting that the transition is primarily set by $S_N$ and its fluctuations. It is worth noting that $|d(S_N)|$ might be measured via protocols already used by recent experiments, particularly in the case of bosonic systems. Recent debates about whether MBL transitions are true phase transitions or crossovers prompt us to seek such an experimentally measurable quantity that can contribute to addressing these concerns.

Secondly, the paper contrasts the distinctive behaviors occurring in the high- and low-energy-density regimes of the phase diagram for the dBH model. We find that: i) The thermalizing states predominantly occupy the low-energy spectrum, as evidenced by energy density distribution, localization lengths, and gap ratios. ii) The low-energy spectrum closely resembles the fermionic (or spin) models, where the transition from the thermalizing state to the scatter MBL phase is visible. We analyze the characteristics of these low-energy localized eigenstates, and their dynamics are also summarized. iii) Bose statistics leads to the formation of distinct clusters in the high-energy spectrum, resulting in the unique cluster MBL phase that significantly differs from the scattering MBL phase in low energies, as reflected through the eigenstate and dynamics properties. iv) Traditional measures such as localization length and gap ratios cannot directly reveal the differences in the localization phenomena between the high- and low-energy sections. We define a new quantity, $|d(S_{\mathrm{v}\mathrm{N}}^n)|$, from the internal structure of $S_{\mathrm{v}\mathrm{N}}$, which effectively captures the distinctions between phase regions in high and low energies.

The two quantities we defined based on the decomposition of entropy by symmetry, $|d(S_N)|$ and $|d(S_{\mathrm{v}\mathrm{N}}^n)|$, prove to be useful in analyzing critical points and distinguishing between different phase regions. This shows that symmetry plays a key role in MBL phases and phase transitions. Currently, a grand portion of the MBL research has not given due attention to this aspect, and the physical quantities employed to analyze MBL phases and transitions have not properly embedded the symmetry information, particularly in those contentious issues.

As can be noticed, there needs a further study on whether the MBL regions are separated by phase transitions or crossovers. Moreover, the bimodal compressibility distribution structure seen in Appendix B seems to predict the existence of a critical regime arising from the Griffiths effect [42, 56, 68–70]. Different implementations of the randomized samples, e.g., quasiperiodic fields, may give rise to qualitatively different results [58, 71, 72]. Although the numerical resource to precisely settle the localization problem for the dBH model is far beyond the current computational capacity, one might still imagine to resort to the available finite-size numerics to search for the emergent new physics when the filling factor and the Hubbard interaction are tuned within the wider parameter ranges.

References

[1] P. W. Anderson. Absence of diffusion in certain random lattices. Phys. Rev., 1958, 109(5): 1492

[2] E. Abrahams, P. W. Anderson, D. C. Licciardello, T. V. Ramakrishnan. Scaling theory of localization: Absence of quantum diffusion in two dimensions. Phys. Rev. Lett., 1979, 42(10): 673

[3] D. Vollhardt, P. Wölfle. Diagrammatic, self-consistent treatment of the Anderson localization problem in d ≤ 2 dimensions. Phys. Rev. B, 1980, 22(10): 4666

[4] S. John. Electromagnetic absorption in a disordered medium near a photon mobility edge. Phys. Rev. Lett., 1984, 53(22): 2169

[5] K. Arya, Z. B. Su, J. L. Birman. Localization of the surface plasmon polariton caused by random roughness and its role in surface-enhanced optical phenomena. Phys. Rev. Lett., 1985, 54(14): 1559

[6] Q. J. Chu, Z. Q. Zhang. Localization of phonons in mixed crystals. Phys. Rev. B, 1988, 38(7): 4906

[7] R. Nandkishore, D. A. Huse. Many-body localization and thermalization in quantum statistical mechanics. Annu. Rev. Condens. Matter Phys., 2015, 6(1): 15

[8] D. M. Basko, I. L. Aleiner, B. L. Altshuler. Metal–insulator transition in a weakly interacting many electron system with localized single-particle states. Ann. Phys., 2006, 321(5): 1126

[9] V. Oganesyan, D. A. Huse. Localization of interacting fermions at high temperature. Phys. Rev. B, 2007, 75(15): 155111

[10] M. Žnidarič, T. Prosen, P. Prelovšek. Many-body localization in the Heisenberg XXZ magnet in a random field. Phys. Rev. B, 2008, 77(6): 064426

[11] J. M. Deutsch. Quantum statistical mechanics in a closed system. Phys. Rev. A, 1991, 43(4): 2046

[12] M. Srednicki. Chaos and quantum thermalization. Phys. Rev. E, 1994, 50(2): 888

[13] M. Serbyn, Z. Papić, D. A. Abanin. Local conservation laws and the structure of the many-body localized states. Phys. Rev. Lett., 2013, 111(12): 127201

[14] A. Chandran, I. H. Kim, G. Vidal, D. A. Abanin. Constructing local integrals of motion in the many-body localized phase. Phys. Rev. B, 2015, 91(8): 085425

[15] S. D. Geraedts, R. N. Bhatt, R. Nandkishore. Emergent local integrals of motion without a complete set of localized eigenstates. Phys. Rev. B, 2017, 95(6): 064204

[16] L. D’Alessio, Y. Kafri, A. Polkovnikov, M. Rigol. From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics. Adv. Phys., 2016, 65(3): 239

[17] A. Lazarides, A. Das, R. Moessner. Fate of many-body localization under periodic driving. Phys. Rev. Lett., 2015, 115(3): 030402

[18] P. Ponte, Z. Papić, F. Huveneers, D. A. Abanin. Many-body localization in periodically driven systems. Phys. Rev. Lett., 2015, 114(14): 140401

[19] D. V. Else, B. Bauer, C. Nayak. Floquet time crystals. Phys. Rev. Lett., 2016, 117(9): 090402

[20] N. Y. Yao, A. C. Potter, I. D. Potirniche, A. Vishwanath. Discrete time crystals: Rigidity, criticality, and realizations. Phys. Rev. Lett., 2017, 118(3): 030401

[21] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, C. Monroe. Observation of a discrete time crystal. Nature, 2017, 543(7644): 217

[22] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, M. D. Lukin. Observation of discrete time-crystalline order in a disordered dipolar many-body system. Nature, 2017, 543(7644): 221

[23] J. A. Kjäll, J. H. Bardarson, F. Pollmann. Many-body localization in a disordered quantum Ising chain. Phys. Rev. Lett., 2014, 113(10): 107204

[24] M. Brenes, M. Dalmonte, M. Heyl, A. Scardicchio. Many-body localization dynamics from gauge invariance. Phys. Rev. Lett., 2018, 120(3): 030601

[25] E. Levi, M. Heyl, I. Lesanovsky, J. P. Garrahan. Robustness of many-body localization in the presence of dissipation. Phys. Rev. Lett., 2016, 116(23): 237203

[26] M. H. Fischer, M. Maksymenko, E. Altman. Dynamics of a many-body-localized system coupled to a bath. Phys. Rev. Lett., 2016, 116(16): 160401

[27] H. P. Lüschen, P. Bordia, S. S. Hodgman, M. Schreiber, S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, U. Schneider. Signatures of many-body localization in a controlled open quantum system. Phys. Rev. X, 2017, 7(1): 011034

[28] J. Ren, Q. Li, W. Li, Z. Cai, X. Wang. Noise-driven universal dynamics towards an infinite temperature state. Phys. Rev. Lett., 2020, 124(13): 130602

[29] M. Žnidarič, A. Scardicchio, V. K. Varma. Diffusive and subdiffusive spin transport in the ergodic phase of a many-body localizable system. Phys. Rev. Lett., 2016, 117(4): 040601

[30] S. P. Kelly, R. Nandkishore, J. Marino. Exploring many-body localization in quantum systems coupled to an environment via Wegner‒Wilson flows. Nucl. Phys. B, 2020, 951: 114886

[31] C. Chamon, A. Hamma, E. R. Mucciolo. Emergent irreversibility and entanglement spectrum statistics. Phys. Rev. Lett., 2014, 112(24): 240501

[32] C. R. Laumann, A. Pal, A. Scardicchio. Many-body mobility edge in a mean-field quantum spin glass. Phys. Rev. Lett., 2014, 113(20): 200405

[33] D. J. Luitz, N. Laflorencie, F. Alet. Many-body localization edge in the random-field Heisenberg chain. Phys. Rev. B, 2015, 91: 081103(R)

[34] I. Mondragon-Shem, A. Pal, T. L. Hughes, C. R. Laumann. Many-body mobility edge due to symmetry-constrained dynamics and strong interactions. Phys. Rev. B, 2015, 92(6): 064203

[35] R. Mondaini, Z. Cai. Many-body self-localization in a translation-invariant Hamiltonian. Phys. Rev. B, 2017, 96(3): 035153

[36] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, I. Bloch. Observation of many-body localization of interacting fermions in a quasirandom optical lattice. Science, 2015, 349(6250): 842

[37] J. H. Bardarson, F. Pollmann, J. E. Moore. Unbounded growth of entanglement in models of many-body localization. Phys. Rev. Lett., 2012, 109(1): 017202

[38] M. Serbyn, Z. Papić, D. A. Abanin. Universal slow growth of entanglement in interacting strongly disordered systems. Phys. Rev. Lett., 2013, 110(26): 260601

[39] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, M. Greiner. Probing entanglement in a many-body-localized system. Science, 2019, 364(6437): 256

[40] J. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, C. Gross. Exploring the many-body localization transition in two dimensions. Science, 2016, 352(6293): 1547

[41] Q. Guo, C. Cheng, Z. H. Sun, Z. Song, H. Li, Z. Wang, W. Ren, H. Dong, D. Zheng, Y. R. Zhang, R. Mondaini, H. Fan, H. Wang. Observation of energy-resolved many-body localization. Nat. Phys., 2021, 17(2): 234

[42] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai, J. Léonard, M. Greiner. Quantum critical behaviour at the many-body localization transition. Nature, 2019, 573(7774): 385

[43] J. Léonard, S. Kim, M. Rispoli, A. Lukin, R. Schittko, J. Kwan, E. Demler, D. Sels, M. Greiner. Probing the onset of quantum avalanches in a many-body localized system. Nat. Phys., 2023, 19(4): 481

[44] J. M. Zhang, R. X. Dong. Exact diagonalization: The Bose–Hubbard model as an example. Eur. J. Phys., 2010, 31(3): 591

[45] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, C. Hubig. Time-evolution methods for matrix-product states. Ann. Phys., 2019, 411: 167998

[46] G. Vidal. Efficient simulation of one-dimensional quantum many-body systems. Phys. Rev. Lett., 2004, 93(4): 040502

[47] M. P. Kennett. Out-of-equilibrium dynamics of the Bose–Hubbard model. ISRN Cond. Matter Phys., 2013, 393616: 39

[48] P. Sengupta, S. Haas. Quantum glass phases in the disordered Bose–Hubbard model. Phys. Rev. Lett., 2007, 99(5): 050403

[49] V. Gurarie, L. Pollet, N. V. Prokof’ev, B. V. Svistunov, M. Troyer. Phase diagram of the disordered Bose–Hubbard model. Phys. Rev. B, 2009, 80(21): 214519

[50] L. Pollet, N. V. Prokof’ev, B. V. Svistunov, M. Troyer. Absence of a direct superfluid to Mott insulator transition in disordered Bose systems. Phys. Rev. Lett., 2009, 103(14): 140402

[51] M. Gerster, M. Rizzi, F. Tschirsich, P. Silvi, R. Fazio, S. Montangero. Superfluid density and quasi-long-range order in the one-dimensional disordered Bose–Hubbard model. New J. Phys., 2016, 18(1): 015015

[52] S. Hu, Y. Wen, Y. Yu, B. Normand, X. Wang. Quantized squeezing and even‒odd asymmetry of trapped bosons. Phys. Rev. A, 2009, 80(6): 063624

[53] J. Šuntajs, J. Bonča, T. Prosen, L. Vidmar. Quantum chaos challenges many-body localization. Phys. Rev. E, 2020, 102(6): 062144

[54] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, D. A. Huse. Avalanches and many-body resonances in many-body localized systems. Phys. Rev. B, 2022, 105(17): 174205

[55] B. Widom. Surface tension and molecular correlations near the critical point. J. Chem. Phys., 1965, 43(11): 3892

[56] T. Vojta. Phases and phase transitions in disordered quantum systems. AIP Conf. Proc., 2013, 1550: 188

[57] T. Orell, A. A. Michailidis, M. Serbyn, M. Silveri. Probing the many-body localization phase transition with superconducting circuits. Phys. Rev. B, 2019, 100(13): 134504

[58] S. X. Zhang, H. Yao. Universal properties of many body localization transitions in quasiperiodic systems. Phys. Rev. Lett., 2018, 121(20): 206601

[59] V. Khemani, D. N. Sheng, D. A. Huse. Two universality classes for the many-body localization transition. Phys. Rev. Lett., 2017, 119(7): 075702

[60] J. Šuntajs, J. Bonča, T. Prosen, L. Vidmar. Ergodicity breaking transition in finite disordered spin chains. Phys. Rev. B, 2020, 102(6): 064207

[61] P. Sierant, D. Delande, J. Zakrzewski. Many-body localization for randomly interacting bosons. Acta Phys. Pol. A, 2017, 132(6): 1707

[62] P. Sierant, D. Delande, J. Zakrzewski. Many-body localization due to random interactions. Phys. Rev. A, 2017, 95: 021601(R)

[63] T. Orell, A. A. Michailidis, M. Serbyn, M. Silveri. Probing the many-body localization phase transition with superconducting circuits. Phys. Rev. B, 2019, 100(13): 134504

[64] Y. Y. Atas, E. Bogomolny, O. Giraud, G. Roux. Distribution of the ratio of consecutive level spacings in random matrix ensembles. Phys. Rev. Lett., 2013, 110(8): 084101

[65] J.ChenC. ChenX.Wang, Energy- and symmetry-resolved entanglement dynamics in disordered Bose−Hubbard chain, arXiv: 2303.14825 (2023)

[66] D. J. Luitz, Y. B. Lev. Absence of slow particle transport in the many-body localized phase. Phys. Rev. B, 2020, 102: 100202(R)

[67] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, J. Sirker. Evidence for unbounded growth of the number entropy in many-body localized phases. Phys. Rev. Lett., 2020, 124(24): 243601

[68] D. J. Luitz, N. Laflorencie, F. Alet. Extended slow dynamical regime close to the many-body localization transition. Phys. Rev. B, 2016, 93: 060201(R)

[69] T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. Das Sarma, I. Bloch, M. Aidelsburger. Observation of many-body localization in a one-dimensional system with a single-particle mobility edge. Phys. Rev. Lett., 2019, 122(17): 170403

[70] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse, M. Knap. Rareregion effects and dynamics near the many-body localization transition. Ann. Phys., 2017, 529(7): 1600326

[71] V. Khemani, D. N. Sheng, D. A. Huse. Two universality classes for the many-body localization transition. Phys. Rev. Lett., 2017, 119(7): 075702

[72] V. Khemani, S. P. Lim, D. N. Sheng, D. A. Huse. Critical properties of the many-body localization transition. Phys. Rev. X, 2017, 7(2): 021013

[73] S. Schierenberg, F. Bruckmann, T. Wettig. Wigner surmise for mixed symmetry classes in random matrix theory. Phys. Rev. E, 2012, 85(6): 061130

[74] S. D. Geraedts, R. Nandkishore, N. Regnault. Many-body localization and thermalization: Insights from the entanglement spectrum. Phys. Rev. B, 2016, 93(17): 174202

Acknowledgements

We thank Zi Cai, Changfeng Chen, Shijie Hu, Haiqing Lin, Mingpu Qin, Jie Ren, and Wei Su for fruitful discussions. This work is supported by the Ministry of Science and Technology of China (Grant No. 2016YFA0300500) and the National Natural Science Foundation of China (No. 11974244). C. C. acknowledges the support from the SJTU start-up fund and the sponsorship from Yangyang Development fund.

RIGHTS & PERMISSIONS

2024 Higher Education Press