Phase diagram, correlation gap, and critical properties of the Coulomb glass
We investigate the lattice Coulomb glass model in three dimensions via Monte Carlo simulations. No evidence for an equilibrium glass phase is found down to very low temperatures, although the correlation length increases rapidly near . A charge-ordered phase (COP) exists at low disorder. The transition to this phase is consistent with the Random Field Ising universality class, which shows that the interaction is effectively screened at moderate temperature. For large disorder, the single-particle density of states near the Coulomb gap satisfies the scaling relation with in agreement with the prediction of Efros and Shklovskii. For decreasing disorder, a crossover to a larger effective exponent occurs due to the proximity of the COP.
In disordered insulators, the localized electrons cannot screen effectively the Coulomb interaction at low temperature. Therefore, many-electron correlations are important in this regime. The long-range repulsion induces a soft “Coulomb gap” in the single-particle density of states (DOS). Efros and Shklovskii (ES) ES () argued that the gap has a universal form near the chemical potential , with in dimensions, and that a saturated bound modifies the variable-range hopping resistivity from Mott’s law to . Both the existence of the gap and the crossover to at low temperature have been confirmed experimentally and in numerical simulations crossover (), but the validity of for has yet to be firmly established. Pseudo ground-state numerical calculations gave li_phillips (), Möbius et al. (1992); sarvestani (); overlin (), surer (), while finite- simulations obtain between and li_phillips (); sarvestani (); overlin () from the filling of the gap as levin (); vjs (); hunt ().
It was also suggested long ago davies () that disordered insulators enter a glass state at low temperature. Ample experimental and numerical evidence of glassy nonequilibrium effects in these systems has been obtained since see (). However, it remains unclear whether these effects are purely dynamical or reflect an underlying transition to an equilibrium glass phase (GP), and whether there is a link between glassiness and the Coulomb gap. Some evidence for a sharp equilibrium transition to a GP was found in simulations of localized charges with random positions grannan (); overlin (); diaz () but not in the presence of on-site disorder diaz (); surer (). In the latter case, the transition would not break any symmetry of the Hamiltonian, similarly to the long-debated Almeida-Thouless transition in spin glasses AT (). These issues have been brought again to the fore by recent mean-field studies pastor (); Pankov and Dobrosavljević (2005); Müller and Ioffe (2004); mueller () which predict a “replica symmetry broken” equilibrium GP below a critical temperature in the presence of on-site disorder. In this GP correlations remain critical, which leads to , and both and the gap width scale as for and large disorder strength mueller ().
In this Letter, we investigate these predictions via extensive Monte Carlo (MC) simulations of the Coulomb glass lattice model with on-site disorder efros (). In addition, we study in detail the transition from the fluid to the charge-ordered phase (COP). For , there is good numerical evidence for an Ising-like transition moebius (). For , mean-field theory predicts a stable COP for Pankov and Dobrosavljević (2005); malik (). Beyond mean field there is some numerical evidence that the COP survives small positional disorder overlin (); surer () and on-site disorder moebius_talk (), but neither the phase diagram nor the critical properties have been investigated. Our results are as follows: (i) A COP exists below the (approximate) phase boundary in Fig. 1. (ii) The fluid-COP transition is consistent with the Random Field Ising model (RFIM) universality class, which shows that the interaction is effectively screened near the transition. (iii) No GP is found for well below the mean-field Pankov and Dobrosavljević (2005), in agreement with the results of Ref. surer () but at lower and in a wider range for . (iv) Due to the long-range interaction, the glass correlation length increases rapidly and possibly diverges as . (v) For large , the DOS scales as near the gap, with a saturated exponent . (vi) As decreases, scaling breaks down above the COP, and an effective power law holds with , in contrast with Ref.surer (). A more extended account will appear later gp ().
Model and simulation – We study the Hamiltonian
where are the occupation numbers for the sites of a hypercubic lattice () with , and is the distance from to . The filling factor is (which gives ). The random on-site energies are independen and Gaussian-distributed with zero mean and variance unity. Energies and temperatures will be in units of and lengths in units of the lattice spacing .
We carry out canonical MC sampling along the paths ABC and ADE in Fig. 1 and at constant . We consider an infinite sphere of periodic images of a central cell and sum over all interactions with the Ewald method with a dipole surface term leeuw (). To reach low temperatures, we use the exchange MC algorithm exchange (). For each realization (sample) , we simulate identical replicas with different () along the simulation path. Every Metropolis steps for single-electron hops, replicas at adjacent , swap their configurations with probability , where , , and , which preserves detailed balance. The simulation time is chosen so that averages over the intervals and agree within the statistical errors, and that the identity , valid for Gaussian disorder, is satisfied. Here, and are the thermal and sample averages and are two independently simulated replicas with the same () method ().
Charge ordering – Fig.2 (top inset) shows the COP order parameter along the paths AB, BC, and DE, where and (we introduce the Ising variables ). The sharp increase demonstrates a transition to a COP. To determine the transition temperature , we measure the finite-size correlation length (CL) cooper ()
where and . Along BC, the data for for different cross (Fig.2, main panel), which signals ballesteros2 () a transition at . We observe similar crossings along AB and DE (not shown) at , in excellent agreement with Refs.overlin (); moebius (), and . The curve interpolates these three points and gives the approximate fluid-COP phase boundary in Fig. 1.
Critical behavior – Since at the fluid-COP transition has a positive specific-heat exponent moebius (), disorder is relevant and the transition will be governed by a random fixed point which, by analogy with the RFIM moebius (), we expect to be at fisher (). Assuming that the transition is second order (indeed the distribution of is unimodal at all for a predominant, and increasing with , fraction of the samples gp ()) we obtain the critical exponents in Table I. and were estimated with the quotient method ballesteros () for the observables and respectively [the quotient estimates from and agree within the errors], while was obtained by fitting to the height of the peak of the susceptibility (data not shown). The peak height for the specific heat increases slowly with (Fig. 2, bottom inset), which suggests either or a logarithmic divergence (). We could not directly estimate in a reliable way, but we obtain from the modified hyperscaling relation fisher () , assuming and using . As shown in Table I, the exponents agree fairly well with the known values for the RFIM middleton (), which suggests that the interaction is effectively short-range near the phase boundary.
Glass phase – Several works have searched for a GP by measuring the parameter davies () or higher cumulants of the overlap between two replicas diaz (). We measure instead the glass CL obtained from Eq. (2) by replacing with the “spin-glass” correlation function . In the fluid phase we have for , where is the bulk CL, thus for and for . In a “many-state” GP mueller (), tends to a constant for large , thus we have . Hence the existence of a GP will be signaled by the crossing of for different near ballesteros2 ().
As shown in Fig.3(a) for , we observe no crossing down to the lowest equilibrated temperature and well below the mean-field glass transition ( for Pankov and Dobrosavljević (2005)). Similar results were found in Ref.surer () for and . We also exclude that a GP occurs at along BC and DE by comparing the crossing temperatures for (not shown) and : For all pairs they differ by less than . Together with the results at constant , this indicates that no GP exists above the dashed line in Fig. 1.
Fig. 3(b) shows that is nearly independent of at large and apart from small finite-size effects, thus is a good estimate of for . At lower we observe , which shows that becomes larger than , and possibly diverges as . Indeed, can be fitted for by a power law with [Fig. 3(b)], which also gives a satisfactory finite-size scaling [Fig. 3(c)]. A divergence would be a nontrivial prediction, since for the ground state of a periodic sample is unique. Because and are rather small, however, we cannot rule out neither that stays finite at , nor an exponent . An interesting question is whether there is a link with the divergence of the screening length found in mean-field theory mueller () and from simple arguments hunt (); tf (). We tested that the CL obtained by replacing with in Eq. (2), remains smaller than unity at all , which suggests that the correlated regions are disordered. Finally, we simulated the short-range RFIM choosing so that the value of is close to the Coulomb glass value at , and found as , which suggests that the large CL is due to the long-range interaction.
Coulomb gap – The single-particle DOS is defined as where is the cost of adding an electron at site of the central cell while leaving the periodic images unchanged. We compute the infinite sum with the Ewald method. Because of the dipole term leeuw (), the DOS has no hard gap dipole (); gp (), unlike for a finite, nonperiodic system. The finite-size effects due to the energy scale turn out to be significant for with , while those due to the sample fluctuations of (of order ) were drastically reduced by shifting the DOS before averaging over the samples shift (). In the gap region , which is our only focus here, one expects the scaling for , with constant as and as . Fig. 4(a,b,c) show scaling plots with for and . For scaling is excellent even for this moderate size, with small deviations for due to finite-size effects. For the data are well fitted by with , which is close to the self-consistent prediction efros () (while Ref.mueller () finds ). As shown in Fig. 4(b) (inset), the finite-size scaling ansatz , which should hold for [with for large ], is also well satisfied with , . Our final estimate is , which provides strong support for a saturated ES bound.
For decreasing , we observe increasingly stronger deviations from the scaling. A fit at low gives an effective exponent for and for [Fig. 4(c)]. We interpret this as a crossover due to the vicinity of the fluid-COP boundary, below which the DOS has a hard gap at . The crossover is apparent in Fig. 4(d): Since for , the plateau for supports , while for decreasing the exponent increases to . The dependence is consistent with the scaling (not shown) with extracted from Fig.4(d) for each value of . Our results differ markedly from Ref.surer (), which reports for the same model at , [however, is much larger than our data, and increases with ]. A similar crossover in the DOS was reported for , where the COP occurs at pikus ().
In conclusion, we presented evidence that no equilibrium glass phase exists in the Coulomb glass, but a saturated ES bound holds. The long-range part of the interaction appears to be irrelevant as to the equilibrium thermodynamics, except for a possible diverging correlation length at , which calls for further investigation.
We thank M. Müller, V. Dobrosavljevic, A.L. Efros, H. Katzgraber, A. Möbius, and G. Zimanyi for discussions. This work is supported by the Generalitat de Catalunya and the Ministerio de Ciencia e Innovación (FIS-2006-13321-C02-01, AP2007-01005). The computations were performed on the BSC-RES node at Universidad de Cantabria and the Albeniz cluster at UB. MP thanks the Aspen Center for Physics for hospitality.
- (1) A. L. Efros and B. I. Shklovskii, J. Phys. C 8, L49 (1975).
- (2) See e.g. A. Möbius, J. Phys. C 18, 4639 (1985); A.G. Zabrodskii, Phil. Mag. B 81, 1131 (2001).
- (3) Q. Li and P. Phillips, Phys. Rev. B 49, 10269 (1994).
- Möbius et al. (1992) A. Möbius, M. Richter, and B. Drittler, Phys. Rev. B 45, 11568 (1992).
- (5) M. Sarvestani, M. Schreiber, and T. Vojta, Phys. Rev. B 52, R3820 (1995).
- (6) M.H. Overlin, L.A. Wong, and C. C. Yu, Phys. Rev. B 70, 214203 (2004).
- (7) B. Surer, H. G. Katzgraber, G. T. Zimanyi, B. A. Allgood, and G. Blatter, Phys. Rev. Lett. 102, 067205(2009).
- (8) E. I. Levin, V. L. Nguyen, B. I. Shklovskii, and A. L. Efros, Sov. Phys. JETP 66, 842 (1987).
- (9) T. Vojta, W. John, and M. Schreiber, J. Phys. Condens. Matter 5, 4989 (1993).
- (10) A. Hunt, Philos. Mag. Lett. 62, 371 (1990).
- (11) J. H. Davies, P. A. Lee, and T. M. Rice, Phys. Rev. B 29, 4260 (1984).
- (12) See e.g. Ref.mueller () and references therein.
- (13) M. Müller and S. Pankov, Phys. Rev. B 75, 144201 (2007).
- (14) E. R. Grannan and C. C. Yu, Phys. Rev. Lett. 71, 3335 (1993).
- (15) A. Díaz-Sánchez, M. Ortuño, A. Pérez-Garrido, and E. Cuevas, Phys. Stat. Sol. (b) 218, 11 (2000).
- (16) J.R.L. de Almeida and D.J. Thouless, J. Phys. A 11, 983 (1978).
- (17) A. A. Pastor and V. Dobrosavljević, Phys. Rev. Lett. 83, 4642 (1999).
- Pankov and Dobrosavljević (2005) S. Pankov and V. Dobrosavljević, Phys. Rev. Lett. 94, 046402 (2005).
- Müller and Ioffe (2004) M. Müller and L. B. Ioffe, Phys. Rev. Lett. 93, 256403 (2004).
- (20) A. L. Efros, J. Phys. C 9, 2021 (1976).
- (21) A. Möbius and U. K. Rößler, arxiv:0904.3723 (2009).
- (22) V. Malik and D. Kumar, Phys. Rev. B 76, 125207 (2007).
- (23) A. Möbius, talk given at TIDS11 (2005).
- (24) M. Goethe and M. Palassini, in preparation.
- (25) S. W. de Leeuw, J. W. Perram, E. R. Smith, Proc. R. Soc. London, Ser. A 373, 27 (1980).
- (26) K. Hukushima and K. Nemoto, J. Phys. Soc. Japan 65, 1604 (1996).
- (27) The largest size simulated is and the number of pairs is , for (ABC, ADE, constant), respectively. The number of samples is between 100 and 676 for , depending on the path, and larger for . For each triplet , we simulate four independent replicas to obtain an unbiased estimate of . We used MC sweeps per replica for paths ABC and , and for ADE, for a total of about computing hours.
- (28) F. Cooper, B. Freedman, and D. Preston, Nucl. Phys. B 210, 210 (1982).
- (29) H. G. Ballesteros et al., Phys. Rev. B 62, 14237 (2000); M. Palassini and S. Caracciolo, Phys. Rev. Lett. 82, 5128 (1999).
- (30) See e.g. D. S. Fisher, Phys. Rev. Lett. 56, 416 (1986).
- (31) H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, and A. Muñoz Sudupe, Phys. Lett. B 378, 207 (1996).
- (32) A. A. Middleton and D. S. Fisher, Phys. Rev. B 65, 134411 (2002).
- (33) A naive application of the Thomas-Fermi theory gives a screening length diverging as .
- (34) The energy change for an electron hop from to in all images receives a positive contribution .
- (35) For each sample we shift by , where .
- (36) F.G. Pikus and A.L. Efros, Phys. Rev. Lett. 73, 3014 (1994).