Experimental characterization of quantum many-body localization transition

As strength of disorder enhances beyond a threshold value in many-body systems, a fundamental transformation happens through which the entire spectrum localizes, a phenomenon known as many-body localization. This has profound implications as it breaks down fundamental principles of statistical mechanics, such as thermalization and ergodicity. Due to the complexity of the problem, the investigation of the many-body localization transition has remained a big challenge. The experimental exploration of the transition point is even more challenging as most of the proposed quantities for studying such effect are practically infeasible. Here, we experimentally implement a scalable protocol for detecting the many-body localization transition point, using the dynamics of a $N=12$ superconducting qubit array. We show that the sensitivity of the dynamics to random samples becomes maximum at the transition point which leaves its fingerprints in all spatial scales. By exploiting three quantities, each with different spatial resolution, we identify the transition point with excellent match between simulation and experiment. In addition, one can detect the evidence of mobility edge through slight variation of the transition point as the initial state varies. The protocol is easily scalable and can be performed across various physical platforms.


INTRODUCTION
Ergodicity and thermalization principles are the foundations of statistical mechanics which imply that a many-body system forgets its local information as it evolves [1,2]. Strikingly, these principles can be violated when the thermalizing dynamics leads to the conservation of local information [3,4]. Many-Body Localization (MBL) is the most recognized phenomenon for this remarkable feature [5][6][7][8]. The MBL is the reminiscent of Anderson localization [9] when the particles of a disordered many-body system interact. The most mysterious feature of the MBL physics is the transition point at which an ergodic phase transforms into a localized one as the disorder strength increases. Despite several theoretical breakthroughs for characterizing the transition point [10][11][12][13][14][15][16][17], its experimental observation is extremely challenging limited to low filling factors [18] or weak interaction [19] regimes.
The MBL transition takes place as the strength of disorder exceeds a threshold value in comparison with the interaction coupling. Unlike quantum phase transition, which only affects the ground state, the MBL transition is more drastic and leaves its impact on the whole spectrum. This makes it difficult to investigate as, for instance, the-oretically it can only be explored via exact diagonalization which restricts us to short chains [20]. Interestingly, each energy eigenstate localizes at a different disorder strength, with the mid-spectrum eigenvectors requiring the maximum disorder. This phenomena, schematically shown in Fig. 1a, is known as mobility edge which has been investigated theoretically [16] and some of its features have been observed experimentally [21,22]. While ergodic and MBL phases are well explored through thermalization studies [2] and investigating local conservation laws [4], several interesting features emerge around the MBL phase transition point which are less understood. This includes scaling properties [10,11,16], anomalous transport [23], critical slowing down [24] and a change in entanglement behavior [17] and energy level statistics [25]. Hence, detection, characterization and understanding the MBL transition point is highly desirable for both fundamental and practical purposes. Most of the quantities, e.g. von Neumann entropy [17], level spacing statistics [25], Schmidt gap [10] and entanglement negativity [11], which have been introduced for identifying the transition point are not experimentally friendly, demanding either costly state tomography protocols [26] or the full knowledge of the energy spectrum [16]. Therefore, most of the experiments [1,18,19,[26][27][28][30][31][32] are performed either in the ergodic phase or deep in the MBL regime, leaving the MBL transition point unexplored.
Here, we propose an experimentally implement a protocol for detecting the MBL transition point using a superconducting transmon qubit array with size N = 12 qubits. We initialize the system in various product FIG. 1. Schematics of the experiment. a, The illustration of the phase diagram of the system across the whole energy spectrum versus disorder strength. Each energy eigenstates localizes at a different disorder strength resulting in many-body mobility edge. b, The twelve transmon qubits, illustrated as red crosses, are arranged in a 1D chain. The direct coupling between them is realized via capacitors. Each qubit has individual Z (yellow) and XY (green) control lines for state manipulation. The twelve readout resonators (blue) are dispersively coupled to their corresponding qubits, and then divided into two groups to couple to the transmission lines (orange). c, The three steps of the experimental procedure are depicted, namely initialization, evolution and readout. For state preparation, the corresponding qubits are excited to |1 by applying X gates. After that, for the state evolution, all the qubits are tunned to their working frequencies ω + h with h being a random number in the interval [−h, h]. After time t, for readout, the qubits are detuned to stop the evolution and then simultaneously measured in σz direction. d, Our protocol can be understood in a simple way. Each circle represents a qubit and the filling color shows the corresponding average population. Namely, empty, full and fractional filling colors stand for n = 0, n = 1 and 0 < n < 1, respectively. The dynamics start with a given initial state. Each column shows the average population of each site at long time dynamics for a specific random instance of qubit frequencies. In the ergodic regime, the system thermalizes with n ∼ 1/2 and thus the variation with respect to different random realizations is small. Deep in the MBL regime the dynamics is almost frozen and n remains close to its initial value leading also to weak dependence to the different random realizations. Remarkably, around the MBL transition, the long time dynamics shows a strong dependence on the random potential. We exploit this feature to identify the MBL transition point.
states and let it evolve under the action of its disordered Hamiltonian until it reaches local equilibrium. Then, we measure three different quantities, namely time auto-correlation, number entropy and Hamming distance, which capture different spatial resolutions varying from single site to a block of finite size and the entire system, respectively. To see the effect of mobility edge, various initial states with different overlap pattern with the eigenstates of the system are considered. Each quantity is measured for several random ensembles. While the averaged quantities vary smoothly across the phase diagram, their standard deviation with respect to different random ensembles peaks at the critical point, revealing the MBL transition.

SUPERCONDUCTING QUBIT ARRAY
We demonstrate our experiment on an array of N =12 superconducting transmon qubits, described by Bose-Hubbard model ( =1) Dynamics of the quantities and equilibrium values. The non-equilibrium dynamics of the system, initialized in a Neel state |0, 1, · · · , 1 , averaged over five random instances, for the three quantities: a, the auto-correlation C(t); b, the number entropy of half-chain S (N/2) (t) ; and c,the Hamming distance D(t ; and f, the Hamming distance Deq . Again, the solid lines are numerical simulations and markers are the experimental data. All the three quantities smoothly changes from the ergodic to the MBL phase without revealing the MBL transition point. The standard deviation, depicted as shadow for numerical simulation and error bars for experimental data, quantifies the uncertainty in our estimations. We note a good agreement between simulation and experiment. where a † ( a ) are the bosonic creation (annihilation) operators at site and n = a † a is the corresponding number operator. The capacitive dipole-dipole interaction leads to the nearest-neighbor hopping J 2π×1.2 MHz. In order to generate disorder in the potential of the system, the qubit frequencies are tuned by DC and pulse signals to be the sum of a constant central frequency ω and a random value h which is drawn from a uniform distribution [−h, h] with h being the disorder strength. The nonlinear on-site interaction U ≈ −22J 1 represents the excess of energy needed for having more than one boson at each site. The schematic of the quantum simulator is shown in Fig. 1b and more details about the device and its parameters are given in the supplementary material [33].

EXPERIMENTAL PROTOCOL
We initialize the system in 10 different product states |Ψ s (0) = |s = |s 1 , s 2 , · · · , s N , (where s = 0, 1 represents the number of bosons at site ) such that the filling factor is f = N =1 n /N = 1/2. The system evolves under the action of the Hamiltonian as |Ψ s (t) = e −ı Ht |Ψ s (0) and then the population configuration is measured for all sites. The protocol is schematically shown in Fig. 1c. We note that, the Hamiltonian commutes with the total number of particles, i.e.  . Averaging over all given initial states leads to an initial state independent probability distribution P (Qeq). By considering different disorder strengths h/J1 = 1 (ergodic), h/J1 = 2 and h/J1 = 3 (near the MBL transition) and h/J1 = 7 (MBL) we estimate the probability distribution of: a, the auto-correlation function P (Ceq); b, the number entropy of half chain P (S (N/2) eq ); and c, the Hamming distance P (Deq). In all the figures, the small rugs are the experimental data, the bars are the histograms and the solid lines are the fitting distributions. In the ergodic regime, the distribution is narrow as the system thermalizes and the final outcomes are determined solely by the thermal state. By increasing the disorder strength h/J1, the outcomes can take a larger range of values as the final results heavily depend on the random potential pattern. Consequently, the probability distribution gets wider. As we go into the MBL regime, the dynamics tends to get frozen, so the final outcomes are mainly determined by the initial state and weakly depend on the random potential pattern. Therefore, the distribution gets narrower again.
The measurement outcomes which do not have the right filling factors are excluded (see [33] for detailed discussions). To ensure we observe the general behavior of the disordered Hamiltonian, we also consider R = 60 distinct random realizations for each initial state. This is accompanied by numerical simulations in which the bosonic modes of each site are truncated to n max = 3.
We first consider the auto-correlation function which we adopt from Ref. [34] and is defined as In the ergodic phase, the local population thermalizes at long times reaching n ∼1/2 which results in C∼0. Deep in the MBL phase the evolution of the system is almost frozen, namely entanglement grows only logarithmically in time [35]. This means that n (t)∼ n (0), which then leads to C∼1.
To go beyond the single site resolution, it is highly desirable to investigate entanglement dynamics in the system. In practice evaluating a true measure of entanglement, such as the von Neumann entropy, is challenging as it demands full quantum state tomography. Here, we instead measure the number entropy [36,37] for a block of size m (for 1 ≤ m ≤ N/2), which is defined as where p n is the probability of finding n bosons in the block of size m. Since the system conserves the total number of bosons, the number entropy S (m) (t) is a lower bound for the von Neumann entropy [37] and recently has attracted a lot of attention [36]. While in the main text we focus on m = N/2, we provide more detailed analysis for other choices of m in the Supplementary Materials [33]. Finally, we consider a global quantity, namely Hamming distance [38], which quantifies how different configurations emerge in the global wave-function of the system making it distinct from the initial one. For a system initialized in |Ψ s (0) = |s , the Hamming distance is define as where s is the configuration that one finds at the output with probability P (s , t) and 0 ≤ d(s , s) ≤ 1 is the classical Hamming distance between the two configurations s and s (i.e. the number of flips that converts s to s divided by N for normalization). In the ergodic phase D(t) is expected to grow in time, ideally reaching 1, due to superposition of multiple configurations in the wavefunction of the system. In contrast, in the MBL phase the Hamming distance remains small as the dynamics is almost frozen and cannot generate many new configurations. The Hamming distance has already been employed to distinguish the ergodic and MBL phases experimen- ; and f, the Hamming distance ∆Deq. The shadow around the theory simulations represents the behavior for various initial states.
Following the footsteps of previous experiments [1,18,19,26,28,[30][31][32], we first plot the three quantities, namely C(t), S (N/2) (t) and D(t), as a function of time for one instance of random realization with different disorder strengths varying from ergodic to MBL phases in Figs. 2a-c. All the three quantities tend to equilibrate after a short transition time and we find good match between the experiments and our numerical simulations, showing that the unitary evolution under the action of the Hamiltonian in Eq. (1) reasonably simulates the behavior of the real quantum device. The scrambling nature of the dynamics in the ergodic phase (i.e. h/J 1 =1 ) makes it very distinct from and an almost frozen evolution in the MBL phase (i.e. h/J 1 =7). In order to better characterize the difference between the ergodic and MBL phases, we focus on the equilibrium values. Let Q (s,r) (t) represent any of the three quantities which depends on time t, a single random instance r and the initial state s. We define the equilib- where t i 's are the measured times in the equilibrium regime and we specifically choose M = 5 points in range 7.9 ≤ Jt ≤ 10.8.
To have a statistical analysis of the behavior at the equilibrium with respect to random realizations for a given disorder strength h, we define the ensemble average and its corresponding standard deviation ∆Q (s) where R = 60 is the total number of random realizations. In order to see the behavior of the system independent of the choice of s, one can average over different initial states to get Q eq = 1 I s Q (s) eq and ∆Q eq = 1 I s ∆Q (s) eq , where I = 10 is the total number of initial states (see [33] for the exact choices of the initial states). In Figs. 2d-f, we plot the average quantities C eq , S (N/2) eq and D eq as a function of the disorder strength h/J 1 . All the quantities clearly show a transition from ergodic to MBL as h/J 1 increase. Although, none of these quantities can directly reveal the transition point from the ergodic to the MBL phase as all of them vary smoothly throughout the phase diagram. However, as we will see later, it is indeed the averaged standard deviation ∆Q eq which is the crucial quantity to characterize the MBL transition point.
It is often more insightful to investigate a probability distribution directly rather than its average value. In our experiment, thanks to fairly large number of random samples (i.e. R = 60), one can estimate the probability distribution of P (Q (s) eq ) at each disorder strength h. To be initial state independent, one can also average over initial states to obtain P (Q eq ) = 1 ) and P (D eq ) are depicted for three different disorder strengths. We notice that, in the ergodic regime the shape of the distribution is almost Gaussian with a small variance. This is due to the scrambling dynamics in the ergodic regime which thermalizes the system locally making it indistinguishable for different random instances (i.e. absence of memory). As the disorder strength increases the probability distribution deviates from Gaussianity and gets wider due to strong dependence on the random ensembles. In fact, near the MBL transition point (h/J 1 ∼ 3) an interplay between thermalization of the ergodic phase and localization of the MBL regime makes the system very sensitive to the random instances resulting in a wide distribution. By further increasing the disorder the system enters the MBL regime where the dynamics is almost frozen (i.e. emergence of memory) and local subsystems remain close to their initial value for all random instances. This makes the role of each random instance irrelevant and thus the distribution gets narrower again. In Fig. 1d, this phenomenon is explained schematically.
Inspired by the probability distribution in Figs. 3a-c, one can exploit the width of the distribution for each initial state |s , quantified through the standard deviation ∆Q and ∆D (s) eq , each for two different initial states |s , as a function of the disorder strength h/J 1 . Strikingly, ∆Q (s) eq shows a very clear peak for all the three quantities, identifying the MBL transition point. Furthermore, each initial state peaks at a slightly different disorder strength h/J 1 . This is a clear evidence of the mobility edge as each initial state has a different overlap pattern with the eigenstates of the Hamiltonian and thus localizes at a different h/J 1 . One can get 10 different transition points corresponding to I = 10 initial states for each of the three quantities (see [33] for the exact values). It is well-known that for any given length, due to the finite size effect and the different convergence rates, distinct quantities may result in slightly different values of the MBL transition point [10]. To extract this point reliably, we adopt a data analysis procedure, based on Bayesian inference (see [33] for details). The transition point is, thus, identified to be h c /J 1 = 2.7 from the experimental data and h c /J 1 = 3.3 from the numerical simulations, with standard deviations determined as 1.1 and 1.2, respectively. The uncertainty comes from the limited number of random samples. The two values of transition points are not only in excellent agreement but also fully consistent with a pure theoretical investigation, based on the von Neumann entropy [39,40], with h c /J 1 ∼2−3 for |U |/J 1 ∼20. The average behavior, which are independent of the initial state, are shown in Figs. 4d-f as a function of h/J 1 . All the three quantities peak at the transition point, matching with the theory prediction. CONCLUSION We have proposed and experimentally realized a protocol for detecting the elusive MBL transition point in a superconducting quantum simulator. The proposed protocol relies on the time evolution of a many-body system under the action of a disordered Hamiltonian. We have focused on the long time evolution where the system is expected to reach an equilibrium. Three quantities, namely auto-correlation, number entropy and Hamming distance, each with different spatial resolution, have been measured to estimate the MBL transition point. The protocol relies on the standard deviation of these quantities with respect to the random samples of the disordered potential. As our results show, across the phase diagram, the sensitivity to random samples is maximum at the MBL transition point resulting in a peak in the standard deviation. The MBL transition point, computed from the three different quantities, are fully consistent with each other and mach well with numerical simulations. In addition, evidence of mobility edge, represented by slightly different transition point for each initial state, can be observed.
A remarkable point about our results is that, neither of the three quantities demands the costly quantum state tomography, needed for measuring the von Neumann entropy [17] or the Schmidt gap [10], thus can be easily extended to larger systems and those platforms which cannot perform tomography measurements. In addition, our protocol is platform independent which provides a clear application for the Noisy Intermediate Scale Quantum (NISQ) simulators to shed light on a difficult problem in many-body physics. In future experiments, by changing the filling factors and developing adjustable couplers [41], we can investigate larger area of the phase diagram and complete the observation of mobility edge for all energy scales.

I. SUPERCONDUCTING QUANTUM PROCESSOR
The experiment is performed on a 12-qubit superconduting quantum processor [S1], arranged in a 1D array, as illustrated in Fig. 1b of the main text. The qubits are Xmon variant [S2] of transmon qubits [S3]. Each nearest-neighboring qubit pairs are coupled via the capacitance between them, yielding an average nearest-neighbor coupling strength of J 1 /2π 11.5 MHz and an average next nearest neighbor coupling of J 2 /2π 1.2 MHz. For each qubit, there are an inductively coupled flux (Z) and a capacitively coupled microwave (XY ) control lines to realize state manipulation. Twelve individual resonators are dispersively coupled to the qubits to realize state readout. The resonators are divided into two groups and each group couples to a transmission line. We use frequency multiplexing technology to simultaneously readout all qubits' states. At the outside of each transmission line, an impedancematched parametric amplifier (IMPA) is used to enhance the readout signal strength. The parameters of the device is listed in Table. S1. Performance of the qubits. f01 and f ah are the idle frequencies and anharmonicities, respectively. T1 and T * 2 are the energy relaxation time and dephasing time, respectively. fr is the frequency of the corresponding readout resonator. f00 and f11 are the probabilities of correctly reading out the qubits' state for |0 and |1 , respectively. IMPA gain is measured as the ratio of the readout signal strengths when IMPAs are turned on and off.

II. EXPERIMENTAL REALIZATION
As shown in Fig. 1c of the main text, the realization of this experiment consists of three steps: (i) state preparation; (ii) system evolution; and (iii) measurement. In state preparation, we apply X gates to corresponding qubits to prepare the initial product state. After that, we detune all qubits to their working points, whose frequencies are randomized distributed in [−h, h] in comparing with the central frequency 4.35 GHz. After the system evolves for a period t, we detune all qubits back to their idle frequencies and perform the simultaneous readout. In our experiment, we only perform the σ z projection measurements. We repeat each cycle for at least 3,000 times to obtain the statistical measurement results. The measurement outcomes which do not have the right filling factors are excluded. In the Supplementary Materials we discuss this in more details. Between the cycles, all qubits are biased at their idle points for 300 µs to initialize all the qubits to their |0 states.

III. CALIBRATION OF FREQUENCY ALIGNMENT
The error of frequency alignment mostly comes from the residual nonlinear Z cross talk [S4]. In our experiment, the central frequency is set at 4.35 GHz. Based on that, the working frequencies are randomized in range [−h, h] in comparing with the central frequency. The calibration of the frequency alignment is thus important. The single round of the calibration consists of two steps. In the first step, we prepare the working frequencies in two different cases. For the positive case, the frequency difference between neighboring qubits is +5 MHz, and the working frequencies for Q 1 to Q 12 range from 4.3225 GHz to 4.3775 GHz. For the negative case, the frequency difference is −5 MHz. For these two cases, we firstly excite one qubit and then detune all qubits to their working frequencies for system evolution. After that, we detune them back to their idle frequencies for measurement. The population propagation is then measured as a function of time. For each case, we excite the twelve qubits sequentially, and then obtain the corresponding time-dependent population distributions. In the second step, we run a Nelder-Mead optimization to find out the best estimation of the frequency differences. We assume that due to the residual Z cross-talk, the frequency differences of each qubit site for the positive and negative cases are the same. Based on that, by numerically simulating the time evolution of the 12-qubit system, we obtain the population distributions as a function of time. We use the square sum of the differences between simulations and experiments as a cost function for optimization. An example of the comparison between the experimental and simulation results is shown in Fig. S1, in which the different behavior for positive and negative cases, caused by the frequency differences, is presented. In the end of this round of calibration, we correct the frequency alignment by adding the differences to the central frequencies. After several rounds of calibration, the maximum frequency difference is reduced from 15.4 MHz to below 4.9 MHz, smaller than 0.43J 1 /2π.

IV. INITIAL STATES
In this section, we show that the behavior of the quantities that we consider in the main text is general for various initial states. However, each initial state has a different overlap pattern with the eigenstates of the Hamiltonian. Each eigenstate localizes at a different disorder strength depending on its energy eigenvalue. Therefore, each initial state may show a different localization point. This is indeed an evidence for the presence of the mobility edge.
Here, we investigate the role of each initial state on Q (s) eq and ∆Q (s) eq for the three quantities, namely the autocorrelation, the number entropy and Hamming distance. In order to ease the visualization, in the following, we only present the data from numerical simulation, since the experimental data follow the same behavior. In Figs. S2a-c (upper panel), we plot the equilibrium auto-correlation C . We note that, the three quantities smoothly change from their expected value in the ergodic phase to the MBL regime with similar qualitative behavior for all the initial states. The standard deviation of equilibrium values of the three quantities with respect to the R = 60 random realizations as a function of disorder strength h/J1 for various initial states |s = |s k is plotted in the lower panels. The lower panel respectively represents the standard deviation of: a,the auto-correlation ∆C ; and c, the Hamming distance ∆D (s) eq . We note that, for all the three quantities the peak of the standard deviation slightly varies for different initial states, witnessing the existence of the mobility edge. expected, all quantities smoothly change from the ergodic to the MBL phase. While, all the initial states qualitatively follow the same behavior there is small difference between them. To further evidence this, we plot the standard deviation ∆Q (s) eq as a function of disorder strength for all the given initial states |s = |s k in Figs. S2a-c (lower panel). As it is clear from the figure, the peak of ∆Q (s) eq takes place at different disorder values for each initial state. This is a strong evidence for the observation of the mobility edge in our system.
To have a better understanding, one may consider the time evolution of a given initial state |s , which in general can be written as where α and | α are the eigenenergies and the eigenstates of the Hamiltonian H for a particular random instance. The overlap of the quantum state of the system with each of the eigenstates | α is time independent and is given by | α |Ψ(t) | 2 = | α |s | 2 . Due to the mobility edge, for any given random instance r each of the eigenstates | α localizes at a different disorder strength h (α,r) c . Therefore, the transition point for the initial state |s and a random disorder instance r is given by One can average over all random realizations as h

V. DATA PROCESSING AND POST SELECTION OF MEASUREMENT RESULTS
The Bose-Hubbar Hamiltonian H, conserves the total number of excitations in the system. All the initial states that we have considered have the filling factor f = 1/2. Nonetheless, due to experimental imperfections some of the measurement outcomes show a different filling factor which can be attributed to either the decay in transmonic qubits or imperfect readout. In addition, due to dephasing the dynamics may not be fully unitary, which can also induce errors in the results. In order to include these effects and possibly compensate them in our numerical simulations, we take two different approaches; (i) post selecting the experimental data, which means excluding the experimental measurement outcomes with wrong filling factors, and use unitary evolution for the numerical simulation; and (ii) keep the raw experimental data and instead use an open quantum system formulation for the numerical simulations. As mentioned before, the data shown in the main text is based on the first approach. Here, we provide a comparison between the two methods. For open quantum system evolution we use a Lindbladian master equation, whereρ is the density matrix of the system, Γ ( ) and γ ( ) are the decay and dephasing rates of the qubit , respectively. Moreover, the Lindblad term is denoted by The values for decay and dephasing rates are taken from the TABLE I   It is worth emphasizing that neither the post selection of experimental data nor the open quantum simulation can fully compensate the imperfections of the experiments. There are several reasons which contribute to this issue. First, the post selection of experimental data only takes into account the loss in the qubits and the first order readout errors (i.e. when only one qubit is flipped) while it cannot compromise the dephasing effects. Second, the Markovian model used for the open quantum system simulation in Eq. (S3) is very simplistic and ignores crosstalks and correlations between different qubits. Third, the experimental measurements of the device parameters, e.g. couplings, dephasing and decay rates etc., are inevitably prone to errors producing uncertainly in the parameter values used in the numerical simulations.
Numerical simulations. All the numerical simulations were performed using the QUTIP Python toolbox [S5]. In particular, for the time evolution we used the QUTIP mesolve master equation solver and the Hamiltonian parameters were taken from Table S1. The local bosonic Hilbert space is truncated at n max = 3. To see the accuracy of this truncation the commutators Ψ s (t)|[ a , a † ]|Ψ s (t) = 1 ± were computed for all evolved states which showed ≤ 10 −5 . In fact, the good agreement between the theory and experiments in our data shows that the assumption of unitary evolution is valid. The Bose-Hubbard model conserves the total number of excitations, that leads to a link between the particle numbers of the two complimentary subsystems. In the main text, we focused on the half-chain number entropy as a tool to characterize the MBL transition. To see the generality of the behavior of this quantity, in this section, we present the number entropy with different subsystem sizes. In Fig. S4a, we plot the equilibrium number entropy S (m),(s) eq as a function of disorder strength h/J 1 , for different subsystem size m and the initial state |s = |s 1 . As the block size m gets larger the number of possible outcomes for the excitation in the block increases and thus the achievable entropy becomes larger. This can be observed in Fig. S4a, when the system is in the ergodic regime, namely small h/J 1 . In the MBL regime, since the dynamics is almost frozen the number entropy is always small for all subsystems sizes. In Fig. S4b, we display the standard deviation of the equilibrium number entropy ∆S (m),(s) eq as a function of disorder strength h/J 1 for different subsystem size m and the initial state |s = |s 1 . By increasing the system size m, the location of the peak of the standard deviation mostly takes place at around h / J 1 ∼ 2 − 3 showing that for detecting the MBL transition using larger subsystem sizes is more reliable.

VII. CONVERGENCE WITH RESPECT TO THE NUMBER OF RANDOM ENSEMBLE
As mentioned before, in our experiment, we employ 10 different initial states and R = 60 random samples. The role of different initial states has already been discussed and in this section we focus on the choice of the number of random realizations R. We compare the behavior of our quantities for the initial state |s 1 when the number of random ensembles is taken to be R = 60 (as performed in our experiment), R = 500 and R = 1000 for which all quantities are expected to converge, . We find that, all the quantities are indeed converged even for R = 60 random realizations. In Fig. S5a-c upper panel (lower panel) we plot, a, the auto-correlation C ; and c, the Hamming distance ∆D (s 1 ) eq . The good agreement between the curves, specially the location at which the standard deviation peaks, shows that the choice of R = 60 random realizations can faithfully captures the statistical behavior with respect to disorder in the system.  For each quantity measured in our experiment , we start with 10 different initial states. The 10 data points obtained from the peaks are within one groups. The three quantities are independent, and the groups for the three quantities are hierarchical. Therefore, to estimate the total mean, which corresponds to the transition point, we use a multilevel modeling approach [S6], which allows us to consider the group level variance and individual level variance.
The data points Y i,j , where i = C, S, D, and j = 1, ..., 10, are assumed to be independently normally distributed within each group i, Y i,j |θ i ∼ N (θ i , δ 2 ), where θ i is the group mean and δ 2 is the individual level variance. The group means θ i are assumed to follow a normal distribution with mean µ and group level variance τ 2 , θ i ∼ N (µ, τ 2 ). Based on this model, we use Bayesian modeling with MCMC algorithms [S6] to estimate group mean µ. Gibbs sampling is used as the MCMC sampler [S6]. For the model used in this study, the form of posterior distribution is already known, and based on that, the conditional posterior distributions for the three parameters are then generated. The start values are sampled from non-informative priors distributions. After assigning starting points, Gibbs sampler randomly draw samples from conditional posterior distributions. At a time, only one component of the parameters is updated. In this study, we have totally three individual chains. The converge is determined by all three chains when between-chain differences are small enough. For each chain, the total iterations is 1 × 10 6 , and we discard first 8 × 10 5 points. Among the last 2 × 10 5 iterations, we retain the samples every 10 iterations to reduce auto-correlation. The final posterior probability distributions for each parameter is then determined with the retained samples. From the mode of the distributions, we obtain the µ and τ , which corresponds to the estimation value of transition point and