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 T=0. 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 g(e,T)=T^\delta f(|e|/T) with \delta = 2.01 +/- 0.05 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.

It was also suggested long ago [11] 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 [12]. 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 [6,14,15] but not in the presence of on-site disorder [7,15]. 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 [16]. These issues have been brought again to the fore by recent mean-field studies [13,17,18,19] which predict a "replica symmetry broken" equilibrium GP below a critical temperature T g in the presence of on-site disorder. In this GP correlations remain critical, which leads to δ = d − 1, and both T g and the gap width ∆ scale as W − 1 2 for d = 3 and large disorder strength W [13]. In this Letter, we investigate these predictions via extensive Monte Carlo (MC) simulations of the Coulomb glass lattice model with on-site disorder [20]. In addition, we study in detail the transition from the fluid to the charge-ordered phase (COP). For W = 0, there is good numerical evidence for an Ising-like transition [21]. For W = 0, mean-field theory predicts a stable COP for d = 3 [18,22]. Beyond mean field there is some numerical evidence that the COP survives small positional disorder [6,7] and on-site disorder [23], 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 T well below the mean-field T g [18], in agreement with the results of Ref. [7] but at lower T and in a wider range for W . (iv) Due to the long-range interaction, the glass correlation length increases rapidly and possibly diverges as T → 0. (v) For large W , the DOS scales as g(ǫ, T ) = T δ f (|ǫ|/T ) near the gap, with a saturated exponent δ ≃ 2.00. (vi) As W decreases, scaling breaks down above the COP, and an effective power law g L (ǫ, T = 0) ∝ |ǫ| δ holds with δ 3, in contrast with Ref. [7]. A more extended account will appear later [24].
Model and simulation -We study the Hamiltonian where n i ∈ {0, 1} are the occupation numbers for the N = L d sites of a hypercubic lattice (d = 3) with N i=1 n i = KN , and r ij is the distance from i to j. The filling factor is K = 1/2 (which gives µ = 0). The random on-site energies ϕ i are independen and Gaussiandistributed with zero mean and variance unity. Energies and temperatures will be in units of e 2 /(κℓ) 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 W = 0.2, 0.5, 1, 2, 4. We consider an infinite sphere of periodic images of a central L 3 cell and sum over all interactions with the Ewald method with a dipole surface term [25]. To reach low temperatures, we use the exchange MC algorithm [26]. For each realization (sample) ϕ = {ϕ i } N i=1 , we simulate identical replicas with different (T, W ) along the simulation path. Every N/2 Metropolis steps for singleelectron hops, replicas at adjacent (T, W ), (T ′ , W ′ ) swap their configurations with probability min (1, p) The simulation time t s is chosen so that averages over the intervals [t s /3, t s ] and [t s /9, t s /3] agree within the statistical errors, and that the identity 2T ] av − 1), valid for Gaussian disorder, is satisfied. Here, · and [·] av are the thermal and sample averages and a, b are two independently simulated replicas with the same (ϕ, T, W ) [27].
Charge ordering - Fig.2 (top inset) shows the COP order parameter M s = [ |m s | ] av along the paths AB, BC, and DE, where m s = N −1 N i=1 σ i and σ i = S i (−1) xi+yi+zi (we introduce the Ising variables S i = 2n i − 1). The sharp increase demonstrates a transition to a COP. To determine the transition temperature T c , we measure the finite-size correlation length (CL) [28] ξ L = 1 2 sin(|k min |/2) where av e ik·rij and k min = (2π/L, 0, 0). Along BC, the data for ξ L (T )/L for different L cross (Fig.2, main panel), which signals [29] a transition at T BC c = 0.0950 (15). We observe similar crossings along AB and DE (not shown) at T AB c = 0.1280 (15), in excellent agreement with Refs. [6,21], and T DE c = 0.031 (2). The curve (T c (W )/T AB c ) 1.60 = 1−(W/0.15) 1.60 interpolates these three points and gives the approximate fluid-COP phase boundary in Fig. 1.
Critical behavior -Since at W = 0 the fluid-COP transition has a positive specific-heat exponent [21], disorder  is relevant and the W = 0 transition will be governed by a random fixed point which, by analogy with the RFIM [21], we expect to be at T = 0 [30]. Assuming that the W = 0 transition is second order (indeed the distribution of m s is unimodal at all T for a predominant, and increasing with L, fraction of the samples [24]) we obtain the critical exponents in Table I. β/ν andγ/ν were estimated with the quotient method [31] for the observables M s andχ L = N [ m s 2 ] av respectively [the quotient estimates from (L, L ′ ) = (6,8), (6,10) and (8,10) agree within the errors], while γ/ν was obtained by fitting aL γ/ν to the height of the peak of the susceptibility N [ m s 2 − |m s | 2 ] av (data not shown). The peak height for the specific heat c L = 1/(N T 2 )[ H 2 − H 2 ] av increases slowly with L (Fig. 2, bottom inset), which suggests either α < 0 or a logarithmic divergence (α = 0). We could not directly estimate ν in a reliable way, but we obtain ν = 1.11(12) from the modified hyperscaling relation [30] (d − θ)ν = 2 − α, assuming α = 0 and using θ =γ/ν − γ/ν = 1.20 (20). As shown in Table I, the exponents agree fairly well with the known values for the RFIM [32], 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 [( n i − 1/2) 2 ] av [11] or higher cumulants of the overlap between two replicas [15]. We measure instead the glass CL ξ G L obtained from Eq. (2) by replacing [ σ i σ j ] av with the "spin-glass" correlation function G(r ij ) = [( S i S j − S i S j ) 2 ] av . In the fluid phase we have G(r) ∼ exp(−r/ξ G ) for ξ G < r ≪ L, where ξ G is the bulk CL, thus ξ G L ∼ ξ G for L ≫ ξ G and ξ G L ∼ L for L ≪ ξ G . In a "many-state" GP [13], G(r) tends to a constant for large r, thus we have ξ G L ∼ L d/2+1 . Hence the existence of a GP will be signaled by the crossing of ξ G L (T )/L for different L near T = T g [29]. As shown in Fig.3(a) for W = 0.5, 1, we observe no crossing down to the lowest equilibrated temperature and well below the mean-field glass transition (T g ≈ 0.037 for W = 0.5 [18]). Similar results were found in Ref. [7] for T ≥ 0.03 and W ≤ 0.4. We also exclude that a GP occurs at T g > T c along BC and DE by comparing the crossing temperatures for ξ G L /L (not shown) and ξ L /L: For all pairs (L, L ′ ) they differ by less than 1%. Together with the results at constant W , this indicates that no GP exists above the dashed line in Fig. 1. Fig. 3(b) shows that ξ G L is nearly independent of L at large T and W = 1 apart from small finite-size effects, thus ξ G L=10 is a good estimate of ξ G for T 0.017. At lower T we observe ξ G L ∼ L, which shows that ξ G becomes larger than ≃ 10, and possibly diverges as T → 0. Indeed, ξ G L=10 (T ) can be fitted for T ∈ [0.017, 0.113] by a power law T −ν ′ with ν ′ = 0.82 [ Fig. 3(b)], which also gives a satisfactory finite-size scaling ξ G L = Lf (T L 1/ν ′ ) [ Fig. 3(c)]. A divergence would be a nontrivial prediction, since for W = 0 the ground state of a periodic sample is unique. Because L and ξ G L are rather small, however, we cannot rule out neither that ξ G stays finite at T = 0, nor an exponent ν ′ = 1. An interesting question is whether there is a link with the T −1 divergence of the screening length found in mean-field theory [13] and from simple arguments [10,33]. We tested that the CL obtained by replacing [ σ i σ j ] av with [ σ i σ j − σ i σ j ] av in Eq. (2), remains smaller than unity at all T , which suggests that the correlated regions are disordered. Finally, we simulated the short-range RFIM choosing W so that the value of N −1 [ i n i ϕ i ] av is close to the Coulomb glass value at W = 1, and found ξ G L ≤ 1 as T → 0, which suggests that the large CL is due to the long-range interaction.
Coulomb gap -The single-particle DOS is defined as is the cost of adding an electron at site i of the central cell while leaving the periodic images unchanged. We compute the infinite sum with the Ewald method. Because of the dipole term [25], the DOS has no hard gap [24,34], unlike for a finite, nonperiodic system. The finite-size effects due to the energy scale L −1 turn out to be significant for |ǫ| ≤ aL −1 with a ≃ 0.3, while those due to the sample fluctuations of µ (of order W/L d/2 ) were drastically reduced by shifting the DOS before averaging over the samples [35]. In the gap region (|ǫ|, T ) ≪ ∆, which is our only focus here, one expects the scaling g L (ǫ, T ) = T δ f (|ǫ|/T ) for aL −1 ≪ (|ǫ|, T ), with f (x) ∼ constant as x → 0 and f (x) ∼ c x δ as x → ∞. Fig. 4(a,b,c) show scaling plots with δ = 2 for L = 10 and W = 4, 2, 0.5. For W = 4 scaling is excellent even for this moderate size, with small deviations for |ǫ| 0.03 due to finite-size effects. For ∆/T ≫ |ǫ|/T 6 the data are well fitted by g 10 (ǫ, T ) = c|ǫ| 2 with c ≃ 1.1, which is close to the self-consistent prediction c = 3/π [20] (while Ref. [13] finds c = 0.2083). As shown in Fig. 4(b) (inset), the finite-size scaling ansatz g L (ǫ, T ) = L −δ h(ǫL), which should hold for T ≪ |ǫ| ≪ aL −1 ≪ ∆ [with h(x) ∼ c|x| δ for large x], is also well satisfied with c = 3/π, δ = 2. Our final estimate is δ = 2.01 ± 0.05, which provides strong support for a saturated ES bound.
For decreasing W , we observe increasingly stronger deviations from the δ = 2 scaling. A fit g 10 (ǫ, T ) = c|ǫ| δ at low T gives an effective exponent δ ≃ 2.3 for W = 2 and δ ≥ 2.8 for W = 0.5 [ 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 T = 0. The crossover is apparent in Fig. 4(d): Since g L (ǫ = 0, T )/T 2 ∝ T δ−2 for aL −1 ≤ T ≤ ∆, the plateau for W = 4 supports δ = 2, while for decreasing W the exponent increases to δ > 3. The L dependence is consistent with the scaling g L (ǫ = 0, T ) = T δ h(T L) (not shown) with δ extracted from Fig.4(d) for each value of W . Our results differ markedly from Ref. [7], which reports δ = 1.83(3) for the same model at W = 0.4, T = 0 [however, g L (ǫ ≃ 0) is much larger than our data, and increases with L]. A similar crossover in the DOS was reported for d = 2, where the COP occurs at W = 0 [36].
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