Taming two interacting particles with disorder

TL;DR

This paper investigates how two interacting particles in a disordered one-dimensional chain exhibit localization properties, revealing how interaction strength and disorder influence localization length and the transition to asymptotic scaling regimes.

Contribution

It provides large-scale numerical analysis of two-particle localization, identifying the conditions for maximum enhancement and the crossover to asymptotic scaling behavior.

Findings

Maximum localization length enhancement occurs at specific interaction strengths and energies.

A crossover to asymptotic scaling occurs for large single-particle localization lengths.

A nonlinear scaling function describes the relationship between two-particle and single-particle localization lengths.

Abstract

We compute the scaling properties of the localization length $\xi_2$ of two interacting particles in a one-dimensional chain with diagonal disorder, and the connectivity properties of the Fock states. We analyze record large system sizes (up to $N=20000$) and disorder strengths (down to $W=0.5$). We vary the energy $E$ and the on-site interaction strength $u$. At a given disorder strength the largest enhancement of $\xi_2$ occurs for $u$ of the order of the single particle band width, and for two-particle states with energies at the center of the spectrum, $E=0$. We observe a crossover in the scaling of $\xi_2$ with the single particle localization length $\xi_1$ into the asymptotic regime for $\xi_1 > 100$ ($W < 1.0$). This happens due to the recovery of translational invariance and momentum conservation rules in the matrix elements of interconnected Fock eigenstates for $u=0$. The entrance into the asymptotic scaling is manifested through a nonlinear scaling function $\xi_2/\xi_1=F(u\xi_1)$.

Equations

$$H_0=\sum _n\left[t(|n⟩⟨n+1|+|n+1⟩⟨n|)+V_n|n⟩⟨n|\right],$$

$$H=H_0\oplus H_0+uP.$$

$$I_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}=\sum _n{\varphi }_{{\mu }^{\prime }}^{\ast }(n){\varphi }_{{\nu }^{\prime }}^{\ast }(n){\varphi }_{\mu }(n){\varphi }_{\nu }(n).$$

$$R_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}\equiv \left|\frac{u{I}_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}}{E_{{\mu }^{\prime }}+E_{{\nu }^{\prime }}-E_{\mu }-E_{\nu }}\right|>1.$$

$$\tilde{G}=\frac{{\tilde{G}}_0}{1-u{\tilde{G}}_0},$$

$$\frac{1}{{\xi }_2}=-\underset{|n-m|\to \infty }{\lim }\frac{\overline{\ln |⟨n,n|\tilde{G}|m,m⟩|}}{|n-m|}$$

$$⟨n,n|G_0(E)|m,m⟩=\sum _{\mu ,\nu }\frac{{\varphi }_{\mu }(n){\varphi }_{\nu }(n){\varphi }_{\mu }^{\ast }(m){\varphi }_{\nu }^{\ast }(m)}{E-E_{\mu }-E_{\nu }}.$$

$$⟨n,n|G_0(E)|m,m⟩=\sum _{\nu }{\varphi }_{\nu }(n)g_{nm}(E-E_{\nu }){\varphi }_{\nu }(m),$$

$$g_{nm}(E)=\sum _{\mu }\frac{{\varphi }_{\mu }(n){\varphi }_{\mu }(m)}{E-E_{\mu }}.$$

$$\text{PN}=\frac{{\xi }_2}{N}\sum _{l=1}^{N/{\xi }_2}P(n\equiv l{\xi }_2),P(n)=\frac{({\sum }_m|{\tilde{G}}_{n,m}|{)}^2}{{\sum }_m|{\tilde{G}}_{n,m}{|}^2}.$$

$$y(x)=ln\overline{⟨n,n|\tilde{G}|m,m⟩}$$

Full text

localization; interaction

I Introduction

The fundamental question about the interplay of Anderson localization Anderson (1958) and many-particle interactions has driven both analytical and numerical studies for decades. Fleishman and Anderson (1980); Shahbazyan and Raikh (1996); Kozub et al. (2000); Wang et al. (2000); Nattermann et al. (2003); Basko et al. (2006) In particular, celebrated many body localization transitions for a macroscopic system with a finite particle density have been predicted using perturbation theory approaches. Basko et al. (2006) It is most natural then to turn attention on a seemingly simple problem of two interacting particles (TIP) in a one-dimensional disordered tight binding chain (with hopping strength $t=1$) which has both ingredients: localization and interaction. While there is little doubt that in one spatial dimension the two particles stay localized for any onsite (and in general short-range) interactions at finite disorder strength, there are conflicting results on how the localization length $\xi_{2}$ of the most extended TIP states scales with the single particle localization length $\xi_{1}$ in the limit of weak disorder $\xi_{1}\gg 1$. Predictions and numerical conclusions for $\xi_{2}\propto\xi_{1}^{\alpha}$ range from $\alpha=1$ Römer and Schreiber (1997); Römer et al. (1999) over $\alpha\approx 1.6$ [10] to $\alpha=2$. Shepelyansky (1994); Imry (1995) Ponomarev et al suggested that the scaling will be modified with a nonuniversal exponent $\alpha(u)$ depending on the strength of interaction $u$, Ponomarev and Silvestrov (1997) and discussed logarithmic corrections $\alpha=a+b\ln(\xi_{1})$. Even with only two particles, the computational task turns difficult since weak disorder values are targeted, and the required system size increases rapidly with decreasing disorder strength. The asymptotic scaling sets in for $\xi_{1}\geq 100$ when momentum conservation correlations and translational invariance begin to be restored in the single particle eigenstates. Krimer and Flach (2015) Yet most numerical studies focused on the more accessible region $\xi_{1}<100$ and therefore yield at best $\xi_{2}\approx 2...3\xi_{1}$. Shepelyansky (1994); Frahm et al. (1995); von Oppen et al. (1996); Frahm (1999); Arias et al. (1999); Waintal et al. (1999); Krimer et al. (2011) A recent Green’s function computation by Ref. Frahm, 2016 entered the scaling regime with reaching $\xi_{1}\approx 400$ and $\xi_{2}\approx 9\xi_{1}$. However, at the chosen interaction strength unavoidable finite size corrections Song and Kim (1997) will bring this number down to $\xi_{2}\approx 6\xi_{1}$ according to our present computations. In addition a completely overlooked impact comes from the interplay of interaction strength $u$ and the eigenstate energy $E$. Also TIP have been studied recently on a three-dimensional lattice in presence of a mobility edge demonstrating a sensitive dependence on the interaction strength. Stellin and Orso (2019)

In this work, we show that in the regime of asymptotic scaling $\xi_{1}\geq 100$ the connectivity between Fock states (non-interacting eigenstates) is strongly selective due to combining energy conservation with emerging momentum conservation. In order to computationally assess $\xi_{2}$ in the asymptotic scaling regime, we extend the projected Green’s function method used in Ref. von Oppen et al., 1996 by adding a finite size scaling and significantly increasing the system size compared to the data presented in the literature by more than one order of magnitude up to $N=20000$, and by systematically varying the energy $E$ and the interaction strength $u$. We show that the largest values of $\xi_{2}$ at given $\xi_{1}$ are obtained for $E=0$ and $u\approx 3$. We report the record values $\xi_{2}\approx 16\xi_{1}$. The entrance into the asymptotic scaling is manifested through a non-linear scaling function $\xi_{2}/\xi_{1}=F(u\xi_{1})$.

II Model

The single particle Anderson Hamiltonian in one space dimension is given by

[TABLE]

where $|n\rangle$ denotes a basis state with one particle located on site $n$, $t$ denotes the nearest neighbor hopping which we fix to unity in all computations $t=1$, and $V_{n}$ is an onsite potential, sampled from a uniform distribution $[-W/2,W/2]$. $W$ characterizes the strength of disorder in the system. The eigenenergies $\{E_{\mu}\}\in[-2t-\frac{W}{2},2t+\frac{W}{2}]$ and the eigenfunctions $\{\phi_{\mu}(n)\equiv|\mu\rangle\}$ of the single particle problem (1) are obtained by diagonalising $H_{0}$. The spectrum is symmetric around $E=0$ and has width $\Delta_{1}=4t+W$. The single particle localization length $\xi_{1}$ controls the exponential decay of an eigenfunction $\phi_{\mu}(|n|\rightarrow\infty)\sim{\rm e}^{-|n|/\xi_{1}}$. For $W<4$ the largest localization length $\xi_{1}$ at $E=0$ is well approximated by $\xi_{1}\approx 100t^{2}/W^{2}$.

Two indistinguishable particles are described with the basis states $|n,m\rangle$, where $n$ and $m\geq n$ stand for their coordinates. The two interacting particles (TIP) Hamiltonian is then defined with the use of the basis states $|n,m\rangle$ as:

[TABLE]

The operator $P\ket{n,m}=\delta_{nm}\ket{n,m}$ projects a two particle state onto the basis states with doubly occupied sites. The parameter $u$ controls the strength of the interaction. Since we address states with the largest localization length in the center of the spectrum with only two particles involved, neither the type of interaction (repulsive with positive $u$ or attractive with negative $u$) nor the particle statistics will influence the quest for the asymptotic scaling. We use bosons for convenience.

III Fock space connectivity

We first evaluate the interaction-induced connectivity in Fock space in order to establish the asymptotic scaling regime parameters, and in order to estimate the relevant energy scales. A Fock state $\left(\mu,\nu\right)\equiv|\mu\rangle\otimes|\nu\rangle$ is an eigenstate of the two particle system for $u=0$ with eigenenergy $E_{\mu,\nu}\equiv E_{\mu}+E_{\nu}$. The interaction induces a matrix element $\langle\mu,\nu|uP|\mu^{\prime},\nu^{\prime}\rangle=uI_{\mu,\nu}^{\mu^{\prime},\nu^{\prime}}$ between two Fock states with the overlap integral

[TABLE]

The energy difference between the chosen two Fock states is given by $\Delta E=E_{\mu^{\prime}}+E_{\nu^{\prime}}-E_{\mu}-E_{\nu}$. A strongly connected pair of Fock states is found if the ratio Krimer and Flach (2015)

[TABLE]

Krimer et al showed in Ref. Krimer and Flach, 2015 that the overlap integrals turn from random like for $W>1$ to selective ones in accord with the restoration of translational invariance and the corresponding momentum conservation for $W<1$, thus invalidating the analytical considerations in Refs. Shepelyansky, 1994; Imry, 1995 in accord with earlier predictions by Ponomarev et al. Ponomarev and Silvestrov (1997)

Here we follow the computational method of Ref. Krimer and Flach, 2015. In a nutshell, we choose a reference single particle state $|\mu_{0}\rangle$ with $E_{\mu_{0}}\approx 0$ and reference Fock state $|\mu_{0},\nu_{0}\equiv\mu_{0}\rangle$. We find all single particle states $\left\{|\nu\rangle\right\}$ whose coordinates $x_{\nu}=\sum_{n}n\phi_{\nu}(n)$ are within the range of one localization length distance from $x_{\mu_{0}}$, i.e. $|x_{\nu}-x_{\mu_{0}}|\leq\xi_{1}$. Each set is sorted with ascending energy $E_{\nu}$ with the convention $\nu_{0}=0$. In the limit of weak disorder the index $\nu$ will be related to a wave number. Each pair $|\mu,\nu\rangle$ corresponds to a new Fock state which may have significant overlap $R>1$ with the reference state. The resulting disorder averaged (5000 realizations) distribution of $R$ is plotted versus $\nu/\xi_{1},\mu/\xi_{1}$ in Fig.1 for three different values of $W=0.5,2,4$. The broad distribution for $W=4$ in Fig. 1(c) is replaced by one major thin resonance line $\nu=\mu$ for $W=0.5$ in Fig. 1(a). The restoring of translational invariance  Frahm (2016); Jacquod et al. (1997) leads to a momentum conservation related selection rule of plane waves with fixed boundary conditions Krimer and Flach (2015) and approximately reads $\nu_{0}+\mu_{0}=\nu+\mu$ which results in $\nu=-\mu$ with our reference state choice.

In addition the largest coupling strength ratio $R$ is obtained for the smallest energy differences $\Delta E$ which imply $E_{\mu}=-E_{\nu}$, resulting again in the condition $\nu=-\mu$ due to the particle-hole symmetry of the single particle spectrum. As a result, for weak disorder $W<1$ the asymptotic connectivity regime of Fock states is observed. In this regime, energy and momentum are approximately conserved, and only a selected set of about $\xi_{1}$ Fock states with $|\nu,-\nu\rangle$ is strongly connected. Their overlap integrals can be approximated replacing $\phi_{\nu}(n)\approx{\rm e}^{i\nu n}/\sqrt{\xi_{1}}$ by plane waves normalized to a box of size $\xi_{1}$. It follows $I\approx 1/\xi_{1}$. The Fock state energies of that group are confined to an interval of width $\Delta_{1}/\xi_{1}$ due to energy conservation, and result in a level spacing $\delta_{2}\approx\Delta_{1}/\xi_{1}^{2}$. We conclude that the selected set of strongly interacting Fock states is characterized by an effective disorder $W_{eff}\approx\Delta_{1}/\xi_{1}^{2}$ and an effective hopping $t_{eff}\approx u/\xi_{1}$. A naive use of the localization length estimate for a corresponding tight binding chain would result in $\xi_{2}\sim 100\frac{u^{2}}{\Delta_{1}^{2}}\xi_{1}^{2}$.

IV Measuring $\xi_{2}$

IV.1 Green function method: Benchmarking the non-interacting case

To proceed we compute the localization length $\xi_{2}$, following Ref. von Oppen et al., 1996 where it was shown that the full Dyson equation for the two-particle Green’s function (GF) $G=(E-H)^{-1}$ can be solved for basis states with double occupancy using the non-interacting $u=0$ GF $G_{0}=(E-H_{0})^{-1}$ as

[TABLE]

where $\tilde{G}=PGP$ and $\tilde{G}_{0}=PG_{0}P$ are projections of the full and non-interacting GFs onto double occupancy sites. This remarkable result allows to assess the wave function of the interacting system on double occupied sites through solving the non-interacting eigenvalue problem, and thus allows to measure the TIP localization length enhancement. The localization length $\xi_{2}$ is then defined as the exponent of the exponential decay of $\tilde{G}$:

[TABLE]

where $\overline{...}$ denotes the disorder average. To obtain the projected GF $\tilde{G}_{0}$ we solve the single particle eigenvalue problem and compute

[TABLE]

The complexity of this expression if $O(N^{4})$ which is prohibitive for large sizes $N$ that we seek to explore. It was shown in Ref. Frahm, 1999 that this complexity reduces to $O(N^{3})$ if one exploites the tridiagonal structure of the single particle Hamiltonian (1). Namely one reorders the summation in Eq. (7) as

[TABLE]

where the single particle Green function $g_{nm}(E-E_{\mu})$ is evaluated using a fast dedicated inversion algorithm for tridiagonal matrices.

The subsequent processing of the data is organised in two steps:

Step 1: disorder averaging and real space fitting. For a given system size $N$ we use GF data for $0.05N\leq n,m\leq 0.95N$ to avoid possible boundary effects and define a function $y(x)\equiv\ln|\langle n,n|\tilde{G}|m,m\rangle|$ where $x\equiv|n-m|$. In method M1 we compute the disorder averaged $\overline{y(x)}$ and then find the best linear fit $a-x/\xi_{2}$ of the TIP localization length $\xi_{2}(N)$. In method M2 we first fit $y(x)$ with a linear function $a-x/\xi$ and compute $\xi_{2}(N)=\overline{\xi}$ as the disorder average of $\xi$. For a given system size $N$, the number of disorder realisations $S$ is chosen such that $NS\approx 10^{6}$ with $N$ going up to $N=20000$.

We benchmark both methods for $u=0$. Figure 2 shows $\xi_{2}(N)$ for $W=0.5,0.75,1.0,1.5,2$ at $u=0$ obtained from both methods M1 and M2. We find that both methods agree once $N\geq N_{*}(W)$ with $N_{*}(W)\approx 6000,2000,1500,750,500$ for $W=0.5,0.75,1,1.5,2$ respectively. This implies a power-law dependence $N_{*}(W)\propto W^{-1.72}$. For larger disorder, $W>2$, the value of $N_{*}(W)$ is effectively zero for any reasonable system size. Importantly, we conclude that we need system size $N>2000$ for disorder strengths $W=0.75,0.5$ for the outcomes of both methods to coincide. These system sizes were not addressed in the literature before.

Step 2: finite size fitting for $N\geq N_{*}(W)$. Using $\xi_{2}(N)=\xi_{2}+a/N$ (see e.g. Ref. Song and Kim, 1997) we extract $\xi_{2}(N\rightarrow\infty)$ using data from the both methods M1, M2. The results are presented in Fig. 3. We observe that the fitting is insensitive to the choice of the methods M1 and M2 as well as to the number of data points (system sizes $N$) used to do the fitting as long as $N\geq N_{*}(W)$. Therefore we use the M1 method in the subsequent analysis for non-zero $u$. Further details about fits and extrapolations are discussed in the Appendix.

According to Ref. von Oppen et al., 1996 $\xi_{2}/\xi_{1}=0.5$ at $u=0$. Our data yield up to $10\%$ error with $\xi_{2}/\xi_{1}\approx 0.56$ for the weakest disorder $W=0.5$. This has to be compared against the $30\%$ error obtained for smaller system sizes in Ref. Frahm, 2016 ($\xi_{2}/\xi_{1}\approx 0.65$ for $W=0.5$).

IV.2 Non-zero interactions

To proceed to nonzero interactions, we compute $\xi_{2}/\xi_{1}$ as a function of $u$ for different values of the energy $E$ and disorder strength $W$ as shown in Fig. 4. For fixed $E$ and $W$ we find that $\xi_{2}$ passes through a maximum $u_{max}$ upon increasing $u$ and decreases with further increasing $u$. This decreasing for $u\geq\Delta_{1}$ is due to doubly occupied site states being tuned out of the spectrum of the remaining two-particle continuum, see Ref. Ponomarev and Silvestrov, 1997 for details. The value of $u_{max}\approx 3t$ is in agreement with the prediction from Ref. Ponomarev and Silvestrov, 1997 as well. In addition we also observe an anomaly Frahm (2016) at $u\approx t$ and $E\approx 1$ for $W\leq 1$ where the localization length $\xi_{2}$ is enhanced as compared to $E=0$. However, the record value of $\xi_{2}$ for each studied disorder strength is found for the band center $E=0$ and the interaction strength $u_{max}\approx 3t$. As shown above, the resonantly interacting Fock state groups contain $\sim\xi_{1}$ member states and are characterized by two energy scales - an effective disorder strength $W_{eff}\approx\Delta_{1}/\xi_{1}^{2}$ and an effective matrix element $t_{eff}\approx u/\xi_{1}$. The ratio of both yields a dimensionless new parameter $t_{eff}/W_{eff}=u\xi_{1}/\Delta_{1}$. In the limit of weak disorder $\Delta_{1}$ turns simply into a constant leaving us with the relevant parameter $u\xi_{1}$. We then plot $\xi_{2}/\xi_{1}\equiv F(u\xi_{1})$ versus $u\xi_{1}$ in Fig. 5(a). For $u\xi_{1}<50$ we find agreement with the data from von Oppen et al. (1996) which indicate $F(x)\sim x$. Lowering $W<1$ leads to an increase of $u\xi_{1}$ and a consequent crossover into the asymptotic scaling regime, which shows a significant slowing down of the increase of the TIP localization length. The corresponding scaling function $F(x)$ turns nonlinear with sublinear growth for large values of $x>50$. The extrapolation of the data from von Oppen et al. (1996) (dashed line) overestimates the TIP localization length by at least a factor of 6 for values $x=u\xi_{1}=1500$. Similarly the data from Ref. Frahm, 2016 overestimate the length by a factor of up to 2 for $x=1500$, probably due to finite size effects which we did take into account. The crossover into the asymptotic regime at $x\approx 50$ is highlighted in Fig. 5(b). Notably the record values of the TIP localization length are obtained for $u\approx 3$. The solid line in Fig. 5(a) connects these data points, and indicates sublinear growth of $\xi_{2}/\xi_{1}$ with $u\xi_{1}$ and a corresponding nonlinear scaling function $F(x)$.

In order to test the single parameter scaling hypothesis, Kramer and MacKinnon (1993) we compute the participation number as

[TABLE]

As before the number of disorder realizations $S$ for a given system size $N$ was fixed such that $NS\approx 10^{6}$. Additionally we perform a finite size scaling to eliminate finite size corrections. The values of the ratio $\text{PN}/\xi_{2}$ for $W=0.5,0.75,1.0,1.5,2.0$ are plotted in Fig. 6. We observe that $1\leq\text{PN}/\xi_{2}\leq 1.5$ which indicates that the extension of the TIP wavefunctions is of the same order as the localization length $\xi_{2}$ which controls their exponential decay, confirming the single parameter scaling.

V Discussion

In conclusion, we showed that two interacting particles in a disordered potential enter an asymptotic scaling regime of their localization length $\xi_{2}$ in units of the single particle localization length $\xi_{1}$ for weak disorder due to the restoring of momentum conservation in the single particle eigenfunctions. The ratio $\xi_{2}/\xi_{1}=F(u\xi_{1})$ grows to record values of $F=15$ for $\xi_{1}=400$ and $u=3$, albeit the growth is much slower than anticipated from earlier numerical studies and reflected in a nonlinear dependence of the scaling function $F(x)$ on its argument. Our findings are supported by the manifestation of resonant couplings between fragile groups of Fock states which have much smaller size ($\xi_{1}$) than originally anticipated ($\xi_{1}^{2}$). Further, the fragility is supported by the presence of particle-hole symmetry due to the bipartite nature of the tight binding chain for a single particle in the limit of weak disorder. This particle-hole symmetry guarantees that resonant groups of Fock states can conserve both momentum and energy. We expect that a violation of particle-hole symmetry e.g. by adding next-to-nearest-neighbour hoppings will reduce the size of the resonantly interacting Fock state groups and further reduce the enhancement factor of $\xi_{2}$ over $\xi_{1}$. The nature of the observed non-linear correction of the scaling function $F(x)$ is an interesting subject of future studies.

Acknowledgements.

We thank Boris Altshuler and Mikhail Fistul for stimulating discussions. This work was supported by the Institute for Basic Science in Korea (IBS-R024-D1). T.E. acknowledges financial support by the Alexander-von-Humboldt foundation through the Feodor-Lynen Research Fellowship Program No. NZL-1007394-FLF-P.

Appendix A Fitting and extrapolating

We explain in this appendix a step by step procedure for finding the localization length $\xi_{2}$ for a specific example with $W=1$, $N=8500$ and interaction strength $u=0.5$. We perform the two steps outlined in the main text in Sec. IV.1. We compute

[TABLE]

averaged over $200$ disorder realisations. By definition (see Eq. (6)), $\xi_{2}(N)$ is extracted by fitting the linear decay of $y(x)$ over the range of $x=|n-m|$, where $n=425(0.05N)$ and $m=8075(0.95N)$. This choice of the range was adapted from the Ref. Frahm, 2016: we discarded the parts close to the boundary to avoid its effect. Figure A.1 shows the average $y(x)$ and the linear fit $y(x)=c-x/\xi_{2}(N)$ for the chosen set of parameters. Similarly $\xi_{2}(N)$ are extracted from the inverse of the slope of the fits obtained for other $N\geq N_{*}(W=1)\approx 1500$.

Next we perform the extrapolation $N\to\infty$ of $\xi_{2}(N)$ to obtain $\xi_{2}(N\to\infty)$. We observe that $\xi_{2}(N)$ show only a mild dependence on $N$ for $N\geq N_{*}(W)$ and rather fluctuates around some average value. Therefore we try both a constant fit $f(x)=c$ and a linear fit $y=m*1/N+c$ from where $\xi_{2}(N\to\infty)=c$. The data and the fits are shown in Fig. A.2. The values of $\xi_{2}(N\to\infty)$ obtained from both fits are reasonably close, always within $5\%$. In the current setting example, $\xi_{2}=272.52$ (linear fit), $267.59$ (constant fit) thereby giving $\xi_{2}=267.59\pm 4.9$. We used the constant fit systematically.