Many-body localization of spinless fermions with attractive interactions in one dimension

We study the finite-energy density phase diagram of spinless fermions with attractive interactions in one dimension in the presence of uncorrelated diagonal disorder. Unlike the case of repulsive interactions, a delocalized Luttinger-liquid phase persists at weak disorder in the ground state, which is a well-known result. We revisit the ground-state phase diagram and show that the recently introduced occupation-spectrum discontinuity computed from the eigenspectrum of one-particle density matrices is noticeably smaller in the Luttinger liquid compared to the localized regions. Moreover, we use the functional renormalization scheme to study the finite-size dependence of the conductance, which resolves the existence of the Luttinger liquid as well and is computationally cheap. Our main results concern the finite-energy density case. Using exact diagonalization and by computing various established measures of the many-body localization-delocalization transition, we argue that the zero-temperature Luttinger liquid smoothly evolves into a finite-energy density ergodic phase without any intermediate phase transition.


Introduction
A much studied model for the interplay between disorder and interactions at finite energy densities are spinless fermions in a one-dimensional lattice with uncorrelated diagonal disorder and nearest-neighbor repulsive interactions. The Hamiltonian, defined on a one-dimensional (1D) lattice of L sites with open boundary conditions, reads The hopping constant is set to unity, t = 1, throughout the paper. The onsite potentials are uniformly and randomly distributed over the interval, ε i ∈ [−W /2, W /2]. The disorder strength is denoted by W and the nearest-neighbor interaction strength by V . The emergent picture, supported mostly by numerical studies [1][2][3][4][5][6][7], is that a many-body localized (MBL) phase is stable at sufficiently strong disorder even in the middle of the many-body energy spectrum. As disorder strength is lowered, an ergodic phase is entered. The ergodic and MBL phase are distinct from each other in many ways: volume versus area-law scaling of the entanglement entropy [8,9], thermalization versus failure of eigenstate thermalization hypothesis and Wigner-Dyson versus Poisson statistics of the level spacing distribution (see [10,11] for a review). The reason for these behaviors is the existence of localized quasi-particles (often called l-bits) in the many-body localized phase [12][13][14][15][16][17] that form an extensive set of conserved quantities (see [10,18,19] for a review). Curiously, the conductivity does not behave according to the naive expectation (i.e., vanishing conductivities in the MBL phase and finite ones in the ergodic phase): in the MBL phase, numerical studies report vanishingly small dc conductivities, consistent with a perfect insulator, while in the ergodic phase, some studies report finite dc-conductivities [20,21], while others provide evidence for subdiffusive dynamics [4,[22][23][24][25] (see also [26] for a critical discussion and [27] for a review).
Most of the literature has, with few exceptions (see, e.g., [6]), formally focussed on repulsive interactions (see the discussion at the end of the introduction), while here we explore the attractive regime. In our work, we emphasize an asymmetry in the phase diagram at weak disorder under sign changes of the interaction. This asymmetry is traced back to the zero-temperature properties. While local repulsive interactions cannot overcome disorder in the ground state which is thus fully localized, the situation is markedly different for attractive interactions. Bosonization studies [28] and a subsequent density matrix renormalization group analysis [29] demonstrated the survival of a Luttinger liquid in the attractive case at half filling (see also [30][31][32][33]). This is not surprising since in general, disorder can even enlarge the region in parameter space for superfluid states at commensurate filling factors, as is The Luttinger liquid first transitions into the localized phase (either at an arbitrarily small energy density or at a finite energy density) via an inverted mobility edge, while ultimately entering into the high-energy density delocalized phase.
particularly well known for the case of bosons with repulsive interactions in the Bose-Hubbard model (see, e.g., [34][35][36]). The reason are the disorder-induced density deviations (or, in other words, fluctuations in the local chemical potential) that drive density locally away from commensurability.
It is therefore an interesting question to ask how that zero-temperature Luttinger liquid connects to the finite-energy density phases [37,38]. Primarily, two scenarios are conceivable (see the sketch shown in Fig. 1). At large energy densities, similar to the case of repulsive interactions, we expect an ergodic phase at weak disorder. This high energy-density delocalized phase can either be directly connected to the zero-temperature Luttinger liquid (the case shown in Fig. 1(a)) or it could be intervened by a localized region (the case shown in Fig. 1(b)). The latter scenario would imply an inverted mobility edge.
Such inverted mobility edges were suggested to exist in bosonic systems without disorder [39], and moreover, even a recent experiment with bosons in a two-dimensional lattice [40] could indicate such a behavior, i.e., a delocalization-localization transition as energy density increases. In that experiment, evidence for a many-body localized regime at high energy densities was found, while (minding details of the specific disorder distribution realized in that experiment), theoretical studies predict a superfluid ground state [41], and thus an extended state, for the values of interaction strength and disorder for which the experiment reports localization. Theoretically, there is evidence for the existence of an inverted mobility edge in the one-dimensional Bose-Hubbard model [42,43]. Ref. [43] considered a two-site model with uncorrelated diagonal disorder while Ref. [42] arrives at the same conclusion yet there, the disorder is in the interactions. Clearly, additional studies are necessary to complete this question for the one-dimensional Bose-Hubbard model with disorder as well.
For our model of spinless fermions, we report evidence that the first scenario is realized, i.e., increasing energy density seems to immediately enlarge the delocalized region and there is no phase transition as energy density increases above the Luttinger-liquid phase. Our results are based on an analysis of the entanglement entropy [8,9], the one-particle density matrix occupation spectrum [5,26], and the level-spacing distribution [1].
An asymmetry between repulsive and attractive interactions has also been found in [6], where the limit of strong interactions in a system of spinless fermions with disorder was considered. In the equivalent language of spin-1/2 degrees of freedom, this corresponds to analyzing the effect of disorder on states with ferromagnetic versus antiferromagnetic order at low energies. In the ground state, both phases become localized upon adding disorder. The strongly attractive side was found to lead to a more stable many-body localized phase at finite energy densities compared to the antiferromagnetic case, which shifts the mobility-edge to higher temperatures.
We should stress that while formally, most studies focussed on repulsive interactions, one can use particle-hole symmetry (which is violated in individual disorder realizations but is restored after disorder averaging) to relate the negative V side to the positive V side (see the discussion in [44]). As a consequence, the low-energy density region at V < 0 maps to the high-energy density at the V > 0 side. Using this argument and by inspection of existing studies of the energy-density versus disorder phase diagram (see, e.g., [3,6,26,44,45]), one can already draw some conclusions on the structure of the phase diagram on the attractive side. The special point of V = 2t is the most studied one, for which there are also full energydensity versus disorder strength phase diagrams available (see, e.g., [26]). For the purpose of clarifying the question of an inverted mobility edge on the attractive side, no conclusive picture arises from the existing results for V = 2t: according to Refs. [29,46], V = −2t sits right at the edge of the Luttinger-liquid phase and any finite disorder drives the system into the localized regime in the ground state (see, e.g., the ground-state phase diagram presented in Figs. 3 and 4(a)). Therefore, in order to answer the question of an inverted mobility edge in this model and to complete the analysis of its full interaction strength -disorder phase diagram, a study of the behavior at weak interaction 0 < V < −2t and disorder strengths 0 < W < 2t is necessary. This paper is organized as follows. In Sec. 2, we define the model and the quantities computed in this work. In Sec. 3, we revisit the ground-state phase diagram and show that the occupation-spectrum discontinuity is smaller in the Luttinger-liquid phase, while finitesize dependencies prohibit an extraction of the phase boundaries from that quantity. We further demonstrate that the zero-temperature disorder-interaction phase diagram can also be obtained from a calculation of the conductance using the functional renormalization group method, with qualitative agreement with other methods. Section 4 contains our main results for the phase diagrams at finite energy density and attractive interactions, obtained from exact diagonalization. We conclude with Sec. 5, where we also discuss open questions.

Model and observables, and numerical methods
We consider the spinless-fermion model defined in Eq. (1). Our choice of units is the one employed in [29], while in most of the current MBL literature, the choice of units derives from the equivalent formulation of Eq. (1) in the spin language [2,3], where then ε ∈ [−W , W ] and t → t /2. Hence, in order to compare with the units used in that literature, our numbers have to be divided by a factor of 4. For instance, the critical value for the delocalizationlocalization transition for a Heisenberg chain (V = 2t in our units) is, at T = ∞, given by W ≈ 3.5t [3], hence in our units, at W ≈ 14t.
We further consider the half-filling subspace N = L/2 and obtain results for both attractive and repulsive interactions.
The target energy density is defined as ε = 2(E − E min )/(E max − E min ). Note that a different convention for ε is used in [3,6,44]. The arithmetic average over disorder realizations is denoted as The observables calculated are the following ones. The one-particle density matrix (OPDM) is defined as where |ψ〉 is the given many-body wavefunction. We here compute the OPDM in many-body eigenstates |ψ〉 = |n〉 with H|n〉 = E n |n〉. By diagonalizing the OPDM, one obtains the occupation spectrum n α and the natural orbitals |φ α 〉, which form a complete basis set of single-particle states for each |n〉. After ordering the occupation spectrum, i.e., n 1 ≥ n 2 ≥ · · · ≥ n L , the discontinuity in the occupation spectrum is defined by ∆n = n L/2−1 − n L/2 since we have N = L/2 particles in the system with L even. In a noninteracting fermionic system, any many-body eigenstate can be written as a Slater determinant of single-particle states. In this case, there will be a set of the natural orbitals with n α = 1, spanning the same vector space as the single-particle eigenstates and therefore, the occupation spectrum is a step function (after reordering the eigenvalues [26]). As a result, a discontinuity ∆n = 1 in the occupation spectrum corresponds to product states in Fock space. One refers to this also as localization in Fock space and a large discontinuity with many eigenvalues being close to one or zero implies only a very weak Fock-space delocalization, reminiscent of a zero-temperature Fermi-liquid [26]. The MBL phase has a nonzero occupation-spectrum discontinuity in our model with repulsive interactions, as was shown in [5,26] and also in a system of hardcore bosons in two dimensions [47]. The properties of OPDMs in systems of one-dimensional hardcore bosons on a lattice and in the presence of disorder were studied in [48,49]. Note that the occupation-spectrum discontinuity is smeared out in extensive superpositions of many-body eigenstates even in the MBL phase, while the distribution as such remains highly nonthermal (see the example of a quantum quench discussed in [50]).
To determine whether particles in the MBL phase are truly localized in real space, we also calculate the inverse participation ratio (IPR) of the natural orbitals, which we define as The inverse participation ratio unveils the localization in real space by taking into account the spatial structure of the natural orbitals. We see that in the limiting cases, IPR = 1/L in the ergodic phase and IPR = 1 in the localized phase. Note that unlike in [5], we do not weigh the natural orbitals by their occupations. The analysis of the statistical properties of the natural orbitals across the entire spectrum for a single disorder realization presented in [26] shows that the natural orbitals are, in first approximation, very similar from many-body to many-body eigenstate, implying that mostly just the occupations change.
The von-Neumann entropy is the bipartite entanglement entropy defined as where ρ A is the reduced density matrix obtained by tracing out half of the system. In the localized regime, the entanglement entropy obeys an area-law scaling, while in the ergodic regime, a volume-law scaling is found [8]. The variance of the entanglement entropy is expected to diverge at the transition [9].
The many-body eigenenergy spectra in the localized and ergodic regimes have different distributions [1]. Therefore, a standard way to study many-body localization is to analyze the different statistics of the adjacent level spacing [1,51]. The adjacent gap ratio is defined as where δ (n) = E n −E n−1 is the difference of adjacent many-body energy levels. Theoretical work predicts that the distribution of level spacings themselves follows a Poisson distribution in the Figure 2: Schematic representation of the flow equations for the self-energy and the effective two-particle scattering (an n-particle vertex has 2n external legs).
localized regime and a Gaussian orthogonal ensemble in the ergodic phase [1]. This implies certain values for the average of the gap ratio in these two cases, namely [r gap ] ≈ 0.3863 and [r gap ] ≈ 0.5307, respectively (see the discussion in [1]).

Density matrix renormalization group
In order to compute the occupation-spectrum discontinuity in the ground state, we use density matrix renormalization group simulations [52,53] in a matrix-product state representation using a single-site algorithm with subspace expansion [54]. We use up to 400 states leading to discarded weights of 10 −14 and an average over 5000 disorder realizations.

Exact diagonalization
In order to map out the phase diagram at finite energy densities, we study the model by exact diagonalization for system sizes L = 12, 14, 16, 18 with at least 3000 realizations for each parameter chosen. The exact diagonalization study of many-body localization is limited by both the large dimension of the subspace with N = L/2, e.g., d im(H 18 ) = 48620 and the requirement for a sufficient number of realizations. For all quantities calculated, except for the adjacent gap ratio, we take only the closest eigenvector for the target energy density in each realization. However, we find it necessary to take the 50 closest eigenvalues in one realization to have acceptable statistics for the adjacent gap ratio. We study the model with up to L = 18 sites at various energy densities ε = {0.025, 0.05, . . . , 1} with the shift-invert spectral transformation H → (H − λI) −1 , where λ is the targeted eigenvalue. The libraries PETSc and SLEPc [55] are used for the problem setup and building the solution. We use an exact shift-and-invert solver. The inverse problem involved is solved by direct linear solver with external package MUMPS [56,57] using parallel Cholesky factorization.

Functional renormalization group
The functional renormalization group (FRG) is one implementation of Wilson's general RG idea for interacting many-particle systems [58]. It can be set up either on the Keldysh contour or the Matsubara axis; since we are interested in ground-state properties, the latter choice is the more convenient one. The starting point of an FRG calculation is to take the noninteracting Green's function G 0 of the system under consideration and to cut it off below an infrared energy scale Λ. In particular, we introduce a multiplicative cutoff in Matsubara frequency space, and consider the flow of many-particle vertex functions (such as the self-energy or the effective interaction) as a function of Λ. The resulting (infinite) set of coupled flow equations can be represented elegantly using Feynman-like diagrams (see Fig. 2). Subsequent re-integration from Λ = ∞ down to the cutoff-free system Λ = 0 amounts to an exact solution of the manyparticle problem. In practice, the infinite hierarchy of flow equations needs to be truncated, rendering the FRG an approximate method which treats interactions perturbatively.
The simplest truncation scheme is to only consider the flow of the single-particle vertex (the self-energy) and to neglect the flow of all higher-order vertex functions. If the two-particle vertex is set to the bare interaction V , the flow of the self-energy associated with Eq. (1) can be expressed simply in terms of effective hopping-matrix elements t Λ l and on-site energies ε Λ l : whereG Λ (iω) denotes the flowing single-particle Matsubara Green's function (an explicit expression can be found in Eq. (8)). At half filling, the initial conditions are given by ε Λ→∞ l = ε l as well as t Λ→∞ l = t. We explicitly allow for a spatial dependence of the two-particle interaction V l in order to model a smooth coupling to free fermion source and drain leads (see below). In the flow equations for ε Λ 1 and ε Λ L , terms that formally involve V 0 and V L are set to zero.
The approximation introduced above is strictly correct only to leading order in the interaction V but contains an infinite resummation of Feynman diagrams (since the self-energy feeds back into its own flow). Similarly, one can obtain a second-order truncation scheme by accounting for the flow of both the two-particle vertex and the self-energy while setting all higher-order vertices to their initial value (zero). In this paper, we partially incorporate second-order contributions by parameterizing the two-particle vertex in terms of flowing effective onsite interactions V Λ l . This is a purely pragmatic approach -the resulting approximation is still strictly controlled only to first order -which improves our results quantitatively but does not change them qualitatively. The corresponding set of flow equations is given by Eq. (7) with V l → V Λ l and complemented by a flow equation for the effective interaction V Λ l , which we do not write down explicitly; it can be found in Ref. [59].
In order to compute the conductance, we couple the system to left and right (source and drain) leads, which, for reasons of simplicity, we model as structureless Fermi liquids (e.g., we take the wide-band limit). Such a free-fermion system can be 'projected out' analytically via equation-of-motion techniques and the calculation ofG Λ (iω) reduces to the inversion of a L × L matrix defined by Due to the tridiagonal structure, this inversion can be carried out with a computational effort scaling linearly with L. The flow equations (7) can be integrated using standard Runge-Kutta routines. Finally, one obtains the conductance g (in units of e 2 /h = 1) from More details on the FRG can be found, e.g., in Refs. [58,59].

Ground-state properties
The existence of a Luttinger-liquid phase in the presence of disorder in the spin-1/2 XXZ chain is a much studied problem (which is equivalent to one of spinless fermions or hardcore bosons). While early work established the existence of such a phase [29][30][31][32][33][60][61][62], there is an ongoing discussion on the nature of the transition between the delocalized superfluid and the localized phase (which, in the language of bosons, is a Bose-glass phase [63]). This question is not at the focus of our work and we refer the reader to the pertinent literature for details [46,[64][65][66][67][68][69][70][71][72].

Occupation-spectrum discontinuity
An early DMRG study used the sensitivity to twisted boundary conditions to locate the transition between localized and extended phase. [29]. Entanglement measures were used to locate this transition in [73,74] and it was pointed out [74] that large system sizes are needed to be in the correct scaling regime, due to the large localization length.
Here, we show the phase diagram from Ref. [29] in the disorder-interaction strength plane in Fig. 3. The figure shows the value of the occupation-spectrum discontinuity [∆n] computed for L = 32 sites using DMRG simulations. Remarkably, the region where [∆n] is small (compared to large W or small |V |/t) coincides with the Luttinger liquid (indicated by the symbols). We attempted an extrapolation of [∆n] in the system size to locate the transition yet obtained no conclusive results, which we attribute to the small system sizes L ≤ 128 considered.

Conductance from FRG
It is known [28] that for weak disorder the localization length diverges as ξ ∼ W −2/(3−2K) as one approaches V = −t from above (where the Luttinger-liquid parameters takes the value K = 3/2). This motivates us to study the problem using a different method -the FRG -whose strengths and shortcomings are orthogonal to those of an exact-diagonalization approach.
The key disadvantage of the FRG is that it is approximate w.r.t. the two-particle interaction V . Within our truncation scheme, all results are strictly correct only to leading order in V ; higher-order contributions are uncontrolled. Thus, the FRG is not a suitable tool to determine the precise position of a phase boundary which is situated away from |V |/t 1. The advantage of the FRG is that one can easily treat large systems of L = 10 6 sites and that leads as well as single-particle disorder can be incorporated exactly. Hence, the FRG can overcome the obstacle of a potentially large localization length [74] and provide further evidence for the existence of the delocalized phase in certain parameter regimes.
We use the FRG to compute the disorder-averaged conductance g in the ground state. To this end, one couples the system to left and right Fermi-liquid leads; in order to avoid backscattering at the interface between the noninteracting leads and the interacting region, we employ a spatially smoothened transition region of about 20 sites. The metallic and localized phases are then characterized by g ∼ const. and g ∼ e −L at large L, respectively. In Fig. 4(a), we show g in the V − W parameter plane for a fixed large system size L = 10 6 , which qualitatively reproduces Fig. 3. The interaction-dependence of the conductance for a fixed disorder strength of W /t = 0.5 is displayed in Fig. 4(b) for three different system sizes ranging from L = 10 4 to L = 10 6 . This illustrates that g becomes independent of L for −2t V −t and thus provides further evidence for the existence of the metallic phase in this regime. Note that our FRG scheme apparently misses the localization at large attractive V ∼ −2t (see below). It is surprising that our simple-minded 'Hartree-Fock-like' FRG approach captures the transition between the delocalized and the metallic phase despite the fact that this transition does not occur for |V |/t 1 where the approximation is controlled. This observation motivates an extension of the FRG to a second-order scheme, which would allow one to access finiteenergy properties for systems of up to L = 10 2 sites. Such an investigation is currently under way; preliminary results show that one can capture the phase-separation transition for large negative V that is missed by the leading-order scheme.

Finite energy densities
We now turn to our main case of interest, finite energy densities ε > 0. Our analysis of the energy-density versus interaction-strength phase diagram is based on calculating the impurity averaged one-particle density matrix occupation-spectrum discontinuity [∆n] [5], the (halfcut) von-Neumann entropy [S vN ] [8], its variance var S vN [9] and the level-spacing distribution [1,2]. Our main goal is to illustrate the quantitative and qualitative differences between attractive and repulsive interactions at weak disorder and low energy densities.

Occupation-spectrum discontinuity
The occupation-spectrum discontinuity [∆n] is plotted in Figs. 5(a)-(d) in the disorder-versusinteraction strength plane for different values of the energy density ε = 0, 0.2, 0.4, 1. Note that the occupation-spectrum discontinuity at V = 2t was already computed for all energy densities in [26]. Figure 5(a) qualitatively reproduces, albeit for much smaller systems, the behavior in the ground state discussed in Sec. 3: in the vicinity of the Luttinger-liquid phase at −2 V /t −1.5, the occupation-spectrum discontinuity is markedly smaller than anywhere else in the phase diagram and the V → −V asymmetry is apparent.
Upon increasing energy density, at both V > 0 and V < 0, a region with small [∆n] emerges at small values of the disorder strength and quickly grows in size. For negative V < 0, this presumably ergodic region appears to be adiabatically connected to the zero energy-density Luttinger liquid, as suggested by the sequence of results for increasing ε presented in Figs. 5(a)-(d). There is also a regime with a reduced occupation-spectrum discontinuity at V ∼ 3t on the repulsive side and weak disorder. This one, however, does not become smaller as system size increases, contrary to the behavior in the Luttinger-liquid regime at negative values of V (data not shown here). We speculate that this reduction in [∆n] for large V and small disorder is inherited from the ground-state degeneracy at W = 0 and V t in the density-wave phase. It is instructive to consider cuts through the ε versus V phase diagrams at constant W . Examples are shown in Figs. 6(a) and (b) for W = 1.2t and W = 10.2t, respectively, for several values of the energy density ε. Obviously, [∆n] = 1 at V = 0, independently of the value of disorder. For ε = 0, the discontinuity is large throughout for all values of V . Increasing ε leads to a quick decay of [∆n] away from V = 0, which is more significant for weak disorder, where we expect an ergodic phase everywhere except for the Anderson insulator at V = 0 and its immediate vicinity. Note the asymmetry with respect to V → −V : at small W /t, the many-body states have a smaller discontinuity for V < 0 compared to V > 0 while the trend is opposite for strong disorder (this applies to small energy densities). While the asymmetry at small W /t is a consequence of the ground-state phase diagram, the (inverted) asymmetry between V and −V at large W /t can be understood based on the arguments given in [6]. Their result is that ferromagnetic states should be more susceptible to localization than antiferromagnetic ones, which is reflected in the dependence of [∆n] on V /t at small energy densities but large W /t (compare Fig. 5(b)). Clearly, due to the small system sizes accessible to exact diagonalization, there is need to check for the robustness of the qualitative trends as L is varied. We present such an analysis in Fig. 7 for parameters in the middle of the Luttinger-liquid phase (V = −1.4t) but for an elevated energy density ε = 0.4 and 2 ≤ W /t ≤ 20. At least for the smallest values of W /t = 2, 4, the data clearly suggest that [∆n] → 0, indicative of a phase with Fock-space delocalization and a continuous occupation spectrum as expected for the ergodic phase. For W /t 10, the data seem to extrapolate to a finite value, as expected for the MBL phase [5]. These finite-size trends clearly establish that there are two different phases separated by a Fock-space delocalization transition.

Natural orbitals and inverse participation ratio
The observation of a large occupation-spectrum discontinuity implies Fock-space localization [5] but not necessarily localization in real space. To demonstrate that in the MBL case, both features come about simultaneously, we can analyze the eigenstates of the OPDM as well, the natural orbitals [5]. Figure 8 shows the full distribution of the inverse participation ratio defined in Eq. (4) in its main panels for (a) the ergodic phase and (b) the MBL phase, both at attractive interactions. Before discussing the IPR, it is illustrative to plot individual natural orbitals, which are shown in the insets: the one in the ergodic phase appears to be extended, while the one for the MBL phase is obviously strongly localized in real space. The full analysis of the distributions P(IPR) of the IPR supports this picture. In the ergodic phase, the typical value of the IPR moves to zero with 1/L as L increases and the distributions become narrower at the same time. Conversely, in the MBL case, the distributions are broad, peaked around a large IPR value and practically L-independent. This is in agreement with the results reported in [5,26] for repulsive interactions and immediately visualizes the localization of the quasi-particles (or l-bits) in real space by approximating their single-particle content via diagonalizing the OPDM [26].

Von-Neumann entropy
As an alternative measure of the transition, we analyze the von-Neumann entropy, which we expect to exhibit a volume law in the ergodic phase but an area law in the MBL phase. We verify this expectation by plotting the half-cut entanglement entropy versus system size in Fig. 9 for two points deep in the ergodic and deep in the MBL region.
The location of the transition can be estimated from the position of the maximum of sample-to-sample fluctuations of the von-Neumann entropy, var S vN [9]. We plot this quantity in Fig. 10 for the same parameters as in Figs. 5(a)-(d). The behavior of var S vN confirms the picture obtained from the occupation-spectrum discontinuity, at least for ε > 0. As energy density increases, the transition between the weak-disorder delocalized and the strong-disorder localized phase moves to larger values of W /t. In the ground state, i.e., at ε = 0, and for small system sizes, only the phase boundary at large negative values of V is resolved by a clear maximum in var S vN , while this is not the case for the transition at small negative values. We believe this to be a consequence of (i) the large single-particle localization length on the attractive side and (ii) the different nature of the Luttinger-liquid to localization transition at ε = 0 (see the discussion in [46,67]). By comparison of Fig. 10(d) and the corresponding data for the occupation-spectrum discontinuity [∆n] shown in Fig. 5(d), one realizes that at a fixed system size var S vN peaks at significantly lower values than where [∆n] exhibits a drop to (approximately) zero. First, this simply reflects that the actual transition point can only be estimated from a finite-size scaling analysis (see the discussion in [75] though). Second, different quantities exhibit different finite-size deviations from large-L behavior. The fact that var S vN has a stronger finite-size dependence with the maximum sitting well below the actual transition point is consistent with the observation of other studies (see, e.g., [76]).

Level-spacing distribution
An example for the W -dependence of the adjacent-gap ratio r gap is presented in Fig. 11

Energy-density versus interaction-strength phase diagram
The main result of our work is an energy-density versus interaction-strength phase diagram at weak disorder. We present such diagrams for W = t in Figs. 12(a)-(c), derived from the occupation-spectrum discontinuity, the half-cut von-Neumann entropy, and the adjacent-gap ratio r gap , respectively. All three quantities give a qualitatively consistent picture: at negative V < 0, the ergodic phase extends down to ε = 0 and thus adiabatically connects to the Luttinger-liquid phase, while for repulsive interactions, there is always a mobility edge as ε increases. Moreover, we can thus rule out the presence of an inverted mobility edge from the data for L ≤ 18 sites in our model. The phase diagram further illustrates the asymmetry under changing the sign of V . We here present results for a fixed system size L = 18, while in principle, one can carry out a finite-size analysis and even attempt a finite-size scaling collapse of the data [3,9]. Such analyses were carried out for V > 0 but the resulting exponents are inconsistent with rigorous bounds [77][78][79]. Hence, we here do not pursue this strategy since our system sizes are not larger than what was used in [3,9]. Moreover, the question of what is the universality class of the delocalization-localization transition in models such as ours is still a topic of ongoing research [9,75,76,[80][81][82][83] By inspection of Figs. 12(a)-(c) one realizes that the occupation-spectrum discontinuity and the adjacent-gap ratio provide the best resolution of the phase diagram as both quantities have small finite-size effects deep in the respective phase. Note that we plot the adjacent-gap ratio only down to ε ≥ 0.025t. The statistics of the adjacent gap ratio does not obey the predictions [1,2] as ε → 0, where finite-size effects are expected to be the largest. The delocalized phase is not that obvious in the ground-state region for this system size, due to the comparably much smaller discontinuities at finite energy densities. Note that (c) does not include the ground-state region and hence does not show the mobility edge on the V > 0 side.

Conclusion
In this work we studied the finite energy-density phase diagram of spinless fermions in a onedimensional lattice with attractive interactions and uncorrelated disorder. Our numerical results illustrate an asymmetry between the V < 0 and V > 0 phase diagrams, rooted in the existence of a delocalized phase in the ground state of the system with attractive interactions and weak disorder [28,29]. This Luttinger-liquid phase connects to a finite-temperature ergodic phase without any intermediate phase transition, as the analysis of the entanglement entropy, the level-spacing distribution and the one-particle density matrix occupation spectrum indicates. This conclusion rules out the existence of an inverted mobility edge in our model, which may still exist in other disordered systems with superfluids present in the ground state, such as the Bose-Hubbard model. In fact, a series of very recent studies [42,43,84] indicates that there is an inverted mobility edge in the one-dimensional Bose-Hubbard model. This different behavior compared to our case of fermions with attractive interactions could be traced back to two observations. First, in the case of bosons, disorder makes the stability region of the superfluid (as a function of W and interaction strength) larger than in the clean case [70]. Second, the Bose-Hubbard model is prone to a sort of dynamical localization even in the absence of disorder: composite objects of multiple bosons themselves can get localized in individual sites and may have very large decay times in the limit of strong interactions [85][86][87]. Increasing energy density has been suggested as one mechanism to induce such heavy objects, which in the presence of disorder would get localized.
We also revisited the ground-state case, with two main observations. First, the Luttingerliquid region features a significantly smaller OPDM occupation-spectrum discontinuity compared to the localized ground-state phases on finite systems, as expected for a Luttinger liquid, which has a vanishing discontinuity in the thermodynamic limit. The localized regions, on the other hand, have a Fock-space localized ground-state wavefunction, just as in the finite-energy density case [5,26]. Second, we used functional renormalization group simulations to compute the conductance at zero temperature, which also resolves the delocalized ground-state region. While the FRG simulations for the conductance are computationally cheap and can be pushed to system sizes as large as L ∼ 10 6 , the phase boundaries, in particular, at large negative interaction strengths V −2t, are quantitatively different from the results of other methods (see, e.g., [46]). Thus, improvements in the FRG scheme are necessary to render the method accurate at larger interactions strengths as well, which poses an interesting direction for further method development.