# Electronic spectrum and hopping conductivity in highly doped lattice systems

## 1 Introduction

The interest to electronic processes in disordered systems was greatly inspired by the fascinating disorder effects in semiconductors, including doped and amorphous ones , and in mesoscopic metallic systems . Unlike the Bloch states in fully periodic systems, the electronic spectrum of disordered systems generally includes both extended and localized states, their coexistence being related to the competition between kinetic and potential energy of Fermi particles. The main consequence of this competition is the possibility for Anderson transition from metallic to insulating state at zero temperature and sufficiently strong disorder . The best studied situation is that of non-interacting electrons (that is, the single-electron approximation) in a certain random field. Efficiency of single-electron theories for each type of electronic states in disordered systems (at the Fermi level $`ϵ_F`$, they are mostly extended in metals and mostly localized in semiconductors) is assured by the presence of a certain small parameter, such as $`\lambda _F/\mathrm{}`$ (where $`\lambda _F`$ is the Fermi wavelength, $`\mathrm{}`$ the mean free path) in metals  and $`na^3`$ in semiconductors  ($`n`$ is the concentration of charge carriers, $`a`$ the lattice parameter).

The following discussion is addressed to the doped semiconducting systems. In traditional semiconductors, typical values of $`na^3`$ do not exceed $`10^4÷10^6`$. This is determined by a very great localization radius $`r_0(20÷50)a`$ of a single localized shallow state and the Mott criterion for metallization in the impurity band: $`nr_0^30.02`$ . At doping levels below this value, the single-electron approach yeilds in a finite Fermi density of (localized) states $`\nu _F=\nu (ϵ_F)`$ and in the Mott law for hopping conductivity vs temperature $`T`$: $`\sigma (T)\{\mathrm{exp}\}(BT^{1/4})\text{ }`$, $`B2.1(r_0^3\nu _F)^{1/4}`$ . However, it was shown by Efros and Shklovskii  that account of Coulomb interactions between such shallow states leads to formation of a “soft gap” in the unperturbed density of states $`\nu _F`$ near the Fermi level so that: $`\nu (ϵϵ_F)^2/e^6`$ until $`|ϵϵ_F|\mathrm{\Delta }=e^3\nu _F^{1/2}/\kappa ^{3/2}`$ (where $`e`$ is the electron charge and $`\kappa `$ the static dielectric constant). Consequently, the Mott law is changed to: $`\sigma (T)\mathrm{exp}(B^{}T^{1/2})`$, $`B^{}e/(\kappa r_0)^{1/2}`$ at sufficiently low temperatures, $`T<T_ce^4r_0\nu _F/\kappa ^2`$. The above theoretic dependencies are in a good agreement with the bulk of experimental data in traditional semiconductors.

A special class of doped materials, displaying semiconducting, metallic, superconducting and various magnetic phases, is comprised by the doped perovskite systems . From the point of view of standard theory of semiconductors, these materials exhibit extremely high values of $`na^30.1÷0.5`$, that is more than three orders of magnitude higher than those observed in traditional semiconductors. So the experimentally observed hopping type of conductance, at sufficiently high temperatures, in perovskite manganites Ln<sub>1-x</sub>A<sub>x</sub>MnO<sub>3</sub> (where Ln is a lanthanide, A an alkali-earth metal) with doping levels $`x0.3`$  is quite a striking fact. It indicates that Fermi states in such materials remain localized even at so heavy doping and, from the before cited Mott criterion, the upper limit for localization radius should be estimated as $`r_00.4a`$. This is an opposite limit to the traditional case of shallow dopants, hence a different evolution of the excitation spectrum can be expected. In particular, the effects of electron-electron interaction can be much more pronounced.

Of course, real perovskite manganites possess many other peculiar properties, as spin-dependent kinetic energy (it is just this dependence that suppresses kinetic energy in the paramagnetic phase) and the related Zener mechanism of double exchange , strong coupling of charge carriers to Jahn-Teller deformations , possible formation of small spin polarons , charge localization  and charge ordering , etc. However all the above mechanisms are usually discussed within uniform (that is, completely ordered) models, while our main focus now is on the effects specific for extremely strong disorder. To this end, we propose a model approximation, starting from a set of strictly localized single-site electronic states, with random energy in each site formed by the superposition of classic Coulomb potentials from other (occupied) sites and from the fixed charged dopants (so that the overall electroneutrality is assured). At the next step, the kinetic energy is considered as a small perturbation, triggering the hops between the nearest neighbor occupied and empty sites (if the necessary energy difference is compensated by phonons, a.c. electric fields, etc.). This model is interesting firstly as an opposite limit to the usual situation when the disorder is treated as a small perturbation of an initially uniform system, and secondly as a new realization of strongly correlated many-body system.

Below in Sec. 2 we define the model parameters and describe the numeric processes to seek the ground state and analyze various types of excitation spectra at zero temperature, as functions of system size $`L`$ and doping concentration $`x`$. Then the size, shape and topology independent behavior of the spectra is rapidly attained with growing $`L`$ (already at $`L=10÷12`$ cell units), and the main findings are:

a non-ergodicity of the full phase space with a definite distribution of local energy minima and practically identical excitation spectra with respect to all typical local minima;

asymmetric deviations from the mean-field $`(ϵϵ_F)^2`$ behavior for the density of single-particle excitations;

vanishing of the density of pair excitations in the limit of low excitation energy;

finiteness of the total density of many-body excitations in the same limit.

The further development of the model in Sec. 3 involves finite temperatures and hopping dynamics through the usual mechanism of deformation potential for electron-phonon coupling. Then the non-ergodic structure of the system phase space at zero temperature leads to its non-Markovian kinetics at finite temperatures. A special numeric algorithm is developed, simulating this statistical process and giving the temperature behavior of electronic specific heat, relaxation times and diffusion coefficient (the latter being related to the conductivity via the Einstein relation). In particular, for the typical concentration value $`x=1/3`$ we found: $`\sigma (T)\mathrm{exp}(bT^\alpha )`$, with $`\alpha 0.87`$, different from the Efros and Shklovskii value 1/2. The conclusions and perspectives for the further studies of the present model are discussed in Sec. 4.

## 2 Description of model. Zero Temperature

The crystalline structure of lanthanum manganite LaMnO<sub>3</sub> is close to ideal perovskite with the (paramagnetic) lattice parameter $`a3.9`$Å(Fig. 1a). Substitution of trivalent La<sup>3+</sup> by divalent alcali-earth ion (say, Ca<sup>2+</sup>) brings an extra hole to the system which resides on a nearby manganese site, changing its state from Mn<sup>3+</sup> to Mn<sup>4+</sup>. For a single dopant in the lattice, an arbitrarily weak tunneling will be sufficient to produce the hole state equally shared between eight Mn sites, nearest neighbors to the dopant Ca. However, for a finite doping, there appear random energy differences between these sites, and if these differences are greater than the tunneling amplitude (the kinetic energy), the hole will mainly occupy the lowest energy site. The localization is also favoured here by the suppression of kinetic energy: firstly due to the presence of intercalating oxygens between manganese sites and secondly due to the incoherence of manganese spins at higher temperatures.

Referring to this situation we consider the model where the dopant ions occupy randomly the central sites in the simple cubic lattice with probability $`x<1`$ (Fig. 1b). Each dopant releases one charge carrier into the crystal and thus acquires a unit

charge $`e`$ of opposite sign. In neglectance of hops between lattice sites, each carrier occupy single lattice site and the total electronic energy includes only Coulomb contributions:

$$E_c=\underset{𝐧}{}c(𝐧)\left[\frac{1}{2}U(𝐧)V(𝐧)\right],$$
(1)   
where $`U(𝐧)=_{𝐧^{}𝐧}c(𝐧)u(|𝐧𝐧^{}|)`$, $`V(𝐧^{})=_{𝐧^{}𝐧}c(𝐧^{})u(|𝐧𝐧^{}\delta |)`$, $`u(r)=e^2/\kappa r`$, $`c(𝐧)`$ is an occupation number for charge carrier in the lattice site $`𝐧=a(n_1,n_2,n_3)`$ (with integer $`n_i`$), while occupation of the dopant site $`𝐧+\delta `$, with $`\delta =a(\frac{1}{2},\frac{1}{2},\frac{1}{2})`$, in the same cell is defined by $`d(𝐧)`$. In what follows we consider fixed random configuration $`d(𝐧)`$ of “frozen” dopants, then the system ground state corresponds to such adjustment of the configuration $`c(𝐧)`$ of carriers that $`E_c`$ is a minimum.

For arbitrary configuration $`c(𝐧)`$, when a carrier is taken off from the site $`𝐧`$ (if occupied) or put into this site (if empty), the full energy respectively decreases or increases by $`ϵ(𝐧)=U(𝐧)V(𝐧)`$, which can be thus associated with single-particle excitations of the Fermi liquid theory. Here both $`U(𝐧)`$ and $`V(𝐧)`$ are random, but only $`V(𝐧)`$, determined by the fixed configuration of dopants, can be considered an usual local random field of Anderson’s model  whereas $`U(𝐧)`$, determined by the variable configuration of carriers, is substantially non-local. Hence each $`ϵ(𝐧)`$ essentially depends on the positions of all other carriers, and the total $`E_c`$ strongly differs from the sum of all $`ϵ(𝐧)`$. In this situation, there is no evidence for unique energy minimum and the structure of phase space can be very complicate. However, the consideration is simplified if one takes in mind that any two configurations in this space, satisfying the same normalization condition, Eq. (2), can be connected by a sequence (a phase trajectory) of single-particle moves between nearest neighbor occupied and empty sites. In fact, the following treatment is limited just to this class of trajectories.

Besides the above indicated single-particle energies $`ϵ(𝐧)`$, the excitation spectrum includes also the so-called pair energies :

$$ϵ(𝐧,𝐧^{})=ϵ(𝐧^{})ϵ(𝐧)u^{}(|𝐧𝐧^{}|),$$

that is the energy change at moving a carrier from the occupied site $`𝐧`$ to empty site $`𝐧^{}`$. The last, “excitonic” term in this expression, accounting for the correlation between different single-site energies, is just responsible for opening of the Coulomb gap in the single-particle spectrum $`\nu _{(sp)}(ϵ(𝐧))`$. If the correlations between different pair energies $`ϵ(𝐧,𝐧^{})`$ and $`ϵ(𝐧^{},𝐧^{\prime \prime })`$ are neglected, which corresponds to the mean-field approximation, a conclusion can be drawn that the pair spectrum $`\nu _p(ϵ(𝐧,𝐧^{}))`$ is ungapped . But, as will be seen from the exact numeric analysis below, in fact the density of such excitations also vanishes at $`ϵ0`$, and this situation may be supposed to exist for higher order excitations as well.

In our numeric procedure we consider lattice samples in the form of finite parallelepipeds and apply the following algorithm. The initial configurations $`d_0(\{𝐧\})`$ and $`c_0(\{𝐧\})`$ are defined by assigning them independently random values 0 or 1 for each site $`𝐧`$, so that the normalization condition holds:

$$\underset{n_1=1}{\overset{L_1}{}}\underset{n_{2=1}}{\overset{L_2}{}}\underset{n_{3=1}}{\overset{L_3}{}}d(n_1+\frac{1}{2},n_2+\frac{1}{2},n_3+\frac{1}{2})=\underset{n_1=1}{\overset{L_1+1}{}}\underset{n_{2=1}}{\overset{L_2+1}{}}\underset{n_{3=1}}{\overset{L_3+1}{}}c(n_1,n_2,n_3)=N$$
(2)   
related to the doping level $`x`$ through $`N=\left[xL_1L_2L_3\right]`$. Then we choose, from all the *nearest neighbor* pairs of occupied sites $`𝐧`$ and empty sites $`𝐧^{}`$ , the pair $`𝐧_\mathrm{𝟎}`$ and $`𝐧_\mathrm{𝟎}^{}`$ corresponding to the minimum value of pair energy: $`ϵ_{min}^{(neighb)}=ϵ(𝐧_\mathrm{𝟎},𝐧_\mathrm{𝟎}^{})`$. If $`ϵ_{min}^{(neighb)}`$ is negative, we change from the configuration $`c_0(\{𝐧\})`$ to a new configuration $`c_1(\{𝐧\})`$, moving the carrier from $`𝐧_\mathrm{𝟎}`$ to $`𝐧_\mathrm{𝟎}^{}`$. This process of single-particle moves is repeated $`m`$ times, until we come to such a configuration $`c_m(\{𝐧\})`$ that $`ϵ_{min}^{(neighb)}`$ is already positive. Then we search for the minimum $`ϵ_{min}^{(all)}`$ of $`ϵ(𝐧,𝐧^{})`$ over *all* occupied $`𝐧`$ and empty $`𝐧^{}`$ in this configuration, and, if it is negative, perform the corresponding move. This procedure is repeated until such a configuration $`c_M(\{𝐧\})`$ is reached that $`ϵ_{min}^{(all)}`$ is positive. Then $`c_M(\{𝐧\})`$ corresponds to a local equilibrium (with respect to single-particle moves) and the respective value of $`E_c=E_{min}[c_0(\{𝐧\})]`$, a functional of the initial configuration, realizes a local minimum of energy.

Next, the system is “shaken up”, that is, for the same initial dopant configuration $`d(\{𝐧\})`$, a new arbitrary initial configuration $`c_0^{}(\{𝐧\})`$ is created, restricted again by the normalization condition, Eq. (2). Then it is found that, with the same equilibration process, some new local equilibrium $`c_M^{}(\{𝐧\})`$ and a new local minimum $`E_{min}^{}[c_0^{}(\{𝐧\})]`$ are obtained. The presence of various local minima is indicative of a non-ergodic structure of the phase space (with a discrete topology restricted to single-particle moves). Since there is a one-to-one correspondence between the initial state and the final state of local equilibrium, the whole phase space of the system gets divided into a number of attraction domains, each corresponding to a definite local minimum. The absolute energy minimum, corresponding to the “true” ground state, can be defined as:

$$E_{min}^{abs}=\underset{c_0(\{𝐧\})}{\mathrm{min}}E_{min}[c_0(\{𝐧\})].$$
(3)   
The consequent “shake ups” and equilibrations define a sort of Monte-Carlo process to approach the ground state and this numeric process can be completed within a reasonable time for not too big system (this is evidenced by the uniqueness of the corresponding configuration $`c_M^{abs}(\{𝐧\})`$).

As a co-product, the process also provides the “spectrum” of local minima $`\nu _{loc}(E)=\delta (EE_{min}[c_0])_{c_0}`$, a new characteristics of the strongly interacting disordered system. At a given dopant configuration $`d(\{𝐧\})`$, the single-particle and the pair spectra, $`\nu _{sp}(ϵ)`$ and $`\nu _p(ϵ)`$, are calculated with respect to each local minimum, including the true ground state. Within accuracy to statistical noise, there is no difference found among the curves taken with respect to different minima. This shows that the “true” ground state is not any outstanding point between other equilibrium points. Finally all the results are averaged over various dopant configurations $`d(\{𝐧\})`$. For numeric simulations of the system, Eqs. (1, 2), we began from cubic samples of increasing size $`L`$, then a size independent behavior corresponding to thermodynamic limit is reached already at $`L10`$, and these averaged characteristics are presented in Figs. 2-5.

The single-particle spectra $`\nu _{sp}(ϵ)`$ at different dopings, shown in Fig. 2a,b, reveal a well-defined Coulomb gap around Fermi energy $`ϵ_F`$, while their asymptotics far from $`ϵ_F`$ is well described by the Gaussian law: $`\nu _{sp}(ϵ)\mathrm{exp}\left[\alpha (ϵϵ_F)^2\right]`$, in agreement with the known results of Lifshitz’s theory of optimal fluctuation for “tail” states in disordered systems .

A new notable feature of these spectra is a pronounced asymmetry between the densities of empty and occupied states near $`ϵ_F`$, which are described by different power laws:

$$\nu _{sp}(ϵ)\{\begin{array}{cc}(ϵϵ_F)^{2+\eta ^{}},& ϵ>ϵ_{F,}\\ (ϵϵ_F)^{2\eta ^{\prime \prime }},& ϵ<ϵ.\end{array}$$
(4)   
Here the asymmetry factors $`\eta ^{}`$ and $`\eta ^{\prime \prime }`$ are not universal, but vary with the doping and tend to zero at $`x1/2`$ (Fig. 3). The latter fact can be easily understood as a consequence of the symmetry between filled and empty sites at this concentration.

Hence the mean-field quadratic law is found to be exact only at half-filling (though the non-universal corrections to it, due to the higher order correlations, can be really small in the case of traditional doped semiconductors with $`r_0a`$).

Another qualitative difference from the mean-field behavior was found in that the density of pair excitations $`\nu _p(ϵ)`$ does not remain constant but tends to zero with $`ϵ0`$ (Fig. 4). Notably, this result is concentration independent. It is of interest for the future studies, to check also the low-energy behavior for the spectra of higher order excitations.

At least, the “spectrum” of local energy minima for the doping level $`x=1/3`$ is shown in Fig. 4. Close to the value of absolute energy minimum $`E_{min}^{abs}`$ (which is proportional to the sample volume), this distribution exhibits quadratic energy dependence (dashed line) and, farther from $`E_{min}^{abs}`$, it falls down by the Gaussian law (dotted line): $`\mathrm{exp}\left[\alpha (EE_{min}^{abs})^2\right]`$, (while the overall distribution width is independent of the sample size and relatively small).

Besides the simplest cubic shape of the samples, we also examined the slabs with $`L_1=L_2=6`$ and $`L_3=30`$. The essential features of single-particle spectrum for them (overall extension and gap asymmetry) at $`x=1/3`$ were found practically identical to those for cubic samples. At least, a modification of the above slab configuration was considered, realizing a “topologically closed” Euclidean bar (Fig. 5) where each vector $`(n_1,n_2,L+n_3)`$ is identified with $`(n_1,n_2,n_3)`$. This form is appropriate for direct simulations of current flow in a closed circuit under external electric field applied along $`𝐞`$<sub>3</sub>. However in this case a special care should be taken to assure the continuity of total energy at transitions of particles through the “topological” interface: $`L1`$, which must be equivalent to any “normal” interface $`n_3n_3+1`$. To this end, the interaction potential $`u(|𝐧𝐦|)`$ in Eq. (1) should be replaced by the modified potential:

$$\stackrel{~}{u}(𝐧,𝐦)=u(|𝐧𝐦|)+u(|𝐧+L𝐞_3𝐦|),$$
(5)   
($`n_3<m_3`$), related to the two “distances” shown in Fig. 5. Evidently, this modified interaction, Eq. (5), which connects each pair of sites in two possible ways, turns to the common Coulomb law in the limit $`L\mathrm{}`$.

Numeric simulations for a “closed” (6,6,30) slab with the modified interaction, Eq. (5), showed no significant difference in the single-particle spectrum, compared to the similar “open” slab and usual interaction $`u(r)`$. Thus, the characteristics of the ground state and excitation spectra found from our simulations are not sensitive to the sample size, shape and topology and should correspond to the true thermodynamic limit of the strongly interacting disordered system, Eq. (1).

## 3 Finite temperatures

At finite temperatures, the system dynamics is determined by the thermally activated hops between nearest neighbor sites. As usually, these hops are considered to be controlled by the electron-phonon interaction in approximation of deformation potential . If we consider only longitudinal phonons with Debye dispersion, the transition rate between the neighboring sites $`𝐧`$ and $`𝐧^{}`$ with the pair energy difference $`ϵ(𝐧,𝐧^{}=ϵ)`$, including the probabilities for both processes with phonon absorption ($`ϵ>0`$) and emission ($`ϵ<0`$), is given by a simple expression:

$$\gamma (ϵ)=\gamma _0\frac{ϵϵ_D\mathrm{sin}(ϵ/ϵ_D)}{\mathrm{exp}(ϵ/T)1}\theta (ϵ_D^2ϵ^2)$$
(6)   
where $`ϵ_D=\mathrm{}s/a`$ is the Debye energy (for sound velocity $`s`$) and $`\gamma _0`$ is some constant proportional to the small tunneling matrix element. The presence of the sine term in Eq. (6) is due to the extremely short range of hops (by one lattice constant) and this factor essentially reduces the transition amplitude at low energies, compared to the usual case of electron-phonon interaction in metals.

If the mean interval of time between two consecutive transitions, in some volume where the correlations are sensible, is much longer than the transition time $`\tau _{tr}`$ itself (this is reasonable for sufficiently small $`\gamma _0`$ and atomically fast $`\tau _{tr}`$), we can consider that only one transition occurs at a time. Also, accordingly to the analysis by Knotek and Pollak , for so strictly localized states we can neglect transitions where more than one electron participate. However, the difficulty with applying Eq. (6) directly to our system consists in the fact that the transition energy for a given pair of sites (in a given transition channel) is not *a priori* defined but depends on the overall configuration and is varied at transitions between other sites. Since those transitions occur at random moments of time, the transition rate for any given channel is itself a random function of time, correlated with all other channels. In mathematical language, this means a realization of non-Markovian branching random process . To manage this problem numerically, we used the following algorithm.

For any initial configuration $`d(\{𝐧\})`$, $`c(\{𝐧\})`$, Eq. (6) defines the transition rates $`\gamma _{𝐧,𝐧^{}}=\gamma [ϵ(𝐧,𝐧^{})]`$ for all appropriate pairs $`𝐧𝐧^{}`$, and these alternatives (the transition channels) can be considered independent. In this approach, only one of the channels can be chosen for each consecutive transition, and this choice is simulated by doing the independent statistical trials for random transition times $`\tau _{𝐧,𝐧^{}}`$, accordingly to the rates $`\gamma _{𝐧,𝐧^{}}`$, and choosing the shortest time from the trial outputs.

In each particular trial, a random number $`\xi `$, $`0<\xi <1`$, is generated, producing the associated random “phase” $`\varphi =\mathrm{ln}\xi `$. It relates to the random output transition time $`\tau `$ at a constant transition rate $`\gamma `$ through: $`\tau =\varphi /\gamma `$. If $`\xi `$ is distributed uniformly: $`P_\xi =\theta (\xi )\theta (1\xi )`$, the related distribution for the phase is: $`P_\varphi =\mathrm{exp}(\varphi )`$, and hence the distribution for transition times: $`P_\tau =\gamma \mathrm{exp}(\gamma \tau )`$. The relation between the random phase and transition time can be generalized to the case when the transition rate $`\gamma `$ is not constant but *a posteriori* definite function of current time $`\gamma =\gamma (t)`$:

$$\varphi =_0^\tau \gamma (t)𝑑t$$
(7)   
(in our specific case this function is stepwise, see Fig. 6). Since the integrand in the r.h.s. of Eq. (7) is nonnegative, there always exists a single finite solution for the transition moment (Fig. 6). This value is determined not only by the trial phase value $`\varphi `$ but also by all the intermediate times $`\tau _i`$ and transition rates $`\gamma _i=\gamma (\tau _i<\tau <\tau _{i+1})`$ expressing complicate correlations between the given transition and the preceding ones. Each time interval $`\mathrm{\Delta }\tau _i=\tau _{i+1}\tau _i`$ is determined by the result of competition between the virtual values for all the channels $`j`$ possible after the moment $`\tau _i`$:

$$\mathrm{\Delta }\tau _i=\underset{j}{\mathrm{min}}(\varphi _i^{(j)}/\gamma _i^{(j)})$$
(8)   
where $`\gamma _i^{(j)}`$ is the transition rate in $`j`$-th channel between the moments $`\tau _i`$ and $`\tau _{i+1}`$, and:

$$\varphi ^{(i)}=\varphi ^{(j)}_{\tau ^{(j)}}^{\tau _i}\gamma ^{(j)}(t)𝑑t$$
(9)   
is the “residual phase at $`\tau _i`$ for $`j`$-th channel. This channel is opened at a certain time moment $`\tau ^{(j)}`$, when it is given the initial random phase value $`\varphi ^{(j)}`$, according to $`P_\varphi `$. Then $`\varphi ^{(j)}`$ is consecutively reduced to $`\varphi _i^{(j)}`$, by Eq. (9), at each intermediate transition, until such a moment $`\tau _i`$ is reached that the virtual value $`\varphi _i^{(j)}/\gamma _i^{(j)}`$ for this channel “wins” the competition, Eq. (8). Then the residual phase value attained at $`\tau _{i+1}=\tau _i+\varphi _i^{(j)}/\gamma _i^{(j)}`$ is just zero, corresponding to an exact solution of Eq. (7) at $`\tau =\tau _{i+1}`$. After a particle has performed a transition through $`j`$-th channel, this channel (together with the whole set of channels $`\overline{j}`$, having common initial site with $`j`$) is considered closed, and a number of new channels is opened for all empty sites, neighbors to the new occupied site. Also, the change of the system configuration produced by the transition in $`j`$-th channel implies all the rates $`\gamma _i^{(j^{})}`$, $`j^{}\overline{j}`$, to be changed for some new values $`\gamma _{i+1}^{(j^{})}`$, and, after opening of new channels with corresponding initial phases and transition rates, the whole process is continued. At each $`i`$-th transition in the system, the values of the energy transfer $`ϵ_i`$ (either positive or negative) and of the time interval $`\mathrm{\Delta }\tau _i`$, past the preceding transition are recorded.

It is important that, unlike a system of independent carriers, the whole passage from opening to closure of a channel is not typically a pair process and it is related to the overall density of many-body states.

Examples of time records for the full energy as functions of full time $`\tau _i=_{i^{}=1}^i\mathrm{\Delta }\tau _i^{}`$, at different values of temperature $`T`$ and concentration $`x=1/3`$ for a cubic sample with $`L=8`$, are shown in Fig. 7. They all demonstrate an initial rapid descent to the regime of dynamical equilibrium, within a certain relaxation time $`t_r`$. This time is almost independent of temperature except at very low temperatures (below $`0.05`$, in our energy units $`e^2/\kappa a`$) when $`t_r`$ grows very rapidly (see inset to Fig. 7).

The latter can serve as an indication of the “freezing” process in the glassy system, though, of course, there is no sharply defined critical temperature in a finite size sample. In the equilibrium regime, both the mean energy $`E`$ and its dispersion $`\delta E=(E^2E^2)^{1/2}`$ are obtained as certain functions of temperature.

Since the volume of our system is kept fixed and the average *in time* of total energy $`E`$ as a function of temperature is known, the specific heat $`C_v`$ can be readily obtained by differentiation: $`C_v=E/T`$ .

The corresponding numeric result is shown in Fig. 8. It displays a linear behavior at low temperatures: $`C_vT`$, characteristic both of the Fermi-liquid systems  and of the glassy systems .

Another distinctive feature of the considered specific heat is a pronounced maximum at $`T2ϵ_D`$, corresponding to saturation of the relaxation channels above

the Debye temperature. Through the relation:

$$S(T)=_0^T\tau ^1C_v(\tau )𝑑\tau $$
(10)   
the thermodynamical entropy $`S(T)`$ is fully determined by the function $`E(T)`$. On the other hand, $`S`$ is also related to $`\delta E`$ :

$$S=ln[\nu _{mb}(E)\delta E],$$
(11)   
which permits to estimate the total density of many-body states $`\nu _{mb}(E)`$, very difficult in other approaches. From the comparison of Eqs. (10) and (11) at $`T0`$, we conclude that $`\nu _{mb}(E)`$ attains *a finite value* near $`E_{min}^{abs}`$, though, for a quantitative accuracy, one perhaps needs more precise and detailed data on $`\delta E(T)`$ and $`E(T)`$ than those in Figs. 7, 8.

Next, the temperature behavior of the hopping conductivity $`\sigma (T)`$, accordingly to the Einstein relation: $`\sigma =ne^2D/(k_BT)`$ , can be estimated from that of the diffusion coefficient $`D`$. Since all the hops have a standard length $`a`$, the diffusion coefficient $`D=a^2/(3\tau _0)`$ is fully determined by the mean lifetime $`\tau _0`$ of a localized state, and the latter is merely the inverse of the average number of hops per one particle per unit time: $`\tau _0=Nlim_i\mathrm{}(\tau _i/i)`$. Then, from the plots of full number of hops $`i`$ vs full time $`\tau _i`$ at different temperatures and the same concentration $`x=1/3`$ (they all turn perfectly linear after the same relaxation time $`t_r`$ as that for energy, Fig. 9a), we deduced the values of $`D`$, proportional to the slopes $`di/d\tau _i`$. Finally, the double logarithmic plot: $`\mathrm{ln}(\mathrm{ln}\sigma _0\mathrm{ln}\sigma )`$ vs $`\mathrm{ln}(1/T)`$, Fig. 9b) (where the fitting parameter $`\sigma _0`$ was adjusted

to get the best linearity), permits to infer the modified hopping conductivity law:

$$\sigma (T)=\sigma _0\mathrm{exp}[(T_0/T)^\alpha ],$$
(12)   
with $`\alpha 0.87`$ and reasonably low $`T_00.1`$. It is of interest to compare these figures with the mean-field values : $`\alpha =1/2`$ and $`T_02.1`$ (the latter is obtained using the estimate for localization radius $`r_00.4a`$ by the Mott criterion at $`x=1/3`$, see Introduction).

## 4 Conclusions

A model of strongly disordered lattice system with long-range Coulomb interactions between localized charge carriers has been considered. The total electronic energy is characterized by the presence of multiple metastable minima (including the true ground state), and different types of excitation spectra over these minima. A numeric procedure, accounting for all many-body correlations in finite size samples, confirms the existence of Coulomb gap in the single-particle spectrum and also provides corrections to the known mean-field theory results, as asymmetry of the gap at non-half-filling, vanishing density of low energy pair excitations, modified temperature exponent for hopping conductivity. The further analysis of this model can involve direct simulations of current flow in a “topologically closed” sample (Sec. 2) and formation of cluster states at finite tunneling amplitudes between nearest neighbour sites with sufficiently small pair energies.

This research was supported by the Portuguese program PRAXIS XXI through the project 2/2.1/FIS/302/94 and under personal Grants BPD 14226/97 (Yu.G.P.) and BM/12717/97 (J. M. V. L.).