Diffusion of excitation and power-law localization in long-range-coupled strongly disordered systems
DDiffusion of excitation and power-law localization in long-range-coupled stronglydisordered systems
K. Kawa ∗ and P. Machnikowski † Department of Theoretical Physics, Wrocław University of Science and Technology, 50-370 Wrocław
We investigate diffusion of excitation in one- and two-dimensional lattices with random on-siteenergies and deterministic long-range couplings (hopping) inversely proportional to the distance.Three regimes of diffusion are observed in strongly disordered systems: ballistic motion at shorttime, standard diffusion for intermediate times, and a stationary phase (saturation) at long times.We propose an analytical solution valid in the strong-coupling regime which explains the observeddynamics and relates the ballistic velocity, diffusion coefficient, and asymptotic diffusion range to thesystem size and disorder strength via simple formulas. We show also that in the long-time asymptoticlimit of diffusion from a single site the occupations form a heavy-tailed power-law distribution.
I. INTRODUCTION
The seminal paper of Anderson [1] anticipated lack ofdiffusion in particular random lattices in three dimen-sions (3D), starting the topic of Anderson localization of(quasi)particles. Although the original discussion con-cerned systems with power-law couplings, much of thesubsequent research dealt with tight-binding-like mod-els with nearest-neighbor coupling and on-site disorder.Within that framework, Mott and Twose [2] proved thelack of diffusion in one dimension (1D). Later, Abra-hams et al. [3] proposed single parameter scaling hy-pothesis and proved the absence of diffusion in two di-mensions (2D), confirming at the same time localizationin one dimension and localization-delocalization phasetransition in three dimensions with respect to the systemsize. Nonetheless, the single parameter scaling hypothe-sis remains an approximate result, demands finite rangehopping, uncorrelated disorder and time-reversal symme-try. It still leaves space for the possibility of diffusion inspecific models, even for dimension less or equal to two.Beside the tight-binding model with on-site disorderand nearest-neighbor inter-site coupling, other theoret-ical models were studied from the point of view of lo-calization and metal-insulator transition [4, 5]. Amongothers, much interest was devoted to models of uncorre-lated diagonal disorder with long range hopping of thepower-law character ∝ /r µ [4, 6]. Such a model rep-resents several physical systems. For instance, Levitov[7] analyzed delocalization of vibrational modes in a 3Dcrystal with dipole interaction ∝ /r . Subsequently,the energy transfer in several systems has been exploredextensively within similar models. In biological light-harvesting systems, energy transport to the reaction cen-ter is mediated by dipole interactions ∝ /r [8–10]. Fur-thermore, many-body systems of nuclear spins with thesame /r interaction have been studied in the context oflocalization [11]. A model with both the the short-range ∗ [email protected] † [email protected] hopping between adjacent sites and long-range dipole-dipole coupling was used for describing energy transferin self-assembled nano-rings [12]. The long-range cou-pling was shown to stabilize the system against disor-der. Moreover, energy transfer has been observed ex-perimentally in the planar quantum dot (QD) ensembles[13]. It has been found that the transfer is even bet-ter for low-density samples i.e. of greater inter-dot dis-tances preventing carriers from tunneling between QDs.Thus, a mechanism different from quantum tunnelingmust be proposed for explaining the energy transfer, anda plausible explanation seems to be the long range dipole-dipole coupling. In fact, an ensemble of QDs seems tobe a coupled system, as its radiation properties can-not be explained as a sum of single emitters [14, 15].The fundamental coupling between the emitters emergeshere from their interaction with a common electromag-netic reservoir, leading to the dispersion-force coupling ∝ cos( kr ) / ( kr ) + sin( kr ) / ( kr ) + cos( kr ) / ( kr ) [16–18],where k is the resonant wave number, which reduces tothe usual /r dipole interaction on short distances butis dominated by the /r term at distances larger than thewave length resonant with the optical transition.A renormalization group analysis by Rodríguez et al. [4], proves the existence of extended states in the modelof uncorrelated diagonal disorder and power-law hopping ∝ /r µ , with an exponent µ greater than the dimensionof the system. Extended states appear in such a systemin the thermodynamic limit in the vicinity of the energyband edge, even in 1D and 2D. Such a model was alsoused for investigating the wave packet propagation in 1D[19]. Time evolution of the wave packet is described hereby its mean square displacement and the participationratio. The localization-delocalization transition occursas a function of the disorder magnitude for < µ < / .For µ > / wave packet tends to localize because of theshort range of hopping decreasing while µ increases. Theparticular case of µ = 1 , relevant to the dispersion forcesat large distances, shows interesting properties already inthe absence of disorder, showing a combination of diffu-sive and super-diffusive transport [20]. However, in disor-dered systems this case has been studied only marginally.It is known that it shows certain criticality features, like a r X i v : . [ c ond - m a t . d i s - nn ] S e p FIG. 1. The system geometry for d = 2 : sites forming aregular lattice on a circular mesa. The initially excited site ismarked by colors. the divergence of the critical disorder strength needed tolocalize the upper band edge [21].In this paper we consider the dynamics of a single exci-tation in arrays of two-level systems with the long-range /r hopping and strong diagonal disorder. We show thata system of finite size in this limit shows three consec-utive phases of excitation transport: a ballistic one, fol-lowed by normal diffusion, and finally saturation of themean-square diffusion range. We analyze also the aver-age distribution of the excitation in the asymptotic (sat-uration) phase and demonstrate a heavy-tail power lawquasi-localization around the initially excited site. Wepoint out that in the strong disorder limit the excita-tion transfer is dominated by the direct coupling withthe initial site, which allows us to reduce the model toa simplified form, which is exactly solvable in a certainrange of parameters. In this way we are able to relatethe transport parameters (ballistic speed, diffusion co-efficient, and asymptotic diffusion range) to the systemsize and disorder strength.The organization of the paper is as follows. In Sec. II,we describe the physical system and theoretical model.In Sec. III we present the results of our numerical simu-lations. The approximate analytical solution for stronglydisordered systems is presented in Sec. IV A, analyzed inSec. IV B and discussed in Sec. IV C. Final conclusionsare presented in Sec. V. II. THE SYSTEM AND THE MODEL
In this section, we describe the physical system understudy and introduce a theoretical model used for numer-ical simulation.We study a d -dimensional ( d = 1 , ) system of N siteson a regular lattice. The system occupies a disk of radius R for d = 2 and a line segment of length R for d = 1 (Fig. 1 shows the geometry for d = 2 ). The number ofsites is related to the dimensionless radius of the system (in units of the lattice constant) by N = ζ d R d , where ζ d is a number that depends on the space dimensionality d .The system is described by the Hamiltonian H = J (cid:88) α (cid:15) α | α (cid:105)(cid:104) α | + (cid:88) αβ V αβ | α (cid:105)(cid:104) β | , (1)where | α (cid:105) represents a basis state localized at the site α , J sets the overall energy scale, J(cid:15) α is the correspondingon-site energy, and JV αβ is the coupling between the sitestates α and β . The dimensionless energies (cid:15) α (in unitsof J ) are uncorrelated normally distributed random vari-ables of zero expected value and standard deviation σ .The coupling V αβ has a long-range character, V αβ = | r α − r β | , for α (cid:54) = β, , for α = β, (2)where r α is the dimensionless position of the site α (inunits of the lattice constant).The central site is initially excited. The diffusion of theexcitation is described by its mean square displacement(MSD) from the origin of the system (cid:10) r ( t ) (cid:11) = (cid:42)(cid:88) α | c α ( t ) | r α (cid:43) , (3)where (cid:104) . . . (cid:105) denotes the average over disorder realiza-tions, c α ( t ) are the coefficients of expansion of the systemstate in the localized basis, | Ψ( t ) (cid:105) = (cid:88) α c α ( t ) | α (cid:105) , (4)and the central site corresponds to the r = . III. SIMULATION RESULTS
In this section we present the results of numerical sim-ulations for the model described in Sect. II in one and twodimensions. The system evolution is found by exact nu-merical diagonalization of the system Hamiltonian. Weinvestigate the character of the excitation diffusion as afunction of the system parameters, i.e., the number ofsites and the magnitude of the disorder characterized byits standard deviation σ .The MSD of the excitation from the central site isshown in Fig. 2 for different system sizes and disorderstrengths. Three consecutive regimes of the excitationtransport can be seen. First, the transport is ballisticwith a certain velocity v , (cid:104) r ( t ) (cid:105) = v t , for < t < t . (5)The velocity depends on the system size, as can be seenin Fig. 2(a,b), but not on the disorder [Fig. 2(c,d)]. At a − − − (a) 1D v t Dt r σ = 1000 N=51N=4001 − − − − − (c) 1D v t Dt r N = 8001 σ = σ = (b) 2D v t Dt r σ = 1000 N=111N=8171 − − (d) 2D v t Dt r N = 4051 σ = σ = h r i − − − h r i time10 − − − − − time − − FIG. 2. The MSD of the excitation from the central siteas a function of time: (a,b) for different sizes of the systembut the same disorder σ = 1000 ; (c,d) for different disordermagnitudes but the same system sizes, as shown. Panels (a,c)and (b,d) show the results for d = 1 and d = 2 , respectively.Dotted lines show the results from the “central atom model”(Sec. IV). All the results are averaged over 1000 repetitions. certain cross-over time t , the normal diffusion with thediffusion coefficient D sets on, (cid:104) r ( t ) (cid:105) = Dt, for t < t < t . (6)As demonstrated by our simulation results in Fig. 2, thediffusion coefficient D grows with the system size and de-creases as the disorder amplitude grows. From continuityrequirement at the crossover one gets t = D/v . (7)Simulations show that this value is size-independent butincreases with the disorder strength. Finally saturationis reached at the second cross-over time t . (cid:104) r ( t ) (cid:105) = r , for t > t . (8)Obviously, the diffusion range must be limited in a finitesystem. However, the values of r presented in Fig. 2,although increasing with the number of nodes, are alwaysconsiderably lower than the system size and decrease asthe disorder grows. Again, the cross-over time is fixed bycontinuity, t = r /D. (9)In order to study the dependence of the dynamicalcharacteristics v , D and r on the system size and dis-order strength, we found the system evolution for a range . . . .
100 1000 t = √ π/σσ = 1000 (c) 1D 100 1000 σ = 1000 (d) 2D . v ∝ √ N D ∝ N σ = 1000 (a) 1D ⇐ ⇒ . σ = 1000 (b) 2D ⇐ ⇒ t N N v D FIG. 3. Size dependence of the dynamical coefficients. (a,b)The ballistic velocity (green circles, right axis) and the dif-fusion coefficient (red squares, left axis) as a function of thenumber of sites. (c,d) The ballistic-to-diffusive crossover timeas a function of the number of sites. Lines show the analyticalresults from the „central atom model” (Sec. IV B). of values of N and σ , and fitted the numerical solutions inthe respective time intervals with the appropriate power-law functions of time according to Eq. (5), Eq. (6), andEq. (8). The dependence of all the dynamical character-istics on the two system parameters turns out to be apower law over at least a decade of parameter variationin each case, with the power-law exponent very close toan integer or a simple fraction.The size dependence of the velocity and diffusion coef-ficient is depicted in Fig. 3(a,b). The velocity increaseswith the number of atoms like N / , while the diffusioncoefficient grows linearly with the number of sites bothin 1D and 2D. As a consequence [Eq. (7)], the first cross-over time is size independent, see Fig. 3(c,d).The dependence of the dynamical parameters on thedisorder strength is shown in Fig. 4(a,b). While thevelocity remains σ -independent, as already concludedabove, the diffusion coefficient is inversely proportionalto σ for any system dimension.Finally, in Fig. 5 and Fig. 6 we analyze the depen-dence of the saturation level r on the system size anddisorder strength, respectively. In 1D, the value of r grows as N for large disorder (Fig. 5(a)), while for mod-erate disorder one observes an approximately power-lawdependence with a lower exponent (Fig. 5(c)). In 2Dthe dependence is r ∝ N / for strongly disorderedsystems (Fig. 5(b)) and for sufficiently short chains atweaker disorder (Fig. 5(d)). This kind of power-law de-pendence in two dimensions means that the saturationlevel grows faster than the system size, hence the excita- − − − − − − − − − t = √ π / σ N = 8001 (c) 1D − − − N = 4051 (d) 2D − v = const . D ∝ / σ N = 8001 (a) 1D ⇐ ⇒ − N = 4051 (b) 2D ⇐ ⇒ t /σ /σ v D FIG. 4. Dependence of the dynamical coefficients on the dis-order strength. (a,b) The ballistic velocity (red circles, leftaxis) and the diffusion coefficient (blue squares, right axis) asa function of the standard deviation of the on-site energies.(c,d) The ballistic-to-diffusive crossover time as a function ofthe standard deviation of the on-site energies. Lines show theanalytical results from the „central atom model”. tion must reach the border of the system for sufficientlylarge number of sites (on the order of σ ). From thatpoint saturation level starts to grow linearly with thesystem size as it is clear from Fig. 5(d). In other words,the saturation level is bounded by the system size. Bothin 1D and 2D, the value of r turns out to be pro-portional to /σ as long as σ (cid:29) , while it starts tooscillate and reaches a constant value as σ → . The lat-ter property simply reflects the uniform spreading of theexcitation across the system, characteristic of an unper-turbed lattice, which roughly sets the upper limit on themean-square displacement. The values corresponding tothe uniform distribution are R / and R / for d = 1 and d = 2 , respectively, and are marked with horizontaldotted lines in Fig. 6. The crossover time t is then pro-portional to N and N / in 1D and 2D respectively, thatis, to the linear size of the system in both cases, and isindependent of disorder, as follows from Eq. (9).By a simple lowest-order “resonance counting” argu-ment, the number of sites resonant to the central (ini-tially occupied) one at a distance x is proportional to /x . For a 1D system, as the chain gets longer, due tooccupation spreading among these resonant sites, the oc-cupation of the distant sites would then grow as ln N , thevalue of r would grow as N , and the long-time occu-pations would be distributed as /x . This prediction for r agrees with the simulation results for very large dis-order but discrepancy is visible at σ = 10 (Fig. 5(c)). Thegrowing occupation of the distant sites should suppress −
100 1000 ∝ N σ = 1000 (a) 1D −
100 1000 ∝ N / σ = 1000 (b) 2D −
10 100 1000 10000 ∝ N σ = 10 (c) 1D
100 1000 10000 ∝ N / ∝ Nσ = 75 (d) 2D r s a t r s a t N N
FIG. 5. The dependence of the diffusion range on the num-ber of sites in one (a,c) and two (b,d) dimensions. Dashedlines show the analytic results from the „central atom model”.Dotted line in (d) is the linear fit to to the further part of theplot. − − N=8001(a) 1D ∝ /σ − − N=4051(b) 2D ∝ /σ r s a t /σ /σ FIG. 6. The dependence of the diffusion range on the stan-dard deviation of the on-site energies in one (a) and two (b)dimensions. Solid lines show the analytical results from the„central atom model”. Dotted lines represent the values cor-responding to uniform distribution over the system. the final occupation of the initial site | c | (the survivalprobability) as − A ln N , until the survival probabilityis reduced enough for the lowest order approach to breakdown. This is indeed confirmed by simulation results pre-sented in Fig. 7(a), where the logarithmic dependence isrepresented by the dashed trend line, although the nar-row range of variability of | c | may be insufficient torigorously prove the subtle logarithmic dependence. Thesaturation at N (cid:38) appears at a rather high value ofthe survival probability. The distribution of occupationsat saturation (asymptotic long-time limit) indeed showsa distinct power-law character over many orders of mag-nitude of the chain length, with an exponent very close to1 for strong disorder (red squares in Fig. 7(b)). This be- . . . . . . (a) σ = 10 10 − − − − −
10 100 1000 (b) N = 8001 | c | N | c j | jσ = 100 σ = 0 . FIG. 7. (a) The final (long-time asymptotic) occupation ofthe initial site in a one-dimensional chain as a function ofthe chain length for σ = 10 . (b) The dependence of theasymptotic occupation on the site index in a one-dimensionalchain for the values of N and σ as shown. The dashed lineshows a logarithmic trend. − − − − − (a) σ = 75 (cid:26) N = 1009 N = 15813 σ = 1000 (cid:26) N = 1009 N = 15813 | c r | r r N = 4053 σ = 1000 N = 4053 σ = 0 . FIG. 8. The dependence of the asymptotic occupation on thedistance from the initially occupied site in a 2-dimensionalsystem: (a) for a fixed σ and two values of N as shown; (b)for a fixed N and two values of σ as shown. havior is characteristic of the strong disorder regime, cor-responding to the power-law dependence of r on σ (leftasymptotics of Fig. 6(a)). It should be contrasted withthe low-disorder limit (right asymptotics in Fig. 6(a)),where the occupations are equally spread all over thechain, as we have already inferred from the asymptoticvalue (blue circles in Fig. 7(b)).A similar power-law localization of the asymptoticstate is observed in a 2D system, as shown in Fig. 8,where we plot the average occupation of a site as a func-tion of its distance from the central site. As can be seen inFig. 8(a) (full symbols), for a moderate disorder strength,the exponent of this power-law dependence increases bymagnitude as the system size grows. The values fromfitting to the power-law part of the data obtained fromsimulations range from − . for N = 197 to − . for N = 32017 (fitting was performed using the central partof the data points, showing the clear power-law depen-dence and may be slightly affected by the choice of thecut-offs). All these values are above − and thereforecorrespond to heavy tail distributions in two dimensionsthat would have a divergent norm if extrapolated to infi-nite size. The values seem to be roughly proportional to ln N over the range of system sizes available for numericalsimulations and we were not able to reliably determinethe asymptotic value of the exponent as N → ∞ .A different situation is observed for strongly disor-dered systems (empty symbols in Fig. 8(a)). Here theslope of the power-law dependence is apparently constantand indeed, the exponents obtained from fitting oscillate(due to inherent randomness and fitting uncertainty) veryclose to the value of − , which precisely corresponds tothe r ∝ N / dependence in Fig. 5(b). This meansthat in this range of system sizes, the asymptotic diffu-sion range grows with the system size by extending the ∝ /r dependence of occupation to larger and larger dis-tances at the cost of the central site occupation, until the ∝ N / dependence breaks down, similar to the situationin Fig. 5(d) but far beyond the range of system sizes ac-cessible in simulations. When comparing the power-lawdependence at the two disorder strengths one arrives atthe somewhat unexpected conclusion that, from the for-mal point of view focused merely on the power-law ex-ponent, localization in moderately disordered systems isstronger (the absolute value of the exponent is larger)than in heavily disordered ones. This is obviously nottrue in terms of the actual values of the occupations thatdecrease as the disorder grows, which simply means thatthe survival probability at the initial site grows with dis-order, as expected. Interestingly, the occupations of themost remote sites are similar for the two different disor-der regimes shown in Fig. 8(a).The trend in the dependence of the average asymp-totic distribution of occupations on the disorder strengthdemonstrated above cannot remain valid towards weakerdisorder strengths, as the occupation should spreadacross all the system in the limit of unperturbed chain.Indeed, as shown in Fig. 8(b), for a very weak disorder,the system tends to a uniform distribution but, rathersurprisingly, even for σ = 0 . a weak localization effectis still visible.In addition to the simulations of the full model dis-cussed so far, we have also studied the dynamics of a“central atom model” in which the central (initially ex-cited) site is coupled to all the other sites in the systemas in the full model, but the other sites are not coupledwith one another, i.e., V αβ = 0 if α (cid:54) = 0 and β (cid:54) = 0 . Theresults are shown in Fig. 2 with dotted lines. One cansee that the system evolution is nearly the same in bothmodels in this parameter range, which means that thedynamics in the strongly disordered case is dominatedby direct jumps to remote places, which is possible dueto the long-range coupling.In the next section we show that the dynamical pa-rameters can be related to the system size and disorderstrength and the analytical relation between the power-law exponents and the system dimension can be found aslong as the “central atom model” is valid. IV. APPROXIMATE ANALYTICAL SOLUTION
In this section we present an approximate analytic ap-proach to the considered problem, which becomes pos-sible within the simplified “central atom model” in thelimit of strong disorder.
A. Solution of the Equation of Motion
The equation of motion for the coefficients of the ex-pansion defined in Eq. (4) has the form i ˙ c α ( t ) = (cid:15) α c α ( t ) + (cid:88) β V αβ c β ( t ) , c α (0) = δ α . (10)We define the Laplace transform f α ( s ) of an amplitude c α ( t ) , f α ( s ) = ∞ (cid:90) e − st c α ( t )d t, (11)where s is a complex variable. Equation (10) in terms ofthe Laplace transform is f α ( s ) = iδ α is − (cid:15) α + (cid:88) β (cid:54) = α V αβ f β ( s ) . (12)Eq. (12) can be iteratively expanded in series dependingonly on the central-site term f , f α ( s ) = V α is − (cid:15) α f ( s ) (13) + (cid:88) β (cid:54) = α,β (cid:54) =0 V αβ is − (cid:15) β V β is − (cid:15) β f ( s ) + . . . , α (cid:54) = 0 and f ( s ) = 1 is − (cid:15) + (cid:88) β (cid:54) =0 V β V β ( is − (cid:15) )( is − (cid:15) β ) (14) + (cid:88) γ (cid:54) = β,γ (cid:54) =0 V β V βγ V γ ( is − (cid:15) )( is − (cid:15) β )( is − (cid:15) γ ) + . . . f ( s ) . The series have the number of terms of the order of N N and the resulting system of equations cannot be solvedin a simple manner. The problem simplifies considerablyin the “central atom approximation”. Then, only the firstterm in Eq. (13) remains, while in Eq. (14) only the firstterm and the first sum survive. After these simplificationone can write f α ( s ) as f α ( s ) = iV β (cid:81) β (cid:54) =0 ,β (cid:54) = α ( is − (cid:15) β ) (cid:81) β ( is − (cid:15) β ) − (cid:80) β (cid:54) =0 (cid:81) γ (cid:54) = β,γ (cid:54) =0 V β V β ( is − (cid:15) γ ) . (15) Next, we transform back to the time domain, by meansof the Mellin’s formula c α ( t ) = 12 πi lim T →∞ γ + iT (cid:90) γ − iT e st f α ( s )d s. (16)To perform the integral we employ residue theorem andobtain c α ( t ) = (cid:88) n b ( α ) n e − iz n t , (17)where z n = is n and s n are the poles of the analytic func-tion given by Eq. (15). All z n must be real. In addition, b ( α ) n = 2 πi Res z n f α ( z ) = V α (cid:81) β (cid:54) =0 ,β (cid:54) = α ( z n − (cid:15) β ) (cid:81) β (cid:54) = n ( z n − z β ) , (18)where Res z n f α ( z ) denotes the residue of f α ( z ) at z n .As long as the coupling is small in comparison to σ/N ,the roots of the denominator of Eq. (15) lie in close vicin-ity to bare energies (cid:15) α , as compared to the typical dis-tance between these roots. Hence, one can associate eachpole with the nearest bare energy, align the numberingand assume z n − (cid:15) β z n − z β ≈ , whenever n (cid:54) = β . With thisapproximation one finds b ( α ) n ≈ V α z − z α , n = 0 , n (cid:54) = α,V α z n − z , n = α (cid:54) = 0 , , otherwise.The amplitudes c α ( t ) then take the form c α ( t ) = b ( α )0 e − iz t + b ( α ) α e − iz α t = e − iz t V α z − z α (cid:16) − e − i ( z − z α ) t (cid:17) and the occupation of the site α becomes | c α ( t ) | = | V α | sin ( δz α t/ δz α / , (19)where δz α = z − z α . Upon inserting this result to theEq. (3) we obtain (cid:10) r ( t ) (cid:11) = (cid:42)(cid:88) r r | V r | n r (cid:88) k =1 sin ( δz k ( r ) t/ δz k ( r ) / (cid:43) , (20)where we decomposed the sum over all the sites into sub-sets of n r sites lying at a distance r from the origin. Allthe sites at a given distance have the same coupling tothe central dot V r = 1 /r . δz k ( r ) are the values of δz k for sites lying at distance r from the origin. δz ( r ) canbe considered random variables with a certain probabil-ity density f r ( δz ) . In the continuum approximation onethen obtains (cid:104) r ( t ) (cid:105) = R (cid:90) ζ d dr d − d r ∞ (cid:90) −∞ f r ( u ) sin ( ut/ u/ d u, (21)where ζ d r d is the number of sites lying inside a d -dimensional sphere.At large δz , the poles are shifted negligibly from thebare energies, hence the distribution function is close tothe on-site energy difference distribution f ∞ ( δ(cid:15) ) (here ∞ refers to infinite distance, hence vanishing coupling, and δ(cid:15) is the bare energy difference), which is a Gaussian dis-tribution with the standard deviation √ σ . However, at δz ∼ V r the pole probability distribution must reflect thelevel repulsion. Its form can be found by noting that apair of sites with a bare energy difference δ(cid:15) > coupledwith a coupling strength V gives rise to a pair of polesseparated by δz = (cid:112) ( δ(cid:15) ) + 4 V . From this, the cumu-lative distribution function for δz follows in the form P (0 < δz k ( r ) < u ) = √ u − V r (cid:90) f ∞ ( x )d x, | u | > V r , (22)and the corresponding probability density is f r ( u ) = (cid:40) f ∞ ( (cid:112) u − V r ) | u | √ u − V r , | u | > V r , , | u | ≤ V r . (23)The above distribution corresponds to the normal dis-tribution of standard deviation σ u = √ σ having a gaparound zero of width V r . The total distribution for asystem of radius R including all possible distances r hasa gap of width V R . B. Regimes of propagation
Eqs. (21) and (23) allow us to explain the three timephases in the excitation evolution.For very short times, t (cid:46) /σ , the function h ( u ) =sin ( ut/ / ( u/ is slowly varying in u and can be ap-proximated by h ( u ) ≈ t over the whole width of thedistribution f r ( u ) , as schematically shown in Fig. 9(a).Eq. (21) then immediately yields (cid:104) r ( t ) (cid:105) = v t , that is,ballistic propagation with the constant velocity v = R (cid:90) ζ d dr d − d r (cid:90) f ∞ ( u )d u = ζ d R d = N. (24)This dependence is shown as dashed lines in Fig. 3(a,b),which perfectly follows the numerical data. FIG. 9. Schematic plots of the probability density f r (red,the same for three cases) and the function h ( u ) (blue) cor-responding to three phases of propagation: (a) ballistic, (b)diffusive, (c) saturation. In the intermediate time range, /σ < t < /V R ,the function h ( u ) probes the central part of the distri-bution but is still insensitive to the narrow central gap(Fig. 9(b)). Then the integral over u in Eq. (21) can beapproximated by ∞ (cid:90) −∞ f r ( u ) sin ( ut/ u/ d u ≈ ∞ (cid:90) −∞ f ∞ ( u )2 πδ ( u )d u = 2 πtf ∞ (0) , from which one finds the diffusive transport (cid:104) r ( t ) (cid:105) = Dt ,with D = R (cid:90) d rζ d dr d − πtf ∞ (0) = √ πNσ . (25)Again, this dependence on N and σ is shown as dashedlines in Fig. 3(a,b) and in Fig. 4(a,b). The agreementwith the numerical data is excellent.Using Eq. (7), the crossover time between ballistic anddiffusive phases is obtained as t = √ π/σ, (26)which agrees with our simulation results, as shown inFig. 4(b), for sufficiently large values of σ . The rangeof the ballistic transport is therefore (cid:104) r ( t ) (cid:105) = v t = N √ π/σ . In 1D it is always much smaller than the systemsize if the disorder is strong.Finally, at t ≈ /V R , the function h ( u ) becomes asnarrow as the the gap in the density function, henceits central part does not contribute, while its oscillat-ing tails are averaged to ˜ h ( u ) = (1 / /u (see Fig. 9).There is no time dependence in this limit, which re-sults in the saturation of (cid:104) r (cid:105) . The saturation level isthen estimated from Eq. (21) as (cid:82) ∞ V r f r ( u )˜ h ( u )d u = √ πr σ exp (cid:2) / ( σ r ) (cid:3) erfc[1 / ( σr )] The resulting saturationsaturation value is r = (cid:104) r ( t ) (cid:105) = R (cid:90) ζ d dr d √ π σ exp (cid:18) σ r (cid:19) erfc (cid:18) σr (cid:19) d r ≈ √ πζ d σ dd + 1 R d +1 (27)where we took into account that V R /σ = 1 /Rσ (cid:28) inthe high disorder energy regime, so the last two terms inthe integral tends to unity (in the zeroth order of Taylorexpansion).The dependence from Eq. (27) is marked by dashedlines in Fig. 5 and Fig. 6 and agrees very well with thesimulation result as long as the size is not too large andthe disorder is strong.The onset of saturation is determined by the continuityof (cid:104) r (cid:105) , Dt = r . Combining Eq. (27) with (25) oneobtains t = 1 πV R , (28)in agreement with the simulations.Within the “central atom model”, the asymptotic oc-cupation of a site at a distance r from the origin is (cid:104)| c r | (cid:105) sat = 2 r ∞ (cid:90) V r f r ( u ) 2 u d u = √ π σr exp (cid:18) σ r (cid:19) erfc (cid:18) σr (cid:19) ≈ √ π σr . (29)This ∝ /r trend obtained from our approximate ana-lytical solution agrees with the numerical data in 1D,shown in Fig. 7 and corresponds to the survival prob-ability − A ln N in the strong disorder regime, whichis consistent with the simulation data in Fig. 7(a). In2D, as we have seen in Fig. 8, the same dependence isobtained for very strongly disordered systems. C. Discussion
As we have seen, the analytical formulas agree very wellwith the simulation results only within a certain limitsof system size and disorder strength. One discrepancyappears when the disorder becomes weak (see Fig. 6).This is obvious, as our “central atom model” is essentiallybased on the assumption that coupling is a perturbationto the on-site energies, which requires a strong disorder.The second discrepancy appears for long chains in Fig. 5.For d ≥ , r ∼ R d +1 , Eq. (27) predicts that the asymp-totic diffusion range grows faster than the system size.This cannot be true for an arbitrary system size and,indeed, the trend in simulations changes at a certain sys-tem size Fig. 5(d), which is not captured by the “centralatom” approximation. We note that the limit of valid-ity of our approximation is σ/N ∼ V R = 1 /R or, using N ∼ R d , R d − ∼ σ . At this limit, (cid:104) r (cid:105) ∼ R , which as-sures consistency of our conclusions. On the other hand,for a one-dimensional chain, r ∼ R , hence the asymp-totic range grows linearly with the system size. Eq. (27)is valid for σ (cid:29) , hence r (cid:28) R and the excitation iseffectively trapped around the original site in the chain. V. CONCLUSIONS
We have studied the diffusion of an initially localizedexcitation in finite lattices with a strong on-site disorderand a long-range coupling (hopping) inversely propor-tional to the distance. We have shown that the diffusionin such a system takes place in three dynamical stages:ballistic transport is followed by normal diffusion andthen by saturation. The numerical findings can be un-derstood with the help of an approximate model whichis valid in the strong disorder regime and allows an an-alytical solution, which relates the dynamical propertiesto the system parameters (size and disorder strength).We have proposed two complementary descriptions oflocalization. The first one emerges from the dynamicsand consists in analyzing the asymptotic range of dif-fusion. The other one consists in studying the spatialdistribution of the average occupations of lattice sites inthe long-time limit and is more directly related to thestructure of energy eigenstates of the disordered system.We have shown, both numerically and analytically, thatthe diffusion range grows proportionally to the systemsize and is always much smaller than the latter in 1Dhence, from this point of view, the excitation remains ef-fectively localized around the initial site. In contrast, in2D systems, the range of diffusion initially grows fasterthan the system size as the latter increases, until at acertain system size the growth slows down so that, insufficiently strongly disordered systems the range againremains much smaller than the system size. While trac-ing the asymptotic diffusion range provides only a sin-gle number characterizing the degree of localization orspreading of the excitation, inspection of the average pro-file of occupations offers a more complete spatial pictureof the localization. We have found out that after a suffi-ciently long time the occupations stabilize into a heavy-tailed power-law distribution both in 1D and 2D that notonly does not have a second moment but would even havea divergent norm when extrapolated to infinite systemsize. Therefore, even if the excitation remains localizedin the sense of showing the diffusion range much smallerthan the system size, it does not have any intrinsic local-ization length and the diffusion reaches (on the average)an arbitrary site of the lattice according to the power-lawdistribution as a function of the distance.
ACKNOWLEDGMENTS
The authors are grateful to Marcin Mierzejewski fordiscussions.Calculations have been partially carried out using re-sources provided by Wroclaw Centre for Networking andSupercomputing (http://wcss.pl), grant No. 203. [1] P. W. Anderson, Phys. Rev. , 1492 (1958).[2] N. Mott and W. Twose, Adv. Phys. , 107 (1961).[3] E. Abrahams, P. W. Anderson, D. C. Licciardello, andT. V. Ramakrishnan, Phys. Rev. Lett. , 673 (1979).[4] A. Rodríguez, V. A. Malyshev, G. Sierra, M. A. Martín-Delgado, J. Rodríguez-Laguna, and F. Domínguez-Adame, Phys. Rev. Lett. , 027404 (2003).[5] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort,M. Zaccanti, G. Modugno, M. Modugno, and M. Ingus-cio, Nature , 895 (2008).[6] F. Domínguez-Adame and V. Malyshev, Am. J. Phys. (2004).[7] L. S. Levitov, Europhys. Lett. , 83 (1989).[8] G. L. Celardo, F. Borgonovi, M. Merkli, V. I.Tsifrinovich, and G. P. Berman, J. Phys. Chem. C ,22105 (2012).[9] M. Sarovar, A. Ishizaki, G. R. Fleming, and K. B. Wha-ley, Nature Phys. , 462 (2010).[10] M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, J. Chem. Phys. , 174106 (2008).[11] G. A. Álvarez, D. Suter, and R. Kaiser, Science , 846 (2015).[12] A. D. Somoza, K. W. Sun, R. A. Molina, and Y. Zhao,Phys. Chem. Chem. Phys. , 25996 (2017).[13] F. V. De Sales, S. W. Da Silva, J. M. Cruz, A. F. Monte,M. A. Soler, P. C. Morais, M. J. Da Silva, and A. A.Quivy, Phys. Rev. B , 1 (2004).[14] M. Scheibner, T. Schmidt, L. Worschech, A. Forchel,G. Bacher, T. Passow, and D. Hommel, Nature Phys. , 106 (2007).[15] M. Kozub, L. Pawicki, and P. Machnikowski, Phys. Rev.B , 121305 (2012).[16] M. J. Stephen, J. Chem. Phys. , 669 (1964).[17] R. H. Lehmberg, Phys. Rev. A , 883 (1970).[18] A. A. Varfolomeev, Zh. Eksp. Teor. Fiz. , 1702 (1970).[19] P. E. De Brito, E. S. Rodrigues, and H. N. Nazareno,Phys. Rev. B , 2 (2004).[20] B. Kloss and Y. Bar Lev, Phys. Rev. A , 1 (2019).[21] F. A. B. F. de Moura, A. V. Malyshev, M. L. Lyra, V. A.Malyshev, and F. Domínguez-Adame, Phys. Rev. B71