\#1

Introduction
============

The dynamics of supercooled liquids approaching the glass transition remains one of the most fundamental problems in condensed matter physics [@review]. Some of the prominent dynamic features in supercooled liquids include the enormous increase in relaxation time scale with lowering temperature, and the nonexponential relaxation in the response to an external perturbation. In addition to these features, an anomaly in transport properties such as breakdown of the Stokes-Einstein (SE) relation in highly supercooled liquids has been observed in experiments [@fujara] and simulations [@mount; @harrowell; @onuki]. Although there exist some theoretical attempts [@stilinger; @tarkiv; @liuopp; @bb; @xia], the underlying microscopic mechanism for the violation of the SE relation is not well understood. Recently, there have been many experimental and simulational studies of supercooled liquids that demonstrate the existence of kinetic heterogeneity which was often invoked to explain the origin of the non-exponential relaxation as well as the breakdown of the SE relation [@sillescu].

In relation to these questions on microscopic slow dynamic features in supercooled liquids, we deemed it worthwhile to investigate whether similar dynamic features can be found in simpler lattice spin systems or lattice gas systems. In this work, we show that the aforementioned features of supercooled liquids, [*i.e.*]{}, slowing-down, non-exponential relaxation and the (analogue) of the breakdown of the SE relation, are also observed in a two-dimensional (2D) lattice Coulomb gas (LCG) system. We also find that the relaxation of the system exhibits a spontaneous appearance of spatial heterogeneity, which we argue is the underlying cause for the non-exponential relaxation and the breakdown of the SE relation.

In recent years, there have been some efforts to find glassy dynamic features in the lattice spin systems with nonrandom interactions [@bulbul]. One of the well-known examples of disorder-free lattice model system is uniformly frustrated XY (UFXY) models in two dimensions, which serve as a model for two dimensional arrays of Josephson junctions under the influence of uniform transverse magnetic fields. Recent work [@kl] has shown that, irrespective of the true nature of the low temperature phase of this system, the equilibrium dynamics of UFXY model in the intermediate range of the temperature for frustration parameter $f$ near $1-g \equiv (3-\sqrt{5})/2 \simeq 0.382$, exhibits a close analogy to that of supercooled liquids. Both spin and vortex dynamics show stretched exponential relaxations with temperature dependent stretched exponents. In order to investigate the dynamics of this system in more detail, we attempted to calculate the self diffusion properties of vortices. However, it turned out to be numerically ambiguous and tricky to trace the trajectories of individual vortices. This is beacause individual vortex around a plaquette is defined in terms of phase angles and one probes the movement of individual vortices not directly but only indirectly through changes of phases, which, at times especially when multi-vortex motion occurs, makes it ambiguous to determine the original position of a vortex corresponding to a new neighboring vortex.

One way to overcome this difficulty was to map the UFXY model onto a LCG via Villain transformation [@villain], where the positive charges in the LCG correspond to the positive current vortices in UFXY models. One can readily probe the diffusive dynamics of charges without ambiguity in the LCG unlike the case of UFXY model. Hence we can investigate both the structural relaxation dynamics and self diffusion dynamics of individual vortices in LCG with relative ease.

With this advantage, we have numerically investigated the equilibrium relaxation dynamics and diffusion characteristics of LCG with charge density factor $f$ near $1-g \simeq 0.382$. We observe that for some range of temperatures above the first order transition, the equilibrium relaxation exhibits slow dynamic features such as stretched exponential relaxation and breakdown of proportionality between diffusive time scale and structural relaxation time scale.

It was a common belief that the 2D UFXY model and the corresponding LCG belong to the same universality class with essentially the same phase transition properties, ground state symmetry, for example. However, recent work on LCG by Gupta, Teitel, and Gingras (GTG)[@gtg] and also another work on UFXY model by Denniston and Tang (DT)[@dt], showed that there exist some difference between the two model systems especially in the case of dense frustration. Both model systems exhibit first order transition but the low temperature vortex configurations in UFXY models are different from the charge configurations of the corresponding LCG for $f$ near $1-g \simeq 0.382$. The underlying cause for this breakdown of Villain approximation in the limit of dense frustration is not known, but probably it is related to application of spin-wave integration to systems having many metastable states with similar energies, that may cause neglect of multi-vortex correlations.

Special interest has been given to the case of $f$ approaching $1-g$ [@halsey; @mychoi]. Consider a system where $f$ equals $p_0 /q_0$ ($p_0$ and $q_0$ are relative primes) which is a rational approximant to $1-g$. Here, in the case of UFXY model, DT argues that the low temperature vortex configuration has lattice periodicity which is of order $q_{0}^{2}$, [*i.e.*]{}, much larger than $q_0$. On the other hand, in the case of LCG, GTG [@gtg] showed, via Monte Carlo (MC) simulations, that the low temperature charge configurations are characterized by arrangements of diagonal stripes that are either completely filled, completely empty, or partially filled with charges that are quite different from those vortex configurations in the corresponding UFXY model. However, GTG did not enumerate the exact patterns of low temperature charge configurations (such as spatial periodicity) for general cases of dense charge filling. In this work, we find that, for the values of $f$ between $1/3$ and $2/5$, there exist a simple regularity in the low temperature charge configuration which consists of periodic arrangements of combinations of two out of three types of striped charge patterns (see Section III).

For wide range of quenching temperatures above the first order transition $T_c$, the equilibrium relaxation continues to slow down with lowering temperature, and the form of the relxations are characterized by the stretched exponential with temperature-dependent exponents. Moreover, we observe that the model exhibits a separation of the two characteristic time scales, i.e., the time scale of single particle diffusion and that of structural relaxation. This feature is quite analogous to the breakdown of the SE relation observed in supercooled liquids. Stretched exponential relaxation is observed to be accompanied by interesting dynamic heterogeneity in the system. It appears that the kinetic heterogeneity is the underlying reason for both the stretched exponetial relaxation and the separation of the relaxation and diffusion time scales.

A convenient measure for dynamic heterogeneity is the so called dynamic cooperativity [@doliwa] of the particle motions. This measures reduction of the effective degrees of freedom. One interesting result from our simulations is that the magnitude of the velocity (or displacement vector) exhibits strong increase in cooperativity of the particle motions. On the other hand, the displacement vector itself shows cooperativity a little smaller than unity due to anti-correlations in the direction of particle motions. This means that the system can be divided into highly mobile regions and relatively inert regions resulting in highly inhomogeneous local mobility distribution. However, there is no macroscopic flow of particles that will generate long range positive correlations between the directions of flows of particles.

When quenched to a temperature below $T_c$, the system is always found to undergo phase ordering via slow coarsening processes. The system therefore does not remain in a supercooled state. Rather it becomes slowly crystalized. It should be emphasized that in this system it is the relaxation for the temperatures [*above*]{} $T_c$ that exhibits slow dynamic behavior which shares some common features with that of supercooled liquids.

Model and Simulation Methods
============================

General 2D LCG [@teitel] is described by the following Hamiltonian that can be mapped from UFXY model by means of Villain transformation [@villain],

$${\cal H}_{CG}={1 \over 2}\sum_{ij}Q_iG(r_{ij})Q_j \label{eq:Hcg}$$

where $r_{ij}$ is the distance between the sites $i$ and $j$, and the magnitude of charge $Q_i$ at site $i$ can take either $1-f$ or $-f$ , where $f$ corresponds to the frustration parameter in the related XY models. Charge neutrality condition $\sum_iQ_i=0$ implies that the number density of the positive charges is equal to $f$. As was mentioned above, we can thus view the system as a lattice gas of $N \cdot f$ charges of unit magnitude upon uniform negative background charges of charge density $-f$ ($N= L^2$ is the total size of the system with the linear dimension $L$). The lattice Green’s function $G(r_{ij})$ solves the equation

$$(\Delta^2 - { 1 \over {\lambda^2}}) G(r_{ij})=-2\pi\delta_{r_{ij},0}$$

where $\Delta^2$ is the discrete lattice Laplacian and $\lambda$ is the screening length which, in normal case of no screening, is set to an infinity. For the case of usual Villain transformation of UFXY model, we have $\lambda = \infty$. But it is included in this equation for generality. Since, in this work, we restrict our attention to only square lattice with periodic boundary conditions, $G(r)$ is given by

$$G(\vec{r})={\pi \over N} \sum_{\vec k \neq 0} { e^{i\bf{k} \cdot r}-1 \over 
{2 - \cos k_x -\cos k_y + 1/\lambda^2} } ,$$

where $\bf{k}$ are the allowed wave vectors with $k_{\mu} = (2 \pi n_{\mu}/ L)$, with $n_{\mu} = 0, 1, \dots, L-1 $. In the case of infinite screening length, for large separation $r$, one gets $G(\vec r)\simeq -\ln r$ [@jrl]. In this work, we consider the limiting case of $ \lambda \rightarrow \infty $ only.

In our MC simulations, the initial disordered random configuration is updated according to the standard Metropolis algorithm by selecting a positive charge at random and moving it over to one of the [*nearest neighbor (NN)*]{} or [*next nearest neighbor (NNN)*]{} sites [@gtg]. We find that this [*NNN*]{} hopping algorithm is particularly effective in terms of simulation time as compared with [*NN*]{} hopping alone, as was emphasized in [@gtg]. Moreover, at low temperature, [*NN*]{} hopping alone presented severe energy barriers to the motions of charges in the case of relatively dense Coulomb gas, [*i.e.*]{}, $f$ approximately larger than $1/3$.

The presented results are averages over $100 \sim 500$ different random initial configurations depending on the temperature. In order to ensure that equilibration is achieved, we calculate the two-time charge density autocorrelation function and locate the waiting time beyond which the autocorrelation function no longer depends on the waiting time. As for the values of the charge density parameter $f$, we chose $f = 55/144 \simeq 0.3819$ that is close to $f=1-g$, and square lattices of linear size $L=36$ is chosen with periodic boundary conditions. This value of $f$ is chosen as a simple rational value that satisfies the two conditions of both being close to $1-g$ and being commensurate with the lattice periodicity $12$ as explained in Section III. We found that qualitative features of relaxation dynamics are the same for other nearby values of the frustration $f$.

Simulation Results and Discussions
==================================

First order transition and low temperature configuration
--------------------------------------------------------

We first discuss the equilibrium phase transition and charge configuration of the system. As was first shown by GTG, we also find that there exist a first order transition in LCG with $f$ near $1-g$. Figure 1 shows temporal snapshots of charge configurations evolving from disordered state into ordered configuration after being quenched to a temperature $T=0.026$. First order nature of the phase transition can easily be confirmed by enumerating the histogram of energies $P(E)$ near the transition temperature [@sw]. $P(E)$ is obtained by counting the occurrences of energies for each of the equally spaced energy bins while performing the equilibrium Monte Carlo simulations (via simple Metropolis algorithm). For a system with first order transition, the energy histogram $P(E)$ becomes bimodal near the transition temperature corresponding to a mixture of ordered state (with lower energy) and a disordered state (with higher energy). The transition temperature $T_c$ can be determined by locating the temperature where the subareas under the two peaks are equal. Figure 2 shows two histograms near the transition temperature, where we could estimate the transition temperature approximately as $T_c \simeq 0.0316$. Since we did not attempt a detaled analysis (including a finite size scaling) of the histogram, we think that this estimate value of the transition temperature should not be considered too seriously for its precision.

We find empirically that there exist a simple regularity in the low temperature charge configuration in LCG (Fig. 3). For the case of values of $f$ in the range $1/3 \leq f \leq 2/5$, it is found that the low temperature configuration becomes quasi-one-dimensional with periodic striped patterns. In the cases of $f=1/3$ and $f=2/5$ the ground state configurations are identical to the low temperature vortex configurations in the UFXY model. However, for values of $f$ in between $1/3$ and $2/5$, the low temperature patterns are found to be, unlike the case of corresponding UFXY model, consisting of periodic arrangements of combinations of two out of three types of striped charge patterns as follows.

First component pattern (type I pattern) is a sequence of three diagonals which are [*empty, filled*]{}, and [*empty*]{} respectively (that may be denoted by [**(010)**]{} in our notation where [**1**]{} refers to a filled diagonal and [**0**]{} refers to an empty diagonal). In other words, it is a pattern with single isolated diagonal filled with charges, that is neighbored by empty diagonals on both sides. Repetition of this pattern alone produces the ground state configuration for the case of $f=1/3$ with spatial periodicity three.

Second component pattern (type II pattern) consists of a sequence of five diagonals that are [*empty, filled, empty, filled*]{}, and [*empty* ]{} respectively, or [**(01010)**]{} in our notation. This may be termed as a double filled diagonal because two filled diagonals are positioned in parallel at second neighbor. This forms the basis of the ground state configuration for the case of $f=2/5$ with lattice periodicity five.

Lastly, the third component pattern (type III pattern) consists of a sequence of seven diagonals that are sequentially [*empty, filled, empty, partially filled, empty, filled*]{}, and [*empty* ]{} [*i.e.*]{}, [**(010p010)**]{} in our notation where [**p**]{} refers to a partially filled diagonal where only part of the diagonal sites are occupied by positive charges. This pattern is essentially a partially filled diagonal enveloped by two filled diagonals on both sides at second neighbor diagonal position, which may be termed as a [*channel*]{} structure. This can form a basis with spatial lattice periodicity seven.

We leave the detailed description of the low temperature charge patterns for the full range of $f$ values between $1/3$ and $2/5$ to the forthcoming publication [@config_cg]. And we describe in this work the low temperature ordered patterns for values of $f$ around $1-g$ only. Near the value of the filling ratio $f=1-g \simeq 0.382$, we find that, among the three patterns above, only two types (type II and type III patterns) participate in the stable charge configurations with the resulting spatial lattice periodicity depending on the combination of the two component patterns.

We find that there exist a value $f=f_c \simeq 0.381$ which separates two regimes with distinct low temperature striped patterns. For values of $f$ in the range $0.36 \lsim f \lsim 0.381$, the stable striped patterns turn out to have periodicity $l_p = 7$ which consists of simple repetitions of channel structures (type III pattern). Note that this periodicity seven refers to the periodicity of the filled diagonals only (neglecting the true periodicity including the charge configurations within the partially filled diagonals).

On the other hand, for values of $f$ in the range $0.381 \lsim f \lsim 0.39$, the stable configuration exhibits a periodicity $l_p=12$, which consists of double filled diagonals (type II) and channels (type III) alternatingly placed. As the value of $f$ continuously increases within the two regimes (in the above), the system in the low temperature stable configuration simply adjusts itself by accomodating the extra charges into the partially filled diagonal channels and thereby changing the charge filling within the channels. The dividing value of $f=f_c \simeq 0.381$ between the two regimes appears to correspond to the value $8/21$ in which case the partially filled diagonals have filling density exactly equal to $2/3$. Our simulations show that the filling density $2/3$ inside the partially filled diagonal plays as a stability limit for the channel structures. Beyond this limit, electrostatic instability probably begins to set in, and rearrangement of the whole charge configuration occurs in order to form a new stable ordered patterns. As was also argued by GTG, in general, at much lower temperature $T_p$ (below $T_c$) the charges within the partially filled channels are expected to exhibit ordering, which would depend sensitively on rationality of the exact filling ratio of charges inside the partially filled diagonals.

An important aspect of our simulations is that one has to choose the lattice size appropriately in order to match the periodicity of the true low temperature configuration in the thermodynamic limit. If, otherwise, one chooses a lattice size that is incommensurate with the periodicity (of striped patterns), then one ends up with defective charge configurations with patches of local ground state configurations. We think that this is probably why GTG got two different equilibrium configurations when two different lattice sizes $L=26$ and $L=52$ are used for $f=5/13$ since these $L$’s turn out to be incommensurate with the true periodicity $l_p =12$.

When the screening length $\lambda$ is finite, then we find the low temperature configuration becomes different from the case of no screening ($\lambda \rightarrow \infty$) in such a way that the partially filled diagonals gets rarer. The influence of the screening effect on the statics and the relaxation dynamics needs further study.

Equilibrium relaxation dynamics
-------------------------------

We now discuss the equilibrium relaxation dynamics of the model above first order transition. In order to probe the structural relaxation of charges, we measured the on-site charge autocorrelation functions,

$$C(t)=\langle \sum_{i=1}^{N} Q_i(0)Q_i(t) \rangle / N,$$

where the bracket $< \cdots >$ represents an average over different random initial configurations.

Shown in Fig. 4a is the on-site charge autocorrelation function $C(t)$ for temperatures from $T=0.1$ down to $T=0.033$. From this figure, we observe a slowing down in the structural relaxation for this temperature range. One can extract a characteristic time scale $\tau(T)$ which, for example, is defined as $C(t=\tau(T))=1/e$ for each temperature $T$. As Fig. 4b clearly shows, the temperature dependence of the relaxation time exhibits a non-Arrhenius behavior. We also checked whether the so-called time-temperature superposition holds for the above autocorrelation functions, which is shown in Fig. 4c. We clearly see that time-temperature superposition is systematically broken by the autocorrelation functions. This is consistent with the fact that the stretched exponents have dependence on temperature as is shown just below.

We find that the relaxation pattern of the correlation function $C(t)$ can be characterized by a power law relaxation $C(t)=1-A t^{b(T)}$ (known as the von Schweidler relaxation) in the early time regime and a stretched exponential relaxation $C(t)=C_0(T) \exp(-A't^{\beta(T)})$ in the late time regime. However, as the temperature gets higher, the regime of validity for early time power law relaxation was significantly reduced and we could better fit the early time relaxation with another stretched exponential form $C(t)=\exp(-A''t^{b'(T)})$. Of course for low temperature regime, we could get $b(T) \simeq b'(T)$.

Fig. 4d shows the temperature dependence of the fitted exponents. We see that non-exponentiality increases as the temperature decreases. These results clearly indicate that the equilibrium relaxation in the 2D LCG above $T_c$ closely resembles the primary relaxation of typical fragile liquids.

One of the main characteristic features of the single particle dynamics is described by the mean square displacement $\langle (\Delta \vec r)^2 \rangle$, which is defined as

$$\langle (\Delta \vec r)^2 \rangle = \langle \sum_{j=1}^{N_Q}
 (\vec{r}_j(t)-\vec{r}_j(0))^2 \rangle / N_Q,$$

where $\vec{r}_j(t)$ is the position vector of the $j$-th charge at time $t$ and $N_Q$ the total number of charges. Figure 5 shows $ \langle (\Delta \vec r)^2 \rangle $ for various temperatures. It exhibits an early time subdiffusive regime and crosses over into late time diffusive regime. Early time subdiffusive behavior is thought to be coming from local frustrated motions of charges before reaching an average displacement of unit lattice spacing. To test the proportionality of the two time scales, the structural relaxation time scale $\tau$ and the diffusion time scale $D^{-1}$, we plot the temperature dependence of the product $4 D\tau$ in Fig. 6. Here, we clearly see that the breakdown of the proportionality between the two time scales is observed for wide range of temperatures below $T=0.1$ and becomes stronger as the temperature is lowered. This separation of the two time scales is due to the weaker temperature dependence of the diffusion coefficient. That is, diffusion is relatively enhanced at lower temperatures. This is quite analogous to the violation of the SE relation ($D=T/a \eta$, where $a$ is a molecular length and $\eta$ is the viscosity of the liquid) observed in experiments on supercooled liquids [@fujara]. Here, we mention that there exists a correlation between the increase of non-exponentiality (as the temperature is lowered) and the increase of the product $4 D\tau$ at low temperatures [@ediger].

If we suppose that there exists a single dominant relaxation mode in the system (and hence one relaxation time scale $\tau$), then we would obtain a simple exponential behavior for the relaxation function $C(t) \sim e^{-t/\tau}$. On the other hand, if the system consists of many regions with different relaxation times, then the relaxation function would be roughly some superposition of exponentials with a broad distribution of relaxation times, which would be in general not expressible in a simple exponential form, but in stretched exponetial form or other more complicated forms.

The fact that there exists a breakdown of proportionality between $\tau$ and $D^{-1}$ can be interpreted in the following way that invokes dynamic heterogeneity. As the temperature is lowered, the system consists of many regions with different relaxation time that comes from different local mobilities. We can easily see that the structural relaxation time is dominated by the least mobile regions, that is, by the regions with the longest relaxation time. On the contrary, the average (long time) diffusion characteristics is dominated by the most mobile regions. In other words, the structural relaxation function and the self-diffusion function, respectively, are probing more or less opposite aspects of the relaxation behavior of the system. For an extreme example, one can imagine a system where half of the whole system is frozen (no motion of the component particles) while the remaining half of the system has finite relaxation time with uniformly distributed mobile particles. For this system, the structural relaxation time $\tau$ would be infinite due to the frozen half of system, but inverse of the average diffusion constant $D^{-1}$ is finite due to the mobile part of the system, leading to an extreme breakdown of SE relation. Above simulation result, thus, can be interpreted as an evidence pointing toward the existence of a kinetic heterogeneity in the relaxation dynamics and mobility of the system.

In fact, the kinetic heterogeneity can be visualized in our system. Typical charge configuration at $T=0.033$, as shown in Fig. 7, exhibits local striped patterns (ordered domains) and interfacial regions due to mismatch between adjacent domains. For a fixed quenching temperature, the average size of these local domains reaches a certain length scale when the system equilibrates. After equilibration, the system structurally rearranges itself going from one configuration to another with local domains of similar length scale. Figure 8 shows the trajectories of moving positive charges over a time interval of 500 MC steps for $T=0.033$ (corresponding to Fig. 7). We can see that there exist local regions with actively moving charges and other regions with relatively immobile charges. Among the active regions, we can find those charges moving along partially filled diagonal channels. We also find some extended interfacial regions where no discernible local order can be identified, that exhibit relatively high mobility. Enhancement of particle diffusion is probably due to the motions of charges along the partially filled diagonals as well as those fluidized motions in the extended interfacial regions. These fastly moving regions in surroundings of very slowly moving regions offer a specific example for spatial heterogeneity in glassy systems [@harrowell; @onuki], which was often thought of as the physical mechanism for breakdown of the SE relation.

One simple way to quantify the degree of dynamic heterogeneity directly from the local motions of particles is to calculate the dynamic cooperativity [@doliwa] for one particle dynamic quantities such as [*e.g.*]{}, displacement vectors $X_i \equiv |\vec{r}_i (t+ \Delta t) - \vec{r}_i (t) | $ between the time $t$ and $t+\Delta t$ for some fixed time interval $\Delta t$. We can also choose $X_i$ to be the vector displacement itself $X_i \equiv \vec{r}_i (t+ \Delta t) - \vec{r}_i (t) $. If there are no correlations between the motions of particles, then the variations of the $X_i$’s will satisfy

$$\sigma [ \sum_{i} X_i ] = \sum_{i} \sigma [X_i ],$$

where $\sigma[x]$ denotes the mean square deviations of the random number $x$, $\sigma [x] \equiv \langle (x-<x>)^2 \rangle$. However, some correlations between the particle motions will increase $\sigma [ \sum_{i} X_i ]$ or anti-correlations will decrease it. Following Doliwa and Heuer, we can define the dynamic cooperativity as $$N_{X}^{coop} \equiv {{\sigma [ \sum_{i} X_i ] } \over 
{ \sum_{i} \sigma [X_i ]}}.$$ In the case of no correlations between the motions of particles, as in (6), we get $N_X^{coop} =1$. If there exist some positive correlated motions between particles, we would get $N_X^{coop} > 1$, while anti-correlations between the motions of particles would result in $N_X^{coop} < 1$. Doliwa and Heuer investigated the dynamic cooperativity of hard sphere systems in 2D and 3D, where they found finite cooperativity ($N_X^{coop} > 1$) for both vector displacement and the scalar magnitude of the displacement, which is consistent with the snapshots of the particle motions in their work. They argue that the dynamic cooperativity measures the total reduction of degrees of freedom due to the correlations. Here we also studied the dynamic cooperativity of the lattice gas particles by calculating $N_{X}^{coop}$ for both the scalar displacement and the vector displacement itself. Interestingly, we found that the scalar displacement exhibited finite dynamic cooperativity (Fig. 9a), while the vector displacement itself showed weak anti-correlations between particles. as shown in Fig. 9b. In the case of scalar displacement, the cooperativity increases at first as the time interval $\Delta t$ increases and reaches its maximum near the $\alpha$-relaxation time scale $\tau$. Then it decrease back to values around unity (corresponding to no correlations) at large $\Delta t$.

Contrasting features of cooperativity for our LCG system and that for the hard sphere systems may be interpreted as follows. In the case of hard sphere systems near the glass transition, the packing density is very high and the inter-particle interaction is a short ranged one. Therefore, the local motions of particles in hard sphere systems are naturally highly correlated in both its direction and magnitude due to the continuity constraint of particles resulting in a large scale flow with directional correlations.

In contrast, in the case of the LCG, the density of particles is relatively low ($f \simeq 0.38 $) as compared with the case of hard sphere systems near the glass transition. In addition to that, charge motions in the LCG is driven by thermal effect. From the snapshots of charge configurations, we see that there exist locally mobile regions as well as locally immobile regions. Locally immobile regions consist of charge configurations that are close to the low temperature striped patterns. Mobile regions, however, consist of charges that are agitated in random directions due to the thermal effect. Thus we do not observe positive dynamic cooperativity in vector displacement, but only the scalar displacement exhibits appreciable positive cooperativity due to the local regions with high mobilities. Hence, heterogeneity still exists in our lattice Coulomb gas in terms of local mobility distribution, but unlike the case of hard sphere systems, there is no appreciable average local flow.

Also, we may look into the nature of the equilibrium dynamics of the system in wave-vector space. Figure 10 shows the structure factor $S(q) \equiv \langle |\rho_q|^2 \rangle $ at equilibrium where $\rho_q \equiv \sum_j \exp [i 
\vec{q} \cdot \vec{r}_j ]/N $ where $q= {{2 \pi }  \over L} n $, $n=1,2, \cdots 2/L $. We see that the structure factor of our LCG shows some similarity to those of dense liquids with first peak corresponding roughly to the inverse of the average distance between charges. Due to the lattice nature of the LCG, the wave vector has cutoff value at $q_{max}=\pi$ as in the figure.

The diffusive properties of the system can be probed by calculating the incoherent scattering function (ISF) $F_S(q,t)$ which is defined as in our model of LCG

$$F_S (q,t) \equiv \langle \sum_{j=1}^{N_Q} \exp i \vec q 
\cdot [\vec r_j(t)- \vec r_j (0)]
\rangle / N_Q,$$

where $\vec r_j(t)$ denotes the position of $j$-th particle on the lattice. Due to the discrete lattice nature of our model sytem, we need to consider the wave-vectors within the first Brillouin zone $q= {{2 \pi }  \over L} n $, $n=0,1,2, \cdots L-1 $. Figure 11 shows the $q$-dependence of $F_S(q,t)$ at temperature $T=0.033$. We find that the long-time behavior of $F_S(q,t)$ also can be fitted to stretched exponential form. For low $q$, the late time $\beta$ exponents were close to one (pure exponential relaxation) but as $q$ increases the exponents decreased down to $\beta \approx 0.73$ for $q = 18 \times 2 \pi /36 $, and $T=0.033$ (Fig. 12). As can be seen from the definition of $F_S(q,t)$, for gaussian distribution for the displacement vector $\Delta \vec{r}_i$, we would get $$F_{G} (q,t) \equiv \langle \exp iq [\Delta r ] \rangle 
= \exp [ -{{q^2} \over 2} \langle (\Delta r )^2 \rangle ].$$ Figure 13 shows that the gaussian approximation is quite good for low $q$. That is, for long distance diffusion, the distribution gets closer to gaussian. However, as $q$ becomes larger, the gaussian approximation gets worse as shown in the figure. Similar features were reported in molecular dynamics simulations on the dynamics of supercooled water [@sciortino].

In summary, we have shown that the 2D LCG with fractional filling of charges exhibits an equilibrium relaxation behavior, above first order melting transition, characterized by two time-regimes of stretched exponetial form with temperature dependent exponents, which is quite similar to the primary relaxation of typical supercooled liquids. We found a strong deviation from proportionality between the diffusive time scale and the structural relaxation time scale resembling the breakdown of SE relation in supercooled liquids. This is accompanied by a characteristic dynamic cooperativity, where the scalar displacement exhibits positive cooperativity while the vector displacement shows anti-correlations leading to the vector cooperativity less than unity. We have identified the microscopic heterogeneous structure which is responsible for this phenomena.

We thank M. D. Ediger, P. Harrowell, K. Kawasaki and S. Teitel for discussions. This work was supported by the Korea Research Foundation Grant (KRF-1999-015-DP0098) (SJL, BK) and (KRF-1998-15-D00089) (JRL).

[99]{} M. D. Ediger, C. A. Angell, and S. R. Nagel, J. Phys. Chem. [**100**]{}, 13200 (1996); most recent developments on the subject can be found in the collection of papers in the conference proceedings such as J. Non-Cryst. Solids [**235-237**]{} (1998), J. Phys. C: Condens. Matter [**10A**]{} (1999), and J. Phys. C: Condens. Matter [**12**]{} (2000).

F. Fujara, B. Geil, H. Sillescu, and G. Fleischer, Z. Phys. B [**88**]{}, 195 (1992); N. Menon, S. R. Nagel, and D. C. Venerus, Phys. Rev. Lett. [**73**]{}, 963 (1994); M. T. Cicerone and M. D. Ediger, J. Chem. Phys. [**104**]{}, 7210 (1996) and references therein.

D. Thirumalai and R. D. Mountain, Phys. Rev. E [**47**]{}, 479 (1993).

D. N. Perera and P. Harrowell, Phys. Rev. Lett. [**81**]{}, 120 (1998).

R. Yamamoto and A. Onuki, Phys. Rev. E [**58**]{}, 3515 (1998); Phys. Rev. Lett. [**81**]{}, 4915 (1998).

J. A. Hodgdon, and F. H. Stillinger, Phys. Rev. E [**48**]{}, 207 (1993); F. H. Stillinger and J. A. Hodgdon, [ibid]{}, [**50**]{}, 2064 (1994).

G. Tarjus and D. Kivelson, J. Chem. Phys. [**103**]{}, 3071 (1995).

C. Z.-W. Liu and I. Oppenheim, Phys. Rev. E [**53**]{}, 799 (1996).

S. Bhattacharyya and B. Bagchi, J. Chem. Phys. [**107**]{}, 5852 (1997).

X. Xia and P. G. Wolynes, cond-mat/0101053.

For recent review on the heterogeneity, see H. Sillescu, J. Non-Cryst. Solids [**243**]{}, 81 (1999) and references therein.

For some of the most recent works, see B. Chakraborty, L. Gu, and H. Yin, J. Phys.: Condens. Matter [**12**]{}, 6487 (2000); A. Lipowski and D. A. Johnson, Phys. Rev. E [**61**]{}, 6375 (2000); M. Swift, H. Bokil, R. D. M. Travasso, and A. J. Bray, Phys. Rev. B [**62**]{}, 11494 (2000).

B. Kim and S. J. Lee, Phys. Rev. Lett. [**78**]{}, 3709 (1997); S. J. Lee and B. Kim, Phys. Rev. E [**60**]{}, 1503 (1999).

J. Villain, J. Phys. (Paris) [**36**]{}, 581 (1975); J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B [**16**]{}, 1217 (1977).

P. Gupta, S. Teitel, and M. J. P. Gingras, Phys. Rev. Lett. [**80**]{}, 105 (1998).

C. Denniston and C. Tang, Phys. Rev. B [**60**]{}, 3163 (1999).

T. C. Halsey, Phys. Rev. Lett. [**55**]{}, 1018 (1985); Physica B [**152**]{}, 22 (1988).

M. Y. Choi and D. Stroud, Phys. Rev. B [**32**]{}, 7532 (1985); Phys. Rev. B [**35**]{}, 7109 (1987); J. S. Chung, M. Y. Choi, and D. Stroud, Phys. Rev. B [**38**]{}, 11476 (1988); S. Y. Park, M. Y. Choi, B. J. Kim, G. S. Jeon, and J. S. Chung, Phys. Rev. Lett. [**85**]{}, 3484 (2000).

B. Doliwa and A. Heuer Phys. Rev. E [**61**]{}, 6898 (2000).

S. Teitel, [*Equilibrium Phase Transitions in Josephson Junction Arrays*]{} in Proceedings of the Sitges Conference on Glassy Systems, E. Rubi, Springer, Berlin (1996); J. P. Straley and G. M. Barnett, Phys. Rev. B [**48**]{}, 3309 (1993).

For further details, see J.-R. Lee and S. Teitel, Phys. Rev. B [**46**]{}, 3247 (1992).

A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. [**61**]{}, 2635 (1988); C. Borgs and R. Kotecký, J. Stat. Phys. [**61**]{}, 79 (1990).

S. J. Lee, B. Kim, and J.-R. Lee unpublished.

M. Cicerone and M. D. Ediger, J. Chem. Phys. [**104**]{}, 7210 (1996).

F. Sciortino, L. Fabbian, S. H. Chen, and P. Tartaglia, Phys. Rev. E [**56**]{}, 5397 (1997).

**FIGURE CAPTIONS**

1.  Snapshots of charge configuration at time steps (a) t=16 MCS, (b) t=4096 MCS, (c) t=65536 MCS, and (d) t=1048576 MCS, for temperature $T = 0.026$ and $f=55/144$, exhibiting coarsening toward an ordered striped state. Positive charges are represented by filled squares.

<!-- -->

1.  Energy histogram near the first order transition temperature (for $T = 0.03165$ and $T=0.0317$).

<!-- -->

1.  Regimes of charge patterns for the range of value of $f$ between $1/3$ and $2/5$. See the text for details.

<!-- -->

1.  \(a) The charge autocorrelation functions for temperatures $T = 0.1$, $0.08$, $0.06$, $0.05$, $0.042$, $0.037$, $0.035$, $0.033$. (b) Arrhenius plot for the relaxation time ($\log (\tau)$ versus $1/T$). (c) Charge autocorrelation functions in (a) replotted in terms of the rescaled time $t/ \tau (T)$ which shows that the time-temperature superposition is broken. (d) Temperature dependence of the $b$ and $\beta$ exponents for charge autocorrelation functions.

<!-- -->

1.  Squared displacement $W(t)$ versus time $t$ for the same temperatures as in Fig. 4a.

<!-- -->

1.  Comparison of the two time scales $D^{-1}$ and $\tau$ ($4D \tau $ versus $T$), which clearly shows that the diffusive time scale increases slowly (as the temperature is lowered) as compared with the structural relaxation time.

<!-- -->

1.  Typical charge configurations at $T=0.033$. Positive charges are represented by filled squares.

<!-- -->

1.  Trajectories of moving positive charges at $T=0.033$ over a time interval of $500$ MC steps. Arrows indicate the directions of single charge motions.

<!-- -->

1.  Dynamic cooperativity for (a) scalar displacement and (b) vector displacement respectively for varying time intervals at various temperatures.

<!-- -->

1.  The structure factor $S(q)$ at $T=0.033$ and $T=0.037$.

<!-- -->

1.  The incoherent intermediate scattering functions at temperature $T = 0.033$ for various wave vectors $q$.

<!-- -->

1.  $q$-dependence of the $b$ and $\beta$ exponents for the intermediate scattering functions at temperature $T = 0.033$.

<!-- -->

1.  Comparison of the Gaussian approximations and the incoherent intermediate scattering functions at temperature $T = 0.033$ for various wave vectors $q$. We can see that the gaussian approximation is worse at large wave vectors.
