First Order Transition in Two Dimensional Coulomb Glass

We have studied the ground states of two dimensional lattice model of Coulomb Glass via Monte Carlo annealing. Our results show a possibility of existence of a critical disorder (Wc) below which the system is in the charge ordered phase and above it the system is in the disordered phase. We have used finite size scaling to calculate Wc = 0.2413, the critical exponent of magnetization {\beta} = 0 indicating discontinuity in magnetization and the critical exponent of correlation length {\nu} = 1.0. The distribution of staggered magnetization for different disorder strengths shows a three peak structure. We thus predict that two dimensional Coulomb Glass shows a first order transition at T=0.


Introduction
Since the introduction of Random field Ising Model, a lot of discussion has been done on the possibility of existence of ferromagnetic ordering below a critical disorder strength (W c ) in two dimensional (2D) RFIM. From Imry-Ma arguments [1], one finds that the lower critical dimension (d c ) below which the ferromagnetic ordered phase gets destroyed is d c ≤ 2, where d = 2 was considered as a limiting case. By doing field theoretical calculations [2,3], later it was suggested that d c = 3. In 1983, Binder [4] gave an argument that the roughening of domain walls would lead to the destruction of ferromagnetic ordering in 2d RFIM. Then an exact result was given by Bricmont and Kupiainen [5] in 1987 where they showed that there exist a ferromagnetic ordering in three dimensional (3d) RFIM. A rigorous theoretical proof was then given by Aizenman and Wehr [6] in 1989 that no long range ordering is possible in 2d RFIM case. Scaling theory for 3d RFIM assuming second order transition at T=0 leads to modified hyperscaling laws [7]. However, in contradiction to all the above arguments Frontera and Vives [8] in 1999 showed numerical signs of transition in 2d RFIM at zero temperature below a critical random field strength. In 2012, Aizenmann [9] too claims existence of phase transition in 2d RFIM. Spasojevic et al [10] have also given some numerical evidences which proves the existence of phase transition in 2d nonequilibrium RFIM at T = 0. Recently, Suman and Mandal [11] have studied some aspects of nonequilibrium 2d RFIM at low temperature. They have commented that 2d RFIM exhibits a phase transition in the disorder parameter even at T > 0. These numerical work suggests a possibility of presence of long range ordering in 2d RFIM at low disorder strengths. If we assume that the above arguments which are claiming that at T = 0 there is a finite disorder strength below which long range ordering in 2d RFIM exist, then the order of transition is another question which is not answered very clearly. This question remains unclear even for 3d RFIM where there is no controversy over the presence of phase transition at zero temperature. Some earlier work suggested a first order transition [12][13][14][15][16][17][18] but there are arguments [19][20][21][22] which favour second order transition in 3d RFIM. Since the value of magnetic component β/ν is very small, it is difficult to determine whether at transition the magnetization vanishes continuously or discontinuously. Numerical study of 3d RFIM at T = 0 claims transitions to be of continuous type [23][24][25][26] but the value of the critical exponent of magnetization (β) is very close to zero in all the papers. We are here interested in determining whether a critical behaviour similar to 2d RFIM exist in the case of 2d Coulomb Glass (CG). Coulomb Glass is defined as a system where all electron states are localized and interactions takes place via long-range Coulomb potential. Since at low temperature, localized electrons are unable to screen the Coulomb interactions effectively, it will be interesting to see how differently these long range interactions effect the phase transition from the RFIM. In our previous work [27] we found the transition from COP to disordered phase was driven by rearrangement of domain wall of the metastable state in charge ordered phase (COP) as W was increased to give disordered phase which indicates phase coexistence around transition. In this paper we have calculated the critical exponents using finite size scaling analysis. We have also investigated the distribution of staggered magnetization which shows a three peak structure characterstic of a first order transition.

Model
We have considered lattice model of CG in 2d as defined by Davies et al [28]. The model is a regular lattice of sites where on each site electron states are localized. The Hamiltonian of this system, defined over a spin configuration {S i = ± 0.5} is given as where the random on-site energies φ i are chosen independently from a probability distribution P(φ i ) which is considered as The unscreened Coulomb interactions is defined as J ij = e 2 /κR ij where κ is the dielectric constant and R ij is the distance between sites i and j which is calculated using periodic boundary conditions. Since the system under consideration possesses particle-hole symmetry hence the chemical potential μ = 0. We have measured energies in units of e 2 /κa where a ≡ 1 is the lattice constant. We have only considered the case of half filling where the number of electrons are half of the total number of sites (N) in the system. We assume that there exist a critical disorder W c below which phase is charge ordered and above which it is disordered [29]. Now for a CG system a COP implies Antiferromagnetic ordering as J ij > 0. So staggered magnetization is the order parameter which is defined as where σ i = (-1) i S i and <...> is the ensemble average performed over different disorder configurations.

Numerical simulation
We have first used Monte Carlo (MC) annealing technique to reach the minimum possible energy state. The initial configuration used for the simulation was completely random with half sites assigned S i = 0.5 and the second half assigned with S i = -0.5 . The annealing was done using the Metropolis algorithm [30]. This algorithm constitutes of a random walk in the entire space of all possible configurations of the system. Since the number of electrons are conserved, one uses Kawasaki Dynamics. A single MC step is performed by randomly choosing two sites i and j, where one site is unoccupied (S i = 0.5) and the second site is occupied (S i = -0.5) for spin-exchange. This exchange is not always executed, which is possible in a simple MC procedure. If the exchange of spins corresponds to relaxation i.e. if the final state of the system yields an energy level which is lower than the initial state, then the exchange is done with a probability 1  ij P (4) On the other hand if the exchange corresponds to a thermal excitation into a state whose energy is higher than the initial state, let us assume by an amount ΔE ij then the spin-exchange is done with a probability given as where ΔE ij = ε j -ε i -1/R ij is calculated using the single particle Hartree energy (ε i ) [31] which is expressed as The Hartree energy of the initial random system is calculated and stored. If the chosen site is the original position of the electron then the Hartree energy remains unchanged but if the electron hops to a new site, then we update all the Hartree energies of the system. When the system is in the low temperature state then the electron chooses its initial state site with a very high probability, so we do not have to recompute our Hartree energies again which speeds up our simulation considerably at low temperature regime. Annealing was done from β = 1 to β = 100 (where β = 1/T) for different system sizes (L=16,32,48). At low temperature we increased the MC steps and our longest run was 5x10 5 for each β above 20. The investigation were carried out at different disorders starting with W=0.0 to W=0.50 for all L. Once annealing was over we then identified the clusters which were created in the minimum energy state using the Hoshen-Kopelman algorithm [32]. The cluster with nearest-neighbour sites that have same σ i was defined as domain. After identifying the domains we then investigated each domain's interactions with all the other domains present in a single disorder realization. We did not find any significant domain-domain interactions present. The domains were then flipped one by one if it decreased the energy of the system. The final state in which no more flipping was possible was then considered as the ground state. By this method we do not claim that we obtained the exact ground state.

Results and discussion
After obtaining ground state, we investigated the average size of the largest and second largest domain at each disorder for the largest system size that we have considered (L=48). Figure 1 shows the volume of the largest and the second largest domains in the ground state. For W<0.20 all the domains got flipped and we found single domain picture which indicates a COP. But on further increasing the disorder strength the size of the second largest domains starts increasing indicative of occurrence of disordered phase. The picture is similar to 2d and 3d RFIM where the ordered state breaks into two large domains [33]. For each disorder strength (W) we have then recorded the average staggered magnetization per spin. The staggered magnetization shows sharp decrease with disorder as shown in figure 2. Calculation was done for different system sizes, L=16, 32, 48 and disorder averaging was done over 800, 400 and 350 realizations of random field respectively. Here the data for |M s | are not rescaled by the factor L -β/ν , which implies that β = 0.
To find the exact value of critical disorder (W c ) and to determine the order of transition we have used standard finite size scaling function for M s which is defined as where β is the exponent of staggered magnetization and ν is defined as correlation length exponent. The scaled data is shown in figure 3, where W c = 0.2413 ; β = 0 ; and 1/ν = 1.0. The value of β = 0 could be indicative of a first order transition but as mentioned in the introduction is not a conclusive proof for a first order transition. In our previous paper [27] we studied the system near the transition region where we found that the domain wall of the metastable state in the charge ordered phase (COP) shifts to give disordered ground state. This indicates that the transition is first order type as the free energy which is equal to the energy at zero temperature has three minimas which are centered around M s = ±0.5, 0. This also shows that for each disorder configuration the staggered magnetization changes discontinuously at W c . Such discontinous jumps for bond energy were also reported for 3d RFIM [18,22,34]. The domains of the metastable state in the COP and the domains in ground state of disordered phase are plotted in figure 4 for L=48.  We have then investigated the distribution of statggered magnetization at different W for L=48. As shown in figure 5, for W<W c , peaks are centered around M s = ± 0.5 and for W>W c around M s = 0 as expected. In the transition region 0.265 < W < 0.285, one gets a three peak structure. The transition region is same as one gets in the staggered magnetization vs W graph in figure 2 for L=48. This three peak picture around the transition region is a further indication of coexisting phases which is a characteristic feature of first order transition. A similar picture for 4d bimodal RFIM was reported by Swift et al [35] but no coexisting phases were found for 3d bimodal and Gaussian RFIM.

Conclusions
To summarize we have investigated the ground state of 2d lattice model of Coulomb Glass. Using finite size scaling analysis, we have shown existence of critical disorder below which the phase is ordered (antiferromagnetic ordering) and above W c = 0.2413 it is disordered. The critical exponent of magnetization, β = 0 indicates that the transition is first order. To verify this point further we have also shown the presence of coexisting phases in two ways. Firstly, the domain in the metastable state of the COP shifts to give ground state in the disordered phase and secondly the three peak picture in the distribution of M s around the transition region characteristic feature of first order transition.