Two-body mobility edge in the Anderson-Hubbard model in three dimensions: Molecular versus scattering states

Most of our quantitative understanding of disorder-induced metal-insulator transitions comes from numerical studies of simple noninteracting tight-binding models, like the Anderson model in three dimensions. An important outstanding problem is the fate of the Anderson transition in the presence of additional Hubbard interactions of strength U between particles. Based on large-scale numerics, we compute the position of the mobility edge for a system of two identical bosons or two fermions with opposite spin components. The resulting phase diagram in the interaction-energy-disorder space possesses a remarkably rich and counterintuitive structure, with multiple metallic and insulating phases. We show that this phenomenon originates from the molecular or scatteringlike nature of the pair states available at given energy E and disorder strength W . The disorder-averaged density of states of the effective model for the pair is also investigated. Finally, we discuss the implications of our results for ongoing research on many-body localization.


I. INTRODUCTION
A central concept in the physics of disordered systems is Anderson localization [1], namely the absence of wave diffusion in certain random media as a result of interference effects between the multiple scattering paths generated by the impurities.To date, this phenomenon has been reported for different kinds of waves, including light waves in diffusive media [2,3] or in disordered photonic crystals [4,5], ultrasound [6], microwaves [7] and atomic matter waves [8,9], to cite a few.
Being an interference effect, Anderson localization crucially depends on the spatial dimension of the system and the underlying symmetries of the associated model, which determines its universality class.In the absence of magnetic fields and spin-orbit couplings, the Hamiltonian of a quantum particle exhibits both time-reversal and spin-rotational symmetries and therefore belongs to the orthogonal class.For an uncorrelated disorder, all wave-functions are then exponentially localized in one and in two dimensions.In three dimensions, however, the energy spectrum contains one or more critical points, called mobility edges, separating localized from extended states.At these points the system undergoes a metalinsulator phase transition, known as Anderson transition [10], which is characterized by universal critical exponents.Mobility edges have been reported [11][12][13] in experiments with non-interacting ultra-cold atoms in threedimensional (3D) speckle potentials.Analogous transition for light waves, despite several claims, have not yet been unambiguously observed, mainly due to the vector character of light [14].
Anderson transitions are difficult to describe analytically and our quantitative understanding relies heavily on numerics.The most studied example of disordered system is a tight-binding model with random on-site energies, known as Anderson model.In first quantization notation, the latter reads where J is the tunneling rate between two nearest neighboring sites n and m, while V n are random variables denoting the local value of the disorder potential.For simplicity, the disorder is assumed to be spatially uncorrelated, V n V n ′ = V 2 n δ nn ′ and to obey a uniform on site distribution, where Θ(x) is the Heaviside function and W is the disorder strength.The position of the mobility edge for the model (1) was first computed in Ref. [15] using transfer matrix techniques.These results, which extended previous work [16] performed for zero energy of the particle, were instrumental to develop approximate semianalytical theories of the Anderson transition, including the self-consistent theory of localization [17][18][19].The Anderson model is currently investigated in three [20] and higher dimensions [21,22] to pinpoint the precise position of the mobility edge and to provide accurate estimates of the universal critical exponents.The same model, but with finite spatial correlations, emerges from the discretization of the Schrödinger equation of a continuum system.In particular, approximating the Laplacian by a second-order finite difference yields Eq. ( 1) with J = 2 /(2m∆ 2 ), where m is the particle mass and ∆ is the lattice spacing.This procedure has recently been applied to obtain precise estimates [23][24][25][26][27][28] for the position of the mobility edge of cold atoms in laser speckle potentials.
A main topic of current research is many-body localization [29][30][31], namely the generalization of Anderson localization to disordered systems of interacting quantum particles.Of particular interest is the existence of many-body mobility edges, namely critical points at finite energy density, separating the many-body localized phase at weak interaction from the metallic, ergodic, phase at strong interaction.Experimental evidence of such critical points has been reported [32][33][34] in experiments with ultra-cold atoms in disordered lattices, implementing either the fermionic or the bosonic Anderson-Hubbard model in various dimensions.Due to the high computational effort, most theoretical studies of systems with a finite density of particles have so far been limited to one-dimensional models [35][36][37][38][39].
A second and complementary approach to interactioninduced Anderson transitions, which is more relevant for the present study, focuses on few-body systems, starting from the solution of the two-particle problem in the presence of disorder.The corresponding Hamiltonian can be written in second quantization as Ĥ = Ĥ0 + Û , where Ĥ0 = Ĥsp ⊗ 1 + 1 ⊗ Ĥsp is the noninteracting part and is the on-site Hubbard interaction of strength U .For 1D systems, the problem of two-particle localization was first addressed by Shepelyansky [40].Using results from random matrix theory, he showed that, in the presence of disorder, two particles coupled via short-range interactions can spread over a distance much larger than the singleparticle localization length, before being ultimately localized.This surprising effect has been confirmed by several numerical studies [41][42][43][44][45][46][47][48][49][50][51] during the last 25 years, although the analytical formula describing the enhancement of the pair localization length at weak disorder is still debated.The localization properties of a system of few (two, three) bosonic atoms in a one-dimensional continuum model laser speckle disorder has been recently addressed [52].
Anderson localization of few interacting photons states in a disordered chain has been discussed theoretically for both linear [47] and nonlinear [53] photonic lattices.Remarkably, signatures of interaction-induced delocalization has been recently observed [54] experimentally in a chain of superconducting qubits simulating the disordered Bose-Hubbard model.In Ref.s [55,56] it was argued that all two-particle states remain localized in one and two dimensions (although the pair localization length can be extremely large), whereas in three dimensions an Anderson transition to a diffusive phase could occur even when all single-particle states are localized.These claims are in stark contrast with subsequent numerical works [57,58], providing evidence of 2D metal-insulator transitions of the pair induced by the Hubbard interactions (although finite-size effects can be an important issue).
Based on large scale numerical simulations, we recently investigated [59] the two-particle problem in three dimensions, focusing on a pair with zero total energy, E = 0. We addressed the localization properties of the system by mapping the original Hamiltonian into an effective single-particle model [see Eq. ( 5) below] describing the center-of-mass motion of the pair, following the lines of Ref. [60].We found that Anderson transitions of the pair The orange and the blue data curves are cuts along the planes E = 0 and E = U .For vanishing interactions, Wc coincides with the single-particle mobility edge calculated in Ref. [15], under the change of variable E = −2ε, where ε is the energy of a single particle.(violet data curve).Panel b) Nature of the pair state as a function of energy and disorder.The two (brown) solid lines define the numerical band edges of the non-interacting two-particle energy spectrum for a given disorder W .They divide the plane in three regions, corresponding to scattering states, attractively and repulsively bound states.For |E| > 12J, the nature of the state changes from molecular to scattering-like as the disorder strength increases (solid vertical arrow), generating multiple Anderson transitions.
are consistent with the orthogonal universality class, although the inclusion of irrelevant variables in the finitesize scaling analysis was crucial to obtain accurate results for the mobility edge.
In Ref. [59] we derived the phase diagram in the interaction-disorder plane for a pair with zero total energy, E = 0.For a given value of the interaction strength U , we found a single critical disorder amplitude W c separating the extended states (W < W c ) from the localized ones (W > W c ).Moreover, we showed that the metalinsulator transition for the pair occurs in a regime where all single-particle states are localized, confirming that interactions favor the delocalization of the pair, irrespective of their attractive or repulsive nature.
The opposite effect, that is interaction-induced localization of the pair, is also possible.Indeed two particles can form attractively or repulsively bound states.For sufficiently strong interactions, so that E ≃ U , these states behave as point-like particles with reduced tunneling rate 2J 2 /|U |.As a consequence, they tend to localize already in the presence of a very weak disorder [60].

Scope of the paper
Building on the results of Ref. [59], in this work we investigate pairs with nonzero total energy and map out the phase boundary between localized and extended states in the interaction-energy-disorder space.This will be done by considering different cuts of the three-dimensional phase diagram along specific planes.Some of these sections are displayed in Fig. 1a.We see that the critical disorder strength along the plane E = U (blue line), exhibits a s-like behavior as a function of the interaction strength, as displayed in Fig. 1b.The double reentrance for intermediate U values signals that in this region the system undergoes three Anderson transitions as W increases, in contrast with the E = 0 case (orange line).
As we shall see, this surprising effect can be explained by the change in the nature of the pair states for increasing disorder.The latter broadens the energy band of a single particle according to −ε be (W ) ≤ ε ≤ ε be (W ), where the (numerical) band edges ±ε be (W ) are computed for a given W value by neglecting Lifshitz tail regions, as explained in Appendix B. As a consequence, the energy spectrum of two noninteracting particles is bound to the interval −2ε be (W ) ≤ ε 1 + ε 2 ≤ 2ε be (W ).The two band edges then divide the energy-disorder plane into three regions, as illustrated in Fig. 1b.If the total energy E of the pair falls inside such interval, the state is scatteringlike, otherwise the state is either attractively bound if E < −2ε be (W ), or repulsively bound if E > 12ε be (W ).
We see from Fig. 1b that for |E| < 12J the pair is described by scattering states for any disorder.The resulting phase diagram for a fixed energy is then fairly similar to the E = 0 case already investigated.For |E| > 12J, however, the nature of the pair states changes from molecular to scattering-like at the disorder threshold W * given by the condition E = ±2ε be (W * ) (vertical arrow).We therefore expect Anderson transitions of molecular states at weak disorder, with W c < W * , followed by a delocalization transition of scattering states at intermediate disorder, with W c > W * .The exploration of these novel metal-insulator transitions of the pair will be the main goal of our paper.
The paper is organized as follows.In Sec.II we review the underlying theoretical approach, which amounts to mapping the two-particle Schrödinger equation into an effective single-particle model with long-range hopping.In Section III we present our numerical results for the two-body mobility edge, while in Section IV we discuss the properties of the (disorder-averaged) density of states of the effective model.Moreover, Section V provides a conclusion and an outlook.In Appendix A we present the numerical method used to derive our results using supercomputers, going from the efficient computation of the matrix K of the effective Hamiltonian to the identification of the critical point of the Anderson transition.In Appendix B we recall the calculation of the single-particle numerical band edge based on the coherent potential approximation.

II. EFFECTIVE SINGLE-PARTICLE MODEL FOR THE PAIR
Hereafter we fix the energy scale by setting J = 1.We address the localization properties of the two-body system via a mapping into an effective single-particle model describing its center-of-mass motion of the pair.The mapping is exact in the subspace of orbitally symmetric wave-functions, describing either bosons or fermions in the spin-singlet state (Hubbard interactions have no effect for identical fermions due to the Pauli principle).
We start by writing the two-particles Schrödinger equation as (E − Ĥ0 )|ψ = Û |ψ , where E is the total energy of the pair.From Eq. (3), we find that the wavefunction obeys the following self-consistent equation where Ĝ(E) = (E Î − Ĥ0 ) −1 is the non-interacting twoparticle Green's function.Eq. ( 4) shows that for contact interactions the wave-function can be completely determined once its diagonal amplitudes f m = m, m|ψ are known.By projecting Eq.( 4) over the state |n, n , we see that these terms obey a closed equation [59][60][61]: where K nm = n, n| Ĝ(E)|m, m .Eq.( 5) can be interpreted as an effective single-particle problem with Hamiltonian matrix K and pseudo-energy λ = 1/U , corresponding to the inverse of the interaction strength.Since K depends explicitly on the total energy, the phase boundary between localized and extended states of the pair will represent a surface in the U − E − W space.
The effective model differs from the Anderson model, Eq. ( 1), in two fundamental aspects.First, the matrix elements of K are unknown and must be calculated numerically.When expressed in terms of the eigen-basis of the single-particle model, Ĥsp |φ r = ε r |φ r , they are given by where φ nr = n|φ r are the amplitudes of the wavefunctions and N is the total number of lattice sites in the grid.Second, the matrix K is fully dense, describing hopping processes between arbitrarily distant sites.The efficient computation of the matrix K, which is the main bottleneck of our approach, is discussed in details in the Appendix A.
For |E| ≫ 1 or W ≫ 1, however, tunneling effects are perturbative and the effective Hamiltonian becomes short-range.To see this, we write the noninteracting twoparticle Green's function as Ĝ(E) = ( Â + T ) −1 , where represents the local part of the Hamiltonian, while accounts for the hopping processes.Next, we expand Ĝ(E) in powers of T , retaining up to second-order terms: Ĝ The second term in the rhs of Eq.( 9) does not contribute to the effective Hamiltonian K, because Â is diagonal in the site basis, whereas T has zero expectation value.The third term contributes through two distinct processes: i) a particle hops from a site to a neighboring one and comes back, while the other does not move; ii) both particles move from one site to the same neighboring site, leading to an effective pair hopping.An explicit calculation yields where δ nm is the Kronecker delta.The rhs of Eq.( 10) defines a tight-binding model for the pair, regarded as a point-like particle.In the absence of tunneling, the matrix K is diagonal, since the two particles can only interact if they share the same lattice site.
The regime |E| ≫ 1, W describes tightly bound states with E ∼ U .In this limit the off-diagonal matrix elements in Eq.( 10) are approximately constant and equal to 1/E 3 .For weak disorder, the effective model ( 5) reduces to [60] 2 showing that tightly bound pairs exhibit a quenched tunneling rate, J b = −2/E, and feel a twice larger disorder strength, W b = 2W (neglecting the small V 2 n correction).Hence we can infer the mobility edge of such states from the known [15] result for the Anderson model.
Eq. ( 10) applies also to the atomic limit, corresponding to W ≫ 1.In this case the short-range nature of the model is ensured by the fact that the amplitudes φ ns of the single-particle wave-functions in Eq.( 6) have support on very few lattice sites.Differently from the molecular regime, the pair tunneling rate cannot be seen as approximately uniform, but depends on the specific values of the disorder potential at the two edges of the bond.In particular, both diagonal and off-diagonal matrix elements of K can take large values when the energy denominators in the rhs of Eq.( 10) become small.

III. TWO-BODY MOBILITY EDGE
Below we mainly investigate two-body systems with total energy E < −12J.The case E > 12J is immedi- The green circles data refer to the pair mobility edge, separating the metallic (M) phase from the insulating (I) one.The green line is a guide to the eye.The dashed lines correspond to the rigorous edges of the interaction band, U = E − W and U = 1/f (E + W ), where f is defined in Eq.( 13); below these lines no pair states exist (grey regions).The dotted lines represent the numerical band edges, which neglect Lifshitz tails, calculated from the disorder-averaged density of states of the effective model.The dot-dashed line corresponds to the mobility edge of point-like molecules, calculated using the numerical data for the single-particle mobility edge from Ref. [15].
ately recovered from our study by using the invariance of the Schrödinger equation under the transformation E → −E, U → −U .The numerical method used to pinpoint the position of the mobility edge of the pair is the same as in Ref. [59] and is briefly discussed in Appendix A.

A. Phase diagrams at fixed energy
We first present our numerical results for a pair with total energy E = −15, focusing initially on the localization properties of the attractively bound states at low disorder.In Fig. 2 we display the calculated boundary between the metallic (M) and the insulating (I) phases (green data points).
In the absence of disorder, the single-particle wavefunctions are plane waves, φ nk = e ikn / √ N , with energy dispersion where k is the lattice momentum.From Eq.( 6) it follows that for E < −12 the solutions of the effective Schrödinger equation ( 5) have the same form, f n = e iQn , where Q is the lattice momentum for the center of mass motion.By direct substitution, one finds For Q = 0 and E < −12, we can calculate the integral in Eq.( 12) analytically, by writing the denominator using the formula 1/x = − +∞ 0 e xt dt, valid for x < 0. This yields λ = f (E), where e Et e 4(cos kx+cos ky +cos kz )t dt with I n (x) being the modified Bessel function of the first kind.For Q = (π, π, π) the integral in Eq.( 12) is also analytical, because ε k = −ε Q−k , and therefore λ = 1/E.Hence for W = 0 molecular states exist for f (E) < λ < 1/E, or equivalently, E < U < 1/f (E).This is evident in Fig. 2 by noticing that 1/f (−15) = −12.995.The dashed curves in Fig. 2 correspond to rigorous band edges of the system, below which no states are allowed, due to energy conservation.To find them, we notice that disorder contributes to the total energy by a term in the interval [−W, W ]. Hence the interaction band of molecular states for E < −12 spreads at most to E − W < U < 1/f (E + W ). Scattering states are instead possible provided that −12 − W < E < 12 + W , independently of the value of the interaction strength.By setting E = −15, this implies that for W ≥ 3 all values of the interaction strength are in principle permitted, whereas for W < 3 only states between the two curves U = E −W and U = 1/f (E + W ) are allowed.
The two dotted lines in Fig. 2 represent the numerical band edges for the pair, calculated from the disorderaveraged density of states of the effective model, as explained in Sec.IV.The regions of the phase diagram between the dotted and the dashed lines correspond to localized states in the Lifshitz-tail regime, where the density of states is very low.
The dot-dashed line in Fig. 2 refers to the mobility edge of the pair regarded as a point-like particle, obeying Eq. (11).It is obtained from the numerical data [15] for the single-particle phase diagram in the ε − W plane, taking into account the rescaled energy ε b = E 2 /U − E − 12/E of the pair as well as the associated hopping rate We see in Fig. 2 that the point-like approximation yields very accurate results for pair states near U = E, but substantially underestimates the size of the metallic phase for weaker interactions.Indeed, such states describe molecules with lower binding energy, so that the corresponding wave-functions can spread over several lattice sites.Fig. 2 shows that the phase diagram is not center-symmetric: the tip is shifted towards the right, showing that weakly bound pairs are more robust against localization than point-like molecules.
We also notice that the point-like approximation misses states at weak interaction, already in the absence of disorder.Indeed, the unperturbed band edges, obtained from the solution of ε b = ±6J b , are given by U = E and U = E 3 /(24 + E 2 ) = −13.55 for E = −15.
We can improve the accuracy of the tight-binding model for pairs by including higher-order tunneling terms in the rhs of Eq.( 8).The third-order term gives zero contribution to the effective Hamiltonian K (like all odd terms), while the fourth-order term gives 3 and a larger pair tunneling rate, J b = −2/E − 120/E 3 = 0.169 for E = −15.Using this last result, the width of the interaction band becomes 12J b = 2.028, in fairly good agreement with our numerics.On the other hand the above fourth-order expansion introduces also second-nearest-neighbor hopping processes, which are not contained in Eq. (11).These and even longer-range hopping terms become more and more important as the energy E increases and the binding energy of the molecule becomes small.
Let us now discuss the localization properties of the pair for stronger disorder.The complete phase diagram for E = −15 is shown in Fig. 3.In Fig. 4 we also display the behavior of the reduced localization length Λ M as a function of the interaction strength, which helps understanding the structure of the phase diagram.The two We see from Fig. 3 that all two-particle states are localized for 1.4 < W < 9.8.In this insulating phase, the region of U values, delimited by the left and right numerical band edges, broadens up as W increases until it covers the entire axis at W = 9.75.Fig. 4 (panels a-c) shows that, for W 9, the two curves for M = 8 and M = 10 tend to further separate out as W increases, as occurs in the single particle problem at strong enough disorder (so that asymptotically Λ M /Λ M ′ = M ′ /M ).This behavior corresponds to localized molecular states.Interestingly, the same panels show that in the Lifshitz tail regions Λ M increases steadily as W increases.
We see from Fig. 4d that for W = 9.4 the two curves show an opposite trend: their distance reduces as W increases, suggesting that the pair has lost its molecular nature, and is better described by a scattering state.This change of behavior should occur when the energy E of the pair falls inside the non-interacting two-particle energy spectrum, as displayed in Fig. 1b.The disorder threshold W * is then given by the condition E = −2ε be (W * ).We compute the single-particle numerical band edge us discussed in Appendix B. The above condition then yields W * = 8.91 for E = −15, confirming the molecule unbinding.Figure 3 shows that, for W = W * (horizontal arrow), the right numerical band edge for the pair crosses the U = 0 axis (corresponding to λ → ∞), as indicated by the star symbol.
We see from Fig. 4e that at W = 9.4 the reduced localization length already possesses a clear absolute mini-mum at U = 0, which then persists for all larger values of the disorder strength, as displayed in the panels f and g of the same figure.This confirms that interactions always favor the delocalization of scattering states.Moreover, by comparing Fig. 4e with Fig. 4f, we see that all scattering states are still localized at W = 9.4, while for W = 10 they are already all extended, except for few states with vanishing interactions.Figure 3 shows indeed that the critical disorder strength is nearly constant, W c ≃ 9.8, with a small bump around U = 0, where W c ≃ 10.5.
The remarkable overlap between the mobility edge and the numerical band edges for strong interactions implies that, in this regime, the rare pair states possess a large mean free path ℓ, as follows from the Ioffe-Regel criterion for the metal-insulator transition, kℓ ∼ 1, k being the effective wave-vector of the pair.
The phase boundary at stronger disorder, where the scattering states ultimately localize, is strongly dependent on the interaction strength, as already observed for the E = 0 case.In particular states with vanishing interaction are the first to localize around W ≃ 14.5, while for |U | 2 the phase transition occurs at much stronger disorder, between W = 23 and W = 24.5.Notice that the metallic phase of scattering states is approximately symmetric under the inversion U → −U .This is also clear from Fig. 4f, showing that the reduced localization length becomes also symmetric under the same tranformation.
Let us now explain how the topology of the phase diagram in the U − W plane is modified by varying the total energy E of the pair.In Fig. 5a we display the results obtained for E = −12.25.In this case the unperturbed band edges are given by U = E and U = 1/f (E) = −8.95.A first striking difference with respect to Fig. 3 is that the two metallic phases of molecular and scattering states are merged together.Interestingly, for U ≥ −8.95 the mobility edge at weak disorder overlaps with the right numerical band edge.In particular the unbinding of molecular states and the subsequent delocalization of scattering states occur almost simultaneously, around W = W * = 2.45.Hence approaching E = −12, all states at low disorder become extended because W * = 0 and the phase diagram for the pair becomes fairly similar to the E = 0 case.
We see from Fig. 5a that the localization of scattering states starts at W ≃ 15.9, for vanishing interactions.Moreover, by comparing it with Fig. 3, we find that the maximum value of the critical disorder strength for the localization of scattering states shifts towards weaker interactions, as the energy E diminishes.
Next, we explore the shape of the phase diagram in the opposite limit, where the energy E of the pair is instead large and negative.In Fig. 5b we show the obtained results for E = −18.In this case the metallic phase of scattering states splits out in two disconnected parts, with support at positive and negative U values, respectively.In other words, there are no metallic pair states for vanishing interactions.With respect to the other values of energy discussed before, finite-size effects are significantly more pronounced and the determination of the phase boundary is therefore less precise.

B. Phase diagram along the E = U plane
We now proceed to discuss the cut of the threedimensional phase diagram of the pair along the E = U plane, which was anticipated in Fig. 1a.The same numerical data are displayed in Fig. 6 (blue circles) together with the previous results for E = 0 (up orange triangles).
While for weak interactions the two data curves re-  Here the two-particle system undergoes three metal insulator transitions as the disorder strength increases, corresponding to localization of molecules, delocalization and subsequent localization of scattering states.These critical points are obtained from Fig. 3 and Fig. 5 (panels a-b) by intersecting the phase boundary with the vertical line at U = E.
It is interesting to note that the critical disorder strengh for the localization of molecules with E = U can be easily computed from the point-like approximation based on Eq. (11).Indeed, from the data of Ref. [15] the single-particle mobility edge at the unperturbed left band edge is W sp c (ε = −6J) ≃ 16J.When expressed in terms of the molecular parameters, this translates into W c ≃ 16/|U |.This is shown in Fig. 6 by the violet double dot-dashed line, which is in very good agreement with our numerics for |U | > 12.

C. Recovering the single-particle mobility edge
A natural question that arises from our discussion is: how the two-body phase diagram in the E − W plane behaves in the limit of vanishing interactions?What is the explicit connexion with the single-particle mobility edge in the ε − W plane? The answer to this question is shown in Fig. 7, where the data symbols correspond to .For E = −18 no critical point is found.The result for E = 0 obtained in Ref. [59] is also shown.The continuos violet line is a guide to the eye connecting the numerical data for the single-particle phase boundary extracted from Ref. [15], upon the change of variable E = 2ε, ε being the single-particle energy.The dashed lines correspond to the rigorous band edges W = −12 ± E, while the dotted lines refer to the numerical band edges of the pair for U = 0 (displayed as solid lines in Fig. 1b).
the critical points at vanishing interactions obtained for E = −15 and E = −12.25 (vertical dashed lines) from the numerical data of Fig. 3 and Fig. 5a (we recall that for E = −18 there are no transitions as U → 0).The corresponding result for E = 0 has also been added.The continuous violet line in Fig. 7 is a guide to eye of the numerical data for the single-particle mobility edge obtained in Ref. [15], expressed in terms of the pair energy E = 2ε.We see that for vanishing interactions, our numerical results for the two-particle mobility edge are fully consistent with the single-particle counterpart (within the numerical accuracy).
Notice that the rigorous and the numerical band edges also agree with the single-particle picture.For instance, the rigorous band edges of the pair for U → 0 are given by the equations −12−W ≤ E ≤ 12+W , which is equivalent to −6 − W/2 ≤ ε ≤ 6 + W/2.The numerical band edge at W = W * , corresponding to the crossing from molecular to scattering states, is fixed by the condition E = ±2ε be (W ), yielding ε = ±ε be (W ), as expected.

IV. DENSITY OF STATES OF THE EFFECTIVE MODEL
While the computation of the transmission amplitude is performed in bar-shaped grids, the DOS can be calcu-lated more accurately using cubic lattices, with L = M , assuming periodic boundary conditions along the three directions.To this end, we compute the matrix K of the effective model with the help of the Woobury matrix identity, as discussed in Appendix A.
The disorder-averaged density of states (hereafter call DOS), expressed as a function of the inverse interaction strength λ = 1/U , is defined as where λ r are the eigenvalues of the kernel K for a given disorder realization and the bar indicates the average over the different disorder realizations.Although this quantity does not show any singular behavior at the critical point of the Anderson transition, it provides useful informations on the distribution of the (pseudo)energy levels which can help us understanding the two-particle phase diagram.
We evaluate the DOS numerically by partitioning the interval [λ min , λ max ], where it is significantly different from zero, into N b bins of equal width ∆λ = |λ max − λ min |/N b .The number of bins used for the evaluation is chosen of the order of the square root of the number of data points per disorder realization, N b ∼ √ N .Let λ j = λ min + ∆λ(j − 1) label the points of the grid, with j = 1, .., N b and let N tr be the total number of disorder realizations considered (in our case N tr = 200).For each bin j and for each disorder realization r, with r = 1, .., N tr , we count the relative number of occurrences p r j , corresponding to the ratio between the number of eigenvalues of the matrix K falling inside the bin and the total number N of eigenvalues.The corresponding value of the DOS is calculated as where the factor ∆λ in the rhs ensures the correct normalization condition, In Fig. 8 we display the DOS of a pair with total energy E = −15 for increasing values of the disorder strength (panel a-f).The vertical arrows mark the position of the numerical band edges, signaling the crossing to the Lifshitz-tail regions.In this work we assume that a given bin j belongs to the Lifshitz tail if the corresponding value of the DOS satisfies where C is a constant of order unity, which for definiteness we choose equal to C = 1/2.The numerical band edges are then obtained as the borders of the region of the λ-spectrum, where Eq.( 16) is satisfied.We have checked that, for the single-particle Anderson model, this working procedure yields results which are consistent with the prediction based on the coherent potential approximation [19].
For very weak disorder (Fig. 8a), the DOS is nonzero only in a narrow region around λ = 1/E = 0.0667, as expected for a tightly bound state.For fixed W , the DOS broadens as the energy E increases, because molecules are less bound, as shown in Fig. 9a for W = 1.The DOS also broadens as the disorder becomes stronger Fig. 8b.This effect is clearly visible in the phase diagram of Fig. 3, where the dotted lines represent the numerical band edge expressed in terms of the interaction strength U = 1/λ.For instance, for E = −15 and W = 7, we see from Fig. 8b that the Lifshitz tails region is given by λ < −0.149 and λ > −0.0485, which translates to −20.64 < U < −6.72.
Entering the scattering regime, at W = W * = 8.91, the support of the DOS becomes unbound, due to the presence of a long-range tail, as shown in Fig. 8c.A power law fit to the tail reveals that the DOS decays algebraically as λ −2 , as displayed in the same panel with the dashed line.This asymptotic behavior signals that the DOS, expressed in terms of the interaction strength as ρK (U ) = ρ K (λ)λ 2 , becomes non zero in the noninteracting limit, ρK (0) = 0; it is therefore a specific feature of the scattering states of the pair.
For stronger disorder, states with repulsive interactions (λ > 0) become also available, as shown in Fig. 8d for W = 24.Differently from the behavior of the reduced localization length (see Fig. 4), the DOS remains strongly asymmetric under a parity transformation λ → −λ, even for rather large values of the disorder strength.This feature can be better understood starting from the atomic limit, where tunneling terms in Eq.( 10) can be neglected, so that the matrix K becomes diagonal, K nm = δ nm /(E−V n ).An explicit calculation then yields where Θ is the unit step function.Eq.( 17) confirms that the DOS behaves as λ −2 , but states with small λ are forbidden due to the energy conservation, yielding W > |E − 1/λ|.An explicit comparison of Eq.( 17) with the full numerical computation of the DOS is shown in Fig. 9b for W = 20 and for three different values of the total energy E of the pair.The dotted lines correspond to the estimate from Eq. (17).We see that for almost all values of λ < 0 the DOS is almost independent of the energy, as expected.The agreement is less good in the strongly interacting regime, corresponding to vanishing λ.Here tunneling effects are important and lead to a finite value of the DOS, ρ K (0) > 0. In contrast, the power-law tails are rather insensitive to such effects, since hopping can always be regarded as perturbative for λ → ∞.
From Fig. 9b we further notice that the DOS becomes more symmetric as the total energy E increases.A full symmetry, however, is recovered only asymptotically for E = 0 [59].

V. CONCLUSION AND OUTLOOK
In this work we have investigated the localization properties of two identical bosons or two fermions with opposite spin moving in a disordered three-dimensional lattice and subject to on-site interactions.The two-body Anderson-Hubbard model provides the simplest example of Anderson transitions in three-dimensional interacting quantum systems.Our theoretical approach is based on an exact mapping of the original Hamiltonian into an effective single-particle model with long-range hopping, whose critical properties are investigated numerical via large-scale simulations (approximately 1.5 million hours of CPU time in state-of-the-art supercomputers).
We found that the two-particle phase diagram in the interaction-energy-disorder space presents an incredibly rich structure characterized by multiple metallic and insulating phases.We showed that this effect originates from the change in the nature of pair states, from molecular to scattering-like, as the disorder strength increases.Our work provides a general framework to study the mobility edge of molecules of arbitrary size, going beyond the point-like approximation holding in the strongly interacting regime.In particular, it allows to describe the behavior of the pair near the dissociation threshold and its subsequent delocalization as a scattering state.Some of our results can readily be tested in current experiments [33] simulating the three-dimensional fermionic Anderson-Hubbard model with atomic gases, by using ultra-diluite samples.These include the observation of interaction-induced delocalization of pairs in regimes where all single-particle states are localized as well as the localization of either attractively or repulsively bound states at low disorder.We hope that our work will help bridging together the field of few-body Anderson localization with its many-body counterpart, at finite particle density.
Finally, our approach can be adapted to investigate the transport properties of other kinds of two-particle systems subject to quenched randomness, like Cooper pairs in strongly disordered atomic gases [62] or superconductors [63][64][65].Investigations of the out-of-equilibrium properties of a Fermi gas undergoing the BCS-BEC crossover in the presence of a random potential [66] are already under way [67][68][69][70].

ACKNOWLEDGEMENTS
We acknowledge fruitful discussions with D. Delande, K. Frahm and S. Skipetrov.
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 665850.This work was granted access to the HPC resources of CINES (Centre Informatique National de l'Enseignement Supérieur) under the allocations 2018-A0040507629 and 2019-A0060507629 supplied by GENCI (Grand Equipement National de Calcul Intensif).In this subsection we outline the numerical procedure followed to efficiently compute the entries of the effective Hamiltonian matrix K for the pair.
We consider a grid of length L and squared transverse section of length M , with L ≤ M .We evaluate the effective Hamiltonian from Eq.( 18), by writing the matrix elements as [43] where G sp (ε) = (εI d − H sp ) −1 is the resolvent of the Anderson model, I d is the identity matrix.Eq. ( 18) shows that the evaluation of the effective Hamiltonian K requires N inversions of N × N matrices, implying that the computational complexity is O(N 4 ).
Fortunately, we can accelerate the calculation of the resolvent exploiting specific properties of the single-particle Hamiltonian, H sp .In the presence of open boundary conditions along the longitudinal direction, the latter possesses a block-tridiagonal structure, each block corresponding to a transverse section of the bar.As a consequence, the resolvent can be written as where A i are M 2 × M 2 symmetric matrices defined by with H sp i being the the Hamiltonian matrix of the i−th block, while 1 and 0 are the identity and the zero matrices, respectively.
Matrices as in Eq.( 19) can be efficiently inverted using recursive techniques [71].To do so, we first compute a sequence of symmetric matrices S i , with i = 1, .., L − 1, using the recurrence relation starting from S L−1 = A −1 L .Let D i and C ij be, respectively, the diagonal and off-diagonal blocks of the matrix G sp that we want to compute (since G sp is symmetric, we can restrict to i > j).These matrices can be determined using the coupled recursive relations for i = 2, .., L − 1, starting from D 1 = A 1 − S 1 .
Using the above procedure, the computational complexity to find G sp reduces to O(L 2 M 6 ), so that the overall cost to evaluate the full kernel K scales with the system size as L 3 M 8 .
Let us now consider the case of periodic boundary conditions along the longitudinal direction.In this case the matrix to invert differs from the rhs of Eq.( 19) by two non vanishing block entries, G sp−1 1L = G sp−1 L1 = 1.Although such a matrix is no longer block-tridiagonal, it can still be inverted efficiently.To see this, we write it as G sp−1 = (B + U t V ) −1 , where B is a block-tridiagonal matrix obtained from the rhs of Eq.( 19) under the change are M 2 × N matrices.After computing the inverse of B using the above procedure, we determine the resolvent G sp via the Woodbury matrix identity: The second term in the rhs of Eq. ( 24) can be calculated using M 2 N 2 elementary operations, which corresponds to the same computational complexity O(M 6 L 2 ) of inverting B. This is consistent with our numerical experiments, showing that, in the presence of periodic boundary conditions along the bar, the time needed to evaluate the effective Hamiltonian approximately doubles.

B. Computation of the critical point
The method followed to extract the position of the mobility edge has been presented in details in [59]; here we briefly outline the main steps.We consider a bar shaped grid, with fixed length L = 150 and transverse size between M = 8 and M = 15, so that L ≫ M .
The logarithm of the transmission amplitude, evaluated at a position n z along the bar, is defined as [16]: where G p (λ) = (λI −K) −1 is the resolvent of the effective model, m ⊥ = (m x , m y ) and n ⊥ = (n x , n y ).We compute the matrix K of the effective Hamiltonian as described in Appendix A. In order to minimize finite-size effects on the transmission amplitude, the boundary conditions on the single-particle Hamiltonian H sp are chosen periodic in the orthogonal directions and open along the transmission axis.For each disorder realization, we evaluate F (n z ) at regular intervals along the bar and apply a linear fit to the data, f f it (n z ) = pn z + q.The Lyapunov exponent is then given by γ M = −p/2, where p is the averaged value of the slope.The critical point W = W c of the metal-insulator transition can be identified by studying the behavior of the reduced localization length Λ M = 1/(γ M M ) for increasing values of the transverse size of the bar.In the metallic phase, Λ M increases as M increases, whereas in the insulating phase it shows an opposite trend.At the critical point Λ M converges to a constant Λ c of order unity, depending on the universality class and the choice of the boundary conditions.In Ref. [59] we show that our numerical results for E = 0 are consistent with the orthogonal universality class, where Λ c = 0.576.This is reasonable, since the effective Hamiltonian K inherits from H sp both the time-reversal and the spin rotational symmetries.
Finite-size effects, drifting the position of the critical point, are however not negligible in our numerics.For this reason, the inclusion of the leading irrelevant variable in the one-parameter scaling ansatz is essential to correctly extrapolate the position of the critical point [59].

VI. APPENDIX B: SINGLE-PARTICLE NUMERICAL BAND EDGE
Neglecting Lifshitz tails, the numerical band edge ε be (W ) for the Anderson model, Eq. ( 1), can be accurately estimated via the coherent potential approximation (CPA) as done in Ref. [19].Here we review the main steps for completeness.We begin by expressing the diagonal term G(ε) = n|(εI − Ĥsp ) −1 |n of the disorder-averaged (translationally invariant) singleparticle Green's function as G(ε) = G 0 (ε − Σ), where is the disorder-free counterpart and Σ is the self-energy.The latter can be found by solving the (self-consistent) CPA equation By substituting the box random potential distribution (2) in Eq.( 27) and performing the integration over the disorder amplitude, we end up with the following equation whose solution yields the self-energy as a function of the single-particle energy and the disorder strength, Σ = Σ(ε, W ). The multi-dimensional integration in Eq.( 26) can be performed analytically following Ref.[72] leading to G 0 (ε) = P (6/ε)/ε, where

FIG. 1 .
FIG.1.Panel a) Critical disorder strength Wc for pair localization as a function of the Hubbard interaction U and the pair total energy E. The orange and the blue data curves are cuts along the planes E = 0 and E = U .For vanishing interactions, Wc coincides with the single-particle mobility edge calculated in Ref.[15], under the change of variable E = −2ε, where ε is the energy of a single particle.(violet data curve).Panel b) Nature of the pair state as a function of energy and disorder.The two (brown) solid lines define the numerical band edges of the non-interacting two-particle energy spectrum for a given disorder W .They divide the plane in three regions, corresponding to scattering states, attractively and repulsively bound states.For |E| > 12J, the nature of the state changes from molecular to scattering-like as the disorder strength increases (solid vertical arrow), generating multiple Anderson transitions.

FIG. 2 .
FIG.2.Zoom of the interaction-disorder phase diagram for a pair with total energy E = −15, in the regime of weak disorder.The green circles data refer to the pair mobility edge, separating the metallic (M) phase from the insulating (I) one.The green line is a guide to the eye.The dashed lines correspond to the rigorous edges of the interaction band, U = E − W and U = 1/f (E + W ), where f is defined in Eq.(13); below these lines no pair states exist (grey regions).The dotted lines represent the numerical band edges, which neglect Lifshitz tails, calculated from the disorder-averaged density of states of the effective model.The dot-dashed line corresponds to the mobility edge of point-like molecules, calculated using the numerical data for the single-particle mobility edge from Ref.[15].

FIG. 3 .
FIG.3.Complete phase diagram in the interaction-disorder plane for a pair with total energy E = −15.The phase boundaries between metallic and insulating phases are displayed by the green symbols.The dashed and dotted lines correspond to the rigorous and the numerical band edges, respectively.The arrow indicates the disorder threshold W * = 8.91, where the nature of the pair wave-function changes from molecular to scattering-like.At this point the right numerical band edge crosses the U = 0 axis, as indicated by the star symbol.For 3 < W < 9.8 (horizontal dashed lines), the pair displays large Lifshitz tail regions.The remaining notation is the same as in Fig.2.

FIG. 5 .
FIG. 5. Topological changes in the phase diagram of the pair for varying energy.Panel a): phase diagram for E = −12.25 showing the two-body mobility edge (orange up triangles) together with the rigorous (dashed-lines) as well as the numerical (dotted lines) band edges.For U ≥ −8.95 the phase boundary at weak disorder basically superposes with the right numerical band edge.The crossing from molecular to scattering states occurs at W = W * = 2.45, as indicated by the star symbol.Panel b): analogous study for E = −18.The two-body mobility edge is displayed by the violet diamonds symbols.The disorder threshold for molecular unbinding is W * = 12.79.This value is slightly smaller than the prediction W * = 13.26 based on the coherent potential approximation, due to finite-size effects.

FIG. 6 .
FIG.6.Phase diagram in the interaction-disorder plane for a pair with total energy E = U (blue circles data).The orange triangles data refer to the phase boundary at E = 0, calculated in Ref.[59].The dashed line corresponds to the molecular result, Wc ≃ 16.0/|U |, obtained by treating the pair as a point-like particle obeying an effective Anderson model, see Eq.(11).The diagram holds for both attractive and repulsive interactions.
main very close (by construction), their behavior in the strongly interacting regime is completely different.For E = U we see that the phase boundary displays a double reentrant (s-like) behavior in the interval 12 < |U | 19.

FIG. 7 .
FIG. 7. Comparison between two-body and single-particle mobility edges for vanishing interactions.The red square symbols denote the two-body data calculated for total energy E = −12.25 (orange dashed line) and E = −15 (green dashed line).For E = −18 no critical point is found.The result for E = 0 obtained in Ref.[59] is also shown.The continuos violet line is a guide to the eye connecting the numerical data for the single-particle phase boundary extracted from Ref.[15], upon the change of variable E = 2ε, ε being the single-particle energy.The dashed lines correspond to the rigorous band edges W = −12 ± E, while the dotted lines refer to the numerical band edges of the pair for U = 0 (displayed as solid lines in Fig.1b).

FIG. 8 .
FIG.8.Disorder-averaged density of states ρK of the effective Hamiltonian for the pair, see Eq. (14), as a function of λ = 1/U .The four panels correspond to increasing values of the disorder strength, while the total energy is fixed to E = −15.The calculation is done assuming a cubic box of sizes L = M = 24 with periodic boundary conditions.The vertical arrows indicate the positions of the numerical band edges, where the Lifshitz tails regions appear.The dashed line in panel (c) corresponds to a power-law fit of the left tail of the data with ρ fit K (λ) = a0λ a 1 yielding a0 = 0.045 ± 0.02 and a1 = −2.03± 0.14.

FIG. 9 .
FIG.9.Disorder-averaged density of states of the effective Hamiltonian for the pair, see Eq. (14), as a function of λ = 1/U , calculated for three different values of its total energy E = −12.25 (orange line), −15 (green line) and −18 (violet line).The left panel (a) corresponds to W = 1, while the right panel (b) refers to W = 20.The grid used is the same as in Fig.8.
APPENDIX A: NUMERICAL METHODA.Evaluation of the matrix K