text:
# Susceptibilities and screening masses in two flavor QCD

## I Introduction

Experiments at the Brookhaven RHIC are now seeking to establish the formation of quark-gluon plasma. Proposals for experimental signatures of the plasma have to take into account its physical nature. Such information is most reliably obtained in lattice simulations of the hot phase of QCD. At present the phase transition temperature, $`T_c`$, is fairly well known tc ; tc2 , and reasonably reliable information on the equation of state has also been obtained eos . Screening and fluctuations of conserved charges are some of the other relevant properties of the plasma. We report extensive new lattice results for these.

We address two different physics questions in this paper. We are primarily interested in the question of baryon number and electric charge fluctuations in the plasma. These are quantities of direct experimental relevance, as has been realized recently koch ; chi0s . We construct them from measurements of quark number susceptibilities qsdef . In addition, we study screening masses in the plasma through the study of spatial correlators, paying special attention to their tensor structures. This gives a new relation between the susceptibilities and screening masses which are verified by our lattice measurements. In particular, we explain why non-perturbative phenomena in the latter are closely connected with deviations from perturbation theory in the former.

Temperatures just above $`T_c`$ are likely to be the most relevant region for applications in heavy ion physics and cosmology. Lattice measurements of the equation of state deviate significantly from the usual high temperature perturbative expansion. This stimulated many attempts to resum the weak coupling expansion, with varying degrees of success resum . Many of these techniques would have definite predictions for the susceptibilities that we measure. Since our measurements also deviate from the weak coupling expansion, they stand as yet another invitation to resum the weak coupling series.

While the world contains six flavors of quarks, it suffices to consider two light (u and d) and a moderately heavy (s) quark for the physics of the QCD phase transition. Starting from a quenched theory where all quark loop effects are turned off, one can envisage successive approximations by which dynamical quarks are switched on in the sequence of their masses. Such an incremental strategy is actually forced by the cost of lattice simulations of full QCD. Earlier we simulated QCD in the quenched approximation, and reported our measurements of various quark susceptibilities qsus . In the simulations we report here, we go to the next level of approximation by using dynamical u and d quarks with a bare mass of about 17 MeV. Our results indicate that at temperatures above $`T_c`$ the change due to the inclusion of light quarks is small. Unquenching the heavier strange quark may thus be correspondingly less important. Even in the cold phase below $`T_c`$, taking still lighter u and d quarks may be more crucial in future than including dynamical strange quarks.

The plan of this paper is the following— technical material on the simulations, including the lattice scales and details of the data taking procedure, is presented in the next section. The definitions and notation, and details of measurements of the screening masses and quark number susceptibilities are given in the following two sections, in that order. The discussion and summary in the final section is designed to facilitate use by those who are mainly interested in applications of our results or wish to have an overview of future directions in lattice measurements of susceptibilities.

## II The simulations and scales

### II.1 Run parameters

We have simulated QCD with two flavors of light dynamical staggered quarks ($`N_f=2`$) with the R-algorithm hybridr  on $`N_t\times N_x^2\times N_z`$ lattices with $`N_t=4`$. For investigations of quark number susceptibilities, we chose spatially symmetric lattices with $`N_x=N_z`$. On the other hand, screening correlators were easier to measure on elongated lattices with $`N_z>N_x`$. Most of our simulations consisted of generating dynamical configurations with two flavors ($`N_f=2`$) of dynamical staggered quarks with the sea quark mass, $`m`$, held fixed in physical units at $`m/T_c=0.1`$ as the temperature was varied. A small set of runs with $`m/T_c=0.075`$ was performed to check the magnitude of sea quark mass effects. The run parameters and statistics collected are shown in Table 1. We emphasize that holding the quark mass fixed is important for physics applications, whereas the older measurements fixed $`m/T`$ instead.

After a few tuning runs which reproduced the results of previous studies hybridr ; nf2 , we chose to run with parameters similar to those used in previous works with two flavors of staggered quarks, i.e., trajectory length of 1 molecular dynamics (MD) time unit, integrated in 100 steps of 0.01 MD time units each, and a conjugate gradient stopping criterion of $`10^5\sqrt{V}`$ on the modulus of the remainder vector ($`V=N_x^2N_zN_t`$). It should be noted that very similar run parameters are often used in $`N_f=4`$ Hybrid Monte Carlo (HMC) runs, where a Metropolis choice ensures that the correct weight for a configuration is obtained even for a finite number of time steps. For $`N_f=2`$ a bias-free estimate of the weight is not guaranteed. In order to estimate systematic errors from this source we also made a test run with half the step size but for the same trajectory length.

Part of the purpose of this work is a systematic comparison of screening masses with different numbers of dynamical quarks. Details of our $`N_f=4`$ simulations are available elsewhere b1 ; dimred . Similar details of runs in quenched QCD have been presented in qsus . However, we also generated a few sets of new quenched configurations. Details of the algorithm and data taking remain as in qsus . The run parameters and statistics for these quenched runs are listed in Table 2.

### II.2 Setting the scale

Relative temperature scales, i.e., $`T/T_c`$, can be deduced from previous estimates of critical couplings for various values of $`N_t`$ at the quark masses of interest to us nf2 . Thus, the temperature of a run with $`N_t=4`$ at $`\beta _c(N_t^{})`$ is $`T/T_c=N_t^{}/4`$. Since $`T_c`$ is known to be $`167\pm 17`$ MeV in physical units tc2 , this allows us to deduce the value of lattice spacing, $`a`$, in physical units using the relation $`T=1/(N_ta)`$. For those values of $`T`$ which cannot be obtained by simulations on lattices with different $`N_t`$ we use the two-loop beta-function of QCD to deduce the lattice spacing, and hence the temperature. A recent global analysis of data assures us that this is possible tc2 . The values of the inverse lattice spacings appropriate to our measurements are given in Table 1. The sea quark masses that we use are $`m=17\pm 2`$ MeV ($`m/T_c=0.1`$) and $`12\pm 1`$ MeV ($`m/T_c=0.075`$).

Another way of setting the lattice spacing is to use previous lattice measurements of the $`\rho `$ meson mass with the same quark mass and coupling rho . This gives us a lattice spacing $`1/a=580\pm 20`$ MeV at $`T_c`$, and hence the dynamical sea quark mass comes out to be $`m=14.5\pm 0.6`$ MeV ($`m/T_c=0.1`$) and $`m=10.9\pm 0.4`$ ($`m/T_c=0.075`$). This consistency between the two estimates is, of course, not unexpected, since it is well known that $`T_c/m_\rho `$ is reasonably independent of lattice spacing nf2 .

The strange quark mass can be chosen in two different ways. One is to take the strange quark mass to be 75-170 MeV pdg . This gives $`m_s/T_c=0.41`$$`1.0`$ 111In quenched QCD, where $`T_c280\pm 10`$ MeV, $`m_s/T_c=0.25`$$`0.60`$.. Alternatively, we can use the ratio $`m_s/m_d=17`$$`25`$ pdg  to get $`m_s/T_c=1.7`$$`2.5`$. Clearly, in the limit when the light quark masses are more realistic, the two procedures should give similar values. However, the heavier the strange quark mass is, the more strongly will it be subject to lattice artifacts. In view of this, we have chosen to work with $`m_s/T_c=0.75`$, corresponding to $`m_s=125`$ MeV. A realistic ratio of the strange and light quark masses will thus need future simulations with lighter u and d quark masses. Furthermore, the fact that the strange quark is reasonably heavy in lattice units indicates that the unquenching effects on it would be small, at least near $`T_c`$.

### II.3 Thermalisation and autocorrelations

Due to the fact that each configuration takes a long time to generate, a major issue in any Fermion simulation is to identify when a run has thermalised. We monitored thermalisation through measurements of plaquettes, Wilson line, the chiral condensate and the number of conjugate gradient iterations needed. Since plaquette measurements are not very noisy, these were a good estimator of the approach to equilibrium. At temperatures above $`T_c`$, the difference of purely spatial and mixed time-space plaquettes is nonzero, and it turns out to develop late during equilibration. This was one of the best criteria we observed for closeness of approach to equilibrium. However, even more efficient was a scatter plot of the stochastic estimator of the condensate ($`\overline{\psi }\psi `$) against the number of conjugate gradient iterations ($`N_{cg}`$). As shown in Figure 1, a very definite directional movement occurs away from equilibrium, but as equilibrium is approached the movement becomes a random walk within a well defined area in the plane of $`\overline{\psi }\psi `$ and $`N_{cg}`$. Monitoring this scatter plot thus offered a possibility of tuning the algorithm to achieve fast thermalisation.

We did this by starting each run from an ordered configuration with small trajectory lengths ($`t0.05`$ MD units). As the simulation proceeded, the plot of $`\overline{\psi }\psi `$ against $`N_{cg}`$ began to slow. When this happened we increased the trajectory length in small steps until it reached unity or resulted in further motion. A post-facto check of the tuning was obtained after subsequent runs verified that the system continued to perform a random walk without directed movement. This method was developed in our earlier $`N_f=4`$ HMC runs and worked well also for these $`N_f=2`$ HMD runs. Using this technique we often managed to achieve thermalisation in 20–30 units of MD time.

A second issue in the control of computer time is the matter of how often measurements of an observable can be made so that two successive measurements are effectively decorrelated. We have estimated autocorrelation times from a variety of gauge and fermion observables and found them to be small. Local operators such as plaquettes seem to be essentially decorrelated in about 2 trajectories. Even long distance observables, such as the pion screening correlator at the longest possible distance, are decorrelated in about 2–3 trajectories. In view of these results, we have taken data for screening masses on every second trajectory, and for susceptibilities every tenth trajectory. This is justified, a posteriori, by the fact that none of the screening masses is found to be very small.

## III Screening correlators and masses

### III.1 Symmetries and transfer matrices

Screening involves the transfer matrix in a spatial direction of the Euclidean finite temperature lattice. Thus, screening correlators and masses are classified according to the symmetries of this transfer matrix, Screening masses are obtained from the exponential decay of the screening correlators

$$C(z)=\frac{1}{N_x^2N_t}\underset{x,y,t}{}M(0)M^{}(x,y,z,t).$$
(1)
The operators $`M(x,y,z,t)`$ are made from combinations of quark and gluon fields which transform according to an irreducible representation (irrep) of the symmetry group, and the angular brackets denote averaging over the correct thermal distribution of fields. We shall also have occasion to use the meson susceptibilities

$$\chi =\underset{z}{}C(z)=\frac{N_z}{V}\underset{x,y,z,t}{}M(0)M^{}(x,y,z,t).$$
(2)
Detailed discussions of the isometries of the lattice and irreps of the transfer matrix can be found for gluon operators in bernd ; saumen  and staggered quark operators in d4h . Here we only mention the features used in this work.

For operators made entirely out of gluon fields, “glueball” operators in the usual shorthand, the relevant symmetry is that of a slice of the lattice. For screening correlators such a slice includes two spatial directions and the Euclidean time direction. Since rotations of the spatial directions into Euclidean time are disallowed, the full cubic group, $`O_h`$, breaks down to the dihedral group, $`D_4^h`$. For bilinear operators constructed from staggered quark fields, the “meson” operators, the symmetries of the transfer matrix are more complicated due to mixing of spin and flavor components and staggering of the quark fields. However, this transfer matrix also carries representations of $`D_4^h`$ as shown in d4h . $`D_4^h`$ has eight one-dimensional and two two-dimensional irreps.

In the continuum, the $`T=0`$ symmetry group is that of $`O(3)`$ rotations of the slice, and breaks down for $`T>0`$ to the symmetry group of the cylinder, $`𝒞=O(2)\times Z_2`$. The lattice groups are subgroups of these continuum groups. The real irreps of $`𝒞`$ are easily obtained. There are two one-dimensional irreps $`0_+`$ and $`0_{}`$, which come from the $`J_z=0`$ components of the even and odd spin representations of $`O(3)`$. There is a countable infinity of two-dimensional real irreps which span the spaces of $`J_z=\pm M`$ for any $`J>0`$. These also carry irreps of the remaining $`Z_2`$ subgroup which corresponds to reflections, $`tt`$, in Euclidean time.

In the continuum limit the scalar irrep $`A_1^+`$ and a non-trivial irrep $`A_2^{}`$ of $`D_4^h`$ both collapse into the $`0_+`$ irrep of the cylinder group; hence these screening masses must be degenerate in this limit. The $`0_{}`$ gets both the $`A_1^{}`$ and $`A_2^+`$, and hence this pair of screening masses must also become degenerate. The four remaining one-dimensional irreps $`B_{1,2}^\pm `$ collapse into the $`J_z=\pm 2`$ irrep of $`O(2)`$, and the two-dimensional irreps $`E^\pm `$ become the $`J_z=\pm 1`$ of $`O(2)`$. All these patterns of degeneracies have been verified in the glue sector of the $`SU(2)`$ pure gauge theory saumen .

These seemingly esoteric group theoretic constructions actually have very simple physical applications. It has long been appreciated that in screening phenomena, particles of different spin may mix, and that different polarization components may have different dispersion relations chin . These correspond, respectively, to mixing of equal $`J_z`$ for different $`J`$ and the distinction between different $`J_z`$ for the same $`J`$. While we discuss only physics for $`T>0`$ and zero chemical potential here, it should be noted that the essential group theory lies in the inequivalence of the spatial and temporal directions. This is also true at finite chemical potential, and hence this group theory is also relevant to such future lattice computations.

### III.2 Results

We have investigated the usual “local” staggered meson operators, i.e., those in which the quark and the anti-quark both land on the same lattice site (see Table 3). These are familiar from their $`T=0`$ characters. The scalar $`\mathrm{𝟏}^{++}`$, S, and the pseudo-scalar $`\mathrm{𝟏}^+`$, PS, both belong to the $`A_1^+`$ of $`D_4^h`$. The three components of the vector $`\mathrm{𝟑}_{}^{\prime \prime \prime \prime }{}_{}{}^{++}`$, V, and the axial-vector $`\mathrm{𝟑}_{}^{\prime \prime \prime \prime }{}_{}{}^{+}`$, AV, can be further reduced. These give two different copies of $`A_1^+`$, and also a $`B_1^+`$. All the $`A_1^+`$ irreps have been measured extensively before.

We have also made the first ever $`T>0`$ investigation of a class of “non-local” staggered mesons, one in which the quark and the anti-quark are separated by one link (see Table 3). These are made gauge invariant by including any product of links that starts on one of these sites and ends on another. The various operators lie in the $`A_2^{}`$ and $`E^{}`$ irreps of $`D_4^h`$.

Chiral symmetry restoration gives the following relations between the staggered local meson correlators—

$$C_{PS}(z)=(1)^zC_S(z)\mathrm{and}C_{AV}(z)=(1)^zC_V(z).$$
(3)
We found that these relations are satisfied for almost all the correlators to great accuracy at all $`T`$. The only small discrepancy is for the S/PS correlators at $`T=1.5T_c`$, where the relation is violated at the 68% confidence level (but not at the 95% confidence level). In order to rule out the finite MD-time step errors as the source of this oddity, we ran a second simulation at this coupling with half the time step. The results remained unchanged. Since our error estimates are stable, we performed another run with $`N_z=24`$. On this longer lattice the S and PS agreed at the 68% confidence level.

For staggered Fermions, correlators are usually fitted to a form such as

$$C(z)=A_1\mathrm{cosh}[M_1(zN_z/2)]+(1)^zA_2\mathrm{cosh}[M_2(zN_z/2)].$$
(4)
The correlator identities in (3) imply that the combinations $`C_V+C_{AV}`$ and $`C_V+(1)^zC_{AV}`$ project out the mass eigenstates in the chiral symmetric phase. We have also made single mass fits to these projections to check for stability of the fitted parameters. In addition we have extracted screening masses by constructing local masses with these projections. The local mass, $`m(z)`$, is defined by the solution of

$$\frac{\mathrm{cosh}[m(z)(z1N_z/2)]}{\mathrm{cosh}[m(z)(z+1N_z/2)]}=\frac{C(z1)}{C(z+1)},$$
(5)
where the correlation functions are estimated by a jack-knife procedure. In order to take care of covariance between various measurements, we take the number of jack-knife bins to be much smaller than the number of measurements.

In the $`N_f=4`$ theory the $`B_1`$ correlators are identically zero b1 . A statistical test of this hypothesis is made through the usual

$$\chi ^2=\underset{ij}{}(y_ih_i)\left(\sigma ^1\right)_{ij}(y_jh_j),$$
(6)
where $`y_i`$ are the measurements of the correlators at different $`z`$, $`\sigma `$ is the covariance of these measurements, and the hypothesis being tested is that $`h_i=0`$. With the data we have gathered in $`N_f=2`$ QCD, we found that

$$\chi ^2=\{\begin{array}{cc}15.8\hfill & \text{(}1.5T_c\text{)}\hfill \\ 13.2\hfill & \text{(}2.0T_c\text{)}\hfill \\ 24.1\hfill & \text{(}3.0T_c\text{)}\hfill \end{array}$$
(7)
Since $`N_z=16`$ in all these measurements, the number of degrees of freedom is 15. Thus, no alternative hypothesis can be supported by our data at the 99% confidence limit for $`T2T_c`$. For $`T=3T_c`$, where there might be a signal, a one-mass fit by the first term of eq. (4) gave a coefficient about $`3\sigma `$ away from zero, and a screening mass $`M/T=4.5\pm 0.2`$. In the quenched theory, on the other hand, the situation was more like the $`N_f=4`$ theory. At all $`T`$ the $`B_1^+`$ correlator was compatible with zero. The fits also gave similar agreement with a vanishing correlator. More detailed studies of the $`B_1`$ correlator with high statistics at several temperatures above $`2T_c`$ might be of interest in the future.

The correlation functions in the remaining irreps are non-trivial. The screening masses obtained are listed in Table 4, where we have also collected measurements made with quenched and $`N_f=4`$ dynamical staggered Fermions on $`N_t=4`$ lattices. Note that the two sets of $`A_1^+`$ correlators, the S/PS and the V/AV, give different screening masses. The V/AV screening mass, $`M_V`$, is consistent with free field theory. At all $`T`$, the $`A_2^{}`$ screening mass is close to the S/PS screening mass, $`M_S`$, but due to the large errors, is also consistent with $`M_V`$. The $`E^{}`$ screening mass is harder to determine because the correlator is very noisy; the central value is consistent with free field theory, but has very large errors.

The two parity projected local V and AV correlators are identical. This implies that not only the screening masses, but also the mixture of states excited by the operator are identical. In other words, the $`a_1`$ and $`b_1`$ screening states mix for $`T>T_c`$. Similar equalities were found in the parity projections of the $`A_2^{}`$ (and $`E^{}`$) correlators coming from the V and AV non-local meson operators, showing that the $`\rho `$ and $`\omega `$ screening states also mix. The non-local S and PS correlators similarly bear evidence for the mixing of $`\pi `$ and $`a_0`$. In addition, the equality of the $`A_2^{}`$ screening mass determined from all four $`A_2^{}`$ correlators at our disposal also argues for a mixing of $`\pi `$ and the $`J_z=0`$ component of the $`\rho /\omega `$. Thus both the phenomena expected in screening correlators are seen— mixing of states of different $`J`$ and different components of the same $`J`$ having different dispersion relations.

## IV Quark number susceptibilities

### IV.1 Definitions

We define the partition function

$$Z=𝒟U\mathrm{e}^{S_g}detM(m_u,\mu _u)detM(m_d,\mu _d)detM(m_s,\mu _s),$$
(8)
where $`S_g`$ is the gluon action of interest and $`M`$ denotes appropriate lattice Dirac operators. The chemical potentials for each flavor can be combined into the singlet, triplet and octet $`SU(3)`$ chemical potentials—

$$\mu _0=\mu _u+\mu _d+\mu _s,\mu _3=\mu _u\mu _d,\mathrm{and}\mu _8=\mu _u+\mu _d2\mu _s.$$
(9)
Our notational convention is that indices such as $`f`$ and $`f^{}`$ run over the flavors $`u`$, $`d`$ and $`s`$, and indices such as $`i`$ and $`j`$ run over the $`SU(3)`$ diagonal generators 0, 3 and 8.

The quark number densities are

$$n_f\left(\frac{T}{V_3}\right)\frac{\mathrm{ln}Z}{\mu _f}=\left(\frac{T}{V_3}\right)\mathrm{tr}M_f^1M_f^{},$$
(10)
where $`M_f^{}=M_f/\mu _f`$, $`V_3=N_x^3a^3`$, and $`T=1/N_ta`$. Conversion to the other basis simply follows using the definitions in eq. (9) and the chain rule of differentiation. Note that $`n_0=(n_u+n_d+n_s)/3`$ and $`n_3=(n_un_d)/2`$ are the baryon number and isospin densities, respectively. These are densities of (conserved) charges and not of quark number. The quark number susceptibilities are the second derivative of the free energy with respect to the chemical potentials

$$\chi _{ff^{}}\frac{n_f}{\mu _f^{}}=\left(\frac{T}{V_3}\right)\left[\frac{1}{Z}\frac{^2Z}{\mu _f\mu _f^{}}\frac{1}{Z}\frac{Z}{\mu _f}\frac{1}{Z}\frac{Z}{\mu _f^{}}\right].$$
(11)
To lighten the notation, we shall put only one subscript on the diagonal parts of $`\chi `$.

We are interested in evaluating the susceptibilities for $`m_u=m_d<m_s`$ at the point $`\mu _f=0`$ for all $`f`$, yielding much simplification. For example, each $`n_f`$ vanishes, a fact that we utilize as a check on our numerical evaluation. Moreover, the product of the single derivative terms in eq. (11) vanishes, since each is proportional to a number density. Finally, since the masses are degenerate, $`M(m_u,0)=M(m_d,0)`$ for each configuration of gauge links under study.

Flavor off-diagonal susceptibilities such as

$$\chi _{ud}=\left(\frac{T}{V_3}\right)\mathrm{tr}M_u^1M_u^{}\mathrm{tr}M_d^1M_d^{}$$
(12)
are given entirely in terms of the expectation values of disconnected loops. Such quantities are numerically hard to compute. We discuss their evaluation later. Since $`M_u=M_d`$, we obtain $`\chi _{us}=\chi _{ds}`$ with each defined by an obvious generalization of the formula above. Of the flavor diagonal susceptibilities we shall use

$$\chi _s=\left(\frac{T}{V_3}\right)\left[\left(\mathrm{tr}M_s^1M_s^{}\right)^2+\mathrm{tr}\left(M_s^1M_s^{\prime \prime }M_s^1M_s^{}M_s^1M_s^{}\right)\right].$$
(13)
$`\chi _u=\chi _d`$ are given by an obvious generalization of this formula. Numerically, the simplest quantity to evaluate is the diagonal iso-vector susceptibility

$$\chi _3=\frac{1}{2}\left(\frac{T}{V_3}\right)\mathrm{tr}\left(M_u^1M_u^{\prime \prime }M_u^1M_u^{}M_u^1M_u^{}\right).$$
(14)
Two more susceptibilities are of interest. These are the baryon number and electric charge susceptibilities,

$$\chi _0=\frac{1}{9}\left(4\chi _3+\chi _s+4\chi _{ud}+4\chi _{us}\right)\mathrm{and}\chi _q=\frac{1}{9}\left(10\chi _3+\chi _s+\chi _{ud}2\chi _{us}\right).$$
(15)
Note that $`\chi _0`$ is the baryon number susceptibility for three flavors of quarks. As a result, this expression differs from the iso-singlet quark number susceptibility for two flavors, defined in qsdef , both in overall normalization and by terms containing strangeness.

Note that quark masses appear in two places: first in the determinant in (8) which defines the weights for the averaging, and second in the trace of the Dirac operators which define the susceptibilities. The first is the dynamical sea quark mass $`m`$, and the second is the valence quark mass $`m_v`$. In principle these can be different. We work with light sea quark masses for the u and d flavors, as described earlier, and use the quenched approximation for the strange quark, i.e., set $`detM_s=1`$. We have earlier reported measurements of these susceptibilities when all the flavors are quenched qsus . For staggered quarks each trace in eqs. (1214) comes with a factor of $`1/4`$ to compensate for flavor doubling. We differ from the conventions of qsdef  by this overall factor.

These susceptibilities are easy to compute in free field theory i.e., for an ideal Fermi gas. As expected, $`n_f=0`$ for $`\mu _f=0`$. All the off-diagonal susceptibilities are also zero. For $`\chi _3`$ we obtain

$$\chi _{FFT}=\frac{1}{4N_tN_x^3}\underset{p}{}\frac{\mathrm{sin}^2p_0\mathrm{cos}^2p_0}{\left[m^2+_\nu \mathrm{sin}^2p_\nu \right]^2},$$
(16)
where the spectrum of momenta is $`p_0=(2\pi /N_t)(n_0+1/2)`$, for $`n_0=0,\mathrm{},N_t1`$, and $`p_\nu =(2\pi /N_x)n_\nu `$, where $`n_\nu =0,\mathrm{},N_x1`$ for $`\nu 0`$. $`\chi _s`$ is also given by the same expression for the appropriate quark mass. $`\chi _0^{FFT}`$ and $`\chi _q^{FFT}`$ can then be obtained using eq. (15).

### IV.2 Optimizing the measurements

The traces required in measurements of the quark number susceptibility were estimated using a well-known stochastic method. Such stochastic estimators for the traces of $`N\times N`$ matrices can be constructed from ensembles of $`N`$-dimensional complex vectors, $`R`$, of which each component, $`z_\alpha `$, is drawn from a Gaussian distribution of unit variance. From the identity

$$\frac{1}{2}z_\alpha ^{}z_\beta \underset{\mu =1}{\overset{N}{}}\mathrm{e}^{|z_\mu |^2/2}\left(\frac{d^2z_\alpha }{2\pi }\right)=\delta _{\alpha \beta },$$
(17)
(here $`d^2z`$ means $`rdrd\theta `$ and the complex number $`z=r\mathrm{exp}i\theta `$) we construct the estimator by multiplying both sides by a matrix element $`A_{\alpha \beta }`$ and summing over the indices. The right hand side gives $`\mathrm{tr}A`$. The integral on the left hand side has an obvious Monte Carlo estimator, giving the relation

$$\mathrm{tr}A=\frac{1}{2N_v}\underset{i=1}{\overset{N_v}{}}R_i^{}AR_i=\overline{R^{}AR},$$
(18)
where $`N_v`$ is the number of vectors used in making the estimate. Note that eq. (17) implies that the factor 2 above is part of the weight in the average over the ensemble of random complex vectors.

The stochastic estimator for $`(\mathrm{tr}A)^2`$ is a small modification of that given in qsdef . We draw $`L`$ sets of independent random Gaussian vectors as before and construct the estimator

$$(\mathrm{tr}A)^2=\frac{2}{L(L1)}\underset{i>j=1}{\overset{L}{}}\left(\overline{R_i^{}AR_i}\right)\left(\overline{R_j^{}AR_j}\right).$$
(19)
If we draw $`N_v`$ vectors in the evaluation of the single trace, then it is simplest to divide these into $`L`$ sets of $`N_v/L`$ vectors for the estimate of each ensemble average in the sum above 222We thank K. Kanaya for a discussion on this estimator..

In qsdef  the sum above was taken over all pairs, and the diagonal term was removed by subtracting a term equal to a stochastic estimator for $`\mathrm{tr}(A^2)`$. In qsus  where that method was used, we showed that reasonably large values, $`N_v100`$ were necessary to avoid a systematic bias in $`\chi _0\chi _3`$. However, the estimate is facilitated by excluding the diagonal term from the sum above. For example, we found that the error in $`\mathrm{tr}A`$ scales as $`1/\sqrt{N_v}`$, whereas that in $`(\mathrm{tr}A)^2`$ scales as $`1/N_v`$, for fixed $`L`$. The dependence on $`L`$ was marginal, as long $`N_v/L>5`$.

Since the simulation time increases at least linearly in volume, we need to make a choice of volume which reduces computation time without introducing large artifacts into the measurements. We have made a detailed study of the volume dependence of susceptibility measurements at $`T=1.5T_c`$. As the data in Table 5 clearly show, there is no significant volume dependence in the values of $`\chi _3`$ and $`\chi _{ud}`$ in going from $`4\times 8^3`$ to $`4\times 16^3`$ lattices.

In view of this, we have chosen to make the remaining measurements on $`4\times 12^3`$ lattices. Close to $`T_c`$ we would expect that the measurements are strongly volume dependent. Based on our estimates of the screening masses we find that even at our smallest temperature of $`1.25T_c`$ the lattice is more than 10 times larger than the longest correlation length. This assures us that finite lattice size effects are negligible at the lower end of our temperature range. On the other end of the scale, by the simple expedient of not going beyond $`3T_c`$, we avoid the finite volume effects such as the onset of spatial deconfinement.

Finally, one must optimize the conjugate gradient stopping criterion, in which the norm of the residual vector is required to be smaller than $`ϵ\sqrt{V}`$. We were pleasantly surprised to find that, in thermalized test configurations, results for $`\chi _3`$ and $`\chi _{ud}`$ changed by less than 1 part in $`10^4`$ on changing $`ϵ`$ in the range from $`0.01`$ to $`10^6`$. Increasing $`ϵ`$ by one order of magnitude meant a decrease of CPU time by 20–25%. These numbers were reproduced for test configurations on three lattice volumes, and a variety of $`T`$. In actual measurements we chose $`ϵ`$ to be $`10^5`$ for $`T1.25T_c`$ and $`10^3`$ otherwise.

### IV.3 Results

Our main results on the $`T`$ and $`m_v`$ dependence of $`\chi _3`$ and $`\chi _{ud}`$ are collected in Table 6 for two flavors of dynamical Fermions with $`m/T_c=0.1`$. The results of a simulation at $`T=1.5T_c`$ with a smaller $`m/T_c=0.075`$ are collected in Table 7. Comparing the two, we see no statistically significant change in $`\chi _3`$ over this range of masses. While $`\chi _3`$ is expected to be strongly dependent on $`m`$ for $`TT_c`$ in the chiral limit, our results indicate that the dependence of $`\chi `$ on the sea quark mass away from the critical region is small. This suggests that the values of $`\chi _3`$ or $`\chi _{ud}`$ in the chiral limit for $`T1.25T_c`$ are within errors those in Table 6. Table 8 shows that all the susceptibilities indeed vanish at $`T=0`$ as expected.

A comparison of the dynamical and quenched theories is displayed in Figure 2 in terms of the ratio $`\chi _3/\chi _{FFT}`$. While the general trend of the data are similar in the two cases, some quantitative differences are visible. The most important of these is in the asymptotic value of $`\chi /\chi _{FFT}`$. For small $`m_v`$ the ratio is 0.85 at $`T/T_c`$ = 3 in the quenched theory, whereas it is 0.88 for $`N_f=2`$. For larger $`m_v`$, $`\chi _3`$ is almost unchanged on unquenching the Fermions. $`\chi _s`$ changes by less than 3% when the light fermions are unquenched, leading us to believe that the unquenching of the strange quark will not affect this quantity significantly.

These observations are not explained in continuum perturbation theory at high temperature, which yields

$$\chi /\chi _{FFT}=12\left(\frac{\alpha _S}{\pi }\right)+8\sqrt{1+\frac{N_f}{6}}\left(\frac{\alpha _S}{\pi }\right)^{3/2},$$
(20)
when the plasmon term is resummed pert . Here $`\alpha _S`$ is the strong coupling at a scale appropriate to the temperature $`T`$. It is easy to see that this ratio is never less than $`0.981`$ for $`N_f=0`$ ($`0.986`$ for $`N_f=2`$). This minimum is reached when $`\alpha _S=\pi /6(6+N_f)`$. A recent analysis showed that $`T_c/\mathrm{\Lambda }_{\overline{MS}}=1.15\pm 0.05`$ for $`N_f=0`$ and $`0.49\pm 0.02`$ for $`N_f=2`$ tc2 . Taking the scale for the strong coupling to be $`2\pi T`$, this means that for $`N_f=2`$ the minimum occurs at $`T3300\mathrm{\Lambda }_{\overline{MS}}`$ (and $`T110\mathrm{\Lambda }_{\overline{MS}}`$ for $`N_f=0`$). Thus in both cases, the weak coupling estimates decrease as a function of temperature in the range studied here, in contrast to our lattice results, which increase. It would be interesting to check whether a full $`𝒪(\alpha _S^{5/2})`$ computation for $`\chi `$ approaches the lattice results, and whether any resummations of the perturbation theory do better resum ; toni .

The off-diagonal susceptibility, $`\chi _{ud}`$, vanishes for $`T>T_c`$ and also at $`T=0`$. As shown in Table 6, the errors on $`\chi _{ud}/T^2`$ for $`T>T_c`$ are of the order of $`10^6`$, and hence this is the precision within which this quantity can be said to vanish for all $`m_v`$. Thus, $`\chi _{ud}`$ and other non-diagonal parts of the flavor space susceptibilities can be totally neglected in constructing other susceptibilities. It is interesting to note that $`\chi _{ud}`$ is zero in an ideal gas, and perturbative contributions start at order $`\alpha _S^2`$. Taking the scale to be $`2\pi T`$, as before, and the same values of $`T_c/\mathrm{\Lambda }_{\overline{M}S}`$, we might then expect $`\chi _{ud}`$ to be non-zero at the level of $`0.04`$$`0.1`$ (for $`N_f=2`$) in the temperature range we studied. The substantially smaller, indeed vanishing, value thus accentuates the puzzle about perturbation theory 333In toni , it is also shown that $`\chi _{ud}`$ starts at order $`\alpha _S^3`$ with a small coefficient, bringing the resummed prediction closer to the data..

The only place where $`\chi _{ud}`$ is definitely non-zero is for $`0<T<T_c`$. A re-analysis of our data in the quenched theory at $`T=0.75T_c`$ qsus  showed that $`\chi _{ud}`$ differs from zero by over four standard deviations. As displayed in Figure 3, $`\chi _{ud}`$ scales inversely with the square of the pion screening mass, $`M_\pi `$, at this temperature. Since it is known that the hadron masses are independent of temperature from $`T=0`$ to some point rather close to $`T_c`$ lowt , the only possible temperature dependence of $`\chi _{ud}`$, for such $`T`$, would be in the value of the constant $`\chi _{ud}M_\pi ^2`$. Note that the spatial lattice size is such that $`M_\pi L10`$, and finite size effects are negligible. In order to ensure this on $`4\times 12^3`$ lattices, we are constrained to use reasonably heavy quarks. The pion screening mass is then rather large. In future we plan to study $`\chi _{ud}`$ for $`T<T_c`$ in more detail by pushing towards lower pion mass with larger lattices, covering a bigger range of $`T`$ and using dynamical quarks.

### IV.4 Relation to screening masses

The physics of a plasma lies in screening. However, it is not always clear how screening appears in different guises. Interestingly it can be shown that quark number susceptibility is directly related to screening, and therefore the non-perturbative physics in the two are likely to be the same. It is possible to present a deductive argument from first principles and use the data to illustrate it. Instead, we first show how the data prefers a relation between $`\chi `$ and the S/PS screening mass $`M_S`$, and then give the group theoretical argument that this is not a numerical coincidence.

The variation of $`\chi _3`$ with $`m_v`$ at fixed $`T`$ is shown in Figure 4. For large enough $`m_v`$, there is an exponential fall of $`\chi _3`$, which can be qualitatively understood as the effect of Boltzmann factors. It is less obvious that there should be a threshold $`m^{}(T)`$ below which $`\chi _3`$ is almost independent of $`m_v`$. The data indicates that $`m^{}(T)T`$. Such a threshold cannot be derived in a weak coupling expansion. Furthermore, this threshold of constancy is exactly what prevents $`\chi _3/\chi _{FFT}`$ from reaching unity as $`m_v0`$, since Figure 4 is just another representation of the data in Figure 2.

The only other observed non-perturbative effect in the quark sector of high-temperature QCD is in the screening mass in the S/PS channel, $`M_S`$. That these two effects are related is shown in the second panel of Figure 4, where we plot $`\chi _3`$ against $`M_S/T`$. A simple exponential relation between them is seen. The threshold we saw before disappears into the relation between $`m_v`$ and $`M_S`$. The remaining question is why the two should be related at all.

The answer is simple, and relies on the group theoretical classification of screening correlators and masses according to the symmetries of the $`T>0`$ transfer matrix in the spatial direction d4h . First, from eq. (14) it can be seen that $`\chi _3`$ is the thermal expectation value of a product of two quark propagators sandwiching $`\gamma _0U`$ and summed over all distances. As a result, this quantity can be related to a susceptibility constructed from a one-link separated meson screening correlator. In the $`T=0`$ notation, this turns out to be a component of the one-link V $`\mathrm{𝟑}^+`$. For $`T>0`$ this component reduces to the irrep $`A_2^{}`$ of the appropriate symmetry group $`D_4^h`$ (see Table 3). Other correlators which lie in the same irrep are the one-link separated S and PS $`\mathrm{𝟑}_{}^{\prime \prime }{}_{}{}^{\pm }`$ and the a specific component of the one-link AV $`\mathrm{𝟑}^{}`$. In the continuum limit the spectrum of $`A_2^{}`$ screening masses is degenerate with that of the $`A_1^+`$, since they come from the same irrep of the continuum $`O(2)`$ symmetry saumen . The smallest $`A_1^+`$ screening mass comes from the S/PS local propagators which are used to extract $`M_S`$. This, is the reason for the close relation between $`\chi _3`$ and $`M_S`$ shown in Figure 4.

Further numerical evidence in favor of this group theoretical argument is the similarity in the relation between $`\chi _\pi (M_S)`$ and $`\chi _3(M_S)`$, as shown in Figure 4. The near equality of the slopes lends support to our earlier observation that the $`A_2^{}`$ screening mass seems to be equal to $`M_S`$. We expect that for smaller lattice spacing or on using a Fermion action which restores the full flavor symmetry, the two screening masses should become equal. In this limit, the lines would be exactly parallel and the difference in values of $`\chi _\pi `$ and $`\chi _3`$ at the same $`M_S`$ would only reflect different operator overlaps with the same state. This, in fact, is a prediction which should be tested in future.

One further piece of information fits neatly into this group theoretic framework. $`\chi _{ud}`$ has the same symmetry as $`\chi _3`$ and therefore it is also related to the S/PS correlators at all non-zero temperatures. Since it vanishes for $`T>T_c`$, the relation is trivial there. However, for non-zero $`T<T_c`$, we have observed that $`\chi _{ud}=Z/M_\pi ^2`$. This relation would be as mysterious as that between $`\chi _3`$ and $`M_S`$ in the absence of the argument given above.

It is interesting to note that the preceding arguments rest entirely on the group theory of screening correlators and masses, i.e., on the spatial direction transfer matrix of the finite temperature system. In terms of the Euclidean time direction transfer matrix, the group theory is the same as at $`T=0`$ and there is no reason for the one-link separated vector temporal correlator to be related to the S/PS. Thus, the physics of quark number susceptibility is related to screening correlators and not to temporal correlators.

## V Discussion and Summary

### V.1 Event-to-event fluctuations

Since $`\chi _q`$ is related to charge fluctuations, it can be observed experimentally in heavy-ion collisions koch . The baryon susceptibility, $`\chi _0`$, and the strange quark susceptibility, $`\chi _s`$, can also be used to measure fluctuations of the corresponding quantities chi0s . Also, $`\chi _s`$ is related by a fluctuation-dissipation theorem to the total amount of strangeness produced in equilibrium, and hence can be used to probe departures from equilibrium in heavy-ion collisions. Given the relevance of these quantities to the experimental search for the quark-gluon plasma, we have collected our results in an easily usable form in Table 9.

The behavior of both $`\chi _q`$ and $`\chi _0`$ above $`T_c`$ are essentially controlled by $`\chi _3`$. It is clear from eq. (15), that for $`TT_c`$, where $`\chi _3\chi _s\chi _{ud}\chi _{us}`$, we must have $`\chi _q/\chi _02`$. Closer to $`T_c`$, but still in the hot phase, if $`\chi _s`$ is also neglected, we would have $`\chi _q/\chi _05/2`$. In fact, we find that this ratio decreases from $`2.13\pm 0.08`$ at $`1.25T_c`$ to $`2.02\pm 0.01`$ at $`3T_c`$. This shows that $`\chi _s`$ is significant at all temperatures in the plasma, but certainly becomes closer to $`\chi _3`$ with increasing $`T`$ (Figure 2). In any case, charge fluctuations are twice as large as baryon number fluctuations. Interestingly, we find that $`\chi _q/\chi _s`$ also falls from $`0.95\pm 0.03`$ at $`T=1.25T_c`$ to $`0.696\pm 0.004`$ at $`T=3T_c`$ (it is 2/3 in the $`T\mathrm{}`$ limit); therefore strangeness fluctuations become more significant than charge fluctuations with increasing $`T`$.

In the low-temperature phase, $`0<T<T_c`$, fluctuations are dominated by the flavor off-diagonal contribution. We have displayed evidence (Figure 3) that $`\chi _{ud}1/m_\pi ^2`$, i.e., fluctuations are due to the lightest particle of given flavor. Although the value of $`\chi _{ud}`$ is extremely small, the measure of fluctuation is the ratio of $`\chi `$ and the entropy density koch . While there are no reliable measurements of the entropy for $`T<T_c`$, it is known to be small. It is presently an open question whether the ratio is smaller for $`T<T_c`$ or on the other side of the QCD phase transition.

Assuming that the inverse relation between $`\chi _{ud}`$ and the pion mass generalizes to other flavor off-diagonal susceptibilities, eq. (15) can be used to predict that $`\chi _q/\chi _01/4`$ for $`T<T_c`$. Corrections to this number are then given in terms of $`(m_\pi /m_K)^2`$ and $`(m_\pi /m_\eta )^2`$. These might further lower the ratio by 15–20%. As a result, in the low-temperature phase we would obtain the hierarchy of fluctuations— $`\chi _0>\chi _q>\chi _s`$, in total contrast to the inverted hierarchy $`\chi _0<\chi _q<\chi _s`$ that we have measured above $`T_c`$.

### V.2 Screening and susceptibilities

In Euclidean quantum field theories at $`T=0`$, transfer matrices in all directions are isomorphic, and one, for example, can focus on the time-direction transfer matrix— $`𝐓=\mathrm{exp}(Ha)`$ (here $`H`$ is the Hamiltonian and $`a`$ the lattice spacing). The eigenvalues of $`𝐓`$ determine the hadron spectrum. The transfer matrix and, hence, its eigenvalues do not change with temperature. Specifically, the symmetry of the transfer matrix remains the rotational $`O(3)`$ symmetry which is used to classify particle states at $`T=0`$.

At finite temperature or chemical potential all directions are not equivalent. As a result, $`𝐓`$ is not the same as any spatial direction transfer matrix. These have a totally different symmetry, $`𝒞`$, that of the cylinder, as discussed extensively in the literature bernd ; saumen ; d4h . If the two lowest eigenvalues belonging to the scalar of $`𝒞`$ become degenerate, then a phase transition occurs. In general, the eigenvalues of this transfer matrix determine screening masses in equilibrium. Breaking of $`O(3)`$ to $`𝒞`$ implies strange phenomena like the ‘mixing’ of different spins, or different states of the same spin taking on independent dispersion relations chin . While these phenomena are strange, they are not new— the group theory is familiar from the textbook examples of the comparison between the hydrogen atom and the hydrogen molecular ion $`H_2^+`$.

Lattice studies of screening correlators have hovered on the verge of this phenomenology. If it has not received widespread attention in the past, that is merely due to the practical difficulties of measuring some of the non-trivial irreps of $`𝒞`$, as we show for the first time in Table 4. The equality of parity projected correlators was used in Section III.2 to point out that the $`\rho `$ and $`\omega `$ mix in screening, as do the $`a_1`$ and $`b_1`$. The equality of screening masses also gives some evidence for the mixing of $`\pi `$ with the $`J_z=0`$ component of the $`\rho /\omega `$, causing a splitting in the screening masses of different components of the latter. We argued in Section IV.4 that the observed relation between $`\chi _3`$ and the screening mass $`M_S`$ or between $`\chi _{ud}`$ and $`M_\pi `$ comes from precisely the same physics as the mixing of different spins at finite temperature. The same argument also implies that since fluctuations are observable, screening phenomena are physical.

### V.3 Future directions

Several directions for future work are clear, and have been discussed in the body of the paper. Here we collect what seems to us the most important and fruitful possibilities.

One direction for future work is to examine as many of the non-local screening correlators as possible, in order to gather further information on all the quantum numbers that screening correlators can come in d4h . As we have seen, this requires large lattices and immense statistics, and may well profit from the use of new noise reduction techniques and improved actions which do not change the symmetries (or the positivity) of the transfer matrix.

The ratio $`\chi _3/\chi _{FFT}`$ is not explained in perturbation theory. We noted this in a quenched computation qsus , and have verified it in the dynamical QCD computation here. The small difference (roughly 3%) between the quenched and dynamical results is also significant. Explanation of these results stand as invitations to those who resum the continuum high-temperature perturbation theory. Future lattice computations will need to push towards the continuum limit. Work in this direction is in progress, and will be reported soon.

As explained before, $`\chi _{ud}`$ vanishes in an ideal gas of quarks, but in an interacting theory would take on a value of order $`\alpha _S^3`$. This is several orders of magnitude larger than the largest result that our measurements can tolerate. This observation also is an invitation to perturbation theory.

The off-diagonal susceptibilities are non-vanishing only in the range $`0<T<T_c`$. This region of temperature is hard to study, since all the complications of $`T=0`$ QCD remain, and none of the simplifications of $`T>T_c`$ QCD start. On the other hand it is important input to studies of fluctuations in heavy ion collisions. We have reported a first quantitative observation here. Future work will push towards more realistic quark masses, taking the thermodynamic limit, and examining a larger range of temperatures.


format: text

text:
# Direct observation of the charging of a 2D electron gas through an incompressible strip in the quantum Hall regime

## Abstract

Using charge accumulation imaging, we measure the charge flow across an incompressible strip and follow its evolution with magnetic field. The strip runs parallel to the edge of a gate deposited on the sample and forms at positions where an exact number of integer Landau levels is filled. An RC model of charging fits the data well and enables us to determine the longitudinal resistance of the strip. Surprisingly, we find that the strip becomes more resistive as its width decreases.

The edge state picture of transport in a 2-dimensional electron gas in a strong magnetic field appears to explain several key properties of the 2DEG in the quantum Hall effect (QHE) regime. The intersection between the constant Fermi energy and the Landau-level energy near the edge of the 2DEG determines the spatial location of the edge states. In a self-consistent electrostatic picture, edge channels are separated by incompressible strips of integer Landau level filling. Resulting models involving alternating compressible and incompressible strips have received considerable theoretical and experimental attention . In particular, groups have employed various local techniques  to resolve the structure of these strips . Subsurface Charge Accumulation (SCA) imaging enables measurements of both the disorder potential and imaging of low compressibility strips in the integer QHE.

In this work, we use SCA to image and characterize the charging properties of the incompressible strips that form in the 2DEG in the presence of a density gradient. The electron gas charges by pumping charge to it across an incompressible strip. We show here that a simple RC charging model fits the data well. Using this model we determine the strip resistance $`R_{strip}`$ as a function of magnetic field and find surprisingly, that the strips become more resistive as their widths decrease.

Our samples are made from a standard GaAs/AlGaAs heterostructure with a carrier concentration of $`1.5\times 10^{11}`$ $`cm^2`$ and transport mobility of $`1.5\times 10^6`$ $`cm^2/V`$ sec. The 2DEG forms at the GaAs-AlGaAs interface 80 nm beneath the surface. We deposited a Cr gate on top of the sample, lithographically patterned in the form of a fingered grating (see schematic in Fig. 1). An ac excitation of 6 mV rms at a frequency of 100 kHz is applied to both the gate and the 2DEG. Due to its capacitance to the grounded body of the scanning microscope, the 2D electron layer charges and discharges in response to the applied signal. By measuring the charge induced on a sharp metal tip scanned $`10`$ nm above the sample surface at a temperature of 0.35 K, we obtain a charging map of the 2DEG. Both the signal in-phase with and lagging $`90^{}`$ from the ac excitation are recorded using a lock-in amplifier.

Figure 1a presents a $`12\times 12`$ $`\mu m`$ SCA image at zero magnetic field. These images result from scanning the tip over an ungated portion of the sample situated between two gate fingers that run parallel to each other. The scan window is positioned so that only one finger appears on the left. In the absence of magnetic field, the 2DEG maintains high conductivity and thus charges uniformly in phase with the applied ac excitation. As a consequence, the signal lagging the excitation (Y-phase) is zero everywhere and the corresponding image (not shown) displays no contrast. The picture changes dramatically in applied magnetic field. At 6.4 Tesla (close to filling factor $`\nu =1`$), charge flowing from the gates, now penetrates the 2DEG only partially (Fig. 1c). In addition, the Y-phase signal in Fig. 1d peaks at the boundary between charging and non charging regions.

In this work, we do not apply any external bias voltage on the gate. However, the electronic density under the metal (chromium) gate is reduced by about $`15\%`$ from its bulk value due to the difference in chemical potential between the metal and GaAs. As a result of this density profile, at high magnetic fields an incompressible strip, running parallel to each finger, forms in regions where electrons completely fill a Landau level. In these areas $`\sigma _{xx}0`$, and the strip prevents charge from flowing freely into the interior of the sample, exactly as depicted in Figure 1c. The regions of the 2DEG with very low conductivity give rise to a phase lagging SCA signal (Fig. 1d). Decreasing the strength of the applied magnetic field reduces the degeneracy of the Landau levels. The strip, as well as the boundary that separates charging from non charging regions, moves closer to the fingers where the density is lower. By following this evolution with field, we determine that the 2DEG electron density is nearly constant between the gate fingers but decreases closer to the gate. We find a density gradient of $`1.2\times 10^{10}cm^2\mu m^1`$ close to the gate. By moving our scan window laterally along the fingers of the gate, we have verified that neighboring strips may close on each other at certain locations, forming closed loops with the 2DEG inside electrically isolated from the rest of the electron layer.

The work-function difference between the Pt-Ir scanning tip and the GaAs surface results in an electric field that perturbs the electronic density of the 2DEG in the vicinity of the tip. As described elsewhere , we perform Kelvin probe measurements and apply a bias voltage to the tip that nulls this electric field. For values away from this nulling voltage, the probe acts as an effective gate on the sample. Fig. 2 shows scans of a 12x12 $`\mu m`$ region at $`6.3`$ Tesla for different tip biases. With the local density under the tip reduced, the condition for integer filling of the $`\nu =1`$ Landau level is satisfied at a region with higher intrinsic electron density. As a result, the location of the observed strip changes and shifts towards the interior of the imaged region. The resulting density suppression, though local, extends to $`1\mu m`$ and connects regions that would otherwise be separated by incompressible fluid.

The shapes of areas with high out-of-phase signal in Fig. 2d resemble sets of arcs. The centers of the arcs are located roughly along the strip position measured at nulling voltage in Fig. 2a suggesting the following non-local gating mechanism: bringing the biased tip on-top of a high density region that does not charge at nulling voltage may reduce the density and permit charging. In this case, there must exist a path between the tip location and the gate, where the filling factor is everywhere less than the nearest integer. The density perturbation induced by the tip decays sharply in the lateral direction. Crucial points along this path are situated at the strip location as measured at nulling. Indeed, in many cases the centers of the arcs such as seen in Fig. 2d, point to the location of features in the unperturbed images (see for example the feature in Fig. 2a marked by “x”). We postulate that each arc indicates tip locations where the tip influence on the “weak link” in the strip is the same . On one side of the arc, the weak link is open and the area under the tip charges fully. On the other side of the arc, the weak link is closed and the area under the tip becomes effectively disconnected from the gate (Fig. 2b).

We can counter changes in the position of the strip created with the tip bias by properly tuning the magnetic field. At $`\nu =1`$, the incompressible strips appear at the same location if we change both magnetic field and tip voltage at a rate of dV/dB=0.5 V/T. At $`\nu =2`$ we find that dV/dB=1.0 V/T creates the same effect. This demonstrates that density changes are compensated by varying the number of electrons that each level can accommodate. The voltage required for compensation should be proportional to the filling factor, in agreement with the experiment. All SCA images and measurements reported in this paper, other than those in Fig. 2 c and d, result from scans taken at nulled tip-sample electric field.

We now describe the results of SCA measurements obtained by scanning the tip along single lines shown schematically in Figure 1. The boundary between the gate finger and the bare sample occurs at $`x1.5\mu m`$ (distance of 1.5 $`\mu m`$ from the leftmost point of the line scan). The magnitude of the resulting signal step reflects the difference in SCA signal measured with the tip positioned above the metal gate ($`0<x<1.5\mu m`$) and above a fully charging region of the 2DEG ($`1.5\mu m<x<3.1\mu m`$). This observed step does not vary with magnetic field (see Fig. 3), indicating no field dependence in the amount of charge flowing in regions of the 2DEG adjacent to the gate.

At zero magnetic field as well as at fields away from the quantum Hall plateaus, the 2DEG charges nearly as a classical 2D metal. This situation changes upon tuning the field close to the integer Hall plateaus. The leftmost sequence (in-phase signal) of line scans in Fig. 3 shows the appearance and evolution of another SCA signal step at fields around 6.4 T. The decrease in signal in the region between neighboring gate fingers results from the formation of incompressible strips at $`\nu =1`$, running parallel to the gate. We point out that this region in the interior, bounded by incompressible strips on both sides, has $`\nu >1`$ and high conductivity ($`\sigma _{xx}`$).

To understand the origin of the $`90^{}`$ lagging SCA signal (curves on the right of Fig. 3), we consider an RC charging model for our 2DEG and assume that the interior of the sample charges through a strip characterized by resistance R while the tip is coupled to the 2DEG via a capacitance C (see inset in Fig. 4). As discussed below, this tip-to-sample capacitance totally dominates over the self-capacitance of the “isolated” $`\nu >1`$ higher conductivity region. The resulting in-phase signal measured at the tip, is proportional to $`1/(1+(\omega RC)^2)`$ while the out-of-phase evolves as $`\omega RC/(1+(\omega RC)^2)`$. The X-signal drops steadily to zero as the resistance of the path leading to the region under the tip increases, while the Y-signal goes through a peak near the roll-off frequency of $`1/(RC)`$. The presence of the out-of-phase signal indicates that there is not enough time to fully charge the 2DEG during one complete cycle of the excitation.

Figure 3 shows that the formation of the strip on the high field side of the Hall plateau (B=6.5 T) begins to impede charging of the interior $`\nu >1`$ region of the sample (Fig. 3, distance $`>4\mu `$). At 6.46 and 6.42 T, we observe a somewhat higher strip resistance, but the interior can still be partially charged during the ac excitation cycle, as evidenced by a nonzero phase lagging signal. With decreasing field, the in-phase and phase lagging signals in the interior $`\nu >1`$ region further drop and eventually saturate beyond B=6.38 T. Further reduction of the field simply pushes the strip closer to the edge of the gate and leaves magnitude of the signals in the interior unchanged. We conclude that the resistance across the strip steadily grows as the strip moves closer to the gate, and at B=6.38 T the strip becomes practically impenetrable to charge flowing from the gate for excitation frequency of 100 kHz. The density gradient is steeper closer to the gate and hence the incompressible strips formed there should be narrower . Therefore, the effective resistance for charging across the strip, counterintuitively, is greater for a narrower strip.

The crude RC model used so far dictates that the measured signal also depends upon the excitation frequency. Lowering the measurement frequency should enhance the charging of the 2DEG through the strip provided that R and C remain constant. We show data consistent with this expectation in the left inset of Fig. 4. A line scan at 6.52T (in a different location from those in Fig. 3), is measured at different frequencies. As expected, the 2DEG charges more fully at lower frequencies, as evidenced by the increased in-phase signal in the $`\nu >1`$ region. The lower trace in that inset corresponds to a magnetic field of 6.44 T where the strip allows essentially no charge to cross it. In contrast to line scans taken at slightly higher fields, this trace remains unchanged even at the lowest experimental frequencies.

We use the corresponding step $`S_0`$ in the in-phase signal between the fully charged and totally non-charged regions that the strip separates (bold trace in the left inset of Fig. 4), to calibrate the sensitivity of our measurement to full charging. We can therefore determine precisely the degree of charging of the 2DEG in the region between the strips at different fields for a range of frequencies and express it as a fraction $`q`$ of $`S_0`$. Figure 4 plots $`q`$ against frequency using data from the curves in the left inset. The smooth curve fits the data points using the RC model discussed earlier and works remarkably well. We have performed extensive measurements at several different locations within the sample and find consistent agreement with the RC model.

We can now extract directly the resistance, R, of the strip by assuming that the value for C is the known geometric capacitance $`C_g`$ between the tip and the 2DEG. We justify this assumption below in an examination of the effects of changing the tip-sample separation. $`C_g`$ is determined by measuring the signal change at zero magnetic field, as we move the tip vertically from ‘infinity’ to within 10 nm from the sample surface, yielding $`C_g`$ 0.5 fF. In addition, since strips form closed loops in our experiment, we are able to follow them to their full extent and hence measure exactly their total length. For the fit in Fig. 4 we obtain $`RC=2\times 10^6`$ s and deduce that at 6.52 T ($`\nu =1`$), $`R_{strip}100`$ M$`\mathrm{\Omega }`$/$`\mu m`$. This resistance increases by 1 order of magnitude to $`1000`$ M$`\mathrm{\Omega }`$/$`\mu m`$ at B=6.48 T. We find that at a field of 6.44 T, the strip resistivity rises beyond our measurement limit of $`10000`$ M$`\mathrm{\Omega }`$/$`\mu m`$. The results reproduce consistently at different locations throughout the sample.

We have also measured the strip resitances at $`\nu =`$ 2 and 4. The strip resistivity is $`10`$ M$`\mathrm{\Omega }`$/$`\mu m`$ at B=3.26 T and $`400`$ M$`\mathrm{\Omega }`$/$`\mu m`$ at 3.24 T (corresponding to $`\nu =2`$). For B=1.60 T ($`\nu =4`$), we get $`R_{strip}1`$ M$`\mathrm{\Omega }`$/$`\mu m`$. This result follows naturally, since at $`\nu =4`$ the energy gap $`\mathrm{}\omega _c`$ is smaller and the strips should be narrower . Temperature dependent measurements of the strip resistivity, could elucidate the mechanism responsible for charge penetrating the incompressible strip (i.e. electron tunneling or hopping). However, these measurements prove difficult with our apparatus.

The effects of changing tip height are depicted in the right inset of Fig. 4. The corresponding line scans are performed at 6.40 T at the same location as the data set in Fig. 3. The capacitance data points are extracted at the same frequencies for 3 different tip-sample separations, $`d=10nm`$, $`100nm`$ and $`300nm`$. Their strong dependence of charging rates on $`d`$ implies that the geometric capacitance $`C_g`$ between the probe and the sample dominates over the self-capacitance of the region of the 2DEG located to the right of the strip. We obtained similar results for lower field values up to the point where the strip becomes totally impenetrable. Indeed, $`C_{self}\epsilon _0a1fF`$, with $`a`$ being a characteristic length scale of the 2DEG region enclosed by the strip ($`10\mu m`$), is about one order of magnitude less than the value calculated for $`C_g`$ in our experiment. Consequently, $`C_g`$ seems to be the relevant parameter in our model and its known value may be used in determining the strip resistance from the product $`RC_g`$. From Fig. 4 we also note that, as expected, the 2DEG charges faster as $`d`$ increases and $`C_g`$ is reduced.

We would like to thank L.S. Levitov for helpful discussions. This work is supported by the Office of Naval Research, the National Science Foundation DMR and the Center for Materials Science and Engineering at MIT.


text:
# Lattice gas automaton approach to “Turbulent Diffusion”

## I Introduction

In an incompressible fluid undergoing two-dimensional flow, the equations describing the trajectories of the fluid elements in terms of the stream function $`\psi `$ take the Hamiltonian form

$`\dot{x}v_x`$ $`=`$ $`_y\psi (x,y;t),`$
$`\dot{y}v_y`$ $`=`$ $`_x\psi (x,y;t),`$ (1)
where $`𝐯`$ is the fluid velocity. In practice the trajectories of the fluid elements are visualized by seeding the fluid with suspended particles which act as passive tracers, whose velocity is given by Eqs.(I). If the flow were stationary, $`\psi _{ss}=\psi (x,y)`$, the system would be conservative and thus integrable, and the trajectories of the tracers would coincide with the curves of the contour plot of the stream function. However the actual flow possesses intrinsic fluctuations,

$`\psi (x,y;t)=\psi _{ss}(x,y)+\delta \psi (x,y;t),`$ (2)
which in general precludes exact analytical solution of Eqs.(I). A consequence of the presence of noise is that the tracers eventually exhibit diffusive behavior. The problem that we address here is the evaluation of the tracer dynamics in the flow (2).

## II Kolmogorov flow

The Navier-Stokes equations for an incompressible fluid subject to an external force $`𝐅`$ read

$`_t𝐯+(𝐯)𝐯`$ $`=`$ $`\rho ^1P+\nu ^2𝐯+𝐅,`$
$`𝐯`$ $`=`$ $`0,`$ (3)
where $`P`$ is the hydrostatic pressure, $`\rho `$ the mass density, and $`\nu `$ the kinematic viscosity. When the flow is two-dimensional and $`𝐅`$ has the form

$`𝐅=\mathrm{\Phi }\mathrm{𝟏}_xcos(\kappa y),\kappa =(2\pi /L),`$ (4)
where L is the extension of the system in the $`y`$-direction, Eqs.(II) describe the periodic Kolmogorov flow zaslavsky . Alternately, with $`\psi `$ defined by Eqs.(I), the Navier-Stokes equations can be recast in the form

$`(_tRe^1^2)\psi +J(\psi ,\mathrm{\Delta }\psi )=Re^1cos(\kappa y);`$ (5)
here $`Re=\psi _0/\nu `$ is the Reynolds number, where $`\psi _0`$ is defined by $`\psi _{ss}=\psi _0cos(\kappa y)`$, and $`J(\psi ,\mathrm{\Delta }\psi )`$ is the Jacobian

$`J(\psi ,\mathrm{\Delta }\psi )=_x\psi _y\mathrm{\Delta }\psi _y\psi _x\mathrm{\Delta }\psi .`$ (6)
By linear stability analysis, one can show that there exists a critical value of the Reynolds number where the system becomes unstable. This is most easily seen by considering a perturbation to the velocity field

$`v_x=v_x^{ss}+\delta v_x,v_y=v_y^{ss}+\delta v_y,`$ (7)
with

$`v_x^{ss}={\displaystyle \frac{\mathrm{\Phi }}{\nu \kappa ^2}}cos(\kappa y),v_y^{ss}=0,`$ (8)
and evaluating the perturbation evolution in response to the imposed forcing. We write the perturbation in terms of its spatial Fourier components as

$`\delta v_x(x,y;t)=e^{ıqx}\varphi (y;t),\delta v_y(x,y;t)=e^{ıqx}\phi (y;t),`$ (9)
with

$`\varphi (y;t)={\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}e^{ıky}\varphi _k(t),\phi (y;t)={\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}e^{ıky}\phi _k(t).`$ (10)
We insert (7) with (9) and (10) into the Navier-Stokes equations which are then Laplace-Fourier transformed to obtain the characteristic equation and the dispersion equation; retaining only the modes $`k=0,k=\kappa `$, and $`k=\kappa `$ (i.e. setting to zero the amplitudes of the modes with $`|k|>\kappa `$), we obtain from the characteristic equation the value of the critical Reynolds number

$`Re_c=\sqrt{2}\nu \kappa ^2{\displaystyle \frac{\kappa ^2+q^2}{(\kappa ^2q^2)^{1/2}}},`$ (11)
and from the dispersion equation the modes

$`\delta v_x(q,t)e^{s_+t},s_+`$ $`=`$ $`\nu \kappa ^2+\nu q^2ϵ,`$ (12)
$`\delta v_y(q,t)e^{s_{}t},s_{}`$ $`=`$ $`\nu q^2ϵ,`$ (13)
with $`ϵ=1(Re/Re_c)^2`$. Equation (12) shows that the mode $`\delta v_y`$ orthogonal to the direction of the forcing (along the $`x`$-axis) exhibits critical slowing down, i.e. $`s_{}0`$ when $`ReRe_c`$. At $`Re=Re_c`$, the Kolmogorov flow becomes unstable, and beyond the bifurcation ($`Re>Re_c`$) goes into a 2-dimensional $`ABC`$ type flow zaslavsky2  with closed streamlines, separatrices, and infinite trajectories between separatrices, as illustrated in Fig.1.

## III Lattice gas automaton dynamics

The flow to be investigated is produced by a lattice gas automaton lga  which we briefly describe. A lattice gas automaton can be viewed as a collection of particles residing in a discrete space, a regular $`d`$-dimensional lattice $``$, where they move at discrete time steps. Associated to each lattice node (with position denoted by $`𝐫`$) there is a finite set of channels (labeled by Latin indices $`i,j,..`$). Each of these channels corresponds to a discrete value of the velocity $`(𝐜_i)`$ that a particle positioned at the specified node may have. An exclusion principle imposes that there be a maximum of one particle per channel. The exclusion principle is important because it allows a symbolic representation of the state of the system in terms of bits, and of its dynamics in terms of operations over sets of bits, which are easily implemented on a computer. The state of the automaton at time $`t`$ is thus described by specifying the configuration on each and every node, i.e. the set of bits $`\{n_i(𝐫,t)\}_{i=1}^b`$ , $`𝐫`$, for an automaton with $`b`$ channels per node. The evolution of the automaton takes place in two stages : propagation and collision, applied sequentially at every time step. Particles are first moved according to their velocity: If channel $`i`$ at node $`𝐫`$ is occupied, the propagation step displaces that particle to channel $`i`$ of node $`𝐫+𝐜_i\mathrm{\Delta }t`$. This updating is done synchronously throughout the lattice. The second stage in the dynamics is a local collision step: As a function of the pre-collisional configuration, a new configuration is chosen by a prescription based on a set of (usually stochastic) rules. This step is crucial to determine the type of physics the automaton will exhibit at the macro- and mesoscopic levels. In particular, the collision step should preserve the quantities that are invariant under the dynamics of the model system. For example, it is possible to construct thermal automata whose rules are such that the number of particles, the momentum and the energy remain unchanged by the collision step GBL . In general, there will be several sets of collision rules consistent with the invariance of the specified constants of motion under the collision step. The choice of a particular set is then dictated by operational convenience, or by the need to explore a particular physical regime. Here we are interested in the two-dimensional hydrodynamic regime and it suffices to implement a simple set of rules governing the automaton dynamics on a triangular lattice FHP . It can be shown, starting from the microscopic equations of the automaton, that, provided the symmetry of the lattice is sufficient, the macroscopic dynamics of the automaton is consistent with the Navier-Stokes equations suarez . It is one of the virtues of the lattice gas automaton that on the sole basis of microscopic rules in accordance with local invariance and symmetry properties, its macroscopic behavior produces correct hydrodynamics.

Furthermore the lattice gas automaton exhibits two important features: (i) It possesses a large number of degrees of freedom; (ii) Its Boolean microscopic nature combined with stochastic micro-dynamics results in intrinsic spontaneous fluctuations, and it has been shown that these fluctuations capture the essentials of actual fluctuations in real fluids grosfils . Therefore the lattice gas automaton can be considered as a reservoir of excitations extending over a wide range of frequencies and wavelengths.

## IV Tracer dynamics

Tracers are implemented as suspended particles in the following way. The tracer, which we denote as a $`t`$-particle, is subjected to the cooperative effects of the fluid particles, and its dynamics results from the local dynamics of the automaton, that is from the combined effects of deterministic advection (due to the non-zero average velocity field resulting from the constraint imposed to the lattice gas) and random advection (due to the automaton intrinsic fluctuations). The $`t`$-particle undergoes displacements according to the average velocity of the fluid particles computed over a local domain of $``$ and over a number ($`\beta `$) of time steps of the automaton ($`\beta \mathrm{\Delta }t`$). In the absence of external force ($`Re=0`$, that is $`\psi _0=0`$ in the notation of Section II), the dynamics of the $`t`$-particle is governed solely by the automaton noise and the tracers undergo random motion leading to diffusive behavior over long distances and long times, i.e. $`\delta x^2(t)=\delta y^2(t)=2D_0t`$, for $`t\mathrm{\Delta }t`$, where $`D_0`$ will be referred to as the coefficient of molecular diffusion, and $`\mathrm{\Delta }t`$ is the elementary time step of the automaton. The value of $`\beta `$ can be varied to modify the effect of noise and thereby to tune the coefficient of molecular diffusion $`D_0`$ of the $`t`$-particle. So the dynamics of the tracers can be represented in a space defined on a plane parallel to the $``$-plane and on a time scale set by a time increment equal to $`\beta \mathrm{\Delta }t`$.

We are now interested in the dynamics of the $`t`$-particles in a fluid subject to an external force. Therefore we impose a bias to the velocity field of the automaton by preparing its initial state by distributing the velocities of the fluid particles on each node of the lattice so that on the average we obtain a periodic velocity profile according to Eqs.(8). With periodic boundary conditions imposed on the automaton universe, the shear triggered by the velocity bias produces a stationary flow in the form of a periodic array of vortices; at a moderately high value of the Reynolds number ($`Re/Re_c3`$) the spatial structure of the streamlines becomes analogous to the topology of the two-dimensional $`ABC`$ flow. In the absence of noise, the $`t`$-particles would follow exactly the streamlines, and their motion would be ballistic. However because intrinsic fluctuations are always present in the automaton (as they are in real fluids), the tracers dynamics can be strongly perturbed, and it is only in the averaged automaton flow that their trajectories reflect the topology of the $`ABC`$ flow, as illustrated in Fig.2.

When we consider the large-scale, long-time limit ($`tL_{}/V_{}`$, where $`L_{}`$ and $`V_{}`$ are the vortex characteristic quantities), we can adopt a Fokker-Planck formulation for the $`t`$-particle distribution function $`f`$, where the average is taken over the ensemble of realizations and, by averaging over a sufficiently large region of space (homogenization hypothesis homo ), we obtain a diffusion-type equation details

$`_t\stackrel{~}{f}=D_{\alpha \beta }^{}_\alpha _\beta \stackrel{~}{f},`$ (14)
where $`𝐃^{}=𝐅(V_{},L_{},D_0)`$ is the effective diffusion coefficient of the $`t`$-particles in the flow. On the basis of a power law assumption, we infer by dimensional analysis that crisanti

$`[D^{}]=[D_0]^\mu [L_{}]^{1\mu }[V_{}]^{1\mu },`$ (15)
indicating that for $`|\mu |<1`$, and $`D_0`$ sufficiently small ($`D_0<L_{}V_{}`$, which is the case in the lattice gas automaton), the elements $`D_{\alpha \alpha }^{}`$ ($`D_{\alpha \beta }^{}=D_{\alpha \alpha }^{}\delta _{\alpha \beta }`$) should be larger than $`D_0`$, with an important quantitative difference between $`D_{\alpha \alpha }^{}`$ with $`\mu >0`$, and $`D_{\beta \beta }^{}`$ with $`\mu <0`$.

We performed lattice gas automaton simulations as described above, and from measurements of the mean-square displacements of the $`t`$-particles $`\delta x^2(t)`$ and $`\delta y^2(t)`$, we find that over sufficiently long times they exhibit diffusive behavior, as shown in Fig.3. Notice that in the absence of external force, the mean-square displacement of the $`t`$-particle (due to the sole effect of the automaton noise) is isotropic and corresponds to molecular diffusion (full line in Fig.3).

The diffusion coefficients

$`D_{}^{}D_{xx}`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{\delta x^2(t)}{2t}},`$
$`D_{}^{}D_{yy}`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{\delta y^2(t)}{2t}}.`$ (16)
can be evaluated quantitatively from the explicit expression of the flow and assuming a Gaussian form for the space and time dependence of the second moment of the noise, which is compatible with the spectrum of the automaton fluctuations. We obtain details

$`D_{}^{}`$ $`=`$ $`{\displaystyle \frac{D_0}{1+\chi }},\chi ={\displaystyle \frac{\overline{v}_x}{(\delta v_x^2)^{\frac{1}{2}}}},`$
$`D_{}^{}`$ $`=`$ $`{\displaystyle \frac{\overline{v}_{x}^{}{}_{}{}^{2}}{2\kappa }}{\displaystyle \frac{1}{D_{}^{}}}{\displaystyle \frac{1+\frac{1}{2}\chi }{1+\chi ^2}}+D_{}^{},`$ (17)
where $`\delta v_x^2`$ denotes the automaton fluctuations, and $`\overline{v}_x`$ the amplitude of the forcing in the velocity field acting on the $`t`$-particle: $`\dot{x}(t)=\delta v_x(t)+\overline{v}_xcos(\kappa y)`$. In the limit $`\overline{v}_x=0`$, one has $`D_{}^{}=D_{}^{}=D_0`$, and for $`\overline{v}_x0`$, it follows from (IV), that $`D_{}^{}>D_{}^{}>D_0`$. We find qualititative agreement between these predictions and our simulation data (see Fig.3).

## V Concluding comments

We have presented an automaton approach to the problem of “turbulent diffusion” in a time-dependent flow with non-trivial average. We have shown the analogy between the stream lines in the averaged flow of the automaton and the phase plane trajectories of the corresponding dynamical system. For the full flow, we obtain agreement between the theoretical value of the diffusion coefficient of the tracers and the corresponding value computed from the lattice gas simulation data. The most important results are: (i) above the critical value of the Reynolds number ($`Re3Re_c`$), the tracers dynamics remains diffusive, and (ii) “super-diffusion” is observed in the direction orthogonal to the direction of the forcing ($`D_{}^{}D_0`$). Higher Reynolds number regimes and flows with non-stationary average are now being investigated.

###### Acknowledgements.

JPB acknowledges support by the Fonds National de la Recherche Scientifique (FNRS, Belgium). EVE was supported by the Association Euratom – Etat Belge (Belgium). DH has benefited from a grant from the Fonds pour la Formation à la Recherche dans l’Industrie et l’Agriculture (FRIA, Belgium).


text:
# Nonequilibrium Phase Transitions in Epidemics and Sandpiles

## I Introduction

Phase transitions arise not only in thermodynamic equilibrium, but in nonequilibrium situations such as the onset of convection , and in models of populations, epidemics, catalysis, cooperative transport, and markets, to cite but a few examples. The simplest models of spatially extended, nonequilibrium systems are lattice Markov processes, or interacting particle systems. There is usually no Hamiltonian; rather, the model is defined by a set of transition rates, that do not satisfy detailed balance with respect to a ‘reasonable’ energy function.

Despite these differences from equilibrium statistical mechanics, phase transitions (which mark a singular dependence of macroscopic properties on control parameters), scaling, and universality, retain their fundamental importance. Such concepts transcend the equilibrium/nonequilibrium boundary. In fact, familiar universality classes, such Ising or Potts, still appear in the nonequilibrium context, along with others not encountered in equilibrium. While new classes are still being identified, and the general principles determing critical behavior far from equilibrium have yet to be specified completely, it has become clear that symmetry and conservation laws play a fundamental role, just as in equilibrium. Some useful reviews of various aspects of nonequilibrium phase transitions can be found in Refs. . Here I discuss some recent work on absorbing-state phase transitions and their relation to self-organized criticality.

## II Absorbing-state phase transitions

Absorbing-state phase transitions arise from a conflict between the spread of activity, and a tendency for this activity to die out . The activity may represent an epidemic , a reaction proceeding on the surface of a catalyst , or spatiotemporal chaos , for example. When the system is translation-invariant, with interactions of finite range, and there are no conserved quantities, the critical behavior is expected to fall generically in the directed-percolation (DP) universality class , as has been verified in studies of a wide variety of systems.

DP and the closely related contact process  exhibit well known scaling properties, for example the stationary order parameter (activity density) grows $`(\lambda \lambda _c)^\beta `$, where $`\lambda `$ is the spreading or reproduction rate for activity. At the critical point, $`\lambda _c`$, starting from a localized active region (e.g., a single active site), the survival probability $`P(t)`$ decays $`t^\delta `$, while the mean number of active sites $`n(t)`$ grows $`t^\eta `$; these scaling laws describe avalanches of activity . For the DP class, a field theory, framed in terms of a nonnegative, scalar order parameter density $`\rho (x,t)`$ yields an upper critical dimension of four, and an $`ϵ`$-expansion for the critical exponents :

$$\frac{\rho }{t}=^2\rho a\rho b\rho ^2+\eta (x,t).$$
(1)
Here $`\eta (x,t)`$ is zero-mean Gaussian noise with autocorrelation

$$\eta (x,t)\eta (x^{},t^{})=\mathrm{\Gamma }\rho (x,t)\delta (xx^{})\delta (tt^{}).$$
(2)
Some other universality classes associated with absorbing-state phase transitions are: compact DP, in which activity is bounded by random walks ; the parity-conserving or directed-Ising class, exemplified by branching-annihilating random walks with an even number of offspring ; and the “conserved nondiffusing field” class . The criteria for a model to belong to the parity-conserving universality class are still not completely clear.

### A Interface representation

A recent development in this area is the application of surface-growth scaling ideas. Consider the contact process, in which, at a given instant, each site of the lattice is either occupied ($`\sigma _i=1`$), or vacant ($`\sigma _i=0`$). Transitions from $`\sigma _i=1`$ to $`\sigma _i=0`$ occur at a rate of unity, independent of the neighboring sites. The reverse can only occur if at least one neighbor is occupied: the transition from $`\sigma _i=0`$ to $`\sigma _i=1`$ occurs at a rate of $`\lambda m`$, where $`m`$ is the fraction of nearest neighbors of site $`i`$ that are occupied; thus the state $`\sigma _i=0`$ for all $`i`$ is absorbing. The stationary order parameter (the activity density, or fraction of occupied sites), is zero for $`\lambda <\lambda _c`$.

We define a set of height variables

$$h_i(t)=_0^t\rho _i(t^{})𝑑t^{}$$
(3)
that represent the cumulative activity at each site. Since the $`h_i`$ are nondecreasing functions of time, we can think of them as describing a a growing or driven interface, and introduce the width $`W`$ via $`W^2(t,L)=\text{var}[h_i(t)]`$, on a lattice of $`L^d`$ sites. The width exhibits Family-Vicsek scaling , but with the important difference that here it saturates only when the process falls into the absorbing state.

At the critical point, the scaling exponents for the width are related to those of DP. For example, the growth exponent (defined via $`Wt^{\beta _W}`$) satisfies $`\beta _W+\delta =1`$, where $`\delta `$ describes the decay of activity at the critical point, starting with all sites active: $`\rho t^\delta `$. This relation follows from a scaling hypothesis for the probability distribution $`P(h;t)`$ of height $`h`$ (at any site):

$$P(h;t)=\frac{1}{\overline{h}(t)}𝒫[h/\overline{h}(t)],$$
(4)
($`𝒫`$ is a scaling function, $`\overline{h}(t)`$ the mean height), that has been verified in simulations . Studies of the Domany-Kinzel model show that the growth exponent exhibits a sharp cusp at the critical point, making it a sensitive indicator of criticality .

## III How sandpiles work

Self-organized criticality (SOC) is usually taken to mean “spontaneous” criticality, or scale-invariance in the absence of identifiable tuning of parameters . In sandpiles, SOC is an absorbing-state phase transition of the same kind as in DP (though with different critical exponents), subject to an external control mechanism . A simple example illustrates how this operates.

Consider a system of activated random walkers (ARW) on a lattice with periodic boundaries in $`d1`$ dimensions. $`N`$ walkers are initially distributed at random over the $`L^d`$ sites, with no restriction on the number of walkers per site. The dynamics, which conserves the number of walkers, is simply that each site with two or more walkers has a rate of 1 to liberate a pair of walkers. Each liberated walker moves independently, without bias, to one of the neighboring sites (in one dimension, from site $`j`$ to $`j+1`$ or $`j1`$ with equal likelihood). Isolated walkers are however paralyzed. The system has two kinds of configurations: active, in which at least one site has $`2`$ walkers, and absorbing, with all walkers immobile. For $`N>L^d`$ only active configurations are possible, and activity continues forever.

For $`NL^d`$ there are both active and absorbing configurations, the latter representing a shrinking fraction of configuration space as the density $`\zeta N/L^d1`$. Given an active initial configuration (a virtual certainty for a random distribution with $`\zeta >0`$ and $`L`$ large), will the system remain active indefinitely, or will it fall into an absorbing configuration? For small $`\zeta `$ it should be easy for the latter to occur, but it seems reasonable that for sufficiently large densities (but $`<1`$), the likelihood of reaching an absorbing configuration becomes so small that the system remains active indefinitely. Mean-field theory and extensive simulations confirm the existence of sustained activity for densities greater than some critical value $`\zeta _c`$, with $`\zeta _c<1`$ . The former predicts $`\zeta _c=1/2`$ while simulation yields 0.9489 for the critical density in one dimension, 0.7170 in two dimensions . Thus the ARW model exhibits a continuous absorbing-state phase transition at a critical density $`\zeta _c`$. The critical exponents are not those of DP, due to a conserved density (i.e., of walkers), which relaxes diffusively in the presence of activity, but is frozen in inactive regions .

Now we make two simple changes in the model:

i) We open the boundaries; walkers can now jump out from the edge.

ii) To compensate for the loss of walkers at the boundaries, we include a source, which inserts a new walker at a random site, if and only if the system has fallen into an absorbing configuration. (Note that the source enjoys global information regarding the state of the system.)

It is easy to see that these innovations amount to a control mechanism that forces the system to its critical point. Suppose that $`\zeta `$ is above the critical value; then there is activity (hence no addition), and loss of walkers when activity reaches the boundary, so that $`d\zeta /dt<0`$. If $`\zeta <\zeta _c`$, on the other hand, there is no activity (and no loss), and the source is actuated, making $`d\zeta /dt>0`$. Thus the system is forced to the critical point, and, quite naturally, exhibits scale invariance in the stationary state. Thus modified, the activated random walkers model is a continuous-time version of Manna’s stochastic sandpile ; the loss and addition mechanisms yield what is commonly called self-organized criticality. Their functioning depends on the conservation of walkers, allowing the control mechanism to be described in terms of insertion and loss of “sand”, without ever making reference to the underlying control parameter $`\zeta `$.

Studies of criticality in sandpiles with a periodic boundaries (so-called fixed-energy sandpiles), have yielded their critical exponents and associated avalanche exponents , and have suggested a field theory for sandpiles . The theory involves two fields, the activity density $`\rho _a`$ and the particle density $`\zeta `$. Upon formal elimination of the latter, one obtains

$`{\displaystyle \frac{\rho _a(𝐱,t)}{t}}`$ $`=`$ $`D_a^2\rho _a(𝐱,t)r(𝐱)\rho _a(𝐱,t)b\rho _a^2(𝐱,t)`$ (5)
$`+`$ $`w\rho _a(𝐱,t){\displaystyle _0^t}𝑑t^{}^2\rho _a(𝐱,t^{})+\sqrt{\rho _a}\eta (𝐱,t).`$ (6)
$`\eta `$ is a Gaussian white noise with autocorrelation $`\eta (𝐱,t)\eta (𝐱^{},t^{})=D\delta (𝐱𝐱^{})\delta (tt^{})`$; $`c,b`$ and $`w`$ are fixed parameters. The theory resembles that for DP, except for the Laplacian memory term, and the spatial dependence of the growth rate $`r`$, inherited from the initial configuration. This theory has yet to yield to renormalization group analysis .

Simulations support the proposal  that stochastic sandpiles fall in the same universality class as that of depinning of a linear elastic interface subject to random pinning forces, at least for $`d2`$. The deterministic Bak-Tang-Wiesenfeld sandpile , by contrast, appears to define a universality class sui generis, and is marked by strong nonergodic effects.

We close this section with some observations on SOC itself. SOC was proposed as a paradigm for explaining the appearance of scale invariance in nature without the need to invoke an outside agent tuning the system to criticality . We have seen, however, that the source \[component ii) of the control scheme\], will require a huge quantity of information ($`L^d`$ bits) in order to decide whether or not to add sand, so that the riddle of scale invariance in nature seems not to have been resolved. A way out of this dilemma is to suppose that the source is does not have access to all this information, but acts on a time scale much longer than the local dynamics. Then we can expect to observe scale invariance over a considerable range of avalanche sizes and durations, large enough, in principle, to agree with observation. A similar consideration applies to strict conservation: admitting a small dissipation rate, certain sandpile models exhibit “quasi-criticality” which may again be sufficient to account for observed scale invariance .

## IV Infinite numbers of absrobing configurations

Infinite numbers of absorbing configurations arise in models of catalysis , and, as noted above, in sandpiles. The simplest model with an infinite number of absorbing configurations is Jensen’s pair contact process (PCP) . In this model, a pair of particles annihilates with probability $`p`$, and creates a new particle (at a vacant neighbor site) with probability $`1p`$. Any configuration devoid of nearest-neighbor pairs is absorbing.

The static critical behavior of the PCP falls in the directed percolation universality class, but in one dimension the spreading exponents $`\delta `$ and $`\eta `$ vary continuously with the particle density $`\varphi `$ in the environment . The spreading exponents take their DP values (e.g., $`\delta =0.16`$), only when the initial density $`\varphi `$ of isolated particles is set to the “natural” value $`\varphi _{nat}0.242`$, i.e., the value generated by the process itself, running at the critical point. (For $`\varphi =0`$, for example, $`\delta 0.27`$.) Numerically, $`p_c`$, $`\delta +\eta `$, and $`z`$ (the exponent governing the spread of the active region: $`R^2t^z`$), appear to be independent of $`\varphi `$. Exponents $`\delta `$, $`\eta `$, and $`z`$ satisfy a generalized hyperscaling relation .

The issue of variable spreading exponents remains unresolved and controversial. While evidence for variable exponents has been found in many simulations , there is as yet no complete theory. A field theory  again involves the order parameter $`\rho `$ coupled to a second field (which may represent the density of isolated particles), that is frozen in regions with $`\rho =0`$:

$$\frac{\rho }{t}=D_\rho ^2\rho a_\rho \rho b_\rho \rho ^2+w_\rho \varphi \rho +\eta _\rho $$
(7)
$$\frac{\varphi }{t}=D_\varphi ^2\rho a_\varphi \rho b_\varphi \rho ^2+w_\varphi \varphi \rho +\eta _\varphi $$
(8)
(Here both noise autocorrelations are $`\rho `$.) Eliminating $`\varphi `$, one obtains an equation for $`\rho `$ with a memory term; the stationary properties have been shown to be the same as DP. While it has not been possible to determinine the spreading exponents analytically, numerical integration of the field theory yields $`\delta `$ and $`\eta `$ that vary with the coupling $`w_\rho `$ between the two fields .

The difficulty of understanding nonuniversal spreading in models with an infinite number of absorbing configurations has motivated the study of simplified descriptions. Noting that anomalous spreading behavior is associated with a long memory of the initial configuration, Grassberger et al. proposed a generalized epidemic process, in which DP (with a unique absorbing configuration), is modified to have long memory . Specifically, the probability $`p^{}`$ for a virgin site to become infected is different from that ($`p`$) for used sites (those that have already been infected and have recovered). These authors find that in one dimension, at the critical point $`p_c`$ (which is independent of $`p^{}`$), the model exhibits variable exponents $`\delta `$ and $`\eta `$ for $`p^{}>p`$, but faster than power-law decay for $`p^{}<p`$, that is, no scaling.

Since the very existence of scaling in a model with a long memory of the initial condition appears to be in doubt, it is helpful to have some simple examples in which a variable exponent can be established analytically. A continuously variable survival probability exponent, $`\delta `$, has been demonstrated for DP confined to a parabola , and for compact DP in a similar geometry . In these cases the space-time boundary is fixed, but recently a similar result was shown for a random walk with a long memory in the form of a movable partial reflector .

Consider an unbiased, discrete-time random walk $`x_t`$ on a one-dimensional lattice, with $`x_0`$ = 1, and the origin absorbing. On visiting a virgin site, the walker is pushed back to its previous position with probability $`r`$ (analogous to $`p^{}<p`$ in the generalized epidemic process). One may think of this as being due to a reflector that sits just beyond the maximal site that has yet been visited. (Initially the reflector is at $`x=2`$.) Asymptotic analysis of the master equation, and exact iteration of transition probabilities for finite times, yield a variable survival exponent $`\delta `$. One may distinguish two kinds of reflector. A soft reflector moves forward one step at every encounter with the walker, even if the latter is reflected back. In this case we find

$$\delta =\frac{1+r}{2}.$$
(9)
A hard reflector, on the other hand, moves only if the walker succeeds in occupying the virgin site, analogous to the spread of activity in the generalized epidemic process. Here the decay exponent can be much larger:

$$\delta =\frac{1}{2(1r)}.$$
(10)
An interesting feature of the hard-reflector case is that for $`r`$ close to 1, decay of $`P(t)`$ looks faster than power-law over many decades: the asymptotic power law only becomes evident at very long times. Variable exponents $`\delta `$ and $`\eta `$ are also found in the case of compact DP with movable reflectors . If and how these findings apply to the more difficult question of spreading in models with an infinite number of absorbing configurations, and in the generalized epidemic process, remains to be seen.

## V Summary

Absorbing-state phase transitions are attracting increasing attention as connections to SOC, spatiotemporal chaos, and interface dynamics are uncovered. At the same time, various puzzles regarding universality and the nature of scaling in models with an infinite number of absorbing configurations remain unsolved, and can be expected stimulate analysis and the search for improved theoretical and computational methods.

Acknowledgments

Much of the work described here grows out of ongoing collaborations with Mikko Alava, Dani ben-Avraham, Miguel-Angel Muñoz, Alessandro Vespignani, and Stefano Zapperi. I have also benefitted from discussions with Márcio Argollo de Menezes, Hugues Chaté, Deepak Dhar, Peter Grassberger, José Guilherme Moreira, Géza Ódor, and Carmen P. C. do Prado. This work was supported by CNPq, Brazil.