Relaxation dynamics of the three-dimensional Coulomb Glass model
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov Relaxation dynamics of the three-dimensional Coulomb Glass model
Preeti Bhandari,
1, 2
Vikas Malik, ∗ Deepak Kumar, † and Moshe Schechter Department of Physical Sciences, Indian Institute of Science Education and Research (IISER) Mohali,Sector 81, S.A.S. Nagar, Manauli P. O. 140306, India Department of Physics, Ben Gurion University of the Negev, Beer Sheva 84105, Israel Department of Physics and Material Science, Jaypee Institute of Information Technology, Uttar Pradesh 201309, India. School of Physical Sciences, Jawaharlal Nehru University, New Delhi – 110067, India. (Dated: November 10, 2020)In this paper, we analyze the dynamics of the Coulomb Glass lattice model in three dimensionsnear a local equilibrium state by using mean-field approximations. We specifically focus on under-standing the role of localization length ( ξ ) and the temperature ( T ) in the regime where the systemis not far from equilibrium. We use the eigenvalue distribution of the dynamical matrix to charac-terize relaxation laws as a function of localization length at low temperatures. The variation of theminimum eigenvalue of the dynamical matrix with temperature and localization length is discussednumerically and analytically. Our results demonstrate the dominant role played by the localizationlength on the relaxation laws. For very small localization lengths we find a crossover from exponen-tial relaxation at long times to a logarithmic decay at intermediate times. No logarithmic decay atthe intermediate times is observed for large localization lengths. PACS numbers: 71.23.Cq, 73.50.-h, 72.20.Ee
I. INTRODUCTION
The term Coulomb Glass (CG) refers to that categoryof disordered insulators that have a sufficiently high dis-order, which leads to localized electronic states coupledwith the Coulomb interactions. The presence of a glassyphase in this model has been predicted theoretically byseveral authors [1–5]. In dimensionless units, the Hamil-tonian for CG lattice model is defined [6] as H{ n i } = N X i =1 ǫ i n i + 12 X i = j e κ | ~r i − ~r j | ( n i − / n j − / ǫ i ’s are the on-site random field energy and theoccupation number n i ∈ { , } . The electrons at site i and j interact via unscreened Coulomb interaction e / ( κ r ij ) where κ is the dielectric constant.Much work has been done to find the ground state ofthe CG model at high disorder. Using mean field ap-proach [2], Monte Carlo simulation [7–9] and other opti-mization [10, 11] approaches it has been found that thereexist many metastable states (pseudo ground states) atlow temperatures. This metastability is responsible forglassy behavior. The density of states (DOS) found in allthese approaches shows a soft gap g ( E ) ∼ E δ around theFermi level [6, 12–18]. The value of δ is very near to thetheoretical prediction of d − d is the dimensionality ofthe system) given by Efros and Shklovskii [6]. RecentlyM¨uller and Ioffe have established a connection betweenthe presence of a glassy phase and the appearance of a ∗ [email protected] † deceased soft gap in three dimensional CG model using locatorapproximation [19]. The formation of a gap in the DOSaffects the conductivity ( σ ) quiet significantly. One cansee that the conductivity changes from the Mott’s law[20, 21] of ln σ ∼ ( T M /T ) / to the Efros-Shklovskii’slaw ln σ ∼ ( T ES /T ) / law [6] at low temperatures.The existence of glass transition in three dimensionalCG has been controversial and is a matter of activeresearch [22–25]. Although, some mean-field analysis,supported by recent numerical analysis [25], do suggestthe presence of a stable glassy phase [19, 26–29]. Non-equilibrium dynamics of structural and spin glasses hasbeen studied using scaling properties of non-stationarycorrelation and response functions [30–35]. Various nu-merical simulations claim that the CG model exhibitsglassy behavior i.e. slow relaxation [36–41], aging [39, 42]and memory effects [43, 44]. Many experimental tech-niques are used to study relaxation in the CG model[42–54] . The basic idea is to introduce a perturbationin the material to push the system out of equilibrium.This leads to an increase in the conductance, whose de-cay with time is then measured. It has been observedthat many materials, amorphous as well as crystalline,show a logarithmic temporal decay in conductance.The study of slow relaxation can be categorizedbroadly into two types of models: ( a ) A quasi-particlemodel, which was proposed by Pollak and Ovadyahu[13, 55]. They considered multi-particle transitions andshowed that the decrease in energy with time is relatedto γ m which is the minimal value of the transition rates γ = τ − exp [ − rξ − EkT ] where r and E are the collectivehopping distance and energy respectively and ξ is thelocalization length. Assuming that the change in con-ductance (∆ G ( t )) and energy are related to each otherlinearly one gets a logarithmic decay in conductance∆ G ( t ) ∝ − ln ( γ m t ) (2)( b ) Second is the local mean-field model, suggested byAmir et al [39, 40, 56]. The dynamics of quite a fewsystems near local stable minima can be described bythe matrix equation dδndt = − Aδn (3)where δn i = n i − f i is the fluctuation of the occupationnumber ( n i ) from its value f i at the local minima. Amir et al have shown [56] that under mean-field approxima-tions and single-particle transitions dynamics, the CGmodel obeys Eq.(3). The regime of low temperatures andsmall localization lengths is considered, and the distribu-tion P ( λ ) ∼ λ is found for the small relaxation rates.This leads to a logarithmic decay of fluctuations in oc-cupation numbers δn ( t ). Assuming that the relaxationof excess conductance ∆ G ( t ) is linear in δn ( t ), one re-covers the logarithmic decay for conductance as givenin Eq.(2). In this approach, the system always remainsnear the local minima, and thus the transition betweendifferent metastable states (multi-particle transitions) iscompletely neglected.Our goal here is to study the relaxation effects in theCoulomb Glass lattice model near a local equilibriumstate by using mean-field approximations. We follow theapproach of Amir et al [56], albeit for a lattice CG model.Within the approach of Amir et al there is disorder insite energies as well as in the position of the sites. Intheir approach [56], and small localization lengths stud-ied, the slow dynamics are mainly due to isolated local-ized states that have a long life-time. However, in thelattice model discussed here, disorder comes only via siteenergies and so the question of isolated states does notcome into the picture. Instead, we find that for smalllocalization lengths, ξ ≪
1, the states near the Fermilevel are very stable and any fluctuations in them re-lax very slowly. The main reason for this slow decayis that the states near the Fermi level are isolated en-ergetically due to the hard gap in the DOS inflicted ontheir near neighbor sites. For all localization lengths andtemperatures, the system always obeys the exponentialrelaxation ( δn ( t ) ∼ exp − ( t/τ max )), at times longer than τ max . The maximum relaxation time ( τ max ) is inverselyproportional to the smallest eigenvalue ( λ min ) of the dy-namical matrix A . Our study shows that λ min dependsupon the localization length as well as temperature.We further find that logarithmic time dependence ofthe relaxation of δn ( t ) at intermediate times is presentonly for small localization lengths, ξ ≪
1, where relax-ation is mainly due to jumps to nearest neighbor sites.The paper is organized as follows. In Sec. II, we haveprovided an overview of our derivation of the linear dy-namical matrix. In Sec. III, we present a detailed discus-sion of our mean-field results obtained numerically andanalytically. And finally in Sec. IV, we provide the con-clusions of our work.
II. DYNAMICS
The most general non-conserved dynamics for the to-tal probability distribution of the spins was developed byGlauber [57]. This was extended to conserved dynamicsby Kawasaki who incorporated the constraint of fixedmagnetization. The Kawasaki formulation [58, 59] ap-plies to CG as the electron number is conserved - whichis equivalent to fixed magnetization. Here we deal withthe probability distribution of P ( n . . . n N ; t ), which in-volves the occupation of all sites in the system. TheKawasaki dynamics holds for the interacting system aswell as for multi-particle dynamics. Since this approachis general, it can be taken beyond mean-field theory.The time evolution of a system can be described usinga generalized master equation [60] ddt P ( { n i } , t ) = − X i = j W i → j P ( { n i } , t )+ X j = i W j → i P ( { n j } , t ) (4)where W i → j denotes the transition rates from state i to j and P ( { n i } , t ) is the probability of finding the systemin state i at time t . The transition rates can be sin-gle or multi-electron transfer. Since we are interestedin Kawasaki dynamics, only transitions that conservethe particle (electron) number will be considered. Usingsingle-particle transitions, the Kawasaki dynamics equa-tion can be rewritten as ddt P ( n . . . n ν ; t ) = − X i = j ω i → j n i (1 − n j ) × P ( . . . , n i . . . n j , . . . ; t ) + X i = j ω j → i ¯ n j (1 − ¯ n i ) × P ( ..., ¯ n i ..., ¯ n j , ... ; t ) (5)where ω i → j is the transition probability from site i to j and ¯ n i = 1 − n i . Now we impose the condition of ”de-tailed balance”, so that the evolution is towards thermalequilibrium. In thermal equilibrium, ω i → j n i (1 − n j ) P eq ( . . . , n i . . . n j , . . . ) = ω j → i ¯ n j (1 − ¯ n i ) P eq ( . . . , ¯ n i . . . , ¯ n j , . . . ) . (6) ω i → j ω j → i = exp [ − βE ( . . . , ¯ n i . . . , ¯ n j , . . . )] exp [ − βE ( . . . , n i . . . , n j , . . . )] . (7)The energy required to transfer an electron from i to j is∆ E ji = E ( . . . , ¯ n i . . . , ¯ n j , . . . )) − E ( . . . , n i . . . , n j , . . . ) , = ǫ j − ǫ i + X m = i K jm n m − X m = j K im n m , = f E ij − f E ji . (8)where f E ij = ǫ j + X m = i K jm n m , = E j − K ji n i (9)and E j is the Hartree energy: E j = ǫ j + P m K jm n m and K jm = r jm is the Coulomb interaction term. Wecan then rewrite Eq.(8) as∆ E ji = E j − E i − K ij ( n i − n j ) , = E j − E i − K ij (10)hence we get ω i → j ω j → i = e − β ∆ E ji . So we choose our transi-tion probability as ω i → j = γ ij τ e β ∆ E ji + 1 (11)where τ is a hopping time scale. With this choice, masterequation takes the form ddt P ( { n l } ; t ) = − X i = j γ ( r ij )2 τ n i (1 − n j ) h f (∆ E ji ) P ( . . . , n i , . . . , n j , . . . ; t ) − f (∆ E ij ) P ( . . . , ¯ n i , . . . , ¯ n j , ... ; t ) i (12)Here f ( E ) = exp ( βE )+1 is the Fermi-Dirac distribution, γ ij = γ ( r ij ) is a factor independent of temperature, butdepends on the distance between sites i and j . For hop-ping electrons γ ij = γ e − r ij /ξ , where γ is a constant.From this, one can derive an equation for time-dependent averages or moments. To connect to the one-particle master equation, we consider N i ( t ) = X { n l } n i P ( n . . . n N ; t ) (13)whose time derivative gives ddt N i ( t ) = − τ X j = k γ ( r jk ) X { n l } n i n j (1 − n k )[ f (∆ E kj ) P ( n . . . n v ; t ) − f (∆ E jk ) P ( n . . . ¯ n j . . . ¯ n k , . . . ; t )] (14)Again, if i = j or i = k , a change of summation variables j ⇋ k makes the two terms cancel. The only survivingterm comes from i = j , Eq.(14) can now be written as ddt N i ( t ) = − τ X k = i γ ( r ik ) X { n l } n i (1 − n k ) h f (∆ E ki ) P ( . . . , n i . . . n k , . . . ; t ) − f (∆ E ik ) P ( . . . , ¯ n i . . . ¯ n k , . . . ; t ) i = − τ X k = i γ ( r ik ) [ h n i (1 − n k ) f (∆ E ki ) i t − h (1 − n i ) n k f (∆ E ik ) i t ] , (15) g ( E ) E T=1.0T=0.33T=0.20T=0.14T=0.10 g ( E ) E T=0.1g(E) = 1.4915 E
FIG. 1: (Colour online) (a) Histogram of the Hartreeenergies E (obtained using Eq. (19) where ǫ i werechosen from a box distribution of width ± W with W = 1) at different temperatures, for L = 16 and µ = 0.(b) Zoom in to low energies of the Hartree energies for T = 0 .
1. Solid line is the best-fit g ( E ) ∝ E . .where h ... i t denotes average at time t . The Eq.(15) isan exact equation. To get a closed set of equations, oneneeds to apply mean-field approximation to Eq.(15). Mean-Field Approximation
The mean-field approximation consists of making theassumption h f ( n ...n N ; t ) i = f ( N ( t ) ...N N ( t )) (16)With this assumption we get ddt N i ( t ) = − τ X k = i γ ( r ik ) [ N i (1 − N k ) f F D ( E k − E i ) − N k (1 − N i ) f F D ( E i − E k )] (17)where E i and E k are the Hartree energies at site i and k respectively. f F D ( E ) = 1 / ( exp [ βE ] + 1) is the FermiDirac distribution. Now let us linearize this equationabout an equilibrium solution: N i ( t ) = f i + δN i (18) E ei = ǫ i + X l K il f l (19)where f i = exp ( βE ei )+1 . Putting Eq.(16) and Eq.(17) intoEq.(15) one gets: ddt δN i = − τ X k = i γ ik " ( f i + δN i )(1 − f k − δN k ) f F D ( E ek − E ei + X l ( K kl − K il ) δN l ) − ( f k + δN k )(1 − f i − δN i ) f F D ( E ei − E ek + X l ( K il − K kl ) δN l ) (20)And the final linear equation using the detailed balanceis: ddt δN i = X k = i " δN i f i (1 − f i ) Γ ik − δN k f k (1 − f k ) Γ ki + 1 T X k = l,i Γ ik ( K kl − K il ) δN l = X l A il δN l (21)where we defineΓ ik = 12 τ γ ( r ik ) f i (1 − f k ) f F D ( E ek − E ei ) (22a)Γ ki = 12 τ γ ( r ki ) f k (1 − f i ) f F D ( E ei − E ek ) (22b) A ii = − X k = i Γ ik f i (1 − f i ) (22c) A il = Γ li f l (1 − f l ) + 1 T X k ( = l = i ) Γ ik ( K kl − K il ) (22d)It is easy to verify that Γ ik = Γ ki . Thus the final linearequation has the same form as the one used by Amir etal [56]. Here A is the linear dynamical matrix governingthe dynamics of the system near equilibrium and Γ ik arethe equilibrium transition rates. The transition rates asdefined in Eq.22(a) and Eq.22(b) can be written asΓ ik = γ exp (cid:18) − r ij ξ (cid:19) exp (cid:18) − T [ | E i | + | E j | + | E i − E j | ] (cid:19) , (23)when the energies | E i | , | E j | and | E i − E j | are greaterthan T . ξ =0.2 P ( λ ) - λ ξ =1.0 P ( λ ) - λ FIG. 2: Distribution of the eigenvalues of thedynamical matrix A obtained by solving Eq.(22), atdifferent temperature ( T = 0 .
33 in blue, T = 0 .
20 in redand T = 0 .
10 in black) for ξ = 0 . ξ = 1 . III. RESULTS AND DISCUSSIONS
In this paper we study a three-dimensional cubic lat-tice of localized states which have random energies andinteract through Coulomb interactions. We model thissystem by a Hamiltonian as defined in Eq.(1). We takethe number of electrons to be half of the total numberof sites in the lattice. All energies are noted in units of e /κa where a is the lattice constant. A. Coulomb Gap
The method-
To calculate the Hartree energy ( E i )given in Eq. (19), we have first calculated the magnetiza-tion, which, approximated within the mean-field theoryis defined as m i = tanh β E i + X k m k r ik ! (24)The above equation was solved self-consistently and thefinal m i ’s were then used to calculate E i ’s using f i =( m i + 1) /
2. We have annealed our data from T = 1 to T = 0 .
1, and the on-site energy ǫ i was chosen randomly P ( λ ) P ( λ ) - λ P ( λ ) ξ = 0.2; T = 0.1 ξ = 0.5; T = 0.1 ξ = 1.0; T = 0.1 FIG. 3: (Colour online) Comparison of the eigenvaluedistribution of the full A-matrix (in black) with theones obtained after neglecting the second term inEq.22(b) (in red) for ξ = 0 . , . , .
0, and T = 0 . ± W/ W = 1and β = 1 /T .It is well established now that in the CG model, a softgap, also called the Coulomb gap, is observed in single-particle DOS at low temperatures. The gap gets filledas the temperature increases. In this paper, the tem-peratures where the soft gap is well established are re-ferred to as low temperatures (i.e. T = 0 . − . g ( E ) ≈ E d − ind-dimensional CG model. In Fig.1(a), one can see forma-tion of a soft gap in the DOS at temperature lower than0.33. We further found that at T = 0 .
1, the DOS can bewell fitted by the relation g ( E ) ∝ E (see Fig.1(b)) assuggested by Efros and Shklovskii. B. Linear Dynamical Matrix
In Fig.2, we show the distribution of the eigenvaluesof a linear dynamical matrix ( A − matrix ) at differenttemperatures and localization lengths. The eigenvalues( λ ) here determine the rate of decay in the system. Withthe decrease in temperature, the shifting of λ towardszero indicates a slowing down of relaxation.For t ≫ λ − min , δn ( t ) behaves as e − λ min t . We now wantto look at the behavior of λ min as a function of tem-perature and localization length. Note that the interac-tion part in the A-matrix (second term in Eq.(22(d)))does not contribute much to the eigenvalue distributionat low temperatures as shown in Fig. 3 for all localizationlengths considered. In-fact for ξ = 1 and ξ = 0 .
5, theeigenvalue distributions (at low T) are mostly determinedby the diagonal part of the A − matrix . Consequently,the lowest eigenvalue of the dynamical matrix A ( λ min )approximately equals the smallest value of A ii (defined A minii ). In Fig.4, we find that A minii ∝ T for large ξ values. We now propose an argument for this behavior: | A ii | T ξ =1.0f(x)=170.898 T | A ii | T ξ = 0.5f(x)=13.9057* T FIG. 4: (Colour online) Lowest value of the diagonalpart of A − matrix at ξ = 0 . , . | A minii | ∝ T relation. A minii ∝ T : Using Eq.(22)c, we calculate the A ii ’s and find thatthey are smallest for sites around the Fermi level ( E i ≃ µ and so f i ≈ . A minii = − X k = i Γ ik = − X r X E electron e − r/ξ e − β | E | F ( r, E ) . (25)Here F ( r, E ) is the probability of finding an electron orhole having the Hartree energy E at a distance r from asite i . Since ξ is large, the electrons will hop to a site soas to minimize the factor ( r/ξ + β | E | ). This means thathops to r >> A minii ∝ − × X r e − r/ξ Z E max g ( E ) e − β | E | dE , (26)where g ( E ) is the density of states (DOS) of single-particle Hartree energies ( E ). As discussed earlier, ourresults (see Fig.1(b)) shows that g ( E ) ∝ E . . Substitut- F ( r = , E k ) k F ( r = , E k ) -2 -1 0 1 2E k (a) (b)(c) (d)T = 1.0 T = 0.33T = 0.10T = 0.20 FIG. 5: (Colour online) (a)-(d) Distribution of theHartree energy on site k , E k (where k are the 6 nearestneighbor sites of i ) at different temperatures, when theHartree energy on site i , E i are chosen from the interval[ − . , . A minii ∝ T . (27)We now look at the behavior of low temperature λ min values at small localization lengths ( ξ = 0 . , . . A ii distribution is different fromthe eigenvalue distribution, but the minimum value re-mains almost the same. More importantly, one shouldnote that the above arguments for A minii for the temper-ature range considered ( T = 0 . − .
2) does not workwell when the localization length is very small. Specif-ically, for small localization lengths the major contribu-tion to A minii comes from the nearest neighbor sites only(i.e. r = 1). So one has to find F ( r = 1 , E ) which isthe two particle nearest neighbor DOS and insert it intoEq.(25). In Fig.5 we show F ( r = 1 , E ) at different tem-peratures for 0 < E i ≤ − .
1. Unlike the full DOS plottedin Fig.1(a), there is a hard gap in F ( r = 1 , E ) for smallenergy electrons at low temperatures. This is not surpris-ing, since if one was working with true ground state, thenthere is a hard gap extending to E ≈ F ( r = 1 , E ).The reason behind it is that the ground state is stableagainst any single electron-hole transition which implies E h − E e − /r eh >
0. This means that | E h | + | E e | > ξ ≪ A minii ≈ − e − /ξ X j e − ∆ E ij /T (28)where ∆ E ij is the energy difference between site i andits nearest neighbors j . In Fig.6, we have shown that A minii ( T ) indeed follows the above relation at small lo-calization lengths. Thus, our analysis of the matrix A ii -28-27.5-27-26.5-26-25.5-25-24.5-24-23.5 0.08 0.1 0.12 0.14 0.16 0.18 0.2 ξ =0.05ln(A ii )=-19.9022 -(0.7714/T)-16-15.5-15-14.5-14-13.5-13-12.5 0.08 0.1 0.12 0.14 0.16 0.18 0.2 l n ( A ii ) T ξ = 0.1ln(A ii )= -9.6141 -(0.6046/T)-10-9.5-9-8.5-8-7.5-7-6.5 0.08 0.1 0.12 0.14 0.16 0.18 0.2 l n ( A ii ) T ξ =0.2ln(A ii )=-3.9218 -(0.5596/T) FIG. 6: (Colour online) Lowest value of the diagonalpart of A − matrix at ξ = 0 . , . , . ln ( | A minii ) | = a + bT relationfor small ξ values.shows that at low temperatures λ min obeys different scal-ing laws for small and large localization lengths.We now look at the behavior of the system for T =0 . t < λ min ( λ > λ min ) for different local-ization lengths. When the localization length is large( ξ = 0 . , .
0) we find that P ( λ ) is almost flat. Thisis shown in Fig.7. This implies an exponential decayfor δn ( t ). δn ( t ) ∼ exp ( − λ t ), where λ is the smallesteigenvalue at which the flat region starts. For small lo-calization lengths ( ξ = 0 . , . P ( λ )vs λ are shown in Fig.8(a-b). For ξ = 0 .
05, one sees sharppeaks at ln ( λ ) = −
20 and − .
3. Since ξ is small, therelaxation is dominated by the nearest neighbor hopping.For all sites for which k = 0 nearest neighbor hops witha decrease in energy are available, A ii is given by A ii = k exp (cid:18) − ξ (cid:19) , (29) ξ =0.5 P ( λ ) - λ ξ =1.0 P ( λ ) - λ FIG. 7: Distribution of the eigenvalues of thedynamical matrix A obtained by solving Eq.(22), forlarge localization lengths at T = 0 . ξ = 0 .
05 we find A ii = exp ( −
20) for k = 1 and A ii = exp ( − .
3) for k = 2, which correspond to the peaks at ln ( λ ) = −
20 and − . ξ = 0 .
06 as shown in Fig.8(b). Soat short times, δn ( t ) will decay according to P p e − λ p t ,where λ p correspond to eigenvalues at which P ( λ ) haspeaks.When nearest neighbor hops, which lead to decreasein energy (∆ E ≤
0) are not possible, one would get atransition to nearest neighbor site with ∆
E >
0. In thiscase A ii can be written as A ii = X j e − r ij /ξ e − ∆ E/T . (30)For the λ − range considered in Fig.8(c-d), we find that A ii ∼ λ . So, using Eq.(30) P ( λ ) = Z ∆ E max ∆ E min δ ( λ − ce − ∆ E/T ) P (∆ E ) d (∆ E ) (31)where c = e − /ξ , ∆ E min ≈ T and ∆ E max ≈ T . Ap-proximating P (∆ E ) (for ∆ E min < ∆ E < ∆ E max ) by auniform distribution on gets P ( λ ) ∼ λ . (32)
10 12 14 16 18 20 22 24-25 -24 -23 -22 -21 -20 -19 -18(a) l n ( P ( λ )) ln( λ ) ξ =0.05 8 10 12 14 16 18 20-22 -21 -20 -19 -18 -17 -16 -15(b) l n ( P ( λ )) ln( λ ) ξ =0.06 19 19.5 20 20.5 21 21.5 22 22.5-25 -24.5 -24 -23.5 -23 -22.5 -22(c) l n ( P ( λ )) ln( λ ) ξ =0.05f(x) = -0.9398x - 1.5195 15.5 16 16.5 17 17.5 18 18.5 19-22 -21.5 -21 -20.5 -20 -19.5 -19 -18.5(d) l n ( P ( λ )) ln( λ ) ξ =0.06f(x) = -0.9288x - 1.4709 FIG. 8: (a)-(b) Log-log plot of the distribution of theeigenvalues of the dynamical matrix A obtained bysolving Eq.(22), for small localization lengths at T = 0 .
1. (c)-(d) For a certain range of λ (see text), ln ( P ( λ )) vs ln ( λ ) has a linear fit.In Fig.8(c-d), we plot P ( λ ) vs λ for the regime wherenearest neighbor activated hoping takes place. In thisregime we find P ( λ ) ∼ /λ . This leads to logarithmictemporal dependence of the relaxation for intermediatetimes. Recently, a crossover from logarithmic time de-pendence to an exponential dependence was shown in anon-equilibrium study [51] of excess conduction ∆ G ( t ) indisordered indium oxide. In this it was shown that as oneapproaches metal insulator transition from the insulatingside, the crossover time becomes smaller. This impliesthat as the disorder in the system decreases and localiza-tion length increases the crossover time to exponentialdecay decreases. Since A minii is equal to λ min , Eq.(28)shows that the λ min increases as localization length in-creases. This implies that τ max = 1 /λ min will decreaseas the localization length increases and crossover fromthe logarithmic behavior to exponential decay happensfaster.At intermediate localization lengths ( ξ = 0 . , .
2) onesees that for large λ ’s, there are peaks corresponding tonext nearest neighbor hops with decrease in energy. For λ = -15 to -13 at ξ = 0 . λ = -9 to -7 at ξ = 0 . ln ( P ( λ )) vs ln ( λ ) has a linear fit but the slope is not equalto −
1. The reason is that for intermediate ξ ’s there iscontribution to λ from next nearest neighbor hops as wellas nearest neighbor hops. In Fig.9(b,e) we have plotted P ( λ ) vs λ for these regions. We get P ( λ ) = a + b/λ withvalue of a ≫ b for both ξ = 0 . ξ = 0 .
2. This formof P ( λ ) implies relaxation behavior of the form a e − λ t + b ln ( t ) where ( λ is the minimum value of λ for the rangeunder consideration). Since a ≫ b , it is quite possiblethat exponential decay will overshadow the logarithmicdecay in relaxation of δn ( t ). l n ( P ( λ )) ln( λ ) ξ =0.1 20000 30000 40000 50000 60000 70000 80000 90000 100000 1 × -6 × -6 × -6 × -6 × -6 × -6 × -6 (b) P ( λ ) - λ ξ =0.1f(x) = 19946.2 + (0.0161/x) 8.5 9 9.5 10 10.5 11 11.5-15.5 -15 -14.5 -14 -13.5 -13(c) l n ( P ( λ )) - . ln( λ ) ξ =0.1f(x) = -1.0x - 4.0327-4-2 0 2 4 6-10 -9 -8 -7 -6 -5 -4 -3(d) l n ( P ( λ )) ln( λ ) ξ =0.2 40 50 60 70 80 90 100 110 0.0005 0.001 0.0015 0.002 0.0025(e) P ( λ ) - λ ξ =0.2f(x) = 38.9002 + (0.0126/x) 2.5 3 3.5 4 4.5 5-9 -8.5 -8 -7.5 -7(f) l n ( P ( λ )) - . ln( λ ) ξ =0.2f(x) = -1.0x - 4.4397 FIG. 9: (a),(d) Log-log plot of the distribution of the eigenvalues of the dynamical matrix A obtained by solvingEq.(22), for intermediate localization lengths at T = 0 .
1. (b),(e) Behavior of P ( λ ) at intermediate times. The solidlines shows that the data fits well with P ( λ ) ∝ a + bλ relation. (c),(f) The fit to the logarithmic distribution aftersubtracting the background constant is close to linear. IV. SUMMARY
We consider here the relaxation properties of the three-dimensional Coulomb Glass lattice model in which allthe electron states are localized and the dynamics occursthrough phonon-assisted hopping among these states.The master equation governing the dynamics of the sys-tem is approximated via mean-field theory.The relaxation law for a range of localization lengthsis studied here. The dependence of the relaxation on thelocalization length can be summarized as follows:(i) For small localization length ξ ≪
1, near neighborhopping is strongly dominant. This results in P ( λ ) ∼ /λ distribution for small λ ’s which leads to a logarithmictemporal dependence of the relaxation at intermediatetimes.(ii) For intermediate ξ values (0.1,0.2) next nearestneighbor (n.n.n) contribution also becomes important inrelaxation. We find that the relaxation is not purely log-arithmic at intermediate times ( P ( λ ) ∝ a + b/λ , a ≫ b ).(iii) For larger values ( ξ = 0 . /λ distribution is seen, consequently no logarithmic tempo-ral dependence.Finally, we looked at relaxation for times t ≫ λ − min .We have found that although the full eigenvalue distri-bution is not much affected by the Coulomb interactionterm in the linear dynamical matrix, one can gain a bet-ter understanding of the behavior of low-temperature dy-namics by looking at the role of the gap in the densityof states in the decay process and the range of hop-ping. The gap in the density of states exists due to the long-range nature of Coulomb interactions and sothe interactions play an important role in the relaxationprocess. For small localization lengths one finds thatthe λ min ∝ e − c/T where c is a constant. This impliesthat time at which exponential decay starts increases ex-ponentially with a decrease in temperature. This mayexplain why the transition from logarithmic decay to ex-ponential decay is not seen in most experiments.Recently, a non-equilibrium study [51] of excess con-duction ∆ G ( t ) in disordered indium oxide showed acrossover from logarithmic time dependence to an expo-nential dependence. The crossover time became smalleras the metal insulator transition was approached fromthe insulating side. This implies that as the disorderin the system decreases and localization length increasesthe crossover time to exponential decay decreases. Inour formalism, the crossover time is 1 /λ min which alsodecreases with increase in localization length.Further work is required to establish results in the casewhere the distance between sites is a continuous variableand not a discrete value (as was the case in the presentwork). ACKNOWLEDGEMENT
PB gratefully acknowledges IISER Mohali and the Is-rael Science Foundation (Grant No. 2300/19) for thefinancial support. Illuminating discussions with A. Amirand Z. Ovadyahu are gratefully acknowledged. We wishto thank NMEICT cloud service provided by BAADALteam, cloud computing platform, IIT Delhi for the com-putational facility and Ben Gurion University of theNegev for access to their HPC resources. [1] J. H. Davies, P. A. Lee and T. M. Rice, Phys. Rev. Lett. , 758 (1982).[2] M. Gr¨unewald, B. Pohlman, L. Schweitzer and D. W¨urtz,J. Phys. C , L1153 (1982).[3] M. Pollak and M. Ortu˜no, Sol. Energy Mater. , 81(1982).[4] M. Pollak, Philos. Mag. B , 265 (1984).[5] E. R. Grannan and C. C. Yu, Phys. Rev. Lett. , 2934(1994).[6] A. L. Efros and B. I. Shklovskii, J. Phys. C: Solid StatePhys. , L49 (1975).[7] A. M¨obius, M. Richter and B. Dritter, Phys. Rev. B. ,11568 (1992).[8] P. Bhandari, V. Malik and S. R. Ahmad, Phys. Rev. B. , 184203 (2017).[9] P. Bhandari and V. Malik, J. Phys.: Condens. Matter , 485402 (2017).[10] S. D. Baranovskii, A. L. Efros, B. L. Gelmont and B. I.Shklovskii, J. Phys. C: Solid State Phys. , 1023 (1979).[11] A. Glatz, V. M. Vinokur, J. Bergli, M. Kirkengen and Y.M. Galperin, J. Stat. Mech. P06006 , (2008).[12] B. I. Shklovskii and A. L. Efros, Electronic Propertiesof Doped Semiconductors, Heidelberg: Springer, Heidel-berg, (1984)[13] M. Pollak, M. Ortu˜no, and A. Frydman, The ElectronGlass, Cambridge University Press, New York, (2013).[14] M. Pollak, Discuss. Faraday Soc. , 13 (1970)[15] G. Srinivasan, Phys. Rev. B , 2581 (1971).[16] J. H. Davies, P. A. Lee, and T. M. Rice Phys. Rev. B ,4260 (1984).[17] J. G. Massey and M. Lee, Phys. Rev. Lett. , 4266(1995).[18] V. Y. Butko, J. F. DiTusa and P. W. Adams, Phys. Rev.Lett. , 1543 (2000).[19] M. M¨uller and L. B. Ioffe, Phys. Rev. Lett. , 256403(2004).[20] N. F. Mott, J. J. Non-Cryst. Solids , 1 (1968).[21] N. F. Mott, Phil. Mag. B , 835 (1969).[22] M. Goethe and M. Palassini, Phys. Rev. Lett. ,045702 (2009).[23] B. Surer, H. G. Katzgraber, T. G. Zimanyi, A. B. Allgoodand G. Blatter, Phys. Rev. Lett. , 067205 (2009).[24] A. M¨obius and M. Richter, Phys. Rev. Lett. , 039701(2010).[25] A. Barzegar, J. C. Anderson, M. Schechter and H. G.Katzgraber, Phys. Rev. B. , 104418 (2019).[26] A. A. Pastor and V. Dobrosavljevi´c, Phys. Rev. Lett. ,4642 (1999).[27] S. Pankov and V. Dobrosavljevi´c, Phys. Rev. Lett. ,046402 (2005).[28] M. M¨uller and S. Pankov, Phys. Rev. B. , 144201(2007).[29] A. J. Bray and M. A. Moore, J. Phys. C: Solid StatePhys. , 2417 (1982).[30] L. C. E. Struik, Physical Aging in Amorphous Polymersand Other Materials (Elsevier, Amsterdam, 1978).[31] J. P. Bouchaud, L. F. Cugliandolo, and J. Kurchan, in
Spin Glasses and Random Fields , edited by A. P. Young(World Scientific, Singapore, 1997) [32] L. F. Cugliandolo, in
Slow Relaxation and Nonequilib-rium Dynamics in Condensed Matter , (Les Houches Ses-sion LXXVII, 1-26 July, 2002).[33] E. Vincent et al., in