# Numerical study of relaxation in electron glasses

## I Introduction

Strongly localized systems are characterized by very slow relaxation rates due to the exponential dependence of the transition rates on hopping length . For a wide range of parameters, the typical times involved are much larger than the experimental times and a glassy behavior is observed. Ben–Chorin et al.  reported on non-ergodic transport in Anderson localized films of indium–oxide and ascribed the phenomena to the hopping transport in non-equilibrium states. Ovadyahu and Pollak  performed further experiments on this system that clearly demostrate the glassy nature of Anderson insulators. Glassy behavior may be obtained independently of the strength of interactions and regardless of their long or short range. In systems with localized states, long hopping lengths result in very long relaxation times. However, it is thought that there are specific features of the glassy relaxation behavior that indeed depend on the type and strength of the interactions involved. If so, relaxation experiments could be an adequate tool for studying the strength of interactions. There has been no systematic study of the effects of interactions on the relaxation properties of strongly localized systems, and in this paper we try to fill this gap as much as possible.

Most properties of systems with localized electronic states strongly depend on interactions. This is especially true for Coulomb glasses where interactions are of a long range character. The non–equilibrium properties of these systems are affected by dynamic correlations in the motion of electrons . One–particle densities of states or excitations are not enough to encompass the whole problem. To deal with such problems, methods were developed  to obtain the low lying states and energies of electron glasses. The states of the system, their energies and the transition rates between them constitute the information needed to compute non-equilibrium properties. We use this information to study energy relaxation for systems with no interactions, with long-range Coulomb interactions and with short-range interactions.

In the next section, we describe the model and the numerical procedure used. In section III, we study the temporal dependence of energy relaxation and, in section IV, we calculate the largest relaxation time $`\tau _2`$ and its dependence on size and temperature. Finally, in section V, we present results about the number of electrons participating in low- energy relaxation processes.

## II Model and numerical procedure

We consider three–dimensional systems in the strongly localized regime, in which quantum overlap energies, $`h`$, arising from tunnelling are much smaller than the other important energies in the problem and are taken into account only to the lowest contributing order, i.e., to zero order for energies and to first order for transition rates. Spin is neglected since exchange energies are proportional to $`t^2`$. We use the standard tight–binding Coulomb gap Hamiltonian :

$$H=\underset{i}{}ϵ_in_i+\underset{i<j}{}n_in_jV_{ij},$$
(1)   
where $`ϵ_i`$ is the random site energy chosen from a box distribution with interval $`[W/2,W/2]`$. For non-interacting systems $`V_{ij}=0`$, while $`V_{ij}=1/r`$ for systems with Coulomb interactions and $`V_{ij}=(0.7/r)^4`$ is the potential chosen for short range interactions. The large value of the Hubbard energy is accounted for by disallowing double occupation of sites.

We study systems with sizes from 248 to 900 sites placed at random (for short range interactions we only consider systems sizes up to 465 sites), but with a minimum separation between them, which we choose to be $`0.5l_0`$ where $`l_0=(4\pi N/3)^{1/3}`$ and $`N`$ is the concentration of sites. We take $`e^2/l_0`$ as our unit of energy and $`l_0`$ as our unit of distance. We choose the number of electrons to be equal to half the number of sites. We use cyclic boundary conditions.

We use two different numerical algorithms to obtain the ground state and the lowest energy many–particle configurations of the systems up to a certain energy. For short range interactions, we employ an algorithm that relaxes the system through certain simultaneous $`n`$–electron transitions . The procedure is repeated for different initial random configurations of the charges until the configuration of lowest energy is found ten times. The configurations thus generated were memorized in terms of site occupation numbers and of energy, whenever this was less than the highest energy configuration in memory storage. We complete the set of low–energy configurations by generating all the states that differ by one– or two–electron transitions from any configuration stored.

For long range interactions, we use an algorithm that consists of finding the low-energy many-particle configurations by means of a three-step algorithm . This comprises local search , thermal cycling  and construction of “neighbouring” states by local re-arrangements of the charges . The efficiency of this algorithm is discussed in Ref. . In the first step, an initial set, $`𝒮`$, of metastable low-energy many-particle states is created. We start from states chosen at random. These states are relaxed by a local search algorithm which ensures stability with respect to excitations from one to four sites. In the second step, this set $`𝒮`$ is improved by means of the thermal cycling method, which combines the Metropolis and local search algorithms. Lastly, the third step completes the set $`𝒮`$ by systematical investigations of the surroundings of the states previously found.

The transition rate $`\omega _{IJ}`$ between configurations $`I`$ and $`J`$ is taken to be

$$\omega _{IJ}=\frac{1}{\tau _0}\mathrm{exp}\left(2r_{ij}/a\right)\mathrm{exp}\left(\frac{E_JE_I}{kT}\right)$$
(2)   
for $`E_J>E_I`$, and without the second exponential for $`E_J<E_I`$. In this equation, $`\tau _0`$ is the inverse phonon frequency, of the order of $`10^{13}`$ s, $`a`$ is the localization radius, which we take equal to $`0.3l_0`$, and $`r_{ij}`$ is the minimized sum of the hopping lengths of the electrons participating in the transition.

The relaxation process is governed by the master equation, which in first order can be written in matrix form as $`𝐩(t+\delta t)=𝐩(t)`$, where $`𝐩`$ is the vector of occupation probabilities in the configuration space, and $``$ the matrix of transition probabilities between states during a time, $`\delta t`$, given by :

$$()_{JI}=\{\begin{array}{cc}\omega _{IJ}\delta t\hfill & \mathrm{for}IJ,\hfill \\ 1_{KI}\omega _{IK}\delta t\hfill & \mathrm{for}I=J.\hfill \end{array}$$
(3)   
We assume that the system initially occupies a set, $`𝒦`$, of $`m`$ configurations with equal probabilities, that is, $`p_K^{(0)}=1/m`$ for $`K𝒦`$, and $`p_L^{(0)}=0`$ for all other $`L`$. The time evolution of $`𝐩`$ is governed by the eigenvalues $`\lambda _i`$ and right eigenvectors $`\stackrel{}{\varphi }_i`$ of $``$. We will assume that the $`\lambda _i`$ are arranged in decreasing order. Rewriting $`𝐩^{(0)}`$ as a linear combination of the $`\stackrel{}{\varphi }_i`$, the probability vector after $`n`$ time steps $`𝐩^{(n)}`$ is given by

$$𝐩^{(n)}=a_1\stackrel{}{\varphi }_1+a_2\stackrel{}{\varphi }_2\lambda _2^n+a_3\stackrel{}{\varphi }_3\lambda _3^n+\mathrm{}$$
(4)   
where $`a_i`$ is the $`i`$-th component of $`𝐩^{(0)}`$ in the basis $`\left\{\stackrel{}{\varphi }_i\right\}`$. At long times (large $`n`$), Eq. (4) approaches equilibrium with time dependences given by $`\lambda _i^n`$. Thus, the relaxation times are given by

$$\tau _i=\frac{1}{|\mathrm{ln}\lambda _i|}$$
(5)   
in units of $`\delta t`$. The final state is $`p_M^{(\mathrm{})}=\mathrm{exp}(E_M/kT)/Z`$ for all $`M`$, where $`E_M`$ is the energy of state $`M`$, and $`Z`$ is the partition function. Clearly $`𝐩^{(\mathrm{})}`$ is a right eigenvector of $``$ with eigenvalue 1, since $`𝐩^{(\mathrm{})}=𝐩^{(\mathrm{})}`$. All the other eigenvalues of $``$ are smaller than 1, since otherwise the system would not tend to the stationary probability distribution. The second largest eigenvalue corresponds to the largest relaxation time of the system. The addition of the other eigenvectors to $`\varphi _1=𝐩^{(\mathrm{})}`$, transfers $`𝐩`$ from high energy states to low energy states at various rates.

We have developed a renormalization method to be able to properly handle the huge range of transition rates involved. Large values of $`\tau _i`$ correspond to $`\lambda _i`$ with values which are very close to unity, Eq. (5), and a direct calculation of $`\tau _i`$, in units of $`\delta t`$, is strongly limited by the numerical precision of the computer. In order to minimize errors, we must choose a $`\delta t`$ which is as large as possible, although this soon yields negative diagonal elements of $``$. We overcome this problem using a renormalization procedure that allows us to increase $`\delta t`$ and to simultaneously keep all terms of $``$ positive. This procedure forms groups of configurations. Each group is made up of configurations connected between themselves by transition rates which are larger than a critical one. The groups are clusters in local equilibrium for times greater than the inverse of the critical transition rate. Firstly, we take a critical transition rate $`\omega _\mathrm{c}`$. Then for each $`\omega _{IJ}`$ larger than $`\omega _\mathrm{c}`$ we define a new equilibrium state, $`M`$, and substitute the original configurations, $`I`$ and $`J`$, by this new state, $`M`$. The transition rates between $`M`$ and any other configuration $`K`$ ($`KI,J`$) are defined as:

$`\omega _{KM}=\omega _{KI}+\omega _{KJ}`$ (6)   
$`\omega _{MK}={\displaystyle \frac{\omega _{IK}}{1+R_M}}+{\displaystyle \frac{\omega _{JK}}{1+R_M^1}}`$ (7)   
where $`R_M`$ is given by:

$$R_M=\frac{\omega _{IJ}}{\omega _{JI}}=\mathrm{exp}\left\{(E_IE_J)/k_\mathrm{B}T\right\}.$$
(8)   
The diagonal matrix elements $`\omega _{MM}`$ are again equal to 1 minus the sum of the non-diagonal elements of the column $`M`$ multiplied by $`\mathrm{\Delta }t`$.

After the matrix $``$ has been renormalized by the previous procedure, we can increase the time scale to a larger interval $`\delta ^{}t=1/\omega _\mathrm{c}`$. With this $`\delta ^{}t`$ we calculate the new elements of $``$. The eigenvalues of the transition matrix will be given now in units of $`\delta ^{}t(>\delta t)`$. We have checked the validity of our renormalization procedure with several samples of small systems where errors are not critical. The method minimizes computer errors in the solution of the eigenproblem as the matrix becomes less ill-conditioned, and allows us to consider large systems, with matrix elements that differ by many orders of magnitude.

## III Temporal dependence

We calculate the temporal dependence of the energy of the system when it relaxes from an initial set of high energy configurations. At very long times, the longest relaxation process involved predominates and we see an exponential relaxation. For shorter times, there is an almost continuous sequence of relaxation times, which gives rise to a power law relaxation $`(EE_{\mathrm{eq}})t^\alpha `$. To obtain the exponent of this law it is convenient to represent the absolute derivative of the energy with respect to time. In Fig. 1 we show $`|\mathrm{d}E/\mathrm{d}t|`$ versus time (in units of $`\tau _0`$) in a double log<sub>10</sub> plot for a sample with Coulomb interactions and 248 sites. The continuous curve corresponds to a temperature $`T=0.004`$, and the dashed curve to $`T=0.005`$. The straight line is a fit to the data in the non–exponential part of both curves, and its slope is equal to $`1.15`$. So the power–law exponent for relaxation is $`\alpha =0.15`$. This exponent is basically independent of temperature for all the systems considered.

We have also studied energy relaxation for systems with short–range interactions and for non–interacting systems. The results for short–range interactions are very similar to those for Coulomb interactions. The power–law exponent is roughly 0.15 and the largest relaxation time is of the same order of magnitude as for Coulomb systems. In Fig. 2 we show $`|\mathrm{d}E/\mathrm{d}t|`$ as a function of time in a double log<sub>10</sub> plot for a non–interacting system with $`N=248`$ sites. The continuous curve is for $`T=0.004`$, and the dashed curve for $`T=0.005`$. The slope of the straight line is again equal to $`1.15`$. There are two differences between the results for interacting and for non–interacting systems. The longest relaxation times are shorter for the latter, and the power–law regime is not very well defined in the absence of interactions. Both figures give the rate of relaxation $`|\mathrm{d}E/\mathrm{d}t|`$ at any time. At very small $`t`$, the interacting systems relax faster than the non–interacting systems. A possible explanation of this is that in the excited state of the interacting systems some electrons get very close to each other. In the initial stages of relaxation these electrons hop away from electrons in the nearest neighbors sites, the whole process being very fast.

Several samples have been checked and in all of them we obtain similar results to Figs. 1 and 2. Two features characterize our relaxation process, the exponent $`\alpha `$ of the power–law regime and the longest relaxation time. The exponents $`\alpha `$ do not appreciably vary from sample to sample, nor with temperature or with the type of interaction. On the other hand, the longest relaxation time drastically changes from sample to sample and with changes in temperature, size and the range of interaction. On average, this time increases with the size of the system and with the strength of the interactions. In the next section we study the longest relaxation time in detail. Now we shall analyze exponent $`\alpha `$.

Temporal relaxation can be described as a sum of parallel exponential relaxation processes, each with its own different relaxation time, $`\tau _i`$. The energy, $`E`$, of the system can be written as a function of time, $`t`$, as follows

$$E(t)=\underset{i>1}{}c_i\mathrm{exp}\left(\frac{t}{\tau _i}\right)+E_{\mathrm{Eq}}$$
(9)   
where $`c_i`$ is the product of the $`i`$-th component of the initial occupation vector, $`a_i`$, and the energy associated to the eigenvector $`\varphi _i`$. This energy is the sum of the components of $`\varphi _i`$ multiplied by the corresponding energies. $`E_{\mathrm{Eq}}`$ is the equilibrium energy, i. e. $`E_{\mathrm{Eq}}=E(t\mathrm{})`$. In Fig. 3 we plot $`(c_i/\tau _i)\mathrm{exp}(t/\tau _i)`$ for the 30 largest eigenvalues of $``$, excluding $`\lambda _1=1`$, as a function of time for a sample with Coulomb interactions and of size $`N=465`$. The solid line represents the temporal derivative of the actual energy as a function of time. This curve is below, but very close to, the envelope of the curves corresponding to the individual relaxation processes. Note how the combination of several simple exponential relaxation processes gives rise to a power law relaxation.

Surprisingly, $`\alpha `$ is fairly independent of temperature, size, type of interaction considered, and localization radius, facts for which we do not have any interpretation. Anyway, the robustness of the exponent could be a signature of self–organized criticality. Similar trends have been found in experimental measurements of the excess conductance of 2D samples excited far from equilibrium . In the absence of magnetic field, the power law exponent of these measurements ranges between 0.27 and 0.29, diminishing with the strength of the magnetic field.

Our results point to the difficulty in extracting information about the effects of interactions from the power law exponent. Nevertheless, the type of interaction significantly affects the longest relaxation times.

## IV Longest relaxation time

We also study the longest relaxation time, $`\tau _2`$, as a function of temperature and the size of the sample for systems with Coulomb interactions, with short-range interactions and for non-interacting systems. In Fig. 4 we plot $`\mathrm{log}_{10}\tau _2`$ versus the inverse of the temperature for the three types of interactions mentioned, Coulomb (solid lines), short-range (dotted-dashed lines) and no-interactions (dashed lines). The number of sites considered are $`N=248`$, 341, 465, 744 and 899, for long range interactions and for non-interacting systems; for short range interactions we did not use the two largest sizes. $`\tau _2`$ increases with sample size, and thus the smallest sample corresponds to the lowest curve, and so on. $``$ denotes averages over site configurations. Fluctuations in $`\tau _2`$ from sample to sample are very large and, as is the case with most properties of disordered systems, one has to average the common logarithm of $`\tau _2`$, rather than $`\tau _2`$ itself. The curves extend over the range of validity of the results. The ’high’ temperature limit $`T_{\mathrm{max}}`$ depends on the energy range $`\mathrm{\Delta }E`$ spanned by the configurations stored. We choose $`T_{\mathrm{max}}=0.1\mathrm{\Delta }E`$. The low temperature limit arises from the discrete nature of the spectrum of configurations and we take it as being equal to the mean energy spacing of the ten lowest energy configurations $`\mathrm{\Delta }ϵ`$.

From Fig. 4 we can conclude that the longest relaxation time depends strongly on the type of interaction. $`\tau _2`$ is one order of magnitude larger for interacting than for non-interacting systems. As we will see, this effect is much larger when extrapolated to macroscopic sizes.

In order to extrapolate the previous results to macroscopic sizes we plotted $`\mathrm{log}_{10}\tau _2`$ as a function of $`L^\beta `$ at a fixed temperature for different values of the exponent $`\beta `$. $`L=N^{1/3}`$ is the length of the side of the system, and $`N`$ is the number of sites. We found that the results for the three types of interactions fit straight lines fairly well when $`\beta =1`$. In Fig. 5 we show $`\mathrm{log}_{10}\tau _2`$ versus $`L^1`$ for systems with Coulomb interactions (dots), short-range interactions (diamonds) and without interactions (squares). The horizontal dashed line represents a macroscopic time, say, one day ($`10^{18}\tau _0`$). The temperature chosen in this plot is $`T=0.0025`$, which is valid for the four sizes employed in both types of interactions. The size of the symbols used roughly corresponds to the standard deviation of $`\mathrm{log}_{10}\tau _2`$. The crossing point of each straight line with the vertical axis is the extrapolation of $`\tau _2`$ to macroscopic sizes. The results are $`\tau _2^{(\mathrm{})}10^{31\pm 1}\tau _0=10^{18\pm 1}`$ s (Coulomb interactions) $`\tau _2^{(\mathrm{})}10^{11\pm 1}`$ s (short-range interactions) and $`\tau _2^{(\mathrm{})}10^{5\pm 1}`$ s (no interactions). It is clear from this figure that the longest relaxation time drastically increases with the strength of interactions, although these results have to be taken with care as they are extracted from a very long extrapolation.

The results presented in Figs. 4 and 5 correspond to a localization radius $`a=0.3l_0`$. For larger values of $`a`$, the relaxation times will decrease, as can be deduced from Eq. (2). We have found empirically that a change in $`a`$ causes a change in $`\tau _2`$ of approximately $`\mathrm{\Delta }\mathrm{log}_{10}\tau _23\mathrm{\Delta }(a^1)`$. The values of $`\tau _2^{(\mathrm{})}`$ are so large for interacting systems that we would expect non–ergodic behaviour for these systems even for much larger localization radii than the one considered here.

## V Variable number relaxation

At zero temperature, the relaxation process is downward in energy and we can assume that the fastest process always dominates, corresponding to a well defined sequence of configurations with decreasing energies. For each transition at $`T=0`$, the shorter the hopping length, the faster the corresponding transition rate. From each configuration, the system chooses the nearest one (in terms of $`r`$) from those with less energy. With this in mind, we have computed for all low-energy configurations the closest one of smaller energy, and have stored the number of electrons $`n`$ participating in the transition.

In Fig. 6 we show the number of electrons, $`n`$, of the fastest transition from an initial configuration as a function of the number of this configuration for a Coulomb interacting sample with 900 sites. At very low energies, the relative importance of many-electron transitions increases. The proportion of transitions with a fixed number of electrons greater than one ($`n>1`$) increases with decreasing energy. Obviously, in the non-interacting case all processes are one-electron transitions.

## VI Conclusions

Our numerical results of relaxation in localized electronic systems show a power law behavior with an exponent close to 0.15 and independent of all the parameters and type of interactions considered. At very long times, we obtain exponential relaxation with a characteristic time that strongly varies with size, localization radius and type of interaction. The extrapolation of this characteristic time to macroscopic sizes predicts values much larger than the typical experimental times, especially for the interacting cases. The strength of interactions in experiments performed on these systems can be deduced from their longer relaxation times.

## VII Acknowledgments

We would like to acknowledge Prof. M. Pollak for useful conversations and a critical reading of the manuscript. We also acknowledge the Dirección General de Investigación Científica y Técnica for financial support, project number PB 96/1118, and for APG’s grant.