Phase Ordering Kinetics of the Asymmetric Coulomb Glass Model

by

Preeti Bhandari1, Vikas Malik2 and Sanjay Puri3

1Department of Physics, Ben Gurion University of the Negev, Beer Sheva 84105,
Israel.

2Department of Physics and Material Science, Jaypee Institute of Information
Technology, Uttar Pradesh 201309, India.

3School of Physical Sciences, Jawaharlal Nehru University, New Delhi – 110067,
India.

###### Abstract

We present results for phase ordering kinetics in the Coulomb glass (CG)
model, which describes electrons on a lattice with unscreened Coulombic
repulsion. The filling factor is denoted by $K\in[0,1]$. For a square lattice
with $K=0.5$ (symmetric CG), the ground state is a checkerboard with
alternating electrons and holes. In this paper, we focus on the asymmetric CG
where $K\lesssim 0.5$, i.e., the ground state is checkerboard-like with excess
holes distributed uniformly. There is no explicit quenched disorder in our
system, though the Coulombic interaction gives rise to frustration. We find
that the evolution morphology is in the same dynamical universality class as
the ordering ferromagnet. Further, the domain growth law is slightly slower
than the Lifshitz-Cahn-Allen law, $L(t)\sim t^{1/2}$, i.e., the growth
exponent is underestimated. We speculate that this could be a signature of
logarithmic growth in the asymptotic regime.

##  1 Introduction

Many systems found in nature are far from equilibrium. We still lack a
complete understanding of non-equilibrium pattern formation in these systems.
An important problem in this context is the temporal evolution of a system
which is rendered thermodynamically unstable by a rapid quench below the
critical temperature $T_{c}$. The system does not instantaneously transform to
its new phase – rather, but small domains of coexisting phases grow with time
to reach the equilibrium state. This evolution is usually referred to as phase
ordering dynamics or domain growth or coarsening [1, 2]. We have a good
understanding of domain growth in pure and isotropic systems. The nature of
evolution towards equilibrium depends on whether or not the order parameter is
conserved, and the nature of defects (interfaces, vortices, monopoles, etc.)
which drive the coarsening process. Typically, non-conserved systems with
short-ranged interactions are characterized by the Lifshitz-Cahn-Allen (LCA)
growth law for the domain size $L$ at time $t$: $L(t)\sim t^{1/2}$ [1, 2].

Recent interest in domain growth problems has focused on experimentally
realistic systems. Most experimental systems have disorder and frustration.
Therefore, it is relevant to study the ordering kinetics of such systems. In a
system with quenched disorder, the dynamics of the system slows down due to
the trapping of domain walls by disordered sites [3, 4, 5]. Another important
system property that affects the growth is the range of the particle-particle
interaction. There have been several analytical [6] and numerical [7, 8]
studies of domain growth in systems with long-range interactions. These have
demonstrated that the growth law is modified from the LCA law if the
interactions are sufficiently long-ranged.

In this paper, we study coarsening in an important model which combines
disorder/frustration with long-range interactions. We focus on the well-known
Coulomb glass (CG) model [9, 10] at half-filling ($K=0.5$), and also at
fillings slightly less than half ($K=0.475,0.45$). In the CG model, the
electronic states are localized due to disorder. Due to the localization of
states, the electron-electron interaction is unscreened. The Hamiltonian for
the CG model is

$\mathcal{H}_{\rm
CG}\\{n_{i}\\}=\sum_{i=1}^{N}(n_{i}-K)\phi_{i}+\frac{1}{2}\sum_{i\neq
j}\frac{e^{2}}{\kappa|\vec{r_{i}}-\vec{r_{j}}|}\left(n_{i}-K\right)\left(n_{j}-K\right),$
(1)

where $N$ is the number of lattice sites, and $n_{i}$ $\in\\{0,1\\}$ is the
occupation number of the electron on site $i$ with position vector
$\vec{r}_{i}$. In Eq. (1), $\phi_{i}$ is the random energy at $i$, and $K$ is
the filling factor. The charge on each electron is equal to $-e$ and a
compensating charge of $+Ke$ is put on each site to maintain charge
neutrality. The second term in Eq. (1) denotes an unscreened Coulombic
interaction between the charges at sites $i$ and $j$ in a medium of dielectric
constant $\kappa$.

At high disorder and low temperature, experiments have shown that this system
exhibits glassy behavior such as aging [11, 12, 13, 14], slow relaxation of
excess conductance [15, 16, 17], and memory effects [18, 19]. Recent
simulations of the $d=3$ CG via population annealing have shown a transition
from a disordered to a glassy phase as the temperature is lowered at high
values of disorder [20]. No glassy phase was found in earlier exchange Monte
Carlo (MC) simulations done by Goethe and Palassini [21] and Surer et al.
[22]. The mean-field study of the model also indicates the presence of a
glassy phase [23, 24, 25].

The CG model is equivalent to an Ising anti-ferromagnet (AF) with spins
interacting via the Coulombic interaction. In $d=3$, equilibrium studies at
half-filling (symmetric CG) have shown that the CG is in the same universality
class as the random field Ising model (RFIM) with ferromagnetic nearest-
neighbor (nn) interactions [21]. The system undergoes a second-order phase
transition from a paramagnetic phase to an AF phase at zero and finite
disorders. At zero disorder, the critical exponents of the CG match those of
the $d=3$ nn Ising model [26, 27]. For disorder amplitudes less than the
critical disorder, the critical exponents are those of the $d=3$ nn RFIM.

In $d=2$, at zero disorder, it has been found that the critical exponents of
the CG match those of the $d=2$ nn Ising model [26, 28]. However, unlike the
$d=2$ nn RFIM, the CG shows a transition from a disordered to a charge-ordered
phase at $T=0$ as the disorder becomes less than a critical value [29].

To the best of our knowledge, there have been far fewer studies of dynamical
properties of the CG. In recent work [30], we investigated the ordering
kinetics of the $d=2$ symmetric ($K=0.5$) CG at zero and small disorder. In
conformity with experiments, the composition of the CG is conserved during the
evolution. However, the relevant order parameter (staggered magnetization) is
a nonconserved variable [32]. For zero disorder, the CG system obeyed the LCA
growth law [30]. Thus, in spite of frustration due to the long-range Coulomb
interaction, the CG kinetics is the same as that of the nn AF. At non-zero
disorder, the system showed a transition from power law to logarithmic growth
similar to that observed in the nn RFIM [33].

In this paper, we focus on the effect of composition asymmetry on the ordering
kinetics of the CG, i.e., we consider filling fractions away from $K=0.5$. We
point out that the Hamiltonian in Eq. (1) is invariant under the
transformation $K\rightarrow 1-K$ and $n_{i}\rightarrow 1-n_{i}$. Thus, the
case with surplus electrons ($K>0.5$) is equivalent to the case with surplus
holes with composition $1-K$. Therefore, without loss of generality, we
consider the case with $K\lesssim 0.5$. For simplicity, we consider the case
with no external disorder. The internal frustration of the system already
yields rich phenomenology, which we will discuss in this paper.

This paper is organized as follows. In Sec. 2, we give details of our
numerical simulations. In Sec. 3, we present detailed numerical results. In
Sec. 4, we conclude this paper with a summary and discussion.

##  2 Simulation Details

Recently, Rademaker et al. [34] used mean-field theory to obtain a finite
temperature phase diagram for the $d=2$ CG model over a range of average
particle density. Using their phase diagram as a reference point, we have
studied ordering kinetics at $K=0.45,0.475,0.50$. The first two compositions
differ only slightly from half-filling, and we assume that the critical
temperature ($T_{c}$) is approximately the same for all three $K$-values. From
the phase diagram proposed by Rademaker et al., one can see that the ground
state has checkerboard order close to half-filling. This implies that one can
still use staggered magnetization as an order parameter for $K=0.45\rightarrow
0.5$. For $K\simeq 0.4$, the phase diagram shows that the ground state has a
stripe order. The ground states become more exotic as $K$ is reduced further.
For this reason, we have chosen $K=0.45$ as the lower limit of asymmetry for
our simulation.

Before proceeding further, it is useful to formulate the CG Hamiltonian in
spin language. We introduce the usual spin notation: $S_{i}=2n_{i}-1$. Then,
Eq. (1) with $\phi_{i}=0$ can be written as

$\mathcal{H}_{\rm CG}\\{S_{i}\\}=\frac{1}{8}\sum_{i\neq
j}\frac{e^{2}}{\kappa|\vec{r_{i}}-\vec{r_{j}}|}(S_{i}+1-2K)(S_{j}+1-2K),\quad
S_{i}=\pm 1.$ (2)

For half-filling ($K=0.5$), the constant factor $1-2K=0$. We would like to
compare domain growth in the $d=2$ CG with a system in which the Coulombic
interaction is screened, i.e., charges (spins) interact via a short-range
potential. The latter is modeled by a nn AF with charge neutrality. The
Hamiltonian can be written as

$\mathcal{H}_{\rm AF}\\{S_{i}\\}=J\sum_{\langle
ij\rangle}(S_{i}+1-2K)(S_{j}+1-2K),\quad S_{i}=\pm 1,$ (3)

where $J>0$ is the strength of the nn interaction, and $\langle ij\rangle$
denotes nn pairs. We will subsequently use the term “CG” to refer to the model
in Eq. (2), and the term “AF” to refer to the model in Eq. (3).

The details of our simulation are as follows. We study coarsening via Monte
Carlo (MC) simulations [35, 36] of the CG and AF with Kawasaki spin-exchange
kinetics. (We will shortly discuss this kinetics in detail.) We use a square
lattice of size $N=A^{2}$, and employ periodic boundary conditions in both
directions. To cope with the long-range Coulomb interactions in Eq. (2), we
have used the Ewald summation technique [37]. The simulation of long-range
interactions is computationally much more demanding than that of nn
interactions. This restricts the size of the CG to $A=512$. The corresponding
size of the AF is $A=2048$.

The initial configuration was prepared by randomly assigning $KN$ sites with
value $S_{i}=+1$, and $(1-K)N$ sites with $S_{i}=-1$. This mimics the
disordered state at high $T$. The system was rapidly quenched at time $t=0$
from the disordered state ($T\gg T_{c}$) to a low-$T$ state ($T\sim T_{c}/2$).
As the number of electrons in the system is conserved, the stochastic kinetics
used to study the ordering process is Kawasaki spin-exchange kinetics [1]. A
pair of opposite spins $S_{n}$ and $S_{m}$ (which are nearest neighbors) are
chosen at random for spin exchange. This exchange is done with probability
$p$, given by the Metropolis algorithm [35, 36]:

$\displaystyle p$ $\displaystyle=$ $\displaystyle 1,\quad\Delta_{nm}\leqslant
0,$ (4) $\displaystyle=$ $\displaystyle
e^{-\beta\Delta_{nm}},\quad\Delta_{nm}>0,$

where $\Delta_{nm}$ denotes the change in energy due to the spin exchange, and
$\beta=1/(k_{B}T)$ (Boltzmann constant $k_{B}=1$). The energy change for the
CG can be computed from Eq. (2) as follows:

$\Delta^{\rm
CG}_{nm}=\frac{1}{4}(S_{m}-S_{n})(\varepsilon_{n}-\varepsilon_{m})-\frac{e^{2}}{\kappa
r_{nm}}.$ (5)

In Eq. (5), $\varepsilon_{n}$ denotes the single-particle Hartree energy:

$\varepsilon_{n}=\sum_{j\neq n}\frac{e^{2}}{\kappa r_{nj}}~{}(S_{j}+1-2K).$
(6)

The $\varepsilon_{i}$’s are computed at the beginning of the simulation and
updated as the spin exchanges take place. The advantage of using the
$\varepsilon_{n}$’s is that they are updated only when a spin exchange attempt
is successful. This standard method for the long-range model significantly
decreases the computational time needed for an update. A single Monte Carlo
step (MCS) corresponds to $N$ attempted updates.

The corresponding expression for $\Delta_{nm}$ for the AF is considerably
simpler, as expected. We use Eq. (3) to obtain

$\Delta^{\rm AF}_{nm}=J(S_{m}-S_{n})\left(\sum_{L_{n}\neq
m}S_{L_{n}}-\sum_{L_{m}\neq n}S_{L_{m}}\right),$ (7)

where $L_{n}$ refers to the nn of site $n$.

All the statistical data presented in this paper has been averaged over 20
independent runs for the CG. The energies and temperatures of the CG system
are calculated in units of $e^{2}/(\kappa a)$, where $a$ is the lattice
constant. The critical temperature for the CG at half-filling [26, 28] is
$T_{c}\simeq 0.1042~{}e^{2}/(\kappa a)$.

For the AF, the statistical results were obtained by averaging over 10
independent runs. In this case, the unit of energy is $J$, and $T_{c}\simeq
2.269~{}J$.

##  3 Detailed Numerical Results

In earlier work, we have established that both the CG and AF in $d=2$ at
$K=0.50$ follow the LCA growth law [30]. We now investigate the role of excess
electrons/holes or asymmetry in ordering kinetics. Without loss of generality,
we consider the case with excess holes ($K<0.5$). As discussed earlier, we
only consider the slightly asymmetric case ($K>0.45$), so that the system
still has a checkerboard-like ground state. The excess holes are distributed
uniformly throughout the system in both the CG and AF. These two situations
are shown schematically in Fig. 1. In the CG, the excess holes form a periodic
structure to minimize the Coulombic energy. No such periodicity is seen in the
AF. The presence of excess holes at the domain interfaces slows down domain
growth for both models. As mentioned earlier, the appropriate order parameter
in this problem is not the composition (which is conserved) but rather the
staggered magnetization (which is not conserved):

$\sigma_{i}=(-1)^{i_{x}+i_{y}}S_{i}.$ (8)

The conserved composition plays the role of an auxiliary slow variable. A
complete description of the kinetics is provided by Model C in the terminology
of Hohenberg-Halperin [31, 32].

In Fig. 2, we present evolution snapshots of $\sigma$-domains at $K=0.45$ for
the CG. The excess holes, which can correspond to either $\sigma_{i}=+1$ or
$\sigma_{i}=-1$ (depending on the location) are seen to be uniformly dispersed
as impurities in the domains of Fig. 2. In Fig. 3, we compare domain
morphologies at different filling factors (at $t=10^{3}$ MCS) for the CG and
AF. For both models, one sees that domain growth slows down as the number of
excess holes or domain impurities increases. (Notice that the snapshots for
$K=0.50$ do not show a perfect checkerboard structure either – this is a
consequence of thermal fluctuations in the equilibrated bulk domains.)

The phase ordering kinetics in a system is usually characterized by the equal-
time correlation function [1]:

$C(\vec{r},t)=\frac{1}{N}\sum_{\vec{R}}\left[\langle\sigma(\vec{R},t)\sigma(\vec{R}+\vec{r},t)\rangle-\langle\sigma(\vec{R},t)\rangle\langle\sigma(\vec{R}+\vec{r},t)\rangle\right],$
(9)

where we have assumed that the system is translationally invariant. This
quantity measures the overlap of staggered magnetization between two spins
separated by $\vec{r}$. Here, $\langle\cdots\rangle$ indicates an average over
independent runs, i.e., with different initial conditions and noise
realizations. If the system is isotropic and is characterized by a single
length scale $L(t)$, $C(\vec{r},t)$ obeys the following dynamical scaling form
[1]:

$C(\vec{r},t)\equiv C(r,t)=f\bigg{(}\frac{r}{L(t)}\bigg{)},$ (10)

where $r$ is the magnitude of $\vec{r}$.

The characteristic domain size $L(t)$ is defined as the distance over which
$C(r,t)$ decays to a fraction (say 0.2) of its maximum value. The scaling
property indicates that the morphology is statistically self-similar in time –
only the scale of the morphology changes. We have confirmed (not shown here)
that dynamical scaling holds for both the CG and AF. In Fig. 4(a), we plot the
scaling functions for the CG with different values of $K$. The data sets are
seen to be in excellent agreement. The solid line denotes the Ohta-Jasnow-
Kawasaki (OJK) function [1]:

$f_{\rm OJK}(r/L)=\frac{2}{\pi}\sin^{-1}\left[\exp(-r^{2}/L^{2})\right].$ (11)

The OJK function is an approximate result for the scaling function of the nn
ferromagnet with nonconserved kinetics. It is obtained by modeling the pattern
dynamics via the motion of domain boundaries or defects [1]. We see that the
OJK function describes the CG ordering kinetics very well, at least for the
small values of asymmetry considered here. In Fig. 4(b), we show the scaling
functions for the AF with small asymmetry. These are also in excellent
agreement with the OJK function.

Next, let us investigate the time-dependence of the domain growth law.
Typically, the evolution of domains in magnetic ordering is curvature-driven,
and is governed by the equation [4, 5]

$\frac{dL}{dt}=\frac{a(L,T)}{L}.$ (12)

Here, $a(L,T)$ is the diffusion constant which, in general, depends upon the
domain size and temperature. Lai et al. (LMV) [4] classified systems with
nonconserved kinetics via the functional dependence of $a(L,T)$. Consider the
situation where coarsening domains of size $L$ encounter energy barriers
$E_{B}(L)$. Then, the growth dynamics is governed by thermally-activated
hopping over these barriers, and the diffusion constant has the form
($k_{B}=1$)

$a(L,T)=a_{0}~{}\exp\left(-\frac{E_{B}(L)}{T}\right),$ (13)

where $a_{0}$ is a constant. For the usual nn Ising model with Glauber spin-
flip kinetics [1], there are no energy barriers ($E_{B}=0$) and $a(L,T)$ is a
constant. Then, Eq. (12) yields the LCA growth law:

$L=(2a_{0}t)^{1/2}.$ (14)

These are designated as Class 1 systems by LMV. For Class 2 systems, the
barrier energy is non-zero but independent of $L$: $E_{B}(L)\equiv E_{B}$. In
that case, one obtains a modified LCA law with a temperature-dependent
prefactor:

$L(t)=\left(2a_{0}e^{-\beta E_{B}}t\right)^{1/2}.$ (15)

LMV and Paul et al. [38] have also considered systems where $E_{B}(L)$ has a
power-law and logarithmic dependence on $L$, respectively. These usually occur
in systems with quenched disorder, where the boundaries of coarsening domains
are trapped by disorder sites [3, 4, 5]. In the case of power-law barriers,
the asymptotic growth law is logarithmic in time [4]. On the other hand, for
logarithmic barriers, the growth law is a power-law in time. However, the
growth exponent depends on the disorder amplitude and temperature.

The plot of $L(t)$ vs. $t$ for the CG is shown in Fig. 5(a). One sees that the
growth is slightly slower than the LCA law for all filling factors considered.
This is clear from the plot of the effective exponent $\theta_{\rm eff}$ vs.
$t$ in Fig. 5(b), where

$\theta_{\rm eff}=\frac{d(\ln L)}{d(\ln t)}.$ (16)

This plot underestimates the LCA exponent ($\theta=1/2$) for the asymmetric
CG. We interpret this as a possible signal of logarithmic growth in the
asymptotic regime. Further, the growth is slower for higher asymmetry. We will
discuss the nature of the asymptotic growth law at length shortly. The
corresponding results for the AF are shown in Figs. 5(c)-(d). These are
consistent with the results for the CG. In earlier work, Das et al. [32] have
obtained similar results for the AF, via both the kinetic spin model and its
phenomenological counterpart (Model C).

We would like to gain a better understanding of the slowing down of the growth
law in the asymmetric case. In this context, following Shore et al. [39], we
investigate the time required for shrinking of a droplet on a square lattice.
Here, the relevant order parameter is the staggered magnetization
$\sigma_{i}$, so we consider the shrinking of a square droplet
($\sigma_{i}=-1$) of size $L_{d}^{2}$ in a background of $\sigma_{i}=+1$. This
is schematically depicted in Fig. 6(a). If $K\lesssim 0.5$, there will be
uniformly scattered impurities (corresponding to excess holes) in the bulk
domains – see Fig. 3. The probability of encountering such an impurity is
$p_{h}=1-2K$. The corresponding shrinking experiment for $K\lesssim 0.5$ is
schematically depicted in Fig. 6(b).

It is convenient to rewrite the first few quadratic terms of $\mathcal{H}_{\rm
CG}$ in Eq. (2) in terms of the staggered magnetization:

$\mathcal{H}_{\rm
CG}\\{\sigma_{i}\\}=-J_{1}\sum_{nn}\sigma_{i}\sigma_{j}+J_{2}\sum_{nnn}\sigma_{i}\sigma_{j}+J_{3}\sum_{nnnn}\sigma_{i}\sigma_{j}+\cdots,\quad\sigma_{i}=\pm
1,$ (17)

where the subscript nnn refers to a sum over next-nearest-neighbor pairs, and
so on. In Eq. (17),

$J_{1}=\frac{e^{2}}{4\kappa a},\quad J_{2}=\frac{e^{2}}{4\sqrt{2}\kappa
a},\quad J_{3}=\frac{e^{2}}{8\kappa a},\quad\cdots.$ (18)

The argument below is based on the first two terms of the right-hand-side
(RHS) of Eq. (17) [39], and illustrates the key features which result in
formation of energy barriers due to excess holes.

First, consider the case $K=0.5$ in Fig. 6(a). We consider the case where
decimation is initiated at the upper left corner. The droplet shrinks by
flipping a corner spin with energy barrier $4J_{2}$. The subsequent flipping
of the edge spins requires zero energy. Thus, the time taken to decimate the
square domain is

$\tau_{d}=t_{0}(L_{d})\exp\left(\frac{4J_{2}}{T}\right).$ (19)

The decimation of droplets is the elementary process that drives the ordering
process. Therefore, we identify $E_{B}(L)$ in Eq. (13) as $4J_{2}$ for the
case $K=0.5$. This yields the LCA growth law in Eq. (15).

Next, consider the case $K\lesssim 0.5$, depicted in Fig. 6(b). Again, the
decimation is initiated at the upper left corner. The corner flip requires
energy $4J_{2}$ for the case shown in Fig. 6(b). However, in this case, we
encounter adjacent impurity spins in the flipping of edge spins. For example,
the impurity shown in Fig. 6(b) yields an additional barrier of $4J_{1}$. In
general, an impurity yields an additional barrier
$\alpha_{1}J_{1}+\alpha_{2}J_{2}$, where the factors $\alpha_{1},\alpha_{2}$
depend on the precise location of the impurity relative to the edge.
Therefore, in the process of peeling off the edge, we encounter the overall
barrier

$E_{B}(L_{d})\simeq 4J_{2}+\gamma(1-2K)L_{d},$ (20)

where the factor $\gamma$ is an average over individual impurity barriers.
Thus, the excess holes give rise to an $L$-dependent barrier in the domain
growth process, yielding a Class 3 system in the LMV classification scheme
[4].

This will modify the asymptotic growth law to a logarithmic form. The relevant
growth equation is obtained from Eqs. (12)-(13) as

$\frac{dL}{dt}=\frac{a_{0}}{L}\exp\left[-4\beta
J_{2}-\beta\gamma(1-2K)L\right].$ (21)

At early times ($t\ll t_{c},L\ll L_{c}$), we can neglect the second term in
the exponent. Then,

$L(t)\simeq\left(2a_{0}e^{-4\beta J_{2}}t\right)^{1/2},\quad t\ll t_{c}.$ (22)

At late times ($t\gg t_{c},L\gg L_{c}$), we approximate $L\simeq L_{c}$ in the
prefactor of the RHS of Eq. (21). This yields the logarithmic growth law:

$L(t)\simeq\frac{1}{\beta\gamma(1-2K)}\ln\left[\frac{a_{0}}{L_{c}}e^{-4\beta
J_{2}}\beta\gamma(1-2K)t\right],\quad t\gg t_{c}.$ (23)

The crossover time and length scales can be obtained by comparison of Eqs.
(22) and (23) – they diverge as $K\rightarrow 1/2^{-}$.

The corresponding results for AF ordering can be obtained by setting $J_{2}=0$
in Eq. (21). A couple of remarks are in order here. Firstly, in the above
discussion, we have only considered the first two terms of $\mathcal{H}_{\rm
CG}$ in Eq. (17). We do not expect the incorporation of higher terms to change
the nature of the crossover ($t^{1/2}\rightarrow$ logarithmic), which is
already captured in the above discussion. However, higher terms will affect
the precise details of the crossover, viz., prefactors of growth laws, etc.
Secondly, for the small asymmetries considered in our simulations, we expect
the crossover to the logarithmic regime to be substantially delayed. The only
signature we see of the asymptotic regime is the under-estimation of the
growth exponent from the LCA value, $\theta=1/2$. Thirdly, it is relevant to
point out that there have been several experimental studies [40, 41] which
have reported logarithmic ordering kinetics in asymmetric mixtures. The above
arguments provide a theoretical basis for understanding the experimental
results.

##  4 Summary and Discussion

Let us conclude this paper with a summary and discussion of our results. We
have studied the domain growth kinetics of a Coulomb glass (CG) model without
any external quenched disorder. The CG is a frustrated system due to the long-
range repulsive Coulombic interaction between electrons. The important
parameter in our simulations is the filling factor $K$ – a system with $N$
sites has $KN$ electrons and $(1-K)N$ holes. For $K=0.5$, there are an equal
number of holes and electrons in the system. Then, the ground state is a
checkerboard where electrons and holes are placed on alternating sites. The CG
with $K>0.5$ is equivalent to the CG with filling fraction $1-K$. Therefore,
without any loss of generality, we consider the CG with $K<0.5$. For
$K\lesssim 0.5$, there is a small excess of holes. In this case, the ground
state is checkerboard-like with excess holes uniformly distributed. For $K\ll
0.5$, more exotic ground states arise, which we do not study in this paper.

We study the phase ordering kinetics in the CG via a kinetic spin model with
long-range antiferromagnetic interactions. For purposes of comparison, we also
study domain growth in the nearest-neighbor antiferromagnet (AF). A suitable
stochastic kinetics is introduced into the system by connecting it with a heat
bath [2]. As the number of electrons is conserved, the appropriate microscopic
kinetics is Kawasaki spin-exchange kinetics, where nearest-neighbor spins are
stochastically interchanged [1]. The relevant order parameter is the staggered
magnetization $\sigma$, which is a nonconserved variable. The kinetics of this
variable is coupled to that of the conserved density field, corresponding to
Model C in the Hohenberg-Halperin framework [31, 32].

Our simulations of the CG and AF with slight asymmetry yield two important
results:
(a) For both models, the domain growth morphology is universal for small
values of asymmetry. The correlation function of the $\sigma$-field is
numerically indistinguishable from that for the symmetric case, $K=0.5$.
Further, the scaling function is in excellent agreement with the Ohta-Jasnow-
Kawasaki (OJK) function, which describes nonconserved ordering dynamics in a
ferromagnet described by Model A [31].
(b) For both models, the domain growth law is somewhat slower than the
Lifshitz-Cahn-Allen (LCA) law: $L(t)\sim t^{1/2}$. Further, the growth law
becomes slower as the asymmetry increases.

To gain a better understanding of the growth law, we study the energy barriers
associated with shrinking of a square droplet in a background of the other
phase. Our results show that the excess holes (impurity spins in the
$\sigma$-domains) present a length-dependent barrier to the peeling process of
droplet edges for both the CG and AF. This barrier is proportional to the
asymmetry ($1-2K$), and is absent at $K=1/2$. Thus, we expect that both the CG
and AF obey the LCA growth law at $K=1/2$. However, for $K\lesssim 1/2$, there
is a crossover to a logarithmic law for $t\gg t_{c}$, where
$t_{c}\rightarrow\infty$ as $(1-2K)\rightarrow 0^{+}$. The values of $1-2K$
studied here are too small for us to explicitly observe the logarithmic
regime. However, we speculate that the under-estimation of the growth exponent
from $\theta=1/2$ in our simulations is a signature of the crossover to the
logarithmic regime. Clearly, much longer simulations with larger systems
(needed to avoid finite-size effects) are required to clearly resolve this
issue. In the context of the long-range interactions considered here, these
simulations would be computationally very demanding.

A promising direction for further research is a study of the ordering dynamics
when the ground states are more exotic, e.g., striped, BCC, FCC, etc. This is
a regime where $K$ is considerably different from 1/2 [34]. Needless to say,
these exotic ground states are of great importance for applications in
materials science and metallurgy.

### Acknowledgments

PB gratefully acknowledges the Kreitman School of Advanced Graduate Studies
for financial support. The HPC facility at Ben Gurion University of the Negev
is gratefully acknowledged for computational resources. VM acknowledges
funding from Science and Engineering Research Board, India under the research
grant CRG/2022/004029.

## References

  * [1] S. Puri and V.K. Wadhawan (eds.), Kinetics of Phase Transitions, CRC Press, Boca Raton (2009).
  * [2] S. Dattagupta and S. Puri, Dissipative Phenomena in Condensed Matter: Some Applications, Springer-Verlag, Berlin (2004).
  * [3] D.A. Huse and C.L. Henley, Phys. Rev. Lett. 54, 2708 (1985).
  * [4] Z.W. Lai, G.F. Mazenko and O.T. Valls, Phys. Rev. B 37, 9481 (1988).
  * [5] S. Puri, Phase Transitions 77, 469 (2004).
  * [6] A.J. Bray and A.D. Rutenberg, Phys. Rev. E 49, R27 (1994).
  * [7] R. Agrawal, F. Corberi, E. Lippiello, P. Politi and S. Puri, Phys. Rev. E 103, 012108 (2021).
  * [8] H. Christiansen, S. Majumder and W. Janke, Phys. Rev. E 103, 052122 (2021).
  * [9] A.L. Efros and B.I. Shklovskii, J. Phys. C 8, L49 (1975).
  * [10] B.I. Shklovskii and A.L. Efros, Electronic Properties of Doped Semiconductors, Springer, Heidelberg (1984).
  * [11] A. Vaknin, Z. Ovadyahu and M. Pollak, Phys. Rev. Lett. 84, 3402 (2000).
  * [12] T. Grenet, Phys. Status Solidi C 1, 9 (2004).
  * [13] Z. Ovadyahu, Phys. Rev. B 73, 214204 (2006).
  * [14] M.T. Shimer, U.C. Tauber and M. Pleimling, Phys. Rev. E 90, 032111 (2014).
  * [15] M. Ben-Chorin, Z. Ovadyahu and M. Pollak, Phys. Rev. B 48, 15025 (1993).
  * [16] G. Martinez-Arizala, C. Christiansen, D.E. Grupp, N. Markovic, A.M. Mack and A.M. Goldman, Phys. Rev. B 57, R670 (1998).
  * [17] T. Grenet, Eur. Phys. J. B 32, 275 (2003).
  * [18] A. Vaknin, Z. Ovadyahu and M. Pollak, Phys. Rev. B 65, 134208 (2002).
  * [19] E. Lebanon and M. Muller, Phys. Rev. B 72, 174202 (2005).
  * [20] A. Barzegar, J.C. Anderson, M. Schechter and H.G. Katzgraber, Phys. Rev. B 100, 104418 (2019).
  * [21] M. Goethe and M. Palassini, Phys. Rev. Lett. 103, 045702 (2009).
  * [22] B. Surer, H.G. Katzgraber, T.G. Zimanyi, A.B. Allgood and G. Blatter, Phys. Rev. Lett. 105, 067205 (2009).
  * [23] A.A. Pastor and V. Dobrosavljevic, Phys. Rev. Lett. 83, 4642 (1999).
  * [24] S. Pankov and V. Dobrosavljevic, Phys. Rev. Lett. 94, 046402 (2005).
  * [25] M. Muller and S. Pankov, Phys. Rev. B 75, 144201 (2007).
  * [26] A. Mobius and U.K. Rossler, Phys. Rev. B 79, 174206 (2009).
  * [27] P. Bhandari and V. Malik, Eur. Phys. J. B 93, 66 (2020).
  * [28] P. Bhandari and V. Malik, Eur. Phys. J. B 92, 147 (2019).
  * [29] P. Bhandari, V. Malik and S.R. Ahmad, Phys. Rev. B 95, 184203 (2017).
  * [30] P. Bhandari, V. Malik and S. Puri, Phys. Rev. E 99, 052113 (2019).
  * [31] P.C. Hohenberg and B.I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * [32] P. Das, T. Saha-Dasgupta and S. Puri, Eur. Phys. J. E 40, 94 (2017).
  * [33] F. Corberi, E. Lippiello, A. Mukherjee, S. Puri and M. Zannetti, Phys. Rev. E 85, 021141 (2012).
  * [34] L. Rademaker, Y. Pramudya, J. Zaanen and V. Dobrosavljevic, Phys. Rev. E 88, 032121 (2013).
  * [35] K. Binder and D.W. Heermann, Monte Carlo Simulation in Statistical Physics: An Introduction, 4th Edition, Springer-Verlag, Berlin (2002).
  * [36] M.E.J. Newman and G.T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford University Press, Oxford (1999).
  * [37] D. Frenkel and B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, San Diego (1996).
  * [38] R. Paul, S. Puri and H. Rieger, Europhys. Lett. 68, 881 (2004); R. Paul, S. Puri and H. Rieger, Phys. Rev. E 71, 061109 (2005).
  * [39] J.D. Shore, M. Holzer and J.P. Sethna, Phys. Rev. B 46, 11376 (1992).
  * [40] R.F. Shannon, C.R. Harkless and S.E. Nagler, Phys. Rev. B 38, 9327 (1988).
  * [41] S.B. Rivers, W.N. Unertl, H.H. Hung and K.S. Liang, Phys. Rev. B 52, 012601 (2005).

Figure 1: Schematic of ground states of the (a) CG and (b) AF with slightly
asymmetric composition: $K\lesssim 0.5$. The circles denote electrons, and the
bars denote holes. The surplus holes are colored red. Figure 2: Domain growth
in the $d=2$ CG at $K=0.45$. We show snapshots at different MCS for a
$512^{2}$ lattice. The order parameter is the staggered magnetization
$\sigma$. The black and white regions correspond to $\sigma>0$ and $\sigma<0$,
respectively. Figure 3: The upper frames show snapshots at $t=10^{3}$ MCS from
domain growth in the $d=2$ CG of size $512^{2}$. We show snapshots at
different filling factors $K$. The order parameter is the staggered
magnetization $\sigma$. The lower frames show corresponding snapshots for the
AF. These frames show a $512^{2}$ corner of a $2048^{2}$ lattice. Figure 4:
(a) Scaled correlation functions, $C(r,t)$ vs. $r/L$, for the CG at $t=10^{3}$
MCS and different $K$. The length scale $L(t)$ is defined as the distance over
which $C(r,t)$ decays to 0.2 of its maximum value. (b) Analogous plot for the
AF. The solid line in both frames denotes the OJK function in Eq. (11). Figure
5: (a) Time-dependence of the characteristic length scale, $L(t)$ vs. $t$ on a
log-log scale, for the CG. We plot data for different $K$-values. The line
labeled $t^{1/2}$ denotes the LCA law. (b) Plot of effective exponent for CG:
$\theta_{\rm eff}$ vs. $t$. The horizontal line denotes the LCA exponent,
$\theta=1/2$. (c) Analogous to (a), but for the AF. (d) Analogous to (b), but
for the AF. Figure 6: Schematic for shrinking of a square droplet (of size
$L_{d}^{2}$) with $\sigma=-1$ in a background of $\sigma=+1$ in the CG at (a)
$K=0.50$ and (b) $K\lesssim 0.50$. We approximate $\mathcal{H}_{\rm CG}$ by
the first two terms in Eq. (17). In (a), there is a barrier for flipping the
corner spin. The edge subsequently peels without any further barriers. In (b),
impurity spins in the background/droplet yield further $L_{d}$-dependent
barriers to the peeling of the edge.