Electron pairing by Coulomb repulsion in narrow band structures

We study analytically and numerically dynamics and eigenstates of two electrons with Coulomb repulsion on a tight-binding lattice in one and two dimensions. The total energy and momentum of electrons are conserved and we show that for a certain momentum range the dynamics is exactly reduced to an evolution in an effective narrow energy band where the energy conservation forces the two electrons to propagate together through the whole system at moderate or even weak repulsion strength. We argue that such a mechanism of electron pair formation by the repulsive Coulomb interaction is rather generic and that it can be at the origin of unconventional superconductivity in twisted bilayer graphene.

Introduction. -The interactions of electrons in narrow band structures play an important role in various physical processes. Often the theoretical analysis is done in the frame work of electrons on a lattice with short range on-site interactions as described in works of . Recently, the interest to narrow band structures with strong Coulomb electronelectron interactions has been inspired by the observation of unconventional superconductivity in magic-angle twisted bilayer graphene (MATBG) [4]. Such structures are characterized by a very high ratio of critical temperature of the superconducting transition to Fermi energy T c /E F [4] and complex phase diagrams of superconducting and insulating phases [5,6]. Experimental results for MATBG clearly show the importance of long range Coulomb electron-electron interactions in these structures [4][5][6][7]. From the theoretical side it has been shown that for small twisted angles the moiré pattern leads to a formation of a superlattice with a unit cell containing more than 10000 atoms that significantly modifies the low-energy structure. Extensive numerical studies by quantum-chemistry methods show the appearance of flat lowest-energy mini bands [8][9][10][11]. These bands are rather narrow and thus the Coulomb interactions play an important role as it was already pointed out in early theoretical studies [12]. The existence of narrow flat bands is also confirmed in recent MATBG experiments [13].
In this Letter we show that the narrow band structure of free electron spectrum leads to a number of unusual properties of their propagation in presence of strong, moderate or even weak Coulomb repulsion between electrons. We discuss these properties on a model of two electrons with Coulomb interaction propagating on a standard tight-binding lattice considered in [1][2][3]. Our results show the appearance of pairing of two electrons induced by a moderate Coulomb repulsion with ballistic pairs propagating over the whole system size. Below we describe the physical properties of such pairs of two interacting particles (TIP).
Quantum tight-binding model of two electrons. -The quantum Hamiltonian of the model in dimensions d = 1 or 2 has the standard form [1][2][3]: where j = (x 1 , x 2 ) (j = (x 1 , x 2 , y 1 , y 2 )) is a multi-index for d = 1 (d = 2); each index variable takes values x 1 , x 2 , y 1 , y 2 ∈ {0, . . . , N − 1} with N being the linear system size with periodic boundary conditions. The first sum in (1) describes the electron hopping between nearby sites on a 1D (or 2D square) lattice with a hopping amplitude taken as energy unit. The second sum in (1) represents a (regularized) Coulomb type long-range interaction with amplitude U and the distance r(j) between two electrons. For 1D we have, due to the periodic boundary conditions, r(j) = ∆x with ∆x = min(∆x, N − ∆x) and relative coordinate ∆x = x 2 − x 1 which is taken modulo N (i.e. ∆x = x 2 − x 1 + N for x 2 − x 1 < 0). For 2D we have r(j) = ∆x 2 + ∆ȳ 2 . Furthermore, we consider symmetric (spatial) wave functions with respect to particle exchange assuming an antisymmetric spin-singlet state (similar results are obtained for antisymmetric wave functions). In absence of interactions at U = 0 the spectrum of free electrons has the standard band structure E = −2 µ=1,2; α cos p µα with µ = 1, 2 being the electron index and α = x in 1D (α ∈ {x, y} in 2D) being the index for each spatial dimension. With periodic boundary conditions each momentum p µα is an integer multiple of 2π/N . Classical dynamics of electron pairs. -The corresponding classical dynamics in 2D is described by the Hamiltonian: H = −2 µ=1,2; α∈{x,y} cos p µα + U C (x 1 , x 2 , y 1 , y 2 ) (2) with U C = U/[1 + (x 1 − x 2 ) 2 + (y 1 − y 2 ) 2 ] and conjugated variables of momentum p µx , p µy and coordinates x µ , y µ (in 1D we have in (2) only p µx and x µ ). In 1D there are two integrals of motion being the total energy E = H and total momentum p + = p 1x + p 2x leading to integrable TIP dynamics. In 2D we have 3 integrals of motion E, p +x , p +y for 4 degrees of freedom and therefore the dynamics of the two electrons is generally chaotic as it is shown in Fig. S1 of [14], Sec. S1.
Writing cos(p 1x ) + cos(p 2x ) = 2 cos(p +x /2) cos[(p 2x − p 1x )/2] (and similarly for y) we see that at given values of p +x , p +y the kinetic energy is bounded by ∆E = 4 α | cos(p +α /2)|. Therefore for TIP states with E > ∆E the two electrons cannot separate and they propagate as one pair. In particular for p +x = p +y = π + δ (with |δ| 1) being close to π there are compact Coulomb electron pairs even for very small interactions U as soon as ∆E ≈ 2d|δ| < U B d with B d = 8d + U being the maximal energy bandwidth in d dimensions. The center of mass velocity of such pairs (in direction α ∈ {x, y}) is v +α = (v 1α + v 2α )/2 = 2 cos(δ/2) sin(p 1α − δ/2) ≈ 2 sin p 1α and it may be close to a maximal velocity v +α = 2. Fig. S1 of [14] clearly confirms the pair formation at small U values and the pair propagation through the whole system. Color plot of wave function amplitudes obtained from the exact quantum time evolution in 1D at times t = 138 ∆t (t = 10 5 ∆t) in left (right) panels (∆t = 1/B1 is an elementary time step of inverse band width B1 = 8 + U ). Top panels show the TIP wave function amplitude |ψ(x1, x2, t)| at x1 and x2 for both axis. Bottom panels show |ψ(p+, ∆x, t)| with ∆x = x2 − x1 (taken modulo N) being the relative coordinate (horizontal axis) and 0 ≤ p+ < 2π being the total momentum (vertical axis). The initial state at t = 0 is a symmetrized state with one electron localized at N/2 and the other one at N/2 + 1, i.e. initial distance ∆x = 1. Other parameters are U = 1 and N = 128; color bar represents (here and in certain subsequent figures) the ratio of the shown quantity (here modulus of wave function amplitude) to its maximal value. Related videos are available at [14,15].
Quantum TIP eigenstates and dynamics in 1D and 2D.
-The factor 2 cos(p + /2) for the kinetic energy (in 1D) at given value of p + can also be obtained from the quantum model. Using a suitable unitary transformation (see [14], Sec. S2) one can transform the quantum Hamiltonian (1) to a block diagonal form where the different blocks on the diagonal correspond, for each value of the conserved quantum number p + , to an effective one-particle tight binding model (in ∆x space) with nearest neighbor hopping matrix element −2 cos(p + /2) and a diagonal potential given by U/(1 + ∆x) (see [14], Sec. S2 for details, especially the boundary conditions and the generalization to the 2D case). Physically, this corresponds to a single particle in ∆x space with a kinetic energy rescaled by 2 cos(p + /2) and moving in a given potential with maximal value U for ∆x being close to 0 or N .
The eigenstates can be efficiently calculated for large system sizes since each block for a given value of p + can be diagonalized individually. Also the quantum time evolution (for = 1) can be efficiently computed (up to N = 1024) by transforming the initial state ψ(x 1 , x 2 ) into block representationψ(p + , ∆x) and computing the time evolution exactly in each p + sector using the exact block eigenstates (details in [14], Sec. S2).
A typical example of the TIP time evolution in 1D is shown in top panels of Fig. 1 for initial electron positions (localized at x µ ≈ N/2) at a distance R = ∆x = 1 for U = 1 and N = 128. We see a wave front of free propagating separated electrons (square at small times in top left panel) and at the same time free propagation of the Coulomb electron pair along the diagonal x 1 = x 2 . At large times the density for separated electrons is uniformly distributed over the whole system while the density for Coulomb pairs is homogeneously distributed over the whole diagonal keeping a relatively small pair size (top right panel). Bottom panels show the block representation |ψ(p + , ∆x, t)| of the same states as in top panels. In this representation the initial state corresponds to two red vertical lines at ∆x = 1 and ∆x = N − 1 for perfect localization at these values and a uniform distribution in p + direction. With increasing time a part of the density stays quite well localized to the initial values with some modest increase of the initial width. The other fraction of the density propagates horizontally and becomes roughly uniform at sufficiently long times. The speed of the horizontal propagation is apparently proportional to cos(p + /2) and for p + ≈ π there is actually a small regime of strong, nearly perfect localization, where the weight of the propagating density is close to 0. In this case the kinetic energy can be considered as a small perturbation of the interaction due to the small ratio 2| cos(p + /2)|/U 1. However even for larger values of this ratio there is still a considerable fraction of the density that stays localized. Fig. 2 shows certain Husimi functions of 1D block eigenstates obtained from the smoothing of the Wigner function on the scale of (see e.g. [16,17]) for two sectors p + ≈ 2π/3 and p + ≈ π. For larger energies there is a clear localization in ∆x at ∆x ≈ 0 and ∆x ≈ N where the interaction potential is maximal and for medium/lower energies there is a near ballistic movement over the full ∆x range. For p + ≈ 2π/3 there is also localization in ∆p with two values close to π in top right panel but the typical localized ∆p value may be different for other eigenstates (not shown) except for the few eigenstates with top energies which are localized in ∆x. For p + ≈ π, where the effective kinetic energy is strongly reduced, there is a larger number of states with localization in ∆x for larger energies and for all cases ∆p appears to be rather delocalized. The states with localization in ∆x contribute mostly in the quantum time evolution shown in Fig. 1 due to the initial localized state at ∆x = 1 thus explaining the significant probability for both particles staying together. Further results of the 1D quantum case are given in Figs. S2-S4 of [14], Sec. S3.
In order to access larger system sizes up to N = 512 we compute for the 2D case the exact quantum time evolution only in individual sectors with p +x = p +y (and exploiting certain symmetries [14]). The initial state is localized at ∆x = ∆ȳ = 1 but also at the chosen value of p +x = p +y . This corresponds to a wave in the center of mass direction with momenta p +x , p +y which is however perfectly localized in the relative coordinate. Fig. 3 clearly shows that there is a considerable probability that both electrons stay together even at very large times. There is however, for p +x,y = 0 or p +x,y ≈ 2π/3, also a certain complementary density for free electron propagation which becomes less visible for larger N = 512. For p +x,y ≈ π there is actually near perfect localization even at N = 128. This confirms the formation of Coulomb electron pairs in the 2D quantum case and as in 1D also the particularly strong localization in relative Color plot of wave function amplitude |ψ(p+x, p+y, ∆x, ∆y)| in block representation obtained from the 2D quantum time evolution at iteration time t = 10 5 ∆t in ∆x-∆y plane for certain p+x, p+y sectors and U = 2 (∆t = 1/B2 = 1/(16 + U )). The initial block state at t = 0 is localized at ∆x = ∆y = 1. All panels show a zoomed region 0 ≤ ∆x, ∆y < 32. The values of p+ = p+x = p+y and N are p+ = 0, N = 128 (top left), p+ = 21π/32 ≈ 2π/3, N = 128 (top right), p+ = 63π/64 ≈ π, N = 128 (bottom left), k = 85π/128 ≈ 2π/3, N = 512 (bottom right). Related videos are available at [14,15].
coordinates for sectors with p +x,y ≈ π We also computed the full space 2D quantum time evolution using the Trotter formula approximation (see e.g. [18] for computational details) with a Trotter time step ∆t = 1/B 2 = 1/(16 + U ) = 1/18 for U = 2 and for N = 128. Here we choose an initial condition with both particles localized at ≈ (N/2, N/2) such hat ∆x = ∆ȳ = 1. The results of Fig. 4 clearly show that there is a significant probability of electron pair formation with a propagation of pairs over the whole system (high probability of ∆x, ∆y at values close to zero in top panels and near the diagonal x 1 ≈ x 2 in bottom panels). At the same time there is also a certain probability to have practically independent electrons with ballistic propagation through the lattice at moderate times (bottom left panel) and an approximate homogeneous distribution over the whole system at large times (bottom right panel).
In order to characterize the pair formation probability we compute the quantum probability w 10 (t) to find both electrons at a finite distance from each other ∆x ≤ ∆R = 10 and ∆ȳ ≤ ∆R. At large times we find that the pair formation probability is w 10 = 0.133 while the probability to have independent electrons is 1 − w 10 .
The time dependence of the pair formation probability w 10 for both types of 2D quantum time evolution is shown in Fig. 5 for U = 2 and N ∈ {128, 512}. For p + ≈ π we have w 10 = 1 with numerical accuracy while for p + ≈ and N = 128, U = 2 (∆t = 1/B2 = 1/(16 + U ) is the Trotter integration time step). Top panels show the zoomed density for 0 ≤ ∆x, ∆y < 32 in the ∆x-∆y plane of relative coordinates obtained from a sum over x1 and y1. Bottom panels show the density in x1-x2 plane obtained from a sum over y1 and y2. The corresponding value of probability near diagonal w10 is w10 = 0.106 (w10 = 0.133) for left (right) panels (see text). Related videos are available at [14, 15]. 2π/3 and p + = 0 this probability decreases with time but saturates at rather high values w 10 ≈ 0.2 and 0.13 respectively. The saturation value w 10 ≈ 0.133 for the full space Trotter formula approximation is very close to the case p + = 0.
The initial state with the pair momentum p + approximately corresponds to the electron Fermi energy E F ≈ −4 cos(p + /2) assuming relatively moderate or weak interaction U B d . Thus the variation of p + in the range 0 ≤ p + ≤ 2π corresponds to the variation of the filling factor ν in the range 0 ≤ ν ≤ 1 with ν = (1 − cos(p + /2))/2. The dependence w 10 (ν) on ν is shown in Fig. 6 for U = 0.5, 2 and N = 256 (data for N = 512 is very close). We see that for small ν 1/2 the electron pair formation probability is relatively small but for 1/3 ≤ ν ≤ 1/2 it takes for U = 2 (U = 0.5) values from w 10 = 0.17 (w 10 = 0.11) to w 10 = 1. The values of w 10 are significantly above the ergodic probability w erg. = (21/N ) 2 ≈ 0.0067 for N = 256 assuming a uniform distribution.
More results for the quantum time evolution and eigenstates properties in 2D are given in [14] (see Figs. S4-S10 of Sec. S4). In particular, there is a clear scaling dependence of w 10 (and other related characteristics) on the ratio 2 cos(p + /2)/U (Figs. S9, S10 of [14]) Discussion. -The presented analytical and numerical analysis definitely shows that a specific energy dispersion law of free electrons in narrow bands leads to a new pos- -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 w 10 t/∆t N=512, p + =0 N=512, p + ≈2π/3 N=512, p + ≈ π N=128, p + ≈2π/3 N=128, Trotter sibility of pair formation induced by Coulomb repulsion between electrons. The pairing of electrons already exists for a relatively moderate, or even small, repulsion strength U . Of course, this analysis is performed in the frame work of only two interacting electrons. However, we conjecture that these Coulomb electron pairs will still exist even at finite electron density ν. Indeed, a propagating pair can feel interactions from other electrons as a certain mean field potential which should not destroy pairs with a relative strong coupling. In presence of an external space inhomogeneous potential the independent electrons can be localized by the potential while the pairs, being protected by their effective coupling energy ∆ c ∼ U/∆R, can remain insensitive to the potential propagating through the whole system (we assume that ∆R is of the order of several lattice units as for w 10 computations). Thus such pairs can form a condensate leading to the emergence of a superconducting state. We argue that the emergence of Coulomb electron pairs in narrow band structures can be at the origin of superconductivity in MATBG experiments.
Acknowledgments. -This work was supported in part by the Programme Investissements d'Avenir ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT (project THETRACOM). This work was granted access to the HPC resources of CALMIP (Toulouse) under the allocation 2019-P0110. Here, we present additional material for the main part of the article.

S1. Classical chaotic dynamics and propagation of two electrons in 2D
We show in Fig. S1 for U = 2 and U = 0.5 two typical trajectories obtained from the canonical equations with respect to the classical 2D Hamiltonian (2) using initial conditions with |x 1 − x 2 | ≈ |y 1 − y 2 | ≈ 1 and conserved total momenta p x1 + p x2 = p y1 + p y2 = π + ∆p close to π. The projection of the trajectory on the (x 1 − x 2 )-(y 1 − y 2 ) plane clearly shows a bounded and chaotic dynamics of two electrons. Obviously both electrons stay close together even though the interaction values are significantly smaller than the total energy width B 2 = 16 of two free electrons. The compact electron pair propagates quasi-ballistically through the whole system which can be seen from the global quasi-linear dependence of y 1 on x 1 and of r = x 2 1 + y 2 1 on t even though there are some loop like structure on short length scales.

S2. Numerical methods for the computation of TIP eigenstates and quantum dynamics in 1D and 2D
In this section we present some technical details about the properties of the quantum Hamiltonian (1) and how to exploit its symmetries in order to efficiently compute eigenfunctions and solve exactly the quantum time evolution of this Hamiltonian. We remind that in (1) the first sum for the kinetic energy is over nearest neighbors j and l such in both multi-indices exactly one of the four indices x 1 , x 2 , y 1 , y 2 differs exactly by ±1 (or ±(N − 1) if one index is 0 and the other N − 1 according to the periodic boundary conditions) while the other three indices are equal.
We also remind the particular notation for the relative coordinate ∆x = (x 2 −x 1 +N ) mod N (and similarly for ∆y), i.e. ∆x = x 2 −x 1 for x 2 ≥ x 1 and ∆x = x 2 −x 1 +N for x 2 < x 1 such that ∆x ∈ {0, . . . N − 1}. However, in order to correctly characterize a physical distance in the space of coordinates with periodic boundary conditions (e.g. for the interaction dependence on distance) it is the quantity ∆x = min(∆x, N − ∆x) (and similarly ∆ȳ) which is relevant since a value of ∆x being close to N : Classical dynamics for 2D case with a classical trajectory obtained by the standard Runge-Kutta 4th order method using a time step of dt = 0.005. The initial condition corresponds to |x1 − x2| ≈ |y1 − y2| ≈ 1 and conserved total momenta being px1 + px2 = py1 + py2 = π + ∆p close to π. Panels of left (right) column correspond to U = 2, ∆p = 0.1 and maximal time tmax = 1500 (U = 0.5, ∆p = 0.02, tmax = 5000). Top panels show the dependence of y1 − y2 on x1 − x2; 2nd row of panels show the propagation of y1 versus x1 and 3rd row of panels provide a zoomed view of the blue rectangles visible in 2nd row panels revealing a loop like structures of the curves; bottom panels show the time dependence of r = x 2 1 + y 2 1 .
corresponds in reality to a physical distance N − ∆x ∆x.

Discussion of 1D case
In order to keep simpler notations we will use in the remainder this section the notation k instead of p + for the total momentum (in 1D) and similarly k 1 , k 2 for p 1 , p 2 which is also more usual for momenta or wave numbers in the quantum case.
In the following we describe some technical details how to exploit the different symmetries for the 1D quantum case. If we denote by |k 1 , k 2 a non-interacting eigenstate with k µ = 2πl µ /N , µ = 1, 2 and l µ ∈ {0, . . . , N − 1} it is obvious that also in the presence of interaction the Hamiltonian (1) only couples states such the total momentum k = k 1 +k 2 is conserved leading to a block diagonal structure when expressing H in the basis of non-interacting eigenstates.
However, for numerical purposes this basis is not very convenient and we prefer to introduce a different basis of states given by: where |x 1 , x 1 + ∆x corresponds to a basis vector in position space |x 1 , x 2 with x 2 = x 1 + ∆x for x 1 , x 2 , ∆x ∈ {0, . . . , N − 1} (the sum x 1 + ∆x is to be taken modulo N ). Furthermore, k = 2πl k /N corresponds to the total discrete momentum with l k ∈ {0, . . . , N − 1}. The quantity (x 1 + ∆x/2) mod N corresponds to the center of mass (note that this quantity may be different from (x 1 + x 2 )/2 due to complications related to the periodic boundary conditions). In the following we call the states |k, ∆x block basis states. For a given state |ψ we introduce in the usual way the wave functions ψ(x 1 , x 2 ) in position space andψ(k, ∆x) in combined (k, ∆x) space (also called block representation) by : One can easily work out that the Hamiltonian has also a block diagonal structure when expressed in terms of the block basis states, i.e. the block matrix elements are given by: with symmetric N × N matricesh (k) , to be called block Hamiltonians, depending on k and with non-vanishing matrix elements: where (S6) applies to 0 ≤ ∆x < N − 1 and in (S7) f k = exp(−kN/2) = (−1) l k with l k = kN/(2π) being the integer index associated to k. U (∆x) = U/[1 + min(∆x, N − ∆x)] is the interaction potential for the 1d case. The matrixh (k) corresponds to a one-dimensional nearest neighbor tight binding Hamiltonian with hopping matrix elements: −2 cos(k/2), potential: U (∆x) and with periodic (anti-periodic) boundary conditions if l k is even (odd). This form is the quantum manifestation of the classical Hamiltonian when rewritten using the total momentum p + = p 1 + p 2 , the relative coordinate ∆x = x 2 − x 1 with associated momentum ∆p = (p 2 − p 1 )/2 : The amplitude of the kinetic energy (hopping matrix element) is proportional to 2 cos(p + /2) (2 cos(k/2)) and becomes very small for p + ≈ π (k ≈ π). At even N there is actually one precise value k = π where the hopping matrix element vanishes. Qualitatively, the eigenvectors ofh (k) at largest energies are expected to be localized at ∆x close to 0 (or N − 1) depending on the ratio C = 2| cos(k/2)|/U . If C 1 there is even a perturbative regime with strong localization and for k being sufficiently close to π this scenario is already possible for modest interaction values such as U = 2 or U = 0.5 which are below the kinetic energy bandwidth 8. Also for C ∼ 1 the eigenfunction amplitudes at small ∆x are enhanced and a state initially localized at small R = ∆x, e.g. R = 1, will retain upon quantum time evolution an enhanced probability to stay close to small ∆x (or small N −∆x) values and propagate only partially to the remaining space.
These expectations are well confirmed by Fig. 2 showing for U = 2, N = 1024, the Husimi functions, defined in classical phase space of the relative coordinate (∆x, ∆p), (see [17] for precise definition and computational details using FFT) of certain eigenvectors ϕ there is either clear localization at ∆x close to 0 or N for larger energies and more or less free propagation for lower energies. If the ratio C = 2| cos(p + /2)|/U is small the range of localized states covers a larger energy interval and the remaining delocalized states feel stronger the shape of the interaction potential as can be seen in bottom-right panel of Fig. 2. Note that in Fig. 2 we have chosen a quite large value of N = 1024 in order to have a nice spatial resolution of the Husimi functions. We have actually also computed certain Husimi functions for even larger values N = 4096 or N = 16384, corresponding to smaller effective values of , which perfectly confirms the above observations but with a reduced width of the classical lines where the Husimi function is maximal.
The block Hamiltonian (S5-S7) obviously commutes with P and using a basis of symmetrized and antisymmetrized states in ∆x space according to the above transformation rule one can transformh (k) in a diagonal 2 × 2 block structure for the two symmetric and antisymmetric sectors. The dimension of the symmetric sector is N S = (N +1+f k )/2 for even N and N S = (N +1)/2 for odd N and the anti-symmetric sector has the dimension N A = N − N S . The dependence of N S on the sign factor f k for even N is due to the fact that the basis state |∆x at ∆x = N/2 either belongs to the symmetric sector for even l k (with k = 2πl k /N ) or to the antisymmetric sector for odd l k while for ∆x = 0 it always belongs to the symmetric sector. In all other cases for 0 < ∆x < N/2 the (anti-)symmetrized basis states are (|∆x + sf k |N − ∆x )/ √ 2 with s = +1 (s = −1). The matrix elements of the (anti-)symmetrized block Hamiltonian are similar to (S5) and (S6) with some slight modifications: (i) the coupling matrix element (S6) implicating the states |0 and |N/2 (for even N ) if they belong to the sector are multiplied with √ 2, (ii) there is no corner matrix element corresponding to (S7) and (iii) for odd N one has to add −2sf k cos(k/2) to the diagonal matrix element (S5) for ∆x = (N − 1)/2.
Our aim is to numerically compute the quantum time evolution of a state |ψ(t) = e −iHt |ψ(0) with initial state |ψ(0) which we choose to be symmetric with respect to particle exchange. An exact way for this is to diagonalize H (eventually in its symmetric sector) and expand |ψ(0) in the basis of (symmetrized) eigenstates. Without any further optimization this requires O(N 6 ) operations for the diagonalization and O(N 4 ) operations for each time value t for which |ψ(t) is computed (the limitation to the symmetrized sector reduces the numerical prefactors by 8 or 4). However, using the block Hamiltonians and the block diagonal structure of (S4) we can simplify the numerical diagonalization significantly by diagonalizing N (or ∼ N/2, see below) N × N matrices (or even N/2 × N/2 matrices in the symmetric sector) which can be done with N O(N 3 ) operations (or even N O(N 2 ) operations when exploiting the tridiagonal structure of the symmetrized block Hamiltonian of each sector). Once the eigenfunctions ϕ (k) l (∆x) ofh (k) are known one can reconstruct the eigenfunctions φ (k) l (x 1 , x 2 ) of the initial Hamiltonian H by the inverse of the transformation (S3): with ∆x = (x 2 − x 1 + N ) mod N . Note that usually the inverse transformation from block to position representation for an arbitrary state would also require a sum over k. However, for eigenstates there is only one k sector with non-vanishing values and the value of k is simply fixed. From the numerical point of view the eigenstate (S9) is not convenient since it is complex (for k = 0 and k = π). Therefore we take the real and imaginary part which provides actually two eigenstates that are linear combinations of the eigenstates (S9) at k and 2π − k. This is related to the time reversal symmetry which provides the identity:h (2π−k) = Th (k) T where T is a diagonal matrix in ∆x space with non-vanishing elements T ∆x,∆x = (−1) ∆x . Due to thish (2π−k) andh (k) have the same eigenvalues and the operator T provides the transformation of the eigenvectors from the latter to the former. Therefore it is only necessary to diagonalize N/2 of the N block Hamiltonians numerically. For k = 0 and k = π (for even N ) the eigenstates (S9) are already real (or purely imaginary) and for k = π the block Hamiltonian is already diagonal since in this case the hopping matrix element (S6) vanishes.
In this way, we can very efficiently construct a basis of real eigenstates of the Hamiltonians which provides a significant acceleration of the time evolution computation. However, for larger values of N we still need a lot of storage for all eigenstates and also the transformations between the initial basis and the eigenstate basis are still quite expensive. Therefore, we apply a further optimization where we transform the initial vector |ψ(0) in block representation using (S3) and then we apply the time evolution individually for each sector of k which is highly efficient and furthermore it only requires to store the "small" block eigenstates ϕ (k) l (∆x) . The state |ψ(t) is transformed back from block to position presentation (in addition both transformations can also be accelerated by FFT). A further advantage is that we can also analyze easily the state |ψ(t) in block representation which is physically very interesting since it shows the (absence or presence of) propagation or localization individually for each k sector corresponding to different hopping ma-trix elements "−2 cos(k/2)" (see bottom panels of Fig. 1). Using this very highly efficient method we have been able to compute the exact full space 1D time evolution up to N = 1024 corresponding to a Hilbert space dimension ≈ 10 6 . We have also verified that the three variants (and further sub-variants with respect to symmetry) of the method provide identical results up to numerical precision for sufficiently small values of N where all methods are possible.

Discussion of 2D case
The block representation (S3) can be generalized to the 2D case providing wave functionsψ(k x , k y , ∆x, ∆y) depending on two total conserved momenta k x , k y in x or y-direction, ∆x and ∆y. In this case the matrix elements in block basis states provide also a block diagonal structure: k x ,k y , ∆x, ∆ỹ| H |k x , k y , ∆x, ∆y = (S10) δk x ,kx δk y ,kyh (kx,ky) (∆x,∆ỹ), (∆x,∆y) with symmetric N 2 × N 2 matricesh (kx,ky) corresponding to a 2D-tight binding model in (∆x, ∆y) space with diagonal matrix elements given by U (∆x, ∆y) = U/(1 + ∆x 2 + ∆ȳ 2 ) and nearest neighbor coupling matrix elements −2 cos(k x /2) (−2 cos(k y /2)) in x (y) direction. The boundary conditions are either periodic or antiperiodic in x (or y) direction according to the parity of the integer number N k x /(2π) (N k y /(2π)).
The two discrete symmetries ∆x ↔ N − ∆x (corresponding to x 1 ↔ x 2 ) and ∆y ↔ N −∆y allow to simplify the diagonalization problem to matrices of size ≈ N 2 /4. (Note that the particle exchange symmetry corresponds to the simultaneous application of both of these symmetries.) If k x = k y one can even exploit a third symmetry with respect to ∆x ↔ ∆y allowing a further reduction of the matrix size to ≈ N 2 /8 ≈ 3.3 × 10 4 for N = 512 which is still accessible to simple full numerical diagonalization. The additional symmetries with respect to the particle exchange symmetry are due to the particular simple form of the initial tight-binding model given as a simple sum of 1D tight-binding models for each direction (for the kinetic energy).
The details with many different cases for the parity of both k x , k y values, and also of N , together with several cases of symmetric or anti-symmetric sectors are quite complicated. For the time evolution of the 2D case we limit ourselves to a single sector with k = k x = k y , even N and also to the totally symmetric case (with respect to the three exchange symmetries mentioned above). In this case each (totally symmetrized) block Hamiltonian has a dimension D k = (N +1+f k )(N +3+f k )/8 ≈ N 2 /8 where f k = (−1) l k is the sign factor associated to k. As initial state we use a totally symmetrized state localized at ∆x = ∆ȳ = R where we choose typically R = 1 but we have also performed certain computations for R = 3 (see e.g. Fig. S8 below). In the original full 4D space of TIP such a state corresponds to a plane wave in the center of mass direction and perfect localization in the relative coordinate. Therefore, the initial spreading along the diagonal (center of mass coordinate) is not visible with such a state (contrary to the full space 1D time evolution where an initial state localized in the center of mass was used as can be seen in top panels of Fig. 1). However the more important spreading in the relative coordinate, which determines the pair formation probability, is of course clearly visible.
In addition to the exact block 2D time evolution we also computed the full space time evolution using the Trotter formula approximation: which is valid for a sufficiently small Trotter time step ∆t and time values t being integer multiples of ∆t. Here H kin represents the kinetic energy part of the Hamiltonian (diagonal in Fourier space) and H pot represents the part with potentials and interactions (diagonal in initial position space). Using a 4D-FFT to transform efficiently between initial position space and Fourier space one can compute the approximate time evolution of an arbitrary given initial state. Further technical details of this method for the case of 2D TIP can be found in [18].
As Trotter time step we choose the inverse bandwidth ∆t = 1/B 2 = 1/(16 + U ) (which is below ∆t = 0.1 used in [18]). Also for the exact time evolution (in 1D and 2D) we measure/present all time values in units of the inverse bandwidth ∆t = 1/B d which corresponds to the smallest time scale of the system. However, for the latter the ratio t/∆t may be arbitrary and is not limited to integer values as for the Trotter formula approximation. Concerning the Trotter formula time evolution we choose a system size of N = 128 (as in [18]) and an initial state being roughly localized at ≈ N/2 for each of the four coordinates with exact initial distance ∆x = ∆ȳ = R and totally symmetrized with respect to the above mentioned three discrete symmetries.
The considered time range for all types of quantum time evolution computations is ∆t × 10 −1 ≤ t ≤ ∆t × 10 6 (∆t ≤ t ≤ ∆t×10 4 ) for the exact 1D/2D (Trotter formula 2D) time evolution using an (approximate) logarithmic scale for the chosen time values. We provide here two example videos and at [15] further videos (see Sec. S5 for details) for several of the densities shown in Figs. 1,3,4 and S7. Color plots of certain symmetric eigenstates ψ(x1, x2) of the 1D Hamiltonian for U = 1 and N = 128 in x1-x2 plane. Top panels correspond to two states in the sector p+ = 0 of total momentum with relative level number (of this sector and with states ordered by increasing energies) being l = 64 (left) and l = 44 (right). Center panels correspond to p+ = 21π/32 ≈ 2π/3 with l = 64 (left) and l = 43 (right). Bottom panels correspond to p+ = 63π/64 ≈ π with l = 62 (left) and l = 43 (right).

S3. Additional data for the 1D quantum case
In Fig. S2 we show some examples of 1D eigenstates of the quantum Hamiltonian (1) obtained from (the real or imaginary part of) Eq. (S9) as explained in the last section from symmetrized block eigenstates (with 65 or 64 levels per symmetrized p + sector). The parameters are U = 1, N = 128 and three values of k = p + = 0, p + ≈ 2π/3 and p + ≈ π. For states with (near) maximal energy (for their respective p + sector) one sees a strong localization on the diagonal indicating a strong localization in the relative coordinate ∆x. For relative level numbers close to 2/3 of the maximal level number the states extend to the full space except for p + ≈ π where a small strip close to the diagonal is excluded. This is a clear effect of energy miss-match and the very strong interaction on the diagonal if compared to the strongly reduced kinetic energy. The (presence or absence of) change of colors in the center of mass direction indicate the periodic oscillations due the wave number p + in this direction. Fig. S3 shows a schematic representation of all symmetrized 1D block eigenstates of certain p + sectors (for U = 1, N = 128). The horizontal axis of each panel corresponds to 0 ≤ ∆x N/2 and the vertical axis to the level number 0 ≤ l 64 (bottom/top for lowest/largest sector energies) of the shown p + sectors. For p + = 0 and p + ≈ 2π/3 there is at top energies a very small number of states with rather strong localization at ∆x close to 0 (values ∆x close to N are not visible due to the symmetrization). Then there is a big majority of states which are completely delocalized on the full ∆x interval with certain wavelength values depending on the level number. At lowest energies there is also a very small number of states with significant localization close at ∆x ≈ N/2 which is the minimum of the 1D interaction potential. For p + ≈ π but different from π the top (bottom) energy ranges of strong localization at ∆x ≈ 0 (∆x ≈ N/2) are strongly enhanced and the intermediate "delocalized" states are quite strongly excluded from small ∆x values due to an energy miss-match caused by the very strong relative interaction (if compared to the reduced kinetic energy). At p + = π exactly there is perfect localization at single ∆x sites for each state since in this case the kinetic energy vanishes exactly and the corresponding block Hamiltonian is already diagonal with eigenvalues given by the interaction values U (∆x).
These observations of Fig. S3 explain quite clearly the behavior of the time evolution of the wave function in block representation visible in the bottom panels of Fig. 1. For sectors with p + ≈ π an initial block state localized at ∆x = 1 will have only contributions from the strongly localized block eigenstates at top energies while other extended block eigenstates have nearly no amplitude at small ∆x and do not contribute in the ini-tial state (when the latter is expanded in terms of block eigenstates). Therefore there is nearly perfect localization for these sectors (with p + ≈ π) also for long times. For other values of p + quite different from π the number of localized top energy states is reduced. However, they still produce an enhanced contribution in the eigenvector expansion of the initial state. Now (most of) the other extended block eigenstates are not excluded/forbidden at small ∆x values and they have also a certain contribution in the eigenvector expansion, but still with smaller coefficients than the localized top energy states. Therefore in the time evolution a quite significant fraction of probability stays localized while the complementary fraction of density propagates through the whole system. To characterize the quantum probability to form compact electron pairs we compute (for the 1D case) the following quantity: where the sum x 1 + ∆x is taken modulo N and ψ(x 1 , x 2 ) represents the wave function (in position representation) of the quantum state for which we want to compute w 10 , typically a state obtained from the exact 1D quantum time evolution. This quantity is just the quantum probability to find both particles in the relative strip ∆x ≤ 10. One can also compute w 10 from the wave functionψ(k, ∆x) in block representation by sums over all k values and ∆x values such that either ∆x ≤ 10 or ∆x ≥ N − 10. Fig. S4 shows the time dependence of w 10 using the same (type of) time evolution data used for Fig. 1 for different parameter combinations with N = 128, 1024, U = 1, 6 for R = 1 or U = 1 for R = 3 (here R is the initial distance of the particles). We see in all cases that w 10 saturates at long times at a finite value which is significantly above its ergodic values w erg. = 21/N (maybe except for N = 128, U = 1 and R = 3) assuming a perfectly uniform wave function amplitude. Increasing R = 1 to R = 3 clearly reduces the pair formation probability while increasing the interaction from U = 1 to U = 6 increases the value of w 10 . For the data at N = 128 there is also a considerable finite size effect.
S4. Additional data for the 2D quantum case Before presenting some additional data for the 2D quantum case, we provide here explicit formulas for the quantities shown in Figs. 4, 5. Assuming that ψ(x 1 , y 1 , x 2 , y 2 ) is the wave function in position space obtained from the Totter formula 2D time evolution top panels of Fig. 4 show the density in ∆x-∆y plane obtained from: |ψ(x 1 , y 1 , x 1 +∆x, y 1 +∆y)| 2 (S13) (with position sums taken modulo N ) and bottom panels show the density in x 1 -x 2 plane obtained from: |ψ(x 1 , y 1 , x 2 , y 2 )| 2 . (S14) The quantum probability w ∆R to find both electrons at ∆x ≤ ∆R and ∆ȳ ≤ ∆R is computed from: w ∆R = ∆x≤∆R,∆ȳ≤∆R ρ rel (∆x, ∆y) (S15) (with ∆x = min(∆x, N − ∆x) and similarly for ∆y). For the case of exact 2D block time evolution w ∆R is computed as in (S15) but with ρ rel (∆x, ∆y) replaced by |ψ(p +x , p +y , ∆x, ∆y)| 2 (at certain given values of p +x , p +y ). Here we mostly use the quantity w 10 corresponding to ∆R = 10 but below we also show some data for w 2 corresponding to ∆R = 2. In Fig. S5, we show for the same cases and raw data (same 2D time evolution wave functions) used for Fig. 4 the time dependence of the inverse participation ratio defined by: (ρ rel,sym is the symmetrized density with respect to the two symmetries ∆x ↔ N − ∆x and ∆y ↔ N − ∆y obtained from ρ rel ) and the typical particle distance: r typ = exp[ ln( ∆x 2 + ∆ȳ 2 + 1) ] − 1 (S17)   where the average · · · is done with respect to the symmetrized density ρ rel,sym . The inverse participation ratio provides roughly the number of sites over which ρ rel,sym (∆x, ∆y) is mostly concentrated or localized. The saturation values at long times well below the ergodic value clearly indicate a considerable pair formation probability confirming the findings of Fig. 5. For p + ≈ π (case of near perfect localization) one can see values slightly above unity (ξ IPR 1.2) for all time scales. Note that for this case w 10 shown in Fig. 5 takes the value 1 with numerical precision (error below 10 −14 ).
The other quantity r typ is mostly dominated by the fraction of density associated to free electron propagation. Therefore the ratio r typ /r typ,erg corresponds to the complementary non-pair formation probability and indeed the quantity 1 − r typ /r typ,erg behaves qualitatively similarly as w 10 . For the case of near perfect localization at p + ≈ π we have for all time scales r typ ≈ √ 2 (corresponding to ∆x = ∆y = 1) with some slight fluctuations ±0.02 at long times.
Globally, we believe that w 10 is more useful and suitable than ξ IPR or r typ to describe the pair formation probability which is the reason why we have chosen to show in the main part of this work the time dependence of w 10 in Fig. 4. Some examples of localized and delocalized 2D block eigenfunctions for the three sectors p + = p +x = p +y = 0, p + ≈ 2π/3 and p + ≈ π and N = 128, U = 2 are presented in Fig. S6. As in 1D top energy eigenstates (with respect to their p + sector) are localized at ∆x ≈ ∆y ≈ 0 and lower energy eigenstates are delocalized. For p + = 0 and p + ≈ 2π/3 the number of localized states is very small and already at sector level numbers slightly below their maximal value delocalization sets in. For p + ≈ π the number of localized states is larger and for the delocalized states the zone close to ∆x ≈ ∆y ≈ 0 is forbidden such that these states do not contribute in the time evolution if the initial state is localized at ∆x = ∆y = R for some small value of R = 1 or R = 3. These findings are very similar to the 1D case (see above discussion of Figs. S3,S4) and help to understand the 2D time evolution results shown in Figs. 3-5.
One can also see, especially for the delocalized eigenstates, a spatial node structure indicating a rough and very approximate angular momentum conservation. However, this is not exact due to the square form of available space, discrete lattice structure and also the cosinus 2D band form of kinetic energy (for the relative coordinates). Furthermore, our limitation to totally symmetrized block eigenstates (with respect to the discrete symmetries explained in Sec. S2) plays a role here.  Fig. 3 with p+ = p+x = p+y = 0 (top left), p+ = 85π/128 ≈ 2π/3 (top right), p+ = 255π/256 ≈ π (center left), all three for N = 512 and t = 10 5 ∆t, and p+ = 63π/64 ≈ π, N = 128, t = 3090 ∆t (center right). Bottom panels correspond to the Trotter formula approximation for N = 128 and t = 9487 ∆t with left panel showing the density in relative coordinates obtained from a sum over x1 and y1 and with right panel showing the density in x1-x2 plane obtained from a sum over y1 and y2. All panels (except bottom right) show to a zoomed region 0 ≤ ∆x, ∆y < 32 in ∆x-∆y plane. At [15] videos corresponding to the cases of all panels of this figure (and also for N = 128 with p+ = 0 and p+ = 21π/32 ≈ 2π/3) are provided.
The 2D pair formation effect still exists at the smaller interaction value U = 0.5 even though it is less pronounced as can be seen in Fig. S7 showing certain wave function densities obtained from both types of 2D time evolution methods. For N = 128 the Trotter formula results indicate the formation of a very slight diagonal in the density ρ XX (x 1 , x 2 ) and a modest localization effect in ρ rel (∆x, ∆y). However, here the density weight corresponding to free electron propagation is quite high. For N = 128 this is also confirmed by the 2D block time evolution results (not shown in Fig. S7 except for one special case) but for N = 512 the density related to free propagation is quite strongly reduced in comparison to N = 128. As can be seen in both top and left center panels of Fig. S7 for N = 512 there is a considerable density for the compact electron pair. However, for p + = 0 (p + ≈ 2π/3) one still clearly see a quite strong (quite modest) density related to free propagation while for p + ≈ π the latter is absent. The right center panel shows the case N = 128 and p + ≈ π (strong localization case at N = 512) at a certain intermediate time when one can briefly see a density related to free propagation which is however smeared out at longer times. The qualitative findings of Fig. S7 for U = 0.5 are confirmed by the time dependence of w 10 shown in the lower panel of Fig. S8 for the cases of Fig. S7. For N = 512 and p + = 0 or p + ≈ 2π/3 the saturation values of w 10 are below 0.1 but still significantly above the ergodic value. For p + ≈ π and N = 512 the saturation value w 10 = 0.99988 is nearly maximal while for N = 128 it is lower but still quite elevated, w 10 ≈ 0.45, showing a certain finite size effect for this case. The saturation of the Trotter formula case at N = 128 is w 10 ≈ 0.05 which is rather low but still roughly twice the ergodic value which is 0.027. The top panel of Fig. S7 shows the time dependence of w 10 for the Trotter formula case and other parameter values U = 0.5, 2, 12 at R = 1 and U = 2 at R = 3 with R being the initial value of ∆x = ∆ȳ. As in 1D (see Fig. S4) stronger (lower) interaction increases (reduces) the pair formation probability and when R is increased the pair formation probability is reduced. In Fig. S9 we show for N = 256 the dependence of the saturation value of w ∆R for ∆R = 2 (top panel) and FIG. S10: Dependence of the 2D inverse participation ratio ξIPR in symmetrized ∆x-∆y plane on the scaling parameter 2 cos(p+/2)/U . The dashed grey line corresponds to the ergodic value ξIPR,erg. = (N/2) 2 ≈ 1.6 × 10 4 assuming a uniform distribution on a square of size N/2 = 128 (due to symmetrization). ξIPR has been computed from the same raw data (time dependent wave functions), using the same parameters and same time average as in Fig. S9. ∆R = 10 (bottom panel), computed as the average over 21 time values t with 10 4 ≤ t/∆t ≤ 10 6 , on the scaling parameter C = 2 cos(p + /2)/U (here we only consider 0 ≤ p + ≤ π such that C ≥ 0). This parameter is the ratio of the hopping matrix element and the interaction amplitude of the effective 2D one-particle tight binding model of the block at p + = p +x = p +y (see Sec. S2 for details). We have also computed the same scaling curves for N = 512 with a reduced density of data points. These curves are identical to (slightly below) the case N = 256 for C 2 (C 2) with a rough factor of 0.5 at C = 8.
The scaling dependence is quite obvious and clearly confirms the role of the total conserved momentum p + for the kinetic energy amplitude of the relative coordinates. At small values of C 1 we have perfect localization (perfect pair formation) with w 2 ≈ w 10 ≈ 1 and then the curves decrease down to a modest local minimum at C ≈ 0.2 which is followed by a local maximum at C ≈ 0.3. Then they continue to decrease down to the value of 2.2 × 10 −3 (1.8 × 10 −2 ) for w 2 (w 10 ) at C = 8. These minimal values are still clearly above the ergodic values 3.8 × 10 −4 (6.7 × 10 −3 ). Fig. S10 shows the same type of scaling curve but for the inverse participation ratio ξ IPR which has a value of ≈ 1 for C 1 for the case of perfect localization (perfect pair formation). Then the curve increases up to a modest local maximum at C ≈ 0.2 followed by a local minimum at C ≈ 0.3. After this it continues to increase until it saturates for C 3 at the value ξ IPR ≈ 6.9 × 10 3 which is clearly below the ergodic value 1.6 × 10 4 .

S5. Videos of TIP time evolution in 2D
Two videos of the 2D time evolution for the cases of (i) Fig. 3 (videoforfig3.avi), top right panel (exact 2D block time evolution of sector p + = 21π/32 ≈ 2π/3) and (ii) Fig. 4 (videoforfig4.avi), bottom panels (density in x 1 -x 2 plane obtained from the Trotter formula 2D time evolution) are provided. For (i) the video is composed of 702 images (25 images per second of video) at time values t 0 = 0 and t i = ∆t × 10 (i−101)/100 for 1 ≤ i ≤ 701 (with ∆t = 1/B 2 = 1/(16 + U )) corresponding to a logarithmic time scale in the range ∆t×10 −1 ≤ t ≤ ∆t×10 6 . For (ii) the video is composed of 464 images at integer multiples of t i = l i ∆t with l i = i for 0 ≤ i ≤ 87 and roughly l i+1 ≈ 10 1/100 l i for 87 ≤ i < 463 (such that l 463 = 10 4 ) corresponding roughly to a logarithmic (linear) time scale for t/∆t > 87 (t/∆t ≤ 87) in the range ∆t ≤ t ≤ ∆t × 10 4 .