Charge ordering in the three-dimensional Coulomb glass at finite temperatures and low disorders
aa r X i v : . [ c ond - m a t . d i s - nn ] M a r Charge ordering in the three-dimensional Coulomb glass at finite temperatures andlow disorders
Preeti Bhandari
Department of Physical Sciences, Indian Institute of Science Education and Research (IISER) Mohali,Sector 81, S.A.S. Nagar, Manauli P. O. 140306, India
Vikas Malik ∗ Department of Physics and Material Science, Jaypee Institute of Information Technology, Uttar Pradesh 201309, India. (Dated: March 6, 2020)In this paper, we have studied the three dimensional Coulomb glass lattice model at half-fillingusing Monte Carlo Simulations. Annealing of the system shows a second-order transition fromparamagnetic to charge-ordered phase for zero as well as small disorders. We have also calculatedthe critical exponents and transition temperature using a finite sizing scaling approach. The MonteCarlo simulation is done using the Metropolis algorithm, which allowed us to study larger latticesizes. The transition temperature and the critical exponents at zero disorder matched the previousstudies within numerical error. We found that the transition temperature of the system decreasedas the disorder is increased. The values of critical exponents α and γ were less and value of ν more than the corresponding zero disorder values. The use of large system sizes led to the correctvariation of critical exponents with the disorder. I. INTRODUCTION
Many different systems, like compensated semiconduc-tors, granular metals, ultra-thin films, etc. are modeledby the Coulomb glass (CG). The term Coulomb glassrefers to systems having localized electrons (AndersonInsulators) interacting via unscreened Coulomb interac-tion. Transport in the CG system is via the hoppingconduction mechanism. The conductivity obeys T / lawgiven by Mott [1, 2] for temperatures where the density ofstates in uniform around the Fermi level. At low temper-atures, there is a soft gap in the density of states whichleads to T / law [3] for conductivity.The disorder manifests itself in the CG system in twoways. In the random site model, the electrons occupyrandom positions in the system. In the lattice model,the electrons occupy lattice sites, but the on-site ener-gies are randomly distributed. The state of the systemis thus characterized by three parameters temperature( T ), the strength of disorder (W), and strength of inter-actions. For the lattice model, the on-site energies aredrawn from Gaussian or box distribution. It was shownby M¨obius et al. [4] that the CG system at zero disor-der, shows a second-order transition from paramagneticto antiferromagnetic phase at finite temperature in twoand three dimensions. The critical exponents found werevery nearly equal to the short range Ising model.Simulations on three dimensional (3 d ) CG, modeledusing a Gaussian distribution [5, 6] have shown a phasediagram similar to the random field Ising model (RFIM).At zero temperature, there exists a critical disorder ( W c )below which one has a charge-ordered (CO) phase, andabove it, one finds a disordered phase. For disorder less ∗ [email protected] than the critical disorder, three dimensional CG modelundergoes a second-order phase transition from the para-magnetic to CO phase at the critical temperature T c .The finite-temperature critical exponents vary with thedisorder. So like RFIM system there exist two fixedpoints ( W = 0, T = T c ) and ( W = W c , T = 0). Boththe work on CG model [5, 6] shows that the critical ex-ponents for W < W c are similar to the ones at zero tem-perature fixed point as seen also in RFIM [7–10]. Soit was concluded that at least at small disorders, bothCG and RFIM belong to the same universality class Oneshould note that in both these works, the on-site ener-gies are modeled by a Gaussian distribution. The nu-merical picture for the Gaussian 3 d model was, to a cer-tain extent, justified using the replica theory formalism[11]. In the two-dimensional lattice CG model, authors[12] have shown that at zero temperature, there exists afirst-order transition from charge-ordered phase to disor-dered phase as the disorder is increased. At small dis-orders ( W < W c ), they found the second-order transi-tion from paramagnetic to charge-ordered phase as thetemperature was decreased [13]. Coarsening dynamicsof two dimensional CG model [14] showed remarkablesimilarity to the dynamics of the RFIM model [15, 16].In all these works, on-site energies were drawn from abox distribution. For disorders greater than critical dis-order, especially at high disorders, many experimentaland theoretical studies have claimed that the CG modelexhibits glassy behavior [21–28]. Mean-field studies [17–20] done on the three-dimensional CG lattice model athigh disorders show a finite temperature glass transitionusing replica symmetry breaking. Numerical studies onthe random site CG model have shown effects like aging,which were explained by the multi valley picture [29–31]. The existence of many minimas (also called pseudoground states) for the same disorder configuration aredue to the competition between disorder and long-rangeCoulomb interactions. In contrast, to these claims, manystudies found no evidence of glass transition for disordersabove W c [5, 32, 33].In this paper we have performed Monte Carlo simula-tion using Metropolis algorithm on the 3d lattice modelat low disorders. Recent studies have used populationannealing [6] and exchange method [5] to study the samemodel. Both these methods are well suited to probe forequilibrium glass transition ( W > W c ). At low disorders,far below the critical disorder Metropolis algorithm sam-pling is an adequate method as seen in our previous workon two-dimensional CG lattice model [13]. Our methodalso allowed us to study the model for much larger systemsizes. This, as our results show gives a more consistentvariation of the critical exponents with disorder.The paper is organized as follows. In Sec.II, we dis-cuss the Hamiltonian of the system and then in Sec.IIIour numerical simulation. The results obtained from thesimulation and its interpretation are presented in Sec.IVand finally the conclusions are provided in Sec.V. II. MODEL
The classical three-dimensional CG lattice model [34,35] can be described by the Hamiltonian H = N X i =1 φ i n i + 12 X i = j e κ | ~r i − ~r j | ( n i − K ) ( n j − K ) (1)where the number occupancy n i ∈ { , } , the random on-site energy φ i are the chosen from a uniform distribution P ( φ ), given by P ( h ) = 1 /W, − W/ ≤ h ≤ W/ , = 0 , otherwise , (2)here W is the strength of disorder. We are concentratingonly on the half-filled case which implies that the fillingfactor K = 1 / J ij = e /κ | ~r i − ~r j | , where κ is the dielectricconstant) using the pseudospin variable S i = n i − /
2. Sothe final Hamiltonian simulated by us can now be writtenas H = N X i =1 φ i S i + 12 X i = j J ij S i S j , S i = ± / . (3) III. NUMERICAL TECHNIQUE
In this paper, we have performed a Monte Carlo simu-lation using simulated annealing on a cubic lattice ( d =3) to find the critical exponents of a CG lattice model.The details of the simulation are as follows: The sys-tem sizes (L) studied in this paper are L = 4 , , , , H = N X i =1 φ i S i + X ≤ i ≤ j ≤ N S i S j erf c ( αr ij ) r ij − α √ π N X i S i + 12 πL X n =0 " n exp (cid:18) − π n L α (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) N X j =1 S j exp (cid:18) − πiL n.r j (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) + 2 π N (cid:12)(cid:12)(cid:12)(cid:12) N X j =1 S j r j (cid:12)(cid:12)(cid:12)(cid:12) (4)where erf c ( x ) = 2 √ π Z ∞ x e − t dt (5)is the complementary error function [38]. Eq.(4) alsohas a regularization parameter α , which must be chosensuch that it maximizes the numerical accuracy. The realspace terms in Eq.(4) are now short-range terms, so onecan easily use a spherical cut-off together with the peri-odic boundary conditions. All the parameters are tunedto ensure a stable convergence. In our simulation, wefound that, α = 5 /L , r c = L/ n c = α r c L/π weresufficient to get the desired results.We found that around W = 0 .
40 the CO state was nolonger the ground state for some configurations for L ≥
8. So we kept our calculations restricted to W = 0 . IV. RESULTS AND DISCUSSIONSA. Staggered magnetization
The phase diagram of the three-dimensional CG model[5, 6, 11] is quite well established now, suggesting aCharge Ordered (CO) phase below a particular value ofdisorder ( W c ). Since the ground state is antiferromag-netic at zero and small disorders, therefore, the staggeredmagnetization, M s (which is the difference between twosublattice magnetizations) is the appropriate order pa-rameter here. M s , can be defined as M s = (cid:20)(cid:28) N N X i =1 σ i (cid:29)(cid:21) , (6)where, staggered spins σ i are related to the Ising spins S i by the following relation: σ i = ( − i x + i y + i z S i . (7) | M s | L = 4L = 6L = 8L = 10L = 12L = 14 χ g C v (a) (b)(c) (d) FIG. 1: (Colour online) (a) Behavior of staggeredmagnetization ( | M s | ) as a function of temperature (T).(b) Dependence of antiferromagnetic susceptibility ( χ )on temperature. (c) Binder ratio ( g ) vs temperature.The crossing of this data for different system sizes givesus a value of the critical temperature, T c = 0 . C v ) vs temperature. All the fourgraphs (a)-(d) are for disorder strength W = 0 . L = 4 (averaged over 600 configurations), L = 6(averaged over 400 configurations), L = 8 (averagedover 300 configurations), L = 10 (averaged over 200configurations), L = 12 (averaged over 200configurations) and L = 14 (averaged over 200configurations).Here, i x , i y , i z denote the x, y, z coordinates of the site i and N = L are the total number of sites in the lattice. InEq.(6), we have used h ... i to denote the thermal averageand [ ... ] for the ensemble average. M s ≈ σ i which finally leads to non-zero M s . One can see this behavior of M s changing from0 to 1 with the decrease in temperature from Fig.(1(a))for W = 0 .
0, Fig.(2(a)) for W = 0 .
10 and in Fig.(3(a))for W = 0 . B. Binder ratio
Another way to monitor phase transition is by usingBinder ratio, g [39], which is defined as g = 12 − [ h M s i ][ h M s i ] ! , (8)where h M s i and h M s i denote the second and fourth mo-ments of the staggered magnetization respectively. In | M s | L = 4L = 6L = 8L = 10L = 12L = 14 χ g C v (a) (b)(c) (d) FIG. 2: (Colour online) (a) Behavior of staggeredmagnetization ( | M s | ) as a function of temperature (T).(b) Dependence of antiferromagnetic susceptibility ( χ )on temperature. (c) Binder ratio ( g ) vs temperature.The crossing of this data for different system sizes givesus a value of the critical temperature, T c = 0 . C v ) vs temperature. All the fourgraphs (a)-(d) are for disorder strength W = 0 .
10 at L = 4 (averaged over 600 configurations), L = 6(averaged over 400 configurations), L = 8 (averagedover 300 configurations), L = 10 (averaged over 200configurations), L = 12 (averaged over 200configurations) and L = 14 (averaged over 200configurations). | M s | L = 4L = 6L = 8L = 10L = 12L = 14 χ g C v (a) (b)(c) (d) FIG. 3: (Colour online) (a) Behavior of staggeredmagnetization ( | M s | ) as a function of temperature (T).(b) Dependence of antiferromagnetic susceptibility ( χ )on temperature. (c) Binder ratio ( g ) vs temperature.The crossing of this data for different system sizes givesus a value of the critical temperature, T c = 0 . C v ) vs temperature. All the fourgraphs (a)-(d) are for disorder strength W = 0 .
20 at L = 4 (averaged over 600 configurations), L = 6(averaged over 400 configurations), L = 8 (averagedover 300 configurations), L = 10 (averaged over 200configurations), L = 12 (averaged over 200configurations) and L = 14 (averaged over 200configurations).the charge-ordered phase (at temperature T < T c ), g ap-proaches the value 1, while it tends to zero at T > T c (paramagnetic phase). When T = T c , g ∗ = g acquires anon-trivial value, the critical Binder ratio. The behav-ior of g vs T are shown in Fig.(1(c)), Fig.(2(c)) and inFig.(3(c)) for W = 0 .
0, 0 .
10 and 0 .
20 respectively. Thecrossing of data at different system sizes gives the valueof T c = 0 . W = 0 . T c = 0 . W = 0 . T c = 0 . W = 0 . C. Susceptibility
Next we have calculated the antiferromagnetic suscep-tibility defined as follows: χ = N [ h M s i − h M s i ] . (9)This can be used as an additional input to analyse thephase transition in the system. In Fig.(1(b)), Fig.(2(b))and Fig.(3(b)) one can see that the peak in χ around T c becomes sharper as L increases. Another importantthing to note is that the distribution of χ shifts to lowertowards lower temperature and also becomes broader aswe increase the disorder from W = 0 . W = 0 . D. Specific Heat
Finally, we calculate the Specific heat, C v as a functionof temperature, using the following definition [40]: C v = 1 N k B T [ h E i − h E i ] (10)where E is the average energy per electron, and k B isthe Boltzmann’s constant. As expected the behavior of C v with variation in temperature and disorder (as shownin Fig.(1(d)), Fig.(2(d)) and Fig.(3(d))) was found to besimilar to what we saw in χ . E. Finite size scaling analysis
For any quantity A ( T, L ), which behaves like A ∞ ∼ t − x in the infinite system when approaching the criticaltemperature, the finite-size scaling function is A ( T, L ) = L x/ν ˜ A ( t L /ν ) (11)where t = T − T c . Using the above standard finite-sizescaling relation [41], we get g ( T, L ) ≈ ˜ g [ L /ν ( T − T c )] (12) χ ( T, L ) ≈ L γ/ν ˜ χ [ L /ν ( T − T c )] (13) C v ( T, L ) ≈ L α/ν ˜ C v [ L /ν ( T − T c )] (14) l og ( χ * ) y = 1.5872 x + 0.8105 log L l og ( C v* ) y = 0.3096 x + 0.1602(a) (b) FIG. 4: (Colour online) (a) A log-log plot of the peakvalue of antiferromagnetic susceptibility as a function of L defined as χ ∗ here vs the system size L . The solid lineis a least-squares straight line fit giving a slope of1 . ± . L defined as C ∗ v here vsthe system size L . The solid line is a least-squaresstraight line fit giving a slope of 0 . ± . L = 4 and 6were not including in the fitting. The exponents werecalculated for disorder strength W = 0 . l og ( χ * ) y = 1.4352 x + 0.9246 l og ( C v* ) y = 0.2019 x + 0.1926(a) (b) FIG. 5: (Colour online) (a) A log-log plot of the peakvalue of antiferromagnetic susceptibility as a function of L defined as χ ∗ here vs the system size L . The solid lineis a least-squares straight line fit giving a slope of1 . ± . L defined as C ∗ v here vsthe system size L . The solid line is a least-squaresstraight line fit giving a slope of 0 . ± . L = 4 and 6were not including in the fitting. The exponents werecalculated for disorder strength W = 0 . -1 0 1L ν (T - T c )01234567 χ / L γ / ν ν (T - T c )00.511.5 C v / L α / ν -1 0 1L ν (T - T c )0.30.40.50.60.70.80.91 g (a) (b) (c) FIG. 6: (Colour online) For W = 0 .
10 using theparameters T c = 0 . ν = 0 .
80: (a) The scalingplot of the antiferromagnetic susceptibility ( χ ) with theparameter γ/ν = 1 . C v ) with the parameter α/ν = 0 . L = 4 isnot including in any of the figures here. -1 -0.5 0 0.5 1L ν (T - T c )0246810 χ / L γ / ν ν (T - T c )00.30.60.91.21.51.8 C v / L α / ν -1 -0.5 0 0.5 1L ν (T - T c )0.40.60.81 g (a) (b) (c) FIG. 7: (Colour online) For W = 0 .
20 using theparameters T c = 0 . ν = 0 .
95: (a) The scalingplot of the antiferromagnetic susceptibility ( χ ) with theparameter γ/ν = 1 . C v ) with the parameter α/ν = 0 . L = 4 isnot including in any of the figures here.Now, the critical exponents α/ν and γ/ν can be extractedby plotting a least-squares straight-line fit of C ∗ v and χ ∗ versus L respectively in a log-log plot (as shown in Fig.(4,5). Here, C ∗ v and χ ∗ are the peak values of C v and χ atthe each system size. The critical exponent of correlationlength ν was used as an adjustable parameter to collapse the data of g (see Fig.(6(c), 7(c)). The value of T c , ν , γ/ν and α/ν at different disorders are summarized inTable.I. For comparison we have also included the criticalexponents at W = 0, which are close to one calculatedby other authors [6], claiming T c = 0 . γ/ν = 2 . ν = 0 .
76 and α/ν = 0 . χ , C v and g are shown in Fig.(6, 7).TABLE I: Critical exponents for the three-dimensionallattice Coulomb glass model. W T c γ/ν ν α/ν ± ± ± ± ± ± V. CONCLUSIONS
We have studied the critical behavior of a three-dimensional Coulomb Glass lattice model at half-fillingusing Monte Carlo simulation via simulated annealing atzero and small disorders. The summary of our results isas follows.(a) We found a finite-temperature phase transitionfrom a high-temperature paramagnetic phase to a low-temperature charge-ordered phase for zero disorder anddisorders less than W c .(b) Finite-size scaling analysis was done to determinethe transition temperature ( T c ), and the critical expo-nents of correlation length ( ν ), susceptibility ( γ/ν ) andspecific heat ( α/ν ). We found that T c , γ/ν and α/ν de-creases as the disorder ( W ) increases but ν increases withincrease in W . The correct variations in the critical ex-ponents with disorder have been achieved because of theuse of larger system sizes. Our estimates of critical ex-ponents at zero disorder are close to the one obtained byAmin et al. [6].(c) We also found that the peak value of susceptibil-ity and specific heat decreases, and the distribution getsbroader with an increase in disorder. ACKNOWLEDGEMENT
We thank late Professor Deepak Kumar for useful dis-cussions on the subject. We wish to thank NMEICTcloud service provided by BAADAL team, cloud com-puting platform, IIT Delhi for the computational facil-ity. PB gratefully acknowledges IISER Mohali for thefinancial support and for the use of High PerformanceComputing Facility. [1] N. F. Mott, J. J. Non-Cryst. Solids , 1 (1968).[2] N. F. Mott, Phil. Mag. B , 835 (1969). [3] A. L. Efros and B. I. Shklovskii, J. Phys. C: Solid StatePhys. , L49 (1975). [4] A. M¨obius, and U. K. R¨ossler, Phys. Rev. B , 174206(2009).[5] M. Goethe, and M. Palassini, Phys. Rev. Lett. ,045702 (2009).[6] A. Barzegar, J. C. Andresen, M. Schechter, and H. G.Katzgraber, Phys. Rev. B , 104418 (2019).[7] A. A. Middleton and D. S. Fisher, Phys. Rev. B ,134411 (2002).[8] N. G. Fytas and V. Martin-Mayor, Phys. Rev. Lett. ,227201 (2013).[9] N. G. Fytas, P. E. Theodorakis, I. Georgiou, and I. Le-lidis, Eur. Phys. J. B , 268 (2013).[10] B. Ahrens, J. Xiao, A. K. Hartmann, and H. G. Katz-graber, Phys. Rev. B , 174408 (2013).[11] V. Malik and D. Kumar, Phys. Rev. B , 125207 (2007).[12] P. Bhandari, V. Malik, and S. R. Ahmad, Phys. Rev. B , 184203 (2017).[13] P. Bhandari and V. Malik, Eur. Phys. J. B , 147(2019).[14] P. Bhandari, V. Malik, and S. Puri, Phys. Rev. E ,052113 (2019).[15] E. Lippiello, A. Mukherjee, S. Puri and M. Zannetti, Eu-rophys. Lett. , 46006 (2010); F. Corberi, E. Lippiello,A. Mukherjee, S. Puri and M. Zannetti, J. Stat. Mech.(2011) P03016.[16] F. Corberi, E. Lippiello, A. Mukherjee, S. Puri and M.Zannetti, Phys. Rev. E, , 144201(2007).[18] A. A. Pastor and V. Dobrosavljevi´c, Phys. Rev. Lett. ,4642 (1999).[19] S. Pankov and V. Dobrosavljevi´c, Phys. Rev. Lett. ,046402 (2005).[20] M. M¨uller and L. B. Ioffe, Phys. Rev. Lett. , 256403(2004).[21] M. Ben-Chorin, Z. Ovadyahu, and M. Pollak, Phys. Rev.B , 15025 (1993).[22] G. Martinez-Arizala, D. E. Grupp, C. Christiansen, A.Mack, N. Markovic, Y. Seguchi, and A. M. Goldman, Phys. Rev. Lett. , 1130 (1997)[23] G. Martinez-Arizala, C. Christiansen, D. E. Grupp, N.Markovic, A. Mack, and A. M. Goldman, Phys. Rev. B , R670 (1998).[24] Z. Ovadyahu and M. Pollak, Phys. Rev. Lett. , 459(1997)[25] A. Vaknin, Z. Ovadyahu, and M. Pollak, Phys. Rev. Lett. , 669 (1998).[26] A. Vaknin, Z. Ovadyahu and M. Pollak, Phys. Rev. Lett. , 3402 (2000).[27] A. Vaknin, Z. Ovadyahu and M. Pollak, Phys. Rev. B. , 134208 (2002).[28] V. Orlyanchik and Z. Ovadyahu, Phys. Rev. Lett. ,066801 (2004).[29] A. Vaknin, Z. Ovadyahu and M. Pollak, Phys. Rev. B. , 6692 (2000).[30] A. B. Kolton, D. R. Grempel and D. Dom´ınguez, Phys.Rev. B. , 024206 (2005).[31] M. Kirkengen and J. Bergli, Phys. Rev. B. , 075205(2009).[32] B. Surer, H. G. Katzgraber, G. T. Zimanyi, B. A. Allgoodand G. Blatter, Phys. Rev. Lett. , 067205 (2009).[33] A. M¨obius and M. Richter, Phys. Rev. Lett. , 039701(2010).[34] B. I. Shklovskii and A. L. Efros, Electronic Propertiesof Doped Semiconductors, Heidelberg: Springer, Heidel-berg, (1984).[35] M. Pollak, M. Ortu˜no, and A. Frydman, The ElectronGlass, Cambridge University Press, New York, (2013).[36] P. P. Ewald, Ann. Phys. , 253 (1921).[37] S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc.R. Soc. A , 27 (1980).[38] M. Abramowitz and I. A. Stegun, Handbook of Mathe-matical Functions with Formulas, Graphs, and Mathe-matical Tables (Dover, New York, 1964).[39] K. Binder, Phys. Rev. Lett. , 693 (1981).[40] M. H. Overlin, L. A. Wong, and C. C. Yu, Phys. Rev. B , 214203 (2004).[41] H. Rieger, Phys. Rev. B52