# Two-body metal-insulator transitions in the Anderson-Hubbard model

F Stellin, G Orso

► To cite this version:

F Stellin, G Orso. Two-body metal-insulator transitions in the Anderson-Hubbard model. International Conference on Computer Simulation in Physics and beyond (CSP 2020), Oct 2020, MOSCOU, Russia. pp.012044, 10.1088/1742-6596/1740/1/012044. [hal-04551555]

HAL Id: hal-04551555

https://hal.science/hal-04551555

Submitted on 6 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Two-body metal-insulator transitions in the Anderson-Hubbard model

To cite this article: F Stellin and G Orso 2021 *J. Phys.: Conf. Ser.* **1740** 012044

View the [article online](#) for updates and enhancements.

You may also like

- Bipolar Room Temperature Ferromagnetic Semiconductor LaMnOP

Eiji Motomitsu, Masahiro Hirano, Hiroshi Yanagi et al.

- [Editorial](#)

- Violation of Leggett inequalities in orbital angular momentum subspaces

J Romero, J Leach, B Jack et al.

<figure><img src="image_1.png" /></figure>

The
Electrochemical
Society

Advancing solid state &
electrochemical science & technology

<figure><img src="image_2.png" /></figure>

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_3.png" /></figure>
PAPER • OPEN ACCESS

# Two-body metal-insulator transitions in the Anderson-Hubbard model

To cite this article: F Stellin and G Orso 2021 *J. Phys.: Conf. Ser.* **1740** 012044

View the [article online](#) for updates and enhancements.

You may also like

- Bipolar Room Temperature Ferromagnetic Semiconductor LaMnOP

Eiji Motomitsu, Masahiro Hirano, Hiroshi Yanagi et al.

- [Editorial](#)

- Violation of Leggett inequalities in orbital angular momentum subspaces

J Romero, J Leach, B Jack et al.

<figure><img src="image_2.png" /></figure>

The
Electrochemical
Society

Advancing solid state &
electrochemical science & technology

<figure><img src="image_3.png" /></figure>

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_4.png" /></figure>
# Two-body metal-insulator transitions in the Anderson-Hubbard model

F Stellin¹ and G Orso²

¹,² Université de Paris, Laboratoire Matériaux et Phénomènes Quantiques, CNRS, F-75013, Paris, France

¹ Université Paris-Saclay, ENS Paris-Saclay, CNRS, Centre Borelli, F-91190, Gif-sur-Yvette, France

E-mail: filippo.stellin@univ-paris-diderot.fr, giuliano.orso@univ-paris-diderot.fr

**Abstract.** We review our recent results on Anderson localization in systems of two interacting particles coupled by contact interactions. Based on an exact mapping to an effective single-particle problem, we numerically investigate the occurrence of metal-insulator phase transitions for the pair in two- (2D) and three-dimensional (3D) disordered lattices. In two dimensions, we find that interactions cause an exponential enhancement of the pair localization length with respect to its single-particle counterpart, but do not induce a delocalization transition. In particular we show that previous claims of 2D interaction-induced Anderson transitions are the results of strong finite-size effects. In three dimensions we find that the pair undergoes a metal-insulator transition belonging to the same (orthogonal) universality class of the noninteracting model. We then explore the phase diagram in the space of energy *E*, disorder *W* and interaction strength *U*, which reveals a rich and counterintuitive structure, endowed with multiple metallic and insulating phases. We point out that this phenomenon originates from the molecular and scattering-like nature of the pair states available at given energy and disorder strength.

## 1. Introduction

The propagation of waves in random media can be suppressed when a certain degree of disorder is exceeded, due to the interference between the multiple scattering paths generated by the impurities. This phenomenon, known as Anderson localization, was initially predicted for spin excitations in semiconductors [1] and has been reported for different kinds of waves, such as light waves [2], ultrasound [3], microwaves [4] and atomic matter waves [5, 6, 7]. Its occurrence is ruled by either spatial dimension and symmetries of the model. For example, when both time-reversal and spin-rotation invariances are preserved, all wave-functions are exponentially localized in one and two spatial dimensions. In three dimensions, the energy spectrum contains both regions of localized and of extended states, which are separated by one or more energy thresholds, dubbed mobility edges [8].

The presence of interactions can substantially modify the one-particle scenario [9], as it was confirmed since the earliest studies on interacting electrons in disordered metals [10], among which the scaling theory proposed by Finkelstein [11]. The puzzling evidence of a metal-insulator transition observed in two-dimensional silicon metal-oxide-semiconductor field-effect transistors [12] has also fostered attention on the field, with further theoretical explanations [13] and numerical proofs [14, 15]. In more recent years, the topic of many-body localization, namely the transposition of the concept of Anderson localization to the Fock space for interacting disordered systems with a finite density of particles, has attracted an increasing interest [16, 17, 18]. While in one-dimensional systems the existence of a many-body mobility edge between non-thermalized states and ergodic ones has been proved for particles subject to short-range interactions, its fate in higher-dimensional systems is still unclear [19, 20].

A complementary viewpoint on interaction-induced delocalization transitions focuses on the transport properties of few-body systems in lattices of infinite size, a problem which is computationally more accessible, especially in dimensions higher than one. This problem was first addressed thirty years ago in the seminal works of Dorokhov [21] and Shepelyansky [22], who showed that, in one dimension, interactions can cause a sizable increase of the pair localization length, with respect to the single-particle one. Since then, the effect has been thoroughly investigated in the presence of long- [14] and short-range interactions [23, 24, 25, 26]. In two dimensions, however, very few numerical studies have been pursued so far [27, 28, 29], as the large computational cost severely limits the system sizes that can be attained. Interestingly, it was claimed [27, 28] that interactions can induce a pair delocalization transition, in

<figure><img src="image_4.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
# Two-body metal-insulator transitions in the Anderson-Hubbard model

F Stellin¹ and G Orso²

¹,² Université de Paris, Laboratoire Matériaux et Phénomènes Quantiques, CNRS, F-75013, Paris, France

¹ Université Paris-Saclay, ENS Paris-Saclay, CNRS, Centre Borelli, F-91190, Gif-sur-Yvette, France

E-mail: filippo.stellin@univ-paris-diderot.fr, giuliano.orso@univ-paris-diderot.fr

**Abstract.** We review our recent results on Anderson localization in systems of two interacting particles coupled by contact interactions. Based on an exact mapping to an effective single-particle problem, we numerically investigate the occurrence of metal-insulator phase transitions for the pair in two- (2D) and three-dimensional (3D) disordered lattices. In two dimensions, we find that interactions cause an exponential enhancement of the pair localization length with respect to its single-particle counterpart, but do not induce a delocalization transition. In particular we show that previous claims of 2D interaction-induced Anderson transitions are the results of strong finite-size effects. In three dimensions we find that the pair undergoes a metal-insulator transition belonging to the same (orthogonal) universality class of the noninteracting model. We then explore the phase diagram in the space of energy *E*, disorder *W* and interaction strength *U*, which reveals a rich and counterintuitive structure, endowed with multiple metallic and insulating phases. We point out that this phenomenon originates from the molecular and scattering-like nature of the pair states available at given energy and disorder strength.

## 1. Introduction

The propagation of waves in random media can be suppressed when a certain degree of disorder is exceeded, due to the interference between the multiple scattering paths generated by the impurities. This phenomenon, known as Anderson localization, was initially predicted for spin excitations in semiconductors [1] and has been reported for different kinds of waves, such as light waves [2], ultrasound [3], microwaves [4] and atomic matter waves [5, 6, 7]. Its occurrence is ruled by either spatial dimension and symmetries of the model. For example, when both time-reversal and spin-rotation invariances are preserved, all wave-functions are exponentially localized in one and two spatial dimensions. In three dimensions, the energy spectrum contains both regions of localized and of extended states, which are separated by one or more energy thresholds, dubbed mobility edges [8].

The presence of interactions can substantially modify the one-particle scenario [9], as it was confirmed since the earliest studies on interacting electrons in disordered metals [10], among which the scaling theory proposed by Finkelstein [11]. The puzzling evidence of a metal-insulator transition observed in two-dimensional silicon metal-oxide-semiconductor field-effect transistors [12] has also fostered attention on the field, with further theoretical explanations [13] and numerical proofs [14, 15]. In more recent years, the topic of many-body localization, namely the transposition of the concept of Anderson localization to the Fock space for interacting disordered systems with a finite density of particles, has attracted an increasing interest [16, 17, 18]. While in one-dimensional systems the existence of a many-body mobility edge between non-thermalized states and ergodic ones has been proved for particles subject to short-range interactions, its fate in higher-dimensional systems is still unclear [19, 20].

A complementary viewpoint on interaction-induced delocalization transitions focuses on the transport properties of few-body systems in lattices of infinite size, a problem which is computationally more accessible, especially in dimensions higher than one. This problem was first addressed thirty years ago in the seminal works of Dorokhov [21] and Shepelyansky [22], who showed that, in one dimension, interactions can cause a sizable increase of the pair localization length, with respect to the single-particle one. Since then, the effect has been thoroughly investigated in the presence of long- [14] and short-range interactions [23, 24, 25, 26]. In two dimensions, however, very few numerical studies have been pursued so far [27, 28, 29], as the large computational cost severely limits the system sizes that can be attained. Interestingly, it was claimed [27, 28] that interactions can induce a pair delocalization transition, in

<figure><img src="image_5.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
contrast to previous theoretical arguments [30] suggesting that the localization length of the two-particle system remains always finite. The exploration of the 3D version of the same problem has remained unattempted so far.

In this proceeding article, we review our recent numerical works [31, 32, 33] on two-body Anderson localization in two and in three spatial dimensions, by highlighting the most important results.

## 2. Model and method

We consider a system of two particles moving in a disordered lattice and coupled by contact (Hubbard) interactions. The two particles can be bosons or fermions in the spin-singlet state, so that the orbital part of the wavefunction is symmetric under particles exchange and is therefore sensitive to the interactions. The Hamiltonian of the two-particle system in the lattice basis {$|i, j\rangle$} is given by:

$$ \hat{H} = -J \sum_{\{i,k\},j} |i,j\rangle\langle k,j| - J \sum_{i,\{j,k\}} |i,j\rangle\langle i,k| + \sum_{i,j} (v_i + v_j)|i,j\rangle\langle i,j| + U \sum_i |i, i\rangle\langle i, i|, \quad (1) $$

where J is the tunnelling amplitude between nearest-neighbour sites, U is the Hubbard interaction strength and $v_i$ represents the value of the disorder potential at site i. As in Anderson's original work [1], we assume that the random energies are uncorrelated, $\overline{v_i v_j} = v_i^2 \delta_{i,j}$ and obey a uniform on-site distribution:

$$ P(v_i) = \frac{1}{W} \Theta\left(\frac{W}{2} - |v_i|\right), \quad (2) $$

where $\Theta(x)$ denotes the Heaviside unit-step function and W is the disorder strength. In the subspace of orbitally symmetric wave-functions, the stationary Schrodinger equation $\hat{H}|\Psi\rangle = E|\Psi\rangle$ can be mapped (exactly) onto a close equation for the diagonal amplitudes $\Psi_n = \langle n, n|\Psi\rangle$ of the wavefunction as in Refs. [34, 23]

$$ \frac{1}{U} \Psi_n = \sum_m K_{nm} \Psi_m, \quad (3) $$

where $K_{nm} = \langle n, n|(E\hat{1} - \hat{H}_0)^{-1}|m, m\rangle$. Here $\hat{H}_0 = \hat{H}(U=0)$ represents the noninteracting two-particle Hamiltonian and E the total energy of the pair. Equation (3) represents an eigenvalue problem and can therefore be interpreted as an effective single-particle Schrodinger equation, describing the centre-of-mass motion of the pair. Notice that, in this effective model, the pseudo-energy is given by the inverse of the interaction strength, $\lambda = 1/U$, while the total energy E is a continuous parameter of the matrix K. Unlike the single-particle Anderson model, where the electron can tunnel only between nearest neighboring sites, the matrix K in Eq. (3) contains hopping processes between arbitrarily distant sites.

As in standard transfer matrix calculations, we address the localization properties of the effective model, Eq. (3) by computing the pair transmission amplitude along bar-shaped d-dimensional lattices with $N = M^{d-1}L$ sites and length $L \gg M$. In order to reduce finite-size effects, we impose periodic boundary conditions in the transverse directions and open boundary conditions along the transmission axis. Labelling as {$|\phi_r\rangle, \varepsilon_r$} the eigenbasis of the one-particle Hamiltonian

$$ \hat{H}^{\text{sp}} = -J \sum_{\{i,k\}} |i\rangle\langle k| + \sum_i v_i |i\rangle\langle i|, \quad (4) $$

the matrix elements of the effective model in Eq. (3) can be written as

$$ K_{mn} = \sum_{r=1}^{N} \phi_{r,n}^{*} \phi_{r,m} G_{m,n}^{sp}(E - \varepsilon_r), \quad (5) $$

where $\phi_{r,n} = \langle n|\phi_r\rangle$ and $G_{m,n}^{sp}(E) = \sum_{s=1}^{N} \phi_{s,m}^{*} \phi_{s,n}(E - \varepsilon_r)^{-1}$ are the entries of the matrix resolvent associated to $\hat{H}^{sp}$. Thanks to the open boundary conditions along the longitudinal direction, the single-particle Hamiltonian possesses a block tridiagonal structure that can be exploited to efficiently compute the Green's function $G^{sp}$ via matrix inversion based on a generalized Thomas algorithm [35]. In this way, the total number of elementary operations needed to compute the matrix K scales as $M^{4(d-1)L^3}$ instead of $M^{4(d-1)L^4}$, as expected from standard algorithms for the inversion of a generic symmetric matrix.

The transmission amplitude of the pair between two transverse sections of the bar, the first situated at the head and the second one at a distance $n_z$ inside the bar, is defined in terms of the resolvent matrix
contrast to previous theoretical arguments [30] suggesting that the localization length of the two-particle system remains always finite. The exploration of the 3D version of the same problem has remained unattempted so far.

In this proceeding article, we review our recent numerical works [31, 32, 33] on two-body Anderson localization in two and in three spatial dimensions, by highlighting the most important results.

## 2. Model and method

We consider a system of two particles moving in a disordered lattice and coupled by contact (Hubbard) interactions. The two particles can be bosons or fermions in the spin-singlet state, so that the orbital part of the wavefunction is symmetric under particles exchange and is therefore sensitive to the interactions. The Hamiltonian of the two-particle system in the lattice basis {$|i, j\rangle$} is given by:

$$ \hat{H} = -J \sum_{\{i,k\},j} |i,j\rangle\langle k,j| - J \sum_{i,\{j,k\}} |i,j\rangle\langle i,k| + \sum_{i,j} (v_i + v_j)|i,j\rangle\langle i,j| + U \sum_i |i, i\rangle\langle i, i|, \quad (1) $$

where $J$ is the tunnelling amplitude between nearest-neighbour sites, $U$ is the Hubbard interaction strength and $v_i$ represents the value of the disorder potential at site $i$. As in Anderson's original work [1], we assume that the random energies are uncorrelated, $\overline{v_i v_j} = \overline{v_i^2} \delta_{i,j}$ and obey a uniform on-site distribution:

$$ P(v_i) = \frac{1}{W} \Theta\left(\frac{W}{2} - |v_i|\right), \quad (2) $$

where $\Theta(x)$ denotes the Heaviside unit-step function and $W$ is the disorder strength. In the subspace of orbitally symmetric wave-functions, the stationary Schrodinger equation $\hat{H}|\Psi\rangle = E|\Psi\rangle$ can be mapped (exactly) onto a close equation for the diagonal amplitudes $\Psi_n = \langle n, n|\Psi\rangle$ of the wavefunction as in Refs. [34, 23]

$$ \frac{1}{U} \Psi_n = \sum_m K_{nm} \Psi_m, \quad (3) $$

where $K_{nm} = \langle n, n|(E\hat{1} - \hat{H}_0)^{-1}|m, m\rangle$. Here $\hat{H}_0 = \hat{H}(U=0)$ represents the noninteracting two-particle Hamiltonian and $E$ the total energy of the pair. Equation (3) represents an eigenvalue problem and can therefore be interpreted as an effective single-particle Schrodinger equation, describing the centre-of-mass motion of the pair. Notice that, in this effective model, the pseudo-energy is given by the inverse of the interaction strength, $\lambda = 1/U$, while the total energy $E$ is a continuous parameter of the matrix $K$. Unlike the single-particle Anderson model, where the electron can tunnel only between nearest neighboring sites, the matrix $K$ in Eq. (3) contains hopping processes between arbitrarily distant sites.

As in standard transfer matrix calculations, we address the localization properties of the effective model, Eq. (3) by computing the pair transmission amplitude along bar-shaped $d$-dimensional lattices with $N = M^{d-1}L$ sites and length $L \gg M$. In order to reduce finite-size effects, we impose periodic boundary conditions in the transverse directions and open boundary conditions along the transmission axis. Labelling as {$|\phi_r\rangle, \epsilon_r$} the eigenbasis of the one-particle Hamiltonian

$$ \hat{H}^{\text{sp}} = -J \sum_{\{i,k\}} |i\rangle\langle k| + \sum_i v_i |i\rangle\langle i|, \quad (4) $$

the matrix elements of the effective model in Eq. (3) can be written as

$$ K_{mn} = \sum_{r=1}^{N} \phi_{r,n}^{*} \phi_{r,m} G_{m,n}^{sp}(E - \epsilon_r), \quad (5) $$

where $\phi_{r,n} = \langle n|\phi_r\rangle$ and $G_{m,n}^{sp}(E) = \sum_{s=1}^{N} \phi_{s,m}^{*} \phi_{s,n}(E - \epsilon_r)^{-1}$ are the entries of the matrix resolvent associated to $\hat{H}^{sp}$. Thanks to the open boundary conditions along the longitudinal direction, the single-particle Hamiltonian possesses a block tridiagonal structure that can be exploited to efficiently compute the Green's function $G^{sp}$ via matrix inversion based on a generalized Thomas algorithm [35]. In this way, the total number of elementary operations needed to compute the matrix $K$ scales as $M^{4(d-1)L^3}$ instead of $M^{4(d-1)L^4}$, as expected from standard algorithms for the inversion of a generic symmetric matrix.

The transmission amplitude of the pair between two transverse sections of the bar, the first situated at the head and the second one at a distance $n_z$ inside the bar, is defined in terms of the resolvent matrix
<figure><img src="image_5.png" /><figcaption>Figure 1. Reduced localization length of the pair plotted as a function of the interaction strength for increasing values of the transverse size M of the grid. The disorder strength is fixed to W = 9 and the pair has zero total energy, E = 0. In Panel (a) the curves for M = {8, 10, 12, 16, 20} and L = 400 are plotted, averaging the Lyapunov exponent over Ntr = {600, 600, 600, 1000} configurations of the random potential respectively. In Panel (b) the same quantities are shown for larger grids with transverse sizes M = {30, 40, 50} and L = 500. The corresponding number of different disorder realizations used is Ntr = {3600, 4400, 2850}, respectively.</figcaption></figure>

$G^{\lambda} = (\lambda 1 - K)^{-1}$ as

$$t_{1n_z} = \sum_{m_{\perp}, n_{\perp}} J^2 |G_{(m_{\perp},1),(n_{\perp},n_z)}^{\lambda}|^2, \quad (6)$$

where the sum is taken over all sites $m_{\perp}, n_{\perp}$ of two corresponding sections. The logarithm of the transmission amplitude (6) is a self-averaging quantity. For each disorder realization, we evaluate it at regular intervals along the bar and apply a linear fit to the data, $\ln t_{1n_z}^{\text{fit}} = pn_z + q$. For a given value of the interaction strength, we evaluate the Lyapunov exponent $\gamma = \gamma(M, U)$ as $\gamma = -\bar{p}/2$, where $\bar{p}$ is the average of the slope over the different disorder realizations. We then infer the transport properties of the system from the behaviour of the reduced localization length, which is defined as $\Lambda_M = (M\gamma)^{-1}$. In the metallic phase $\Lambda_M$ increases as M increases, whereas in the insulating phase it does the opposite. At the critical point, $\Lambda_M$ becomes constant for values of M sufficiently large. Hence the critical point $U = U_c$ of the Anderson transition can be identified by plotting the reduced localization length versus U for different values of the transverse size M and looking at their common crossing points. In the following sections we will consider the hopping amplitude as the energy scale of the system and thus we will set $J = 1$.

### 3. Absence of 2D MITs for the pair

In two-dimensional lattices we have computed the reduced localization length in the region of the parameters space in which a delocalization transition was first predicted [27, 28]. In Ref. [27], the metal-insulator transition was predicted to occur for $E = 0$ and $W = 9.3 \pm 0.5$, based on transmission-amplitude calculations performed on rectangular strips of width $M = 10$ and length $L = 62$ at most.

As we notice from Fig. 3, $\Lambda_M$ possesses a nonmonotonic dependence on U, as previously found in one dimension [25], with a peak near $U = 6$ whereas, at the band centre ($E = 0$), $\Lambda_M$ presents the symmetry $\Lambda_M(-U) = \Lambda_M(U)$. As shown in Panel (a) of Fig. 3, the curves corresponding to different values of M between $M = 8$ and $M = 20$ intersect each other around $U = 1$, suggesting a possible phase transition, as already reported in Refs. [27, 28].

This result clearly conflicts with the general expectation that 2D Anderson transitions are forbidden in models where both time-reversal and spin-rotational symmetries are preserved. For the effective model K, the two invariances are inherited from the single-particle Anderson model. In order to test the stability of the result against finite-size effects, we have performed additional calculation using larger grids, going from $M = 30$ to $M = 50$. The results, plotted in Panel (b) of Fig. 3 show that the crossing points have completely disappeared, suggesting that the pair localizes in an infinite lattice, irrespectively of the specific value of U. We therefore conclude that the former results [27, 28] were plagued by severe finite-size effects, due to the limited system sizes accessible at that time, and no Anderson transition can take
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Reduced localization length of the pair plotted as a function of the interaction strength for increasing values of the transverse size *M* of the grid. The disorder strength is fixed to *W* = 9 and the pair has zero total energy, *E* = 0. In Panel (a) the curves for *M* = {8, 10, 12, 16, 20} and *L* = 400 are plotted, averaging the Lyapunov exponent over *N*<sub>tr</sub> = {600, 600, 600, 1000} configurations of the random potential respectively. In Panel (b) the same quantities are shown for larger grids with transverse sizes *M* = {30, 40, 50} and *L* = 500. The corresponding number of different disorder realizations used is *N*<sub>tr</sub> = {3600, 4400, 2850}, respectively.</figcaption></figure>

$G^\lambda = (\lambda \mathbf{1} - K)^{-1}$ as

$$t_{1n_z} = \sum_{m_{\perp}, n_{\perp}} J^2 |G_{(m_{\perp},1),(n_{\perp},n_z)}^{\lambda}|^2, \qquad (6)$$

where the sum is taken over all sites $m_{\perp}, n_{\perp}$ of two corresponding sections. The logarithm of the transmission amplitude (6) is a self-averaging quantity. For each disorder realization, we evaluate it at regular intervals along the bar and apply a linear fit to the data, $\ln t_{1n_z}^{\text{fit}} = pn_z + q$. For a given value of the interaction strength, we evaluate the Lyapunov exponent $\gamma = \gamma(M, U)$ as $\gamma = -\bar{p}/2$, where $\bar{p}$ is the average of the slope over the different disorder realizations. We then infer the transport properties of the system from the behaviour of the reduced localization length, which is defined as $\Lambda_M = (M\gamma)^{-1}$. In the metallic phase $\Lambda_M$ increases as $M$ increases, whereas in the insulating phase it does the opposite. At the critical point, $\Lambda_M$ becomes constant for values of $M$ sufficiently large. Hence the critical point $U = U_c$ of the Anderson transition can be identified by plotting the reduced localization length versus $U$ for different values of the transverse size $M$ and looking at their common crossing points. In the following sections we will consider the hopping amplitude as the energy scale of the system and thus we will set $J = 1$.

### 3. Absence of 2D MITs for the pair

In two-dimensional lattices we have computed the reduced localization length in the region of the parameters space in which a delocalization transition was first predicted [27, 28]. In Ref. [27], the metal-insulator transition was predicted to occur for $E = 0$ and $W = 9.3 \pm 0.5$, based on transmission-amplitude calculations performed on rectangular strips of width $M = 10$ and length $L = 62$ at most.

As we notice from Fig. 3, $\Lambda_M$ possesses a nonmonotonic dependence on $U$, as previously found in one dimension [25], with a peak near $U = 6$ whereas, at the band centre ($E = 0$), $\Lambda_M$ presents the symmetry $\Lambda_M(-U) = \Lambda_M(U)$. As shown in Panel (a) of Fig. 3, the curves corresponding to different values of $M$ between $M = 8$ and $M = 20$ intersect each other around $U = 1$, suggesting a possible phase transition, as already reported in Refs. [27, 28].

This result clearly conflicts with the general expectation that 2D Anderson transitions are forbidden in models where both time-reversal and spin-rotational symmetries are preserved. For the effective model *K*, the two invariances are inherited from the single-particle Anderson model. In order to test the stability of the result against finite-size effects, we have performed additional calculation using larger grids, going from $M = 30$ to $M = 50$. The results, plotted in Panel (b) of Fig. 3 show that the crossing points have completely disappeared, suggesting that the pair localizes in an infinite lattice, irrespectively of the specific value of $U$. We therefore conclude that the former results [27, 28] were plagued by severe finite-size effects, due to the limited system sizes accessible at that time, and no Anderson transition can take
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Phase boundary between localized (labelled by I) and extended (M) states in the (U, W) plane for a pair with zero total energy, *E* = 0. The dashed horizontal line corresponds to the noninteracting limit, *W* = *W*<sub>c</sub><sup>sp</sup> = 16.54. Inset: disorder-averaged density of states ρ<sub>K</sub> of the effective Hamiltonian of the pair calculated for *W* = 23.5 using a cubic grid of sizes *L* = *M* = 20 with periodic boundary conditions. (b) Phase diagram for a pair with total energy *E* = *U* (blue circles data) [32]. The orange triangles data refer to the phase boundary at *E* = 0, calculated in Ref. [31], which is added for comparison. The dot-dashed purple line at low disorder corresponds to the molecular result, *W*<sub>c</sub> ≃ 16.0/*|*U*|*, obtained by treating the pair as a point-like molecule obeying an effective Anderson model, see Eq. (7). The two diagrams hold for both attractive and repulsive interactions.</figcaption></figure>

place for a pair subject to short-range interactions in 2D disordered lattices. Our work therefore reconciles
the numerics with the one-parameter scaling theory of localization and the concept of universality classes.

Based on a least-square variational procedure [36], in Ref. [33] we have verified that our numerical
results for the reduced localization length obey the one-parameter scaling ansatz Λ<sub>M</sub> = f(ξ/M), where
f(x) is the scaling function and ξ is the localization length of the infinite system, up to a scale factor.
Notice that the scale factor is a pure number, independent of the value of the interaction strength.
For weak to moderate interactions, ξ displays an exponential increase over more than three orders of
magnitude, reaching its maximum value around U ~ 6.

We have also investigated the case of a non-zero energy of the pair. In particular our numerical
results show that Λ<sub>M</sub>(E, U, W) ≤ Λ<sub>M</sub>(0, U, W). As a consequence, the pair localizes with an even shorter
localization length with respect to the E = 0 case [33].

4. Two-body mobility edges in 3D

The scenario changes when three-dimensional systems are examined. We have first determined the critical
point for zero total energy of the pair, *E* = 0, and fixed disorder strength *W* = 23.5. By carrying out the
same analysis presented in the previous section, we identified a critical region where the curves Λ*<sub>*M*</sub>(*U*)
obtained for different *M* cross each other. Since the values of *M* numerically accessible are rather small,
the crossing point drifts to stronger interaction and upwards as the system sizes increase. For the largest
system sizes, the height of the crossing point becomes close to the universal value, Λ*<sub>*c,orth*</sub>* = 0.58* ± 0.09,
expected for the orthogonal universality class, within the accuracy of our numerics. Then we use this
information to pinpoint the position of the critical point by including the leading irrelevant variable
*y*<sub>orth</sub> = 3.3 ± 0.6 [37] in the scaling analysis, thus accounting for the drift of the crossing point. In this
way we were able to identify the critical interaction strength as *U*<sub>c</sub> = 2.16 ± 0.04 [31].

The same procedure has then been systematically applied to explore the phase diagram in the space
of energy, disorder and interaction. We first determined the phase boundary at zero total energy of the
pair, which is displayed in Fig. 4. We see that, for zero total energy of the pair, Anderson transitions occur
for disorder strength *W* > 16.5, where all single-particle states are localized [37]. The critical disorder
strength exhibits a nonmonotonic behaviour as a function of the interaction, reaching its maximum value
at *U* ~ 14. This feature can be understood from the behaviour of the disorder-averaged density of states
ρ<sub>K</sub>(λ) = Σ<sub>r</sub> δ(λ − λ<sub>r</sub>)/N of the effective model, which is plotted in the inset of Panel (a) of Fig. 4.
Notice that for *E* = 0 this quantity is parity-symmetric, (ρ<sub>K</sub>(λ) = ρ<sub>K</sub>(−λ)), like the reduced localization
length. The weakly interacting states, located in the tails of the density of states, where ρ<sub>K</sub> ∝ λ<sup>−2</sup>,
localize at lower disorder than those close to the peaks, occurring at λ = ±*W*<sup>−1</sup>.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** (a) Phase boundary between localized (labelled by I) and extended (M) states in the (U, W) plane for a pair with zero total energy, *E* = 0. The dashed horizontal line corresponds to the noninteracting limit, *W* = *W*<sub>c</sub><sup>sp</sup> = 16.54. Inset: disorder-averaged density of states ρ<sub>K</sub> of the effective Hamiltonian of the pair calculated for *W* = 23.5 using a cubic grid of sizes *L* = *M* = 20 with periodic boundary conditions. (b) Phase diagram for a pair with total energy *E* = *U* (blue circles data) [32]. The orange triangles data refer to the phase boundary at *E* = 0, calculated in Ref. [31], which is added for comparison. The dot-dashed purple line at low disorder corresponds to the molecular result, *W*<sub>c</sub> ≃ 16.0/*|*U*|*, obtained by treating the pair as a point-like molecule obeying an effective Anderson model, see Eq. (7). The two diagrams hold for both attractive and repulsive interactions.</figcaption></figure>

place for a pair subject to short-range interactions in 2D disordered lattices. Our work therefore reconciles
the numerics with the one-parameter scaling theory of localization and the concept of universality classes.

Based on a least-square variational procedure [36], in Ref. [33] we have verified that our numerical
results for the reduced localization length obey the one-parameter scaling ansatz Λ<sub>M</sub> = f(ξ/M), where
f(x) is the scaling function and ξ is the localization length of the infinite system, up to a scale factor.
Notice that the scale factor is a pure number, independent of the value of the interaction strength.
For weak to moderate interactions, ξ displays an exponential increase over more than three orders of
magnitude, reaching its maximum value around U ~ 6.

We have also investigated the case of a non-zero energy of the pair. In particular our numerical
results show that Λ<sub>M</sub>(E, U, W) ≤ Λ<sub>M</sub>(0, U, W). As a consequence, the pair localizes with an even shorter
localization length with respect to the E = 0 case [33].

4. Two-body mobility edges in 3D

The scenario changes when three-dimensional systems are examined. We have first determined the critical
point for zero total energy of the pair, *E* = 0, and fixed disorder strength *W* = 23.5. By carrying out the
same analysis presented in the previous section, we identified a critical region where the curves Λ*<sub>*M*</sub>(*U*)
obtained for different *M* cross each other. Since the values of *M* numerically accessible are rather small,
the crossing point drifts to stronger interaction and upwards as the system sizes increase. For the largest
system sizes, the height of the crossing point becomes close to the universal value, Λ*<sub>*c,orth*</sub>* = 0.58* ± 0.09,
expected for the orthogonal universality class, within the accuracy of our numerics. Then we use this
information to pinpoint the position of the critical point by including the leading irrelevant variable
*y*<sub>orth</sub> = 3.3 ± 0.6 [37] in the scaling analysis, thus accounting for the drift of the crossing point. In this
way we were able to identify the critical interaction strength as *U*<sub>c</sub> = 2.16 ± 0.04 [31].

The same procedure has then been systematically applied to explore the phase diagram in the space
of energy, disorder and interaction. We first determined the phase boundary at zero total energy of the
pair, which is displayed in Fig. 4. We see that, for zero total energy of the pair, Anderson transitions occur
for disorder strength *W* > 16.5, where all single-particle states are localized [37]. The critical disorder
strength exhibits a nonmonotonic behaviour as a function of the interaction, reaching its maximum value
at *U* ~ 14. This feature can be understood from the behaviour of the disorder-averaged density of states
ρ<sub>K</sub>(λ) = Σ<sub>r</sub> δ(λ - λ<sub>r</sub>) / N of the effective model, which is plotted in the inset of Panel (a) of Fig. 4.
Notice that for *E* = 0 this quantity is parity-symmetric, (ρ<sub>K</sub>(λ) = ρ<sub>K</sub>(-λ)), like the reduced localization
length. The weakly interacting states, located in the tails of the density of states, where ρ<sub>K</sub> ∝ λ<sup>-2</sup>,
localize at lower disorder than those close to the peaks, occurring at λ = ±*W*<sup>-1</sup>.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Phase diagrams in the (U, W) plane for a pair with E = -12.25 (a), E = -15 (b) and E = -18 (c). The phase boundaries are denoted by the orange, green and purple triangles respectively. The numerical band edges are represented by the dotted lines, whereas the black dashed lines denote the exact band edges, which delimit the regions in which no two-particle states are present, shaded in light grey. The brown arrows indicate the disorder strength at which the molecular state evolve into a scattering state. For this value of disorder the right numerical band edge crosses the U = 0 axis, as signalled by the brown stars.</figcaption></figure>

Making *J* explicit, in the strongly interacting regime, that corresponds to *E* ~ *U* with |*E*| ~ ~*W*, *J*, the effective Hamiltonian takes the same form of the single-particle Anderson model, but endowed with modified parameters, describing the motion of tightly bound pairs [34, 32]

$$ \hat{K}_{\text{TB}} \Psi(\mathbf{n}) = \frac{2J^2}{E} \sum_{i=1}^{3} \Psi(\mathbf{n} + \mathbf{e}_i) + \frac{2J^2}{E} \sum_{i=1}^{3} \Psi(\mathbf{n} - \mathbf{e}_i) + \frac{2v_{\mathbf{n}}(E + 2v_{\mathbf{n}})}{E} \Psi(\mathbf{n}) = E_{\text{TB}} \Psi(\mathbf{n}), \quad (7) $$

where the tunnelling amplitude is $J_{TB} = -2J^2/E$, the disorder typical value is $W_{TB} = 2W$, and $E_{TB} = E^2(\lambda - 1/E - 12J^2/E^3)$ denotes the effective energy of the molecule. Notice that $|J_{TB}| \sim 2J^2/|U| \gg J$, implying that the molecule localizes already for weak disorder. In particular we can estimate the position of the mobility edge for tightly bound molecules from the corresponding results for the single-particle Anderson model, upon substitution of the bare parameters with the molecular ones.

In the Panel (b) of Fig. 4 we display the phase boundary between localized and extended states calculated along the direction *E* = *U*. For weak interactions, the mobility edge remains very close to the result previously obtained for *E* = 0. In the opposite strongly interacting regime, the system exhibits a single Anderson transition, corresponding to the localization of tightly bound molecules with $W_c \simeq 16/|U|$, as evidenced by the purple dot-dashed curve in Fig. 4 (b). For intermediate values of the interaction strength, $12 \lesssim U \lesssim 19$, the critical disorder strength displays an interesting s-like shape, signalling that the system undergoes three different metal-insulator transitions (instead of one) at given interaction strength.

To the purpose of understanding the emergence of this region, we have built the three diagrams at fixed energy displayed in Fig. 3, corresponding to $E = -12.25, -15, -18$. We notice that, at very low disorder, the pair exhibits a metallic phase corresponding to delocalized molecule, which shrinks in size as *E* becomes large and negative, due to the increased inertia of the bound state. A striking feature of the three phase diagrams is the existence of disorder-induced delocalization transitions, resulting in a metallic phase for intermediate values of *W*. This counterintuitive result reflects a fundamental change in the nature of the pair state as the disorder strength increases. For sufficiently weak disorder and $|E| > 12$, the pair behaves as a molecule and, in particular, its energy falls outside the noninteracting two-particles energy spectrum. Neglecting Lifshitz tails [38], the latter corresponds to the interval $[-2\epsilon_{BE}, 2\epsilon_{BE}]$, where $\pm\epsilon_{BE}(W)$ are the numerical band edges of the single-particle model, which can be easily determined by making use of the coherent potential approximation (CPA) [39, 40]. As *W* increases, the noninteracting energy spectrum broadens and the energy *E* will inevitably fall inside it, transforming the molecule into a scattering state, which is more akin to delocalize. This change of behaviour is highlighted by the brown lateral arrows in the phase diagrams of Fig. 3.

In Fig. 3, the phase boundary for scattering states (at high disorder) does not depend on the sign of the interaction, and, when the latter vanishes, the critical disorder is consistent with the one found in the single-particle case [41]. In accordance with Ioffe-Regel criterion [42], the states in the limit of infinite interactions (in magnitude) must possess a large mean free path, since the phase boundary follows closely the numerical band edge, which is represented by the dotted brown curve and delimits the region
<figure><img src="image_8.png" /><figcaption>Figure 3. Phase diagrams in the (U, W) plane for a pair with E = -12.25 (a), E = -15 (b) and E = -18 (c). The phase boundaries are denoted by the orange, green and purple triangles respectively. The numerical band edges are represented by the dotted lines, whereas the black dashed lines denote the exact band edges, which delimit the regions in which no two-particle states are present, shaded in light grey. The brown arrows indicate the disorder strength at which the molecular state evolve into a scattering state. For this value of disorder the right numerical band edge crosses the U = 0 axis, as signalled by the brown stars.</figcaption></figure>

Making *J* explicit, in the strongly interacting regime, that corresponds to *E* ~ *U* with |*E*| ~ *W*, *J*, the effective Hamiltonian takes the same form of the single-particle Anderson model, but endowed with modified parameters, describing the motion of tightly bound pairs [34, 32]

$$ \hat{K}_{\text{TB}} \Psi(\mathbf{n}) = \frac{2J^2}{E} \sum_{i=1}^{3} \Psi(\mathbf{n} + \mathbf{e}_i) + \frac{2J^2}{E} \sum_{i=1}^{3} \Psi(\mathbf{n} - \mathbf{e}_i) + \frac{2v_{\mathbf{n}}(E + 2v_{\mathbf{n}})}{E} \Psi(\mathbf{n}) = E_{\text{TB}} \Psi(\mathbf{n}), \quad (7) $$

where the tunnelling amplitude is $J_{TB} = -2J^2/E$, the disorder typical value is $W_{TB} = 2W$, and $E_{TB} = E^2(\lambda - 1/E - 12J^2/E^3)$ denotes the effective energy of the molecule. Notice that $|J_{TB}| \sim 2J^2/|U| \gg J$, implying that the molecule localizes already for weak disorder. In particular we can estimate the position of the mobility edge for tightly bound molecules from the corresponding results for the single-particle Anderson model, upon substitution of the bare parameters with the molecular ones.

In the Panel (b) of Fig. 4 we display the phase boundary between localized and extended states calculated along the direction *E* = *U*. For weak interactions, the mobility edge remains very close to the result previously obtained for *E* = 0. In the opposite strongly interacting regime, the system exhibits a single Anderson transition, corresponding to the localization of tightly bound molecules with $W_c \simeq 16/|U|$, as evidenced by the purple dot-dashed curve in Fig. 4 (b). For intermediate values of the interaction strength, $12 \lesssim U \lesssim 19$, the critical disorder strength displays an interesting s-like shape, signalling that the system undergoes three different metal-insulator transitions (instead of one) at given interaction strength.

To the purpose of understanding the emergence of this region, we have built the three diagrams at fixed energy displayed in Fig. 3, corresponding to $E = -12.25, -15, -18$. We notice that, at very low disorder, the pair exhibits a metallic phase corresponding to delocalized molecule, which shrinks in size as *E* becomes large and negative, due to the increased inertia of the bound state. A striking feature of the three phase diagrams is the existence of disorder-induced delocalization transitions, resulting in a metallic phase for intermediate values of *W*. This counterintuitive result reflects a fundamental change in the nature of the pair state as the disorder strength increases. For sufficiently weak disorder and $|E| > 12$, the pair behaves as a molecule and, in particular, its energy falls outside the noninteracting two-particles energy spectrum. Neglecting Lifshitz tails [38], the latter corresponds to the interval $[-2\epsilon_{BE}, 2\epsilon_{BE}]$, where $\pm\epsilon_{BE}(W)$ are the numerical band edges of the single-particle model, which can be easily determined by making use of the coherent potential approximation (CPA) [39, 40]. As *W* increases, the noninteracting energy spectrum broadens and the energy *E* will inevitably fall inside it, transforming the molecule into a scattering state, which is more akin to delocalize. This change of behaviour is highlighted by the brown lateral arrows in the phase diagrams of Fig. 3.

In Fig. 3, the phase boundary for scattering states (at high disorder) does not depend on the sign of the interaction, and, when the latter vanishes, the critical disorder is consistent with the one found in the single-particle case [41]. In accordance with Ioffe-Regel criterion [42], the states in the limit of infinite interactions (in magnitude) must possess a large mean free path, since the phase boundary follows closely the numerical band edge, which is represented by the dotted brown curve and delimits the region
where $\rho_K$ is negligible. The nonmonotonicity of the upper phase boundary of extended scattering states is clearly reminiscent of the one computed for $E = 0$ (see Fig. 4 (a)), but the maxima of $W_c$ generally occur for lower values of the disorder.

From the three panels of Fig. 3, it is possible to observe that the single metallic region for $E = -12.25$ (panel a) splits into two regions for $E = -15$ related to molecular and scattering states (panel b), then the latter break up in turn into two disconnected regions for $E = -18$ (panel c). In that case, there are no delocalized states at vanishing interactions, in accordance with the single-particle case [41].

Finally, the grey regions of the diagrams of Fig. 3 correspond to configurations unaccessible to the system, where the density of state is exactly zero. These regions are limited by the black dashed curves, corresponding to the rigorous band edges, which have been calculated analytically [32].

## 5. Conclusions

Based on large-scale numerics, we have investigated the localization properties of a two-particle system described by the Anderson-Hubbard model. We have first proved that Anderson transitions for the pair do not occur in two dimensions, as opposed to previous claims [27, 28]. In three dimensions, we have also shown that the universality class of the transition is not affected by interactions and that the critical disorder at the band centre is a nonmonotonic function of interaction strength, a feature we have interpreted by analyzing the behaviour of the disorder-averaged density of states. The phase diagram gets particularly rich when the total energy is finite and becomes larger (in modulus) than the bandwidth, corresponding to $|E| > 12$. In this case, the phase diagram exhibits multiple Anderson transitions as the disorder increases, reflecting the underlying change from molecular to scattering states. Our results represent the zero-particle-density limit of the many-body localization problem [17, 18].

## 6. Acknowledgements

We acknowledge D. Basko, N. Cherroret, D. Delande, K. Frahm, C. Monthus, S. Pilati, T. Roscilde and S. Skipetrov for fruitful discussions. 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. The numerical study of the two-body localization problem reviewed in this work has required in total about three millions hours of CPU times. For its realization, access was granted to the HPC resources of CINES (Centre Informatique National de l'Enseignement Supérieur) under the allocations 2016-c2016057629, 2017-A0020507629, 2018-A0040507629, 2019-A0060507629 and 2020-A0080507629 supplied by GENCI (Grand Équipement National de Calcul Intensif).

## 7. References

[1] Anderson P W 1958 *Phys. Rev.* **109** 1492

[2] Schwartz T, Bartal G, Fishman S and Segev B 2007 *Nature (London)* **446** 52–55

[3] Faez S, Strybulevych A, Page J H and Lagendijk A and van Tiggelen B A 2009 *Phys. Rev. Lett.* **103** 155703

[4] Chabanov A A, Stoytchev M and Genack A Z 2000 *Nature (London)* **404** 850–853

[5] Billy J, Josse V, Zuo Z, Bernard A, Hambrecht B, Lugan P, Clément D, Sanchez-Palencia L, Bouyer P and Aspect A 2008 *Nature (London)* **453** 891–894

[6] Roati G, D'Errico C, Fallani L, Fattori M, Fort C, Zaccanti M, Modugno G, Modugno M and Inguscio M 2008 *Nature (London)* **453** 895–898

[7] Lopez M, Clément J F, Szriftgiser P, Garreau J C and Delande D 2012 *Phys. Rev. Lett.* **108** 095701

[8] Mott N F 1967 *Adv. Phys.* **16** 49

[9] Fleishman L and Anderson P W 1980 *Phys. Rev. B* **21**(6) 2366–2377

[10] Giamarchi T and Schulz H J 1988 *Phys. Rev. B* **37**(1) 325–340

[11] Finkelstein A M 1983 *Sov. Phys. J. Exp. Theor. Phys.* **57** 97

[12] Kravchenko S V, Kravchenko G V, Furneaux J E, Pudalov V M and D'Iorio M 1994 *Phys. Rev. B* **50**(11) 8039–8042

[13] Punnoose A and Finkel'stein A M 2005 *Science* **310** 289–291

[14] Shepelyansky D L 2000 *Phys. Rev. B* **61**(7) 4588–4591

[15] Song P H and Shepelyansky D L 2000 *Phys. Rev. B* **61**(23) 15546–15549

[16] Gorny I V, Mirlin A and Polyakov D G 2005 *Phys. Rev. Lett.* **95** 206603

[17] Basko D M, Aleiner I L and Altshuler B L 2006 *Ann. Phys.* **321** 1126–1205

[18] Abanin D A, Altman E, Bloch I and Serbyn M 2019 *Rev. Mod. Phys.* **91** 021001

[19] De Roeck W and Huveneers F 2017 *Phys. Rev. B* **95**(15) 155129

[20] Wahl T B, Pal A and Simon S H 2019 *Nature Physics* **15** 164–169

[21] Dorokhov O 1990 *Sov. Phys. JETP* **71** 360

[22] Shepelyansky D L 1994 *Phys. Rev. Lett.* **73**(19) 2607–2610

[23] von Oppen F, Wettig T and Müller J 1996 *Phys. Rev. Lett.* **76**(3) 491–494
where $\rho_K$ is negligible. The nonmonotonicity of the upper phase boundary of extended scattering states is clearly reminiscent of the one computed for $E = 0$ (see Fig. 4 (a)), but the maxima of $W_c$ generally occur for lower values of the disorder.

From the three panels of Fig. 3, it is possible to observe that the single metallic region for $E = -12.25$ (panel a) splits into two regions for $E = -15$ related to molecular and scattering states (panel b), then the latter break up in turn into two disconnected regions for $E = -18$ (panel c). In that case, there are no delocalized states at vanishing interactions, in accordance with the single-particle case [41].

Finally, the grey regions of the diagrams of Fig. 3 correspond to configurations unaccessible to the system, where the density of state is exactly zero. These regions are limited by the black dashed curves, corresponding to the rigorous band edges, which have been calculated analytically [32].

## 5. Conclusions

Based on large-scale numerics, we have investigated the localization properties of a two-particle system described by the Anderson-Hubbard model. We have first proved that Anderson transitions for the pair do not occur in two dimensions, as opposed to previous claims [27, 28]. In three dimensions, we have also shown that the universality class of the transition is not affected by interactions and that the critical disorder at the band centre is a nonmonotonic function of interaction strength, a feature we have interpreted by analyzing the behaviour of the disorder-averaged density of states. The phase diagram gets particularly rich when the total energy is finite and becomes larger (in modulus) than the bandwidth, corresponding to $|E| > 12$. In this case, the phase diagram exhibits multiple Anderson transitions as the disorder increases, reflecting the underlying change from molecular to scattering states. Our results represent the zero-particle-density limit of the many-body localization problem [17, 18].

## 6. Acknowledgements

We acknowledge D. Basko, N. Cherroret, D. Delande, K. Frahm, C. Monthus, S. Pilati, T. Roscilde and S. Skipetrov for fruitful discussions. 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. The numerical study of the two-body localization problem reviewed in this work has required in total about three millions hours of CPU times. For its realization, access was granted to the HPC resources of CINES (Centre Informatique National de l'Enseignement Supérieur) under the allocations 2016-c2016057629, 2017-A0020507629, 2018-A0040507629, 2019-A0060507629 and 2020-A0080507629 supplied by GENCI (Grand Équipement National de Calcul Intensif).

## 7. References

[1] Anderson P W 1958 *Phys. Rev.* **109** 1492

[2] Schwartz T, Bartal G, Fishman S and Segev B 2007 *Nature (London)* **446** 52–55

[3] Faez S, Strybulevych A, Page J H and Lagendijk A and van Tiggelen B A 2009 *Phys. Rev. Lett.* **103** 155703

[4] Chabanov A A, Stoytchev M and Genack A Z 2000 *Nature (London)* **404** 850–853

[5] Billy J, Josse V, Zuo Z, Bernard A, Hambrecht B, Lugan P, Clément D, Sanchez-Palencia L, Bouyer P and Aspect A 2008 *Nature (London)* **453** 891–894

[6] Roati G, D'Errico C, Fallani L, Fattori M, Fort C, Zaccanti M, Modugno G, Modugno M and Inguscio M 2008 *Nature (London)* **453** 895–898

[7] Lopez M, Clément J F, Szriftgiser P, Garreau J C and Delande D 2012 *Phys. Rev. Lett.* **108** 095701

[8] Mott N F 1967 *Adv. Phys.* **16** 49

[9] Fleishman L and Anderson P W 1980 *Phys. Rev. B* **21**(6) 2366–2377

[10] Giamarchi T and Schulz H J 1988 *Phys. Rev. B* **37**(1) 325–340

[11] Finkelstein A M 1983 *Sov. Phys. J. Exp. Theor. Phys.* **57** 97

[12] Kravchenko S V, Kravchenko G V, Furneaux J E, Pudalov V M and D'Iorio M 1994 *Phys. Rev. B* **50**(11) 8039–8042

[13] Punnoose A and Finkel'stein A M 2005 *Science* **310** 289–291

[14] Shepelyansky D L 2000 *Phys. Rev. B* **61**(7) 4588–4591

[15] Song P H and Shepelyansky D L 2000 *Phys. Rev. B* **61**(23) 15546–15549

[16] Gorny I V, Mirlin A and Polyakov D G 2005 *Phys. Rev. Lett.* **95** 206603

[17] Basko D M, Aleiner I L and Altshuler B L 2006 *Ann. Phys.* **321** 1126–1205

[18] Abanin D A, Altman E, Bloch I and Serbyn M 2019 *Rev. Mod. Phys.* **91** 021001

[19] De Roeck W and Huveneers F 2017 *Phys. Rev. B* **95**(15) 155129

[20] Wahl T B, Pal A and Simon S H 2019 *Nature Physics* **15** 164–169

[21] Dorokhov O 1990 *Sov. Phys. JETP* **71** 360

[22] Shepelyansky D L 1994 *Phys. Rev. Lett.* **73**(19) 2607–2610

[23] von Oppen F, Wettig T and Müller J 1996 *Phys. Rev. Lett.* **76**(3) 491–494
[24] Frahm K M 1999 *Eur. Phys. J. B* **10** 371

[25] Frahm K M 2016 *Eur. Phys. J. B* **89** 115

[26] Thongjaomayum D, Andreanov A, Engl T and Flach S 2019 *Phys. Rev. B* **100**(22) 224203

[27] Ortuño M and Cuevas E 1999 *Europhysics Letters* **46** 224

[28] Roemer R A, Leadbeater M and Schreiber M 1999 *Ann. Phys. (Leipzig)* **8** 675

[29] Chattaraj T 2018 *Condens. Matter* **3**(4) 38

[30] Imry Y 1995 *Europhys. Lett.* **30** 405

[31] Stellin F and Orso G 2019 *Phys. Rev. B* **99**(22) 224209

[32] Stellin F and Orso G 2020 *Phys. Rev. Research* **2**(3) 033501

[33] Stellin F and Orso G 2020 *Phys. Rev. B* **102**(14) 144201

[34] Dufour G and Orso G 2012 *Phys. Rev. Lett.* **109** 155306

[35] Jain J, Li H, Cauley S, Koh C K and Balakrishnan V 2007 *Purdue ECE Technical Reports. Paper 357*

[36] McKinnon A and Kramer B 1983 *Z. Phys. B* **53** 1

[37] Slevin K and Ohtsuki T 2014 *New Journal of Physics* **16** 015012

[38] Lifshitz I 1963 *Sov. Phys. JETP* **17** 1159-1170

[39] Elliott R J, Krumhansl J A and Leath P L 1974 *Rev. Mod. Phys.* **46**(3) 465-543

[40] Vollhardt D and Wölfle P 1992 *Self-Consistent Theory of Anderson Localization Electronic Phase Transitions* ed Hanke, W and Kopaev Yu V (Elsevier) pp 1-78

[41] Bulka B, Schreiber M and Kramer B 1987 *Z. Physik B - Condensed Matter* **66** 21

[42] Ioffe A and Regel A 1960 *Prog. Semicond.* **4** 237
[24] Frahm K M 1999 *Eur. Phys. J. B* **10** 371

[25] Frahm K M 2016 *Eur. Phys. J. B* **89** 115

[26] Thongjaomayum D, Andreanov A, Engl T and Flach S 2019 *Phys. Rev. B* **100**(22) 224203

[27] Ortuño M and Cuevas E 1999 *Europhysics Letters* **46** 224

[28] Roemer R A, Leadbeater M and Schreiber M 1999 *Ann. Phys. (Leipzig)* **8** 675

[29] Chattaraj T 2018 *Condens. Matter* **3**(4) 38

[30] Imry Y 1995 *Europhys. Lett.* **30** 405

[31] Stellin F and Orso G 2019 *Phys. Rev. B* **99**(22) 224209

[32] Stellin F and Orso G 2020 *Phys. Rev. Research* **2**(3) 033501

[33] Stellin F and Orso G 2020 *Phys. Rev. B* **102**(14) 144201

[34] Dufour G and Orso G 2012 *Phys. Rev. Lett.* **109** 155306

[35] Jain J, Li H, Cauley S, Koh C K and Balakrishnan V 2007 *Purdue ECE Technical Reports. Paper 357*

[36] McKinnon A and Kramer B 1983 *Z. Phys. B* **53** 1

[37] Slevin K and Ohtsuki T 2014 *New Journal of Physics* **16** 015012

[38] Lifshitz I 1963 *Sov. Phys. JETP* **17** 1159-1170

[39] Elliott R J, Krumhansl J A and Leath P L 1974 *Rev. Mod. Phys.* **46**(3) 465-543

[40] Vollhardt D and Wölfle P 1992 *Self-Consistent Theory of Anderson Localization Electronic Phase Transitions* ed Hanke, W and Kopaev Yu V (Elsevier) pp 1-78

[41] Bulka B, Schreiber M and Kramer B 1987 *Z. Physik B - Condensed Matter* **66** 21

[42] Ioffe A and Regel A 1960 *Prog. Semicond.* **4** 237