# Coulomb gap in the one-particle density of states in three-dimensional systems with localized electrons

###### Abstract

The one-particle density of states (1P-DOS) in a system with localized electron states vanishes at the Fermi level due to the Coulomb interaction between electrons. Derivation of the Coulomb gap uses stability criteria of the ground state. The simplest criterion is based on the excitonic interaction of an electron and a hole and leads to a quadratic 1P-DOS in the three-dimensional (3D) case. In 3D, higher stability criteria, including two or more electrons, were predicted to exponentially deplete the 1P-DOS at energies close enough to the Fermi level. In this paper we show that there is a range of intermediate energies where this depletion is strongly compensated by the excitonic interaction between single-particle excitations, so that the crossover from quadratic to exponential behavior of the 1P-DOS is retarded. This is one of the reasons why such exponential depletion was never seen in computer simulations.

## I Introduction

The study of localized, interacting electrons in a disordered system has been a source of interesting physics for nearly half a century. The canonical example of such a system is a lightly doped, compensated, -type semiconductor, where electrons become localized on donor sites SE1984 . The study of electron-electron interactions in the localized regime was initiated by Pollak Pollak and Srinivasan Srinivasan . Efros and Shklovskii ES1975 ; Efros later argued that the single particle density of states (1P-DOS) tends to zero at the Fermi level as a result of the long-range part of the Coulomb interaction between electrons. They proposed the following universal soft gap in the 1P-DOS near the Fermi level at temperature , which is called the Coulomb gap and depends only on the dimensionality of the system:

(1) |

and

(2) |

Here the reference point for the energy is the Fermi level and denotes the square of the electron charge divided by the dielectric constant. Eqs. (1) and (2) were derived for the case when the bare DOS, which is the DOS without Coulomb interactions, has a non-zero value at .

The primary method for quantitative, theoretical study of the Coulomb gap is through computer simulations Baran1979 ; Levin ; Lee ; Mobius ; Li ; Zimanyi ; Mobius3 ; Palassini ; Mobius2 , which mostly confirm the above results for the 1P-DOS. The most important experimental manifestation of the Coulomb gap is its effect on the variable range hopping conductivity , which as a consequence of the Coulomb gap obeys the law ES1975 ,

(3) |

observed in dozens of experimental works. Here is a constant and , where is the localization length and is a numerical constant that depends on the dimensionality: SE1984 while Nguyen ; tsig2002 .

Eqs. (1) and (2) are derived from the requirement that the ground state be stable with respect to the transfer of a single electron from an occupied to an empty state. In this introduction we will just mention that this depletion of the bare 1P-DOS is a result of the excitonic attraction of electrons and holes near the Fermi level, which eliminates electrons and holes with small that are close in space, pushing them to higher energies. In other words, one can say that electrons and holes that remain in the ground state and contribute to the 1P-DOS “intimidate” each other.

While Eqs. (1) and (2) reflect the effect of this first-order stability criterion, a more challenging question is how the 1P-DOS is affected by higher-order stability criteria, in which two or more electrons are transferred simultaneously. It has been shown Efros that such criteria do not significantly modify the linear gap in the 2D case [Eq. (1)], but they should modify the quadratic gap in the 3D case [Eq. (2)]. Previous authors Efros ; Baran1980 have suggested that compact dipole excitations intimidate single particle excitations, thereby making the 1P-DOS exponentially small in the limit of small enough energies. It was also understood that this effect does not modify the law of variable-range hopping conductivity [Eq. (3)] because when dipole excitations play an important role hopping conductivity occurs by multi-electron excitations called electronic polarons, for which the DOS always obeys Eq. (2) Efros ; SE1984 ; ES1985 . The depleted 1P-DOS could in principle be seen in tunneling experiments MarkLee .

However, such strong, exponential depletion of the 1P-DOS was never seen in computer experiments Baran1979 ; Levin ; Lee ; Mobius ; Li ; Zimanyi ; Mobius3 ; Palassini ; Mobius2 . For large simulated samples this is the case because computer simulations are not able to reach the ground state of the system or cannot enforce a sufficiently large number of higher-order stability criteria. In this paper we show that even for samples where the ground state can be reached the crossover from a quadratic to an exponentially small 1P-DOS with decreasing energy is strongly retarded. The reason is that the survival probability of single particle excitations in the presence of intimidating dipoles remains a weaker function of than Eq. (2) over a wide range of energies. The crucial idea of this paper is that the bare DOS multiplied by can be considered as an effective new bare DOS. We show that because of the remarkable stability of with respect to changes in the bare DOS, such modification produces only a small effect on the 1P-DOS over a substantial range of intermediate energies.

Finally, at even smaller energies, where the exponential effect of dipoles reduces the modified bare DOS so strongly that the mutual electron-hole intimidation plays a weak role, the 1P-DOS decreases exponentially and is proportional to . At present, computer simulations are not able to find the ground state in systems that are large enough to explore such small energies.

The remainder of this paper is organized as follows. In Sec. II we present the simple model Hamiltonian that is used in most numerical studies of the Coulomb gap and we discuss the role of mutual intimidation of single particle excitations, including its description by the self-consistent equation (SCE) Efros ; ES1985 . In Sec. III we recapitulate old theoretical results related to dipole excitations and their mutual intimidation Baran1980 . We then discuss the depletion of the bare DOS by surviving dipole excitations Efros ; Baran1980 and we use this depleted DOS as a bare DOS to solve the SCE for . In Sec. IV we study analytically the stability of the Coulomb gap for hypothetical cases where the bare DOS follows a power law, and we confirm all qualitative conclusions obtained numerically in Sec. III.

## Ii Excitonic effect on single particle density of states

In order to remind the reader of the derivation of Eqs. (1) and (2), we start from the model Hamiltonian Efros :

(4) |

The electrons described by this Hamiltonian are assumed to occupy sites on a cubic lattice with lattice constant . is the occupation number of site and is the distance between sites and . The quenched random site energies are distributed uniformly within the interval , where is some positive number, so that the bare DOS created by the random site energies is the “box function” . Here, and denotes the Heaviside step function. In order to make the system neutral each site is given a positive background charge . Due to electron-hole symmetry, the chemical potential . The single particle energy at site is given by

(5) |

At zero temperature (in the ground state of the system) all states with are occupied and all states with are empty, since the ground state must be stable with respect to transfer of an electron from an occupied site to infinity or from infinity to an empty site. This condition is only the first stability criterion of the ground state. The second stability criterion considers the transfer of one electron from some site that is occupied in the ground state to a site that is vacant in the ground state. The change in resulting from such a transfer,

(6) |

must be positive for the stability of the ground state. The last term in Eq. (6) reflects the simple fact that the ground state energy includes the potential of the transferred electron, which initially was at site .

One can see the origin of the Coulomb gap from Eq. (6). Since and , the first two terms give a positive contribution to , while the third term is negative. Thus, any two sites with energy close to the Fermi level should be separated in space to keep . To see how this produces the Coulomb gap, consider sites whose energies fall in a narrow band around the Fermi level. According to Eq. (6), any two sites in this band with energies on opposite sides of the Fermi level must be separated by a distance not less than . Therefore the concentration of such sites cannot exceed . Thus, the 1P-DOS must vanish when tends to zero at least as fast as .

In this way we arrive at a simple upper bound for the 1P-DOS. Applying the principle of the micro-canonical distribution to the disordered system we find that all states except the forbidden ones are homogeneously distributed in phase space. Thus, the bound is an exact one if we enforce only the criterion . Applying all other stability criteria can only reduce the 1P-DOS.

The numerical coefficients in Eqs. (1) and (2) were later found by solving the SCE for , which is designed to describe mutual intimidation of single particle excitations. In 3D this equation has the form Efros ; ES1985

(7) |

where is the bare DOS.

The asymptotic solution of Eq. (7) at small energies leads to Eq. (2) Efros . Here is the width of the Coulomb gap, defined as the energy at which Eq. (2) matches , so that . Solution of the 2D version of the SCE at small energies leads to Eq. (1). Remarkably, these asymptotic solutions do not depend on or and in this sense are universal. It should be noted that Eq. (7) does not describe the corresponding increase in high energy states that accompanies the depletion of low energy states. In this sense its solution is not exactly normalized to contain the same total number of states as , but this should not be very important for .

At only moderately small energies, where is not asymptotically small, the integral equation of Eq. (7) has no known analytical solution. Nonetheless, one can find numerical solutions to the SCE by using a simple iteration procedure, wherein an initial guess is made for the function [say, that of Eq. (2)] and then inserted into the right-hand side of Eq. (7) to produce a new estimate for . This iteration can be repeated an arbitrary number of times until a convergent solution is found for the function over the desired range of energies convergencefootnote .

## Iii Dipole excitations and their interaction with single particles excitations

Until now we have taken into account only the conditions that be minimized with respect to the change of one or two occupation numbers. An analysis of the role of stability criteria with respect to the change of three or more occupation numbers shows that in the 2D case the single-electron approach, used above, is good, while in the 3D case the physics is more complicated Efros .

It is convenient to talk about this problem by introducing a dipole branch of excitations. For a given electron transition from a site that is occupied in the ground state to an empty site, the energy is given by Eq. (6). The result of such a transition is the formation of an “electron-hole pair”. If the energy of this transition satisfies , then the pair can be thought of as two independent single-particle excitations: an extra electron on site and an extra hole on site . If , on the other hand, the pair constitutes a strongly bound small dipole ES1975 ; Efros ; Baran1980 . The typical arm of such dipoles is . At small these dipoles are spatially separated from each other and do not participate in dc conductivity. They do contribute, however, to the ac conductivity and the low temperature thermodynamics ES1975 ; Efros ; Baran1980 .

In the first approximation the interaction between these dipole excitations was assumed to be weak ES1975 ; Efros and their density of states (2P-DOS) was predicted to be constant at small energies and given by . A more careful consideration Baran1980 ; ES1985 shows that dipole-dipole interaction leads to some mutual intimidation, resulting in the logarithmic depletion of . Specifically, . At the same time, the remaining soft dipoles have dipole arms that are shorter than , with an average arm . This happens because shorter dipoles interact more weakly and therefore survive with larger probability.

Let us now describe the interaction between dipole excitations and single particle excitations. We first note that dipole excitations are more abundant than single particle excitations due to the strong depletion of the 1P-DOS at the Fermi level. Thus, it is natural to think that after finding the 2P-DOS one can use it to study the role of dipoles in intimidating single-particle excitations. Following this logic Efros Efros noticed that in 3D a dipole with the proper orientation and proximity can intimidate a single particle excitation. Using the constant 2P-DOS , Efros showed that the probability that a single particle will survive this interaction is given by . An improved understanding of the effects of dipole-dipole interactions Baran1980 leads to a logarithmic increase of this probability. A calculation along the lines of Ref. Baran1980, shows that a more accurate expression for the survival probability is

(8) |

where is a numerical constant Baran1980 . Both the logarithmic reduction of and the logarithmic reduction of the size of dipoles discussed above contribute to the logarithmic factor in Eq. (8) that increases the probability . The subscript in emphasizes that it describes the effect of intimidation by a single dipole.

Later in Ref. Baran1980, , the authors realized that collective intimidation by many dipoles surrounding a one particle excitation is stronger than intimidation by a single dipole. This can be seen by considering that around an empty site (a positive single particle excitation) there is a certain probability to find a cloud of many dipole moments oriented away from the hole. Such dipoles can re-polarize (move their electron away from hole) when an electron is brought to the empty site from a low-energy occupied site. This process lowers the system’s total energy. To guarantee that the original state is not vulnerable to such an event, and, therefore, is the true ground state, one has to exclude such a polarization atmosphere.

The probability of a hole or electron surviving in the ground state was calculated in Ref. Baran1980, :

(9) |

Remarkably, the argument of the exponential in is exactly the square of the argument of the exponential of . Thus, in the interesting case where the latter argument is larger than unity, the intimidation of single particle excitations is determined by the collective effect result of Eq. (9).

Eq. (9) is of course derived in the limit that . Attempting to use it at leads to nonphysical divergence, since tends to zero. Instead, the correct behavior for is that it should approach unity at . One can model such behavior while still maintaining the correct small energy asymptotic form in Eq. (9) by introducing a multiplicative “crossover function” into the exponent of Eq. (9):

(10) |

The function should have the properties and . One should also ensure that approaches zero at faster than , so that the divergent behavior of Eq. (9) at is eliminated. One obvious choice for is

(11) |

with . In the numeric results presented in Fig. 1 we use , which produces a smooth crossover from the asymptotic behavior of Eq. (9) to at .

We can now consider how the 1P-DOS is affected by the elimination of single particle states, as given by the small survival probability in Eq. (10). Asymptotically at very small , goes to zero faster than and should play a dominating role in the 1P-DOS. This is why the fact that numerical data for in 3D are close to Eq. (2) looked puzzling SE1984 ; ES1985 .

We would like to argue here that the dominance of Eq. (10) requires very small . One can notice that around the probability still has the strength of the first power of . It reaches the strength of only at . This is seen in Fig. 1, where the probability is plotted by the dashed blue line.

Still, one may think that the effect of mutual intimidation of single particle excitations and the effect of intimidation of single particles by dipole excitations are multiplicative, so that to obtain one should multiply the solution of Eq. (7) for the bare “box” DOS, which we dub (thin red line), by from Eq. (10). This product is shown by the dashed-dotted green curve in Fig. 1.

We argue that instead of directly multiplying the final result of Eq. (7) by , one should think that the role of is to modify the bare DOS . Thus, to find one should replace in Eq. (7) by and solve this new self-consistent equation for . The result of this calculation is shown by the thick black line in Fig. 1, as determined by the numerical iteration procedure described in Sec. II.

One can see that at the 1P-DOS stays much closer to Eq. (7) than the product . Thus, in the range of energies , the solution of Eq. (7) is very stable with respect to modification of the bare density of states. The mechanism of this stability is as follows. When becomes moderately small, the exponential factor on the right-hand side of Eq. (7) grows sharply, supporting an almost unchanged value for . In other words, any attempt at moderate reduction of by dipole intimidation is met by resistance from the compensating effect of weakening of the intimidation by single-particle excitations. This is an extension of the universality of the small energy solution of Eq. (7) with respect to varying the bare “box” density of states SE1984 ; ES1985 .

Returning to Fig. 1, we note that at , where the probability changes substantially faster than , our result for follows the exponentially small . At present, however, this is a purely theoretical range, since computer simulations can find the ground state at only for lattices of size or smaller PalassiniFuture . For such small samples, finite size effects Baran1979 limit the range over which one can determine the 1P-DOS of infinite systems to PalassiniFuture . At these energies the value of in the ground state follows Eq. (2) PalassiniFuture , in agreement with our prediction in Fig. 1 (the thick black curve).

## Iv Stability of the Coulomb gap: Analytical solutions

In the previous section we dealt with the stability of the 1P-DOS Coulomb gap with respect to moderate depletion of the bare DOS by solving the SCE numerically. In this section we would like to illustrate this stability for an analytically solvable model. We consider this model only in its 3D version, but the two dimensional version is similar.

Assume that the bare “box” DOS used above is modified at small energies by the power law factor

(12) |

where . In order to find out how this changes the low energy solution of Eq. (7), we write it in the differential form

(13) |

At small values of and the solution of this equation is

(14) |

where is a numerical factor. Substituting Eq. (14) into Eq. (13) one gets

(15) |

At we are back to the bare “box” and Eq. (14) coincides with Eq. (2). One can see that even at the density of states is still quadratic and the only effect of the depletion of the bare DOS in Eq. (12) is to reduce the numerical coefficient . In this way at small energies the solution is remarkably stable to the strong change in . The reason for this stability is that at the number of available sites in a small energy band near exceeds the number of electrons that is permitted by their interaction. Of course, the electrons cannot arrange themselves as comfortably as in the case of constant , and this fact is reflected in the reduction of the coefficient .

On the other hand, when the number of places for electrons in a small energy band near is less than the number of electrons permitted by their interaction. Therefore, has the same energy dependence as . This can be seen formally from Eq. (7). Substitution of into the integral term with allows one to set inside the integral at . Then at small one gets , where is a constant given by

(16) |

Here we have used the fact that is zero at , so that the integral in Eq. (16) converges.

In the critical case, , one can show that the asymptotic solution at is

(17) |

where is a small correction to . We were not able to find this correction analytically, but numerical solutions to the SCE suggest that at .

The analytical solutions of this section confirm our numerical experience from the previous section. Namely, that when the bare density of states at small energies is larger than the one given by Eq. (2), the solution of Eq. (7) is very close to Eq. (2). On the other hand, when is smaller than Eq. (2) at small energies the density of states is close to .

We are grateful to M. Goethe and M. Palassini for making available to us their yet unpublished data on the 1P-DOS and 2P-DOS ; their data stimulated this work. We also thank them and A. Möbius for their comments on the draft of this work. We acknowledge the hospitality of the KITP, where this work was started. B. S. acknowledges financial support from the NSF.

## References

- (1) B. I. Shklovskii, A. L. Efros, Electronic Properties of Doped Semiconductors, Springer-Verlag, Berlin (1984).
- (2) M. Pollak, Disc. Faraday Soc. 50, 13 (1970).
- (3) G. Srinivasan, Phys. Rev. B 4, 2581 (1971).
- (4) A.L. Efros, B.I. Shklovskii, J. Phys. C 8, L49 (1975)
- (5) A. L. Efros, J. Phys. C 9, 2021 (1976)
- (6) S. D. Baranovskii, A. L. Efros, B. L. Gelmont, B. I. Shklovskii, J. Phys. C 12, 1023 (1979).
- (7) E. I. Levin, V. L. Nguyen, B. I. Shklovskii, A. L. Efros, Zh. Eksp. Teor. Fiz. 92, 1499 (1987) [Sov. Phys. JETP 65, 842 (1987).]
- (8) J. H. Davies, P. A. Lee, T. M. Rice, Phys. Rev. B 29, 4260 (1984).
- (9) A. Möbius, M. Richter, and B. Drittler, Phys. Rev. B 45, 11568 (1992).
- (10) Q. Li and P. Phillips, Phys. Rev. B 49, 10269 (1994).
- (11) B. Surer, H. G. Katzgraber, G. T. Zimanyi, B. A. Allgood, and G. Blatter, Phys. Rev. Lett. 102, 067205 (2009).
- (12) A. Möbius, M. Richter, Phys. Rev. Lett. 105, 039701 (2010).
- (13) M. Goethe and M. Palassini, Phys. Rev. Lett. 103, 045702 (2009).
- (14) A. Möbius, P. Karmann, and M. Schreiber, J. Phys: Conf. Ser. 150, 022057 (2009).
- (15) V. L. Nguyen, Sov. Physics – Semiconductors 18, 335 (1984).
- (16) D. N. Tsigankov and A. L. Efros, Phys. Rev. Lett. 88, 176602 (2002).
- (17) S. D. Baranovskii, B. I. Shklovskii, and A. L. Efros, Sov. Phys. JETP 51, 199 (1980).
- (18) A. L. Efros, B. I. Shklovskii, in A. L. Efros and M. Pollak (eds.), Electron-Electron Interaction in Disordered Systems, Elsevier, Amsterdam, (1985) p.409.
- (19) J. G. Massey and M. Lee, Phys. Rev. Lett. 75, 4266 (1995).
- (20) In fact, the iterative numerical solution to Eq. (7) converges only very slowly at , with the solution oscillating as a function of iteration number. The convergence of the solution can be greatly sped up if one introduces a relaxation procedure, such that each new solution for is averaged with the previous one using a weighted arithmetic or geometric mean. All numeric solutions presented in this paper use this relaxation procedure, and they were found to be completely independent of the nature of the relaxation process and of the choice of initial guess for .
- (21) M. Goethe, M. Palassini, to be published.