Many-body localization in the Fock space of natural orbitals

Like as for the first referee, we would express our gratitude to the second referee for the careful reading of the manuscript and for bringing up valuable suggestions for improvements. A revised version (v2) of the manuscript has been submitted. Please see our replies below.

$\textbf{Point 1}$

The first part of the paper describing the basis idea is clearly written. The construction of the Fock state basis using the natural orbitals could be described in more detail. The interpretation of the data in Fig. 1 is a bit unclear. For n=0, it is clear that the participation is tending to 1 in the localized phase. Why doesn’t the same effect occur for n>0? Why does n=1 curve for P^{(n)} appear qualitatively different from the n=0 curve? Even in the thermal phase n=3,4 curves are around 0.5. I would have expected to be smaller. The natural orbitals presumably in the thermal phase are completely delocalized. Is that true? Since, this is an important aspect of the work and probably the subsequent analysis depends crucially on the behaviour of this quantity, it would be suitable if the the data and its interpretation is clarified substantially. Presumably, further clarity can be achieved by studying the finite size scaling of this quantity.

In the revised manuscript, we have elaborated on the construction of the Fock state basis and the interpretation of $P^{(n)}$ in more detail. As pointed out in the revised manuscript, the scaling of $P^{(n)}$ with $n$ depends on the system size in a non-trivial way, such that no simple finite-size scaling analysis is possible for this quantity. With the aim of giving a flavour of the dependence of $P^{(n)}$ on the system size, we have added a plot of $P^{(n)}$ as a function of $W$ for $L=14$ next to the one for $L=16$ in Figure 1 of the revised manuscript.

The question whether the natural orbitals are delocalized in the thermal phase has been studied in Ref. [23] of the revised manuscript. Even though the natural orbitals are far more extended in the thermal phase compared to the many-body localized phase, it is unclear whether this can be seen as true delocalization. We have added the above statement in the revised manuscript.

$\textbf{Point 2}$

The average PR shown in fig. 2 shows the two different regimes. The exponential growth with system size at low disorder in the thermal and the convergence to 1 for all system sizes clearly show the two regimes quite clearly. The authors use this to conclude that the MBL regime sets in at lower disorder. This is likely correct but the analysis used to reach this conclusion is a bit loose. Scaling the data to extract the critical disorder strength would more thoroughly demonstrate this point. The authors should do this for further clarity.

Please see our reply to point 1 of report 1.

$\textbf{Point 3}$

The authors also use the variance of the participation to find the crossover regime. For system sizes 10 and 12, the peak is imperceptible but becomes stronger for system sizes 14 and 16. The delocalised eigenstates satisfy the eigenstate thermalization hypothesis. What is the scaling of the participation ratio expected in the thermal phase based on that? By using the appropriate scaling the location of the critical regime can be better estimated from this quantity.

Unfortunately, we were not able to use scaling arguments to make a better estimate of the location of the crossover when focusing on the variance of the participation ratio. Figure 5 of the revised manuscript shows that the eigenstates acquire a non-trivial structure in a rather wide region around the crossover. As a consequence, scaling arguments based on the behaviour expected when the eigenstate thermalization hypothesis applies can not be used without inducing large uncertainties.

We would like to thank the referee for the careful reading of the manuscript and for bringing up constructive comments. A revised version (v2) of the manuscript has been submitted. Please see our replies below.

$\textbf{Point 1}$

The authors perform an analysis of statistical errors over the ensemble of 1000 single central eigenstates obtained from independent disorder realizations. For this, they split their ensemble in two halves and take the difference in the results of the two subsamples as an estimate for the uncertainty. While this procedure yields a very rough estimate of the error, it suffers from large uncertainties in the error estimate and also potentially overestimates the uncertainty, since the underlying amount of data relies only on half samples. I suggest that the authors change their analysis to use standard resampling methods such as either the jackknife or bootstrap method to obtain more reliable standard error estimates. The such obtained errorbars should be attached in all figures. I expect that these error estimates may explain the nonmonotonic behavior of the PR and in particular var(PR) in Figs 2 and 3, most pronounced for L=16.

In the revised manuscript, we have added error bars to all plots except the one in the first figure (justified below). The errors are estimated by using the jackknife resampling method. Also, we have indicated the number of samples in the captions of the figures.

An error estimate only has a clear interpretation if the sampled data roughly follows a Gaussian distribution. In view of this, we have shifted the focus from the participation ratio to the logarithm of the participation ratio. We have verified that this does not qualitatively alter any of the results. Figure 1 shows averages taken from various qualitatively different distributions, for which we were not able to find an error estimate with a clear interpretation. As the figure serves only an illustrative purpose, we think this is not a major issue.

$\textbf{Point 2}$

The analysis of the evolution with system sizes should be done in a more systematic way. In particular I think that it would be very interesting to use the data from Fig. 2 and extract the scaling with system size in the two bases. Apparently, the PR grows exponentially with the system size, and one could estimate the corresponding scaling exponent by fitting an ansatz of the form PR(L) = aexp(dL). The exponents d appear to be slightly different in the two bases and it looks like d roughly goes to zero at the MBL transition. In the computational basis, the scaling of log(PR) in the MBL phase has been reported to show a logarithmic growth with system size, which might be a difference to the behavior in the Fock basis.

In the revised manuscript, we have included an analysis of the system size scaling of the data shown in Figure 2 of the revised manuscript. Due to the limited range of system sizes and the rather fuzzy system-size dependence over a wide range of disorder strengths, we were unfortunately not able to obtain reliable estimates of the corresponding scaling coefficients.

$\textbf{Point 3}$

It is interesting but not surprising to see that also the PR in the Fock basis shows a peak in its variance in the ergodic phase close to the MBL transition. This is probably for the same reason as the peak in the eigenstate entanglement entropy variance observed in the same parameter regime. If the authors could present the full histograms of the $\text{PR}$ in the Fock basis, it would be interesting to check how the distribution looks like at different disorder strengths.

We have included histograms of the logarithm of $\text{PR}$ for several disorder strengths in Figure 5 of the revised manuscript.

$\textbf{Point 4}$

If the behavior suggested in 3) is indeed observed, it would be fully consistent with the previously observed phenomenology in the distributions of the entanglement entropy, which exhibit bimodal distributions close to the critical point, pointing to a mixture of localized and delocalized states.

We would like to thank the referee for bringing this to our attention. In the revised manuscript, we have added references to Phys. Rev. B 94, 184202 (2016) and Phys. Rev. X 7, 021013 (2017) in the introduction and in Section 3.