text:
# ADP-01-58/T490 Cooling for instantons and the Wrath of Nahm

## 1 INSTANTONS

In 1975 Belavin et al.  determined that the vacuum of QCD contained non-trivial local minima of the Euclidean action,

$$S_E=\frac{1}{2g^2}d^4x\mathrm{Tr}F_{\mu \nu }F_{\mu \nu }.$$
(1)
Since this action is required to be finite in the case of semi-classical solutions (i.e. solutions obtained by minimizing the action) the definition of the field-strength tensor

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu +[A_\mu ,A_\nu ]$$
(2)
implies that

$$A_\mu U^1_\mu U$$
(3)
as $`|x|`$ goes to infinity. The gauge transformation $`U`$ is interpreted as a mapping from $`S^3SU(2)`$, that is, a mapping from the sphere at spacial infinity to the gauge group $`SU(2)`$ of the transformation $`U`$ ($`SU(2)`$ results can be readily extended to $`SU(3)`$). The mapping $`U`$ can then be shown to have a winding number given by the equation

$$n=\frac{1}{16\pi ^2}d^4x\mathrm{Tr}F_{\mu \nu }\stackrel{~}{F}_{\mu \nu }$$
(4)
and hence we can associate a winding number with each vacuum state. Instantons are interpreted as tunnelling amplitudes between vacuum states with different winding numbers, and their topological charge is given by the difference between the winding numbers of the states that they connect. Further important properties of (anti-)instantons are spherical symmetry in the four dimensions of Euclidean spacetime, topological charge of $`\pm 1`$, and an associated action of $`S_0=8\pi ^2/g^2`$ which is independent of the radius of the instanton.

## 2 CONNECTING THE LATTICE TO THE CONTINUUM

The lattice is a 4-dimensional approximation to Euclidean spacetime. This means that physics on the lattice should ideally approach physics in the continuum if we can take the limit of infinite lattice volume and infinitessimal lattice spacing ($`V\mathrm{}`$ and $`a0`$). Unfortunately reducing the lattice spacing indefinitely is not practical, nor is increasing the total lattice volume, since these approaches to the continuum limit rapidly increase the time taken to perform computations.
We have chosen to remove the non-physical ‘edge’ of the lattice (simulating an infinite volume) by defining periodic untwisted boundary conditions on our lattices. We are therefore performing numerical QCD on a 4-toroidal mesh of points. The discretization errors introduced by the finite lattice spacing are reduced algebraically by Symanzic improvement  and non-classical errors arising from the self-couplings of the gluon fields are dealt with by the use of mean-field improvement .
Since the lattice we use is topologically equivalent to a 4-torus, we may consider the consequences of the Nahm transform  for the stability of the configurations we study. The Nahm transform is a duality mapping which interchanges between an $`SU(N)`$ configuration with topological charge $`Q`$ on the torus and an $`SU(Q)`$ configuration with topological charge $`N`$ on the dual torus. Since there are no instanton solutions in $`U(1)`$, the Nahm transform implies  that there can be no self-dual $`|Q|=1`$ configurations on the torus.

## 3 THE LATTICE ACTION AND TOPOLOGICAL CHARGE

### 3.1 The Wilson Action

The Yang-Mills action on the lattice  is

$$S_{\mathrm{Wil}}=\beta \underset{x}{}\underset{\mu <\nu }{}\left[1\frac{1}{N}\left(Re\mathrm{Tr}W_{\mu \nu }^{(1\times 1)}(x)\right)\right]$$
(5)
where $`W_{\mu \nu }^{(m\times n)}`$ is the $`m\times n`$ Wilson loop in the $`\mu \nu `$ plane. When we expand the plaquette $`W_{\mu \nu }^{(1\times 1)}`$ around the point $`x`$ and set $`\beta =6/g^2`$ it can be easily shown that (5) approaches the continuum Yang-Mills action, plus $`𝒪(a^2)`$ corrections

$$S_{\mathrm{Wil}}\frac{1}{2}d^4x\mathrm{Tr}F_{\mu \nu }F_{\mu \nu }+𝒪(a^2)$$
(6)
in the continuum limit $`a0`$.

### 3.2 Improving the Action

We can eliminate $`𝒪(a^2)`$ errors by constructing the Wilson action from a linear combination of the plaquette and the average of the $`1\times 2`$ and $`2\times 1`$ rectangular Wilson loops. We may of course choose to employ other loops to eliminate higher-order error terms. DeForcrand et al.  have previously constructed an action free from $`𝒪(a^2)`$ and $`𝒪(a^4)`$ errors by using the five planar Wilson loops

* $`L^{(1,1)}=W_{\mu \nu }^{(1\times 1)}`$
* $`L^{(2,2)}=W_{\mu \nu }^{(2\times 2)}`$
* $`L^{(1,2)}=\frac{1}{2}\left(W_{\mu \nu }^{(1\times 2)}+W_{\mu \nu }^{(2\times 1)}\right)`$
* $`L^{(1,3)}=\frac{1}{2}\left(W_{\mu \nu }^{(1\times 3)}+W_{\mu \nu }^{(3\times 1)}\right)`$.
* $`L^{(3,3)}=W_{\mu \nu }^{(3\times 3)}`$

A general improved action can then be written in the form

$`S_{\mathrm{Imp}}`$ $`=`$ $`c_1S(L^{(1,1)})+c_2S(L^{(2,2)})+c_3S(L^{(1,2)})`$ (7)
$`c_4S(L^{(1,3)})+c_5S(L^{(3,3)}).`$
where the $`c_1,\mathrm{},c_5`$ are improvement factors which take the values

$`c_1`$ $`=`$ $`(1955c_5)/9`$
$`c_2`$ $`=`$ $`(164c_5)/9`$
$`c_3`$ $`=`$ $`(640c_564)/45`$
$`c_4`$ $`=`$ $`1/52c_5`$
and $`c_5`$ is a free parameter which we can use to “tune” the action. By setting $`c_5=0`$ we create a 4-loop improved action, denoted 4LIS for short, and by setting $`c_5=1/10`$ we eliminate the contribution of $`S(L^{(1,2)})`$ and $`S(L^{(1,3)})`$, thereby creating a 3-loop improved action, denoted 3LIS for short. For the purposes of investigating a 5LIS we choose to follow DeForcrand et al. and set $`c_5=1/20`$, midway between the 3LIS and 4LIS values. It is important to note that the 3LIS and 4LIS are just special cases of the general 5LIS, and all choices of $`c_5`$ are free from discretization errors up to and including $`𝒪(a^4)`$, however their $`𝒪(a^6)`$ errors will be different, leading to slightly different results.

### 3.3 Topological Charge

On the lattice the topological charge is calculated as the sum of the local topological charge density over all lattice sites,

$$Q=\underset{x}{}q(x)=\frac{1}{32\pi ^2}\underset{x}{}ϵ_{\mu \nu \rho \sigma }\mathrm{Tr}F_{\mu \nu }F_{\rho \sigma }.$$
(8)
If we have a dilute instanton gas with $`n_I`$ instantons and $`n_A`$ anti-instantons the total topological charge of the configuration will be

$$Q=n_In_A.$$
(9)
We hence desire the topological charge calculated on the lattice to approach integer values as the configuration being studied is cooled towards a dilute instanton gas condition. Unfortunately discretization errors tend to lead to non-integer values of $`Q`$. Fortunately, however, we can improve the topological charge.

### 3.4 The Improved Field-Strength Tensor

To improve the topological charge we have chosen to improve the field-strength tensor $`F_{\mu \nu }`$ directly, by analogy with the improvement of the action, and substitute this into Eq. (8), the definition of the topological charge.
Consider the full expansion of the plaquette

$`W_{\mu \nu }^{(1\times 1)}`$ $`=`$ $`e^{ig{\scriptscriptstyle A𝑑x}}`$ (10)
$`=`$ $`1+ig{\displaystyle A𝑑x}{\displaystyle \frac{g^2}{2}}({\displaystyle A𝑑x})^2`$
$`+𝒪(g^3)`$
$`=`$ $`1+ig\left[a^2F_{\mu \nu }+𝒪(a^4)\right]{\displaystyle \frac{g^2a^4}{2}}F_{\mu \nu }^2`$
$`+𝒪(a^6,g^3).`$
We wish to extract the second term (within the square parentheses) on the last line. This may be achieved by making the following construction;

$`W_{\mu \nu }^{(1\times 1)}`$ $`=`$ $`1+ig{\displaystyle A𝑑x}{\displaystyle \frac{g^2}{2}}({\displaystyle A𝑑x})^2`$
$`+𝒪(g^3),`$
$`W_{\mu \nu }^{(1\times 1)}`$ $`=`$ $`1ig{\displaystyle A𝑑x}{\displaystyle \frac{g^2}{2}}({\displaystyle A𝑑x})^2`$ (11)
$`+𝒪(g^3).`$
Writing

$$𝒜=W_{\mu \nu }^{(1\times 1)}W_{\mu \nu }^{(1\times 1)}$$
(12)
we may obtain a term proportional to the field strength

$`{\displaystyle \frac{i}{2}}\left(𝒜{\displaystyle \frac{1}{3}}\mathrm{Tr}𝒜\right)`$ $`=`$ $`g{\displaystyle A𝑑x}+𝒪(g^3)`$
$`=`$ $`ga^2F_{\mu \nu }+𝒪(ga^4).`$
Notice that in order to enforce the tracelessness of the Gell-Mann matrices we have subtracted one-third of the trace of $`𝒜`$.
To construct an improved field strength tensor we utilise the same five planar Wilson loops that are used in the construction of the improved action, but in this case we calculate $`W_{\mu \nu }^{(m\times n)}`$ from the clover average of four $`m\times n`$ Wilson loops in the $`\mu \nu `$ plane. The improved field strength tensor can then be written

$`F_{\mu \nu (\mathrm{Imp})}`$ $`=`$ $`k_1F_{\mu \nu }(L^{(1,1)})+k_2F_{\mu \nu }(L^{(2,2)})`$ (14)
$`+k_3F_{\mu \nu }(L^{(1,2)})+k_4F_{\mu \nu }(L^{(1,3)})`$
$`+k_5F_{\mu \nu }(L^{(3,3)})`$
where the $`k_1,\mathrm{},k_5`$ are improvement factors which take the values

$`k_1`$ $`=`$ $`19/955k_5`$
$`k_2`$ $`=`$ $`1/3616k_5`$
$`k_3`$ $`=`$ $`64k_532/45`$
$`k_4`$ $`=`$ $`1/156k_5`$
and $`k_5`$ is a free parameter, which enables us to create 3-loop, 4-loop, and 5-loop improved definitions of the field-strength, all of which are $`𝒪(a^4)`$-improved.
In addition to constructing $`Q`$, the field strength may be used in (5) to create a “reconstructed action”. Since this action is improved in a different manner to the cooling action, comparison of both actions may serve as a mechanism with which to test the improvement procedures.
In this brief report we will simply “test the water” with the values $`c_5=k_5=1/20`$ (5-loop), $`c_5=k_5=0`$ (4-loop), and $`c_5=1/10`$, $`k_5=1/90`$ (3-loop), and show that we obtain excellent results even though many other 5-loop operators, corresponding to other choices of $`c_5`$ and $`k_5`$, are possible.

## 4 RESULTS AND COMPARISON

### 4.1 Constructing Configurations

The configurations used in our investigations are generated using the Cabibbo-Marinari pseudo-heatbath algorithm . We generate $`SU(3)`$ configurations by looping twice over three diagonal $`SU(2)`$ subgroups. From a cold start (all link values set to zero) we thermalize for 5000 sweeps using an $`𝒪(a^2)`$-improved action and a fixed mean-link value, and then select configurations every 500 sweeps to create an ensemble of configurations with non-trivial topology, which are statistically distinct from each other.
To reveal the semi-classical structure of the configurations we use a cooling algorithm which likewise uses diagonally embedded $`SU(2)`$ subgroups, with appropriate link partitioning .
The results we present are produced from configurations on a $`12^3\times 24`$ lattice with $`\beta =4.60`$, providing a lattice spacing of $`a=0.125`$ fm.

### 4.2 Comparison of Cooling Schemes

In practice, because of the large discretization errors present, a cooling scheme based upon the plaquette action will destroy topological structure over the whole lattice if it proceeds for long enough. However, a cooling scheme using an improved action will stabilise at a state with non-trivial structure. We would anticipate cooling schemes with small discretization errors will reproduce expected continuum behaviour, especially integer values of $`Q`$ and $`S/S_0`$, and remain stable for hundreds of cooling sweeps. When comparing several cooling schemes, the one which fulfills these criteria best will be the preferred scheme for studying gauge field dynamics.

In Figure 1 we can see how the use of a 2-loop, $`𝒪(a^2)`$-improved action in the cooling algorithm stabilises a configuration at a state with $`S/S_06`$, indicating the presence of six instantons and/or anti-instantons, while plaquette-based cooling gradually destroys the structure of the configuration.
In Figure 2 we compare the values of the action obtained with 2-loop, 3LIS, 4LIS, and 5LIS-based cooling on the same configuration used in Figure 1. We can see that the 4LIS curve drops well below the desired integer value of six, indicating the presence of large negative discretization errors. This comparison and similar work with other configurations have lead us to conclude that 3LIS-based cooling gives the most integer-like results in general, and furthermore it is faster than 4LIS or 5LIS-based cooling, since it requires fewer computational steps.

### 4.3 Comparison of Topological Charge Improvement

In Figure 3 we compare the various improved $`Q`$ values obtained as we cool a configuration. The 3-loop improved operator gives the most continuum-like results, coming to within around 1 part in 6000 of an integer, and has the smallest computational cost of all the $`𝒪(a^4)`$-improved operators. Similar comparisons with other configurations have lead us to select the 3-loop operator as our preferred choice for studying the topological charge of cooled configurations.

### 4.4 Investigation of the Nahm Transform

Having established our preferred operators for investigating instanton dynamics, we are now in a position to compare the stability of $`|Q|=1`$ and $`|Q|1`$ configurations, to determine if we can find any evidence that self-dual $`|Q|=1`$ configurations are not permitted on the 4-torus, as implied by the Nahm transform. Figure 4 shows an overlay of the reconstructed action and topological charge of CFG 1, a $`|Q|1`$ configuration, and the corresponding values for CFG 2, a $`|Q|2`$ configuration (shifted down by an increment of one unit so that they overlap with those of CFG 1 on the scale of this diagram). These configurations are separated by 8000 thermalization sweeps. We can clearly see that the $`|Q|2`$ configuration stabilises quickly at a self-dual near-integer value, while the other configuration does not achieve self-duality and eventually destabilises. Investigations of other configurations have shown equivalent behaviour.

## 5 CONCLUSIONS

Highly improved operators enable us to cool configurations to stable non-trivial self-duality. Counting the number of instantons and anti-instantons present by evaluating $`S/S_0`$ and $`Q`$ allows us to detect evidence of continuum-like behaviour. The long-term stability of $`|Q|2`$ configurations indicates that our cooling algorithms are highly reliable, and hence the destabilization of $`|Q|1`$ configurations is due to the consequences of the Nahm transform, indicating that on the lattice, just as in the continuum, a configuration may not simultaneously fulfill the two criteria $`S/S_0=1`$ and $`|Q|=1`$ on the 4-torus.


format: text

text:
# Structural Probe of a Glass Forming Liquid: Generalized Compressibility

## I Introduction

The glass transition is still not well understood despite extensive study. There have been two main theoretical approaches to the problem: dynamic and thermodynamic. Theories in the first category view the glass transition as a kinetic phenomenon characterized by a growing relaxation time and viscosity Ediger96 ; Fourkas97 ; Gotze92 ; Donati98 ; Yamamoto97b . When the relaxation time exceeds the measurement time, particle motion appears to be arrested resulting in the glass transition. One of the most prominent theories espousing this view is the mode coupling theory (MCT) in which ideally the relaxation time diverges at a temperature $`T_C`$ above the experimental glass transition temperature Gotze92 . Interesting and fruitful concepts such as dynamic inhomogeneities Donati98 ; Yamamoto97a ; Tracht98  and the influence of the energy landscape on relaxation processes Goldstein69 ; Sastry98  have resulted from this approach. The thermodynamic viewpoint attributes the glass transition to an underlying phase transition hidden from direct experimental observation by extremely long relaxation times Ediger96 ; Fourkas97 ; Adam65 ; Gibbs58 ; Mezard99 . In most scenarios there is an underlying second order phase transition associated with a growing correlation length which produces diverging relaxation times as well as diverging static susceptiblities Kirkpatrick89 ; Sethna91 ; Kivelson95 ; Ernst91 ; Dasgupta91 ; Menon95 . More recently Mezard and Parisi Mezard99 ; Parisi97b  have argued that the underlying transition is actually a random first order transition signaled by a jump discontinuity in the specific heat.

Experimentally the glass transition is characterized by both kinetic and thermodynamic features. For example in the supercooled liquid kinetic quantities such as the viscosity and relaxation time grow rapidly as the temperature is lowered. When the system falls out of equilibrium below a certain temperature, thermodynamic quantities exhibit features reflecting the glass transition. For example as the system is cooled the specific heat has a step–like form and the dielectric constant exhibits a peak at a frequency dependent temperature.

In an effort to better characterize the glass transition we introduce a novel probe which we call a generalized compressibility Yvon58 . Unlike the specific heat which monitors energy fluctuations, this linear compressibility is a function of the microscopic structure of the system: it depends solely on the positions of the particles and not on their previous history. Since we do not need to compare the system’s state at different times, it is not a dynamic or kinetic quantity. Rather it is a thermodynamic quantity in the sense that it is purely a function of the microstate of the system dictated by its location in phase space. The generalized compressibility is easy to compute numerically, and it is simpler than the dielectric constant which involves both the translation and orientation of electric dipoles. In addition it does not suffer from finite size effects that can often plague measurements of the ordinary compressibility deduced from simulations. The generalized compressibility can be calculated in either the canonical or grand canonical ensembles. In particular it is well defined for a system with fixed volume $`V`$ and particle number $`N`$ in contrast to the ordinary compressibility which is defined for a system that has fluctuations in $`N`$ or $`V`$. The generalized compressibility should be directly measurable experimentally in colloidal suspensions of polystyrene spheres Weeks00  and possibly in other systems as well. In this paper we present measurements of this quantity in a molecular dynamics simulation of a two component system of soft spheres. We find that the linear generalized compressibility drops sharply as the temperature decreases below the glass transition temperature $`T_g`$. The drop becomes more and more abrupt as the measurement time increases. This is consistent with the structural arrest associated with a kinetic transition in which the system falls out of equilibrium. Similar results are seen as the density is increased at fixed temperature.

The paper is organized as follows. Section II describes the molecular dynamics simulations. Section III describes how the relaxation times and mode coupling $`T_C`$ are determined. These are useful for setting the time and temperature scales. Section IV describes our specific heat measurements which show a peak at the glass transition. Section V derives the expressions for the linear and nonlinear generalized compressibilities, and shows our results for these quantities. The linear generalized compressibility shows an abrupt drop at the same temperature and density as the peak in the specific heat. Section VI compares the ordinary isothermal compressibility with our linear generalized compressibility and shows the advantages of the latter. Section VII gives our results for the diffusion constant. Section VIII explains the relation between the linear generalized compressibility and the static structure factor. Finally we summarize our results in Section IX. A brief description of some of these results as well as results for a single component fluid that forms a crystal was reported earlier Carruzzo02 .

## II Molecular Dynamics Simulation

We have performed a molecular dynamics simulation on a glass forming liquid Weber85 ; Kob94  consisting of a 50:50 binary mixture of soft spheres in three dimensions. The two types of spheres, labelled A and B, differ only in their sizes. The interaction between two particles a distance $`r`$ apart is given by $`V_{\alpha \beta }(r)=ϵ[(\sigma _{\alpha \beta }/r)^{12}+X_{\alpha \beta }(r)]`$ where the interaction length $`\sigma _{\alpha \beta }=(\sigma _\alpha +\sigma _\beta )/2`$ with $`\sigma _B/\sigma _A=1.4`$ ($`\alpha `$, $`\beta =`$ A, B). For numerical efficiency, we set the cutoff function $`X_{\alpha \beta }(r)=r/\sigma _{\alpha \beta }\lambda `$ with $`\lambda =13/12^{12/13}`$. The interaction is cutoff at the minimum of the potential $`V_{\alpha \beta }(r)`$. Energy and length are measured in units of $`ϵ`$ and $`\sigma _A`$, respectively. Temperature is given in units of $`ϵ/k_B`$ where $`k_B`$ is Boltzmann’s constant, and time is in units of $`\sigma _A\sqrt{m/ϵ}`$ where $`m`$, the mass of the particles, is set to unity. The equations of motion were integrated using the leapfrog method Rapaport95  with a time step of 0.005. During each run the average density $`\rho _o=N/L^3`$ was fixed, and the temperature was kept constant using a constraint algorithm Rapaport95 . The volume $`V=L^3`$. $`N=N_A+N_B`$ is the total number of particles. The system occupies a cube with dimensions ($`\pm `$ L/2, $`\pm `$ L/2, $`\pm `$ L/2) and periodic boundary conditions.

We have done sweeps of both temperature and density. We fix the parameters so that crystallization is avoided upon cooling or when the density is increased. For the temperature sweeps, we fix the density at $`\rho _o=0.6`$. For $`\sigma _B/\sigma _A=1.4`$, this corresponds to a packing fraction of 1.04. Having a packing fraction larger than 1 means that each particle was interacting with other particles most, if not all, of the time. Our measuring procedure is the following. For runs where we cool the system, we start each run at a high temperature (T=1.5) and lower the temperature in steps of $`\mathrm{\Delta }T=0.05`$. At each temperature we equilibrate for $`10^4`$ molecular dynamics steps (md steps) and then measure the quantities of interest for $`N_\tau `$ additional md steps where $`N_\tau =10^5`$, $`2\times 10^5`$, $`10^6`$, $`3\times 10^6`$, or $`10^7`$. All the particles move at each md step. The results are then averaged over up to 40 different initial conditions (different initial positions and velocities of the spheres). We have done some runs in which we heat the system of particles by starting at our lowest temperature $`T=0.1`$ with a configuration obtained by cooling the system. We then increased the temperature in steps of $`\mathrm{\Delta }T=0.05`$. As before we equilibrate at each temperature for $`10^4`$ time steps and then measure quantities for an additional $`10^6`$ time steps.

We have also done some density sweeps in which we fix the temperature (T=1.0) and systematically change the density. The glass transition occurs as we increase the density. Colloidal experiments often study the glass transition as a function of density. We start each run at a low density ($`\rho =0.4`$) and increase the density in steps of $`\mathrm{\Delta }\rho =0.025`$. At each density we equilibrate for $`10^4`$ md steps and then measure the quantities of interest for $`N_\tau `$ additional md steps.

The glass transition occurs either as the temperature is lowered or as the density is raised. It is worth noting that temperature and density can be combined into a dimensionless parameter $`\mathrm{\Gamma }`$ Barrat90 :

$$\mathrm{\Gamma }=\rho \sigma _{eff}^3/T^{\frac{1}{4}}\sigma _{eff}^3=\underset{\alpha \beta }{}n_\alpha n_\beta \sigma _{\alpha \beta }^3$$
(1)
where $`\sigma _{eff}`$ represents an effective diameter for particles in the mixture. The concentration of each type of particle is given by $`n_A=N_A/N`$ and $`n_B=1n_A`$. For our simulations $`n_A=n_B=0.5`$. $`\mathrm{\Gamma }`$ is the relevant parameter when the particles spend most of their time sampling a nonzero interparticle potential, i.e., for $`\rho ^{1/3}<\sigma _{eff}`$. Thus $`\mathrm{\Gamma }`$ is particularly useful for interparticle interactions which fall off with distance as a power law and do not have a cutoff beyond which the interaction is zero. When cooling from the liquid phase, the glass transition is known to occur around $`\mathrm{\Gamma }=1.45`$ Barrat90 .

We have looked for phase separation of the two types of spheres by examining the distribution of large and small spheres in the neighborhood of large spheres and in the neighborhood of small spheres. We see no evidence for phase separation at either high ($`T=1.5`$) or low ($`T=0.15`$) temperatures at a density of $`\rho =0.6`$.

## III Relaxation Times and Mode Coupling $`T_C`$

As points of reference for the time and temperature scales, it is useful to find the mode coupling $`T_C`$ and the $`\alpha `$ relaxation times. We can find the relaxation times using the intermediate scattering function $`F(\stackrel{}{k},t)`$ which is a useful probe of the structural relaxation. It is the spatial Fourier transform of the van Hove correlation function $`G(\stackrel{}{r},t)`$ and the inverse time transform of the dynamic structure factor $`S(\stackrel{}{k},\omega )`$. There are two different types of intermediate scattering function: the self (incoherent) intermediate scattering function $`F_s(\stackrel{}{k},t)`$ and the full (coherent) intermediate scattering function $`F(\stackrel{}{k},t)`$.

In a computer simulation, the self (incoherent) part of the partial intermediate scattering function $`F_{s,\alpha }(\stackrel{}{k},t)`$ can be calculated directly using Hansen

$$F_{s,\alpha }(\stackrel{}{k},t)=\frac{1}{N_\alpha }\underset{i=1}{\overset{N_\alpha }{}}e^{i\stackrel{}{k}(\stackrel{}{r}_i(t)\stackrel{}{r}_i(0))}$$
(2)
where the subscript $`\alpha `$ refers to the particle type, A or B. $`\stackrel{}{r}_i(t)`$ is the position of particle $`i`$ at time $`t`$, and $`\mathrm{}`$ refers to an average over different configurations. The wave vector $`\stackrel{}{k}=2\pi \stackrel{}{q}/L`$ where $`\stackrel{}{q}`$ is a vector of integers. For an isotropic system $`F_{s,\alpha }(\stackrel{}{k},t)`$ depends only on the magnitude $`k=|\stackrel{}{k}|`$. We will choose $`k=k_{max}`$ where $`k_{max}`$ is the position of the first maximum of the partial static structure factor $`S_\alpha (k)`$. In Figure 1 we show the self intermediate scattering function $`F_{s,B}(k,t)`$ versus time at temperatures below the caging temperature ($`T0.4`$). The caging temperature is the highest temperature at which a plateau is present in the intermediate scattering function versus time. The plateau represents the temporary localization of a particle by a cage of other particles surrounding it.

Mode coupling theory is applicable in the temperature range below the caging temperature and somewhat above the mode coupling $`T_C`$. We define the relaxation time $`\tau _s`$ by $`F_s(k,\tau _s)=1/e`$. We determine the relaxation times for the seven highest temperatures shown in Figure 1 and then fit the temperature dependence of $`\tau _s(T)`$ to the mode coupling form $`\tau _s(T)(TT_C)^\gamma `$ to find $`T_C`$. For the self part of the intermediate scattering function, the actual value of $`\tau _s`$ increases as the magnitude of the wave vector decreases Kob95b . However, the value of $`T_C`$ is independent of $`k`$. $`\tau _s(T)`$ versus temperature and the mode coupling fit are shown in Figure 2. We find the best fit with the mode coupling temperature $`T_C=0.303`$ which corresponds to $`\mathrm{\Gamma }=1.46`$. Note that $`T_C`$ is determined from measurements made at temperatures where the system is equilibrated. Also shown in Figure 2 is the fit to the Vogel–Fulcher form $`\tau _s(T)=A\mathrm{exp}[B/(TT_{VF})]`$ with $`T_{VF}=0.21`$ which corresponds to $`\mathrm{\Gamma }=1.60`$. In doing the Vogel–Fulcher fit, we were able to use a much broader range of temperatures since $`T_{VF}`$ is much lower than the mode coupling $`T_C`$.

The full intermediate scattering function $`F(\stackrel{}{k},t)`$ is given by Hansen

$$F_\alpha (\stackrel{}{k},t)=\frac{1}{N_\alpha }\rho _{\stackrel{}{k},\alpha }(t)\rho _{\stackrel{}{k},\alpha }(0)$$
(3)
where the Fourier transform of the density $`\rho _\stackrel{}{k}(t)=_{i=1}^Ne^{i\stackrel{}{k}\stackrel{}{r}_i(t)}`$. The subscript $`\alpha `$ refers to the particle type, A or B. The longest $`\alpha `$ relaxation time can be determined from the full intermediate scattering function evaluated at $`k=k_{max}`$ Rinaldi01 . We set $`k_{max}=2\pi 8.3666/L`$ ($`L=8`$) which is the location of the first peak in the structure factor for type B particles. We define the $`\alpha `$ relaxation time $`\tau `$ as the time where $`F_\alpha (k_{max},t)`$ decays to $`1/e`$. At a temperature $`T=0.2895858`$ which is just below the mode coupling $`T_C`$, we find that $`F_B(k,t)`$ has fallen to $`1/e`$ at $`\tau =(1.0\pm 0.1)\times 10^6`$ md time steps for a system with 512 particles of which half are type B. This gives us a time scale by which to compare other times such as our run times. This value of $`\tau `$ shows no signs of aging Kob97b ; Kob00  and stays about the same even after 108 time steps. At higher temperatures this relaxation time is much shorter.

The runs used to determine the intermediate scattering function were done in a slightly different way from the other measurements. These runs were performed at a given temperature and density for $`N_\tau `$ md time steps with no change in temperature or density. The runs were started from a configuration that had been equilibrated at that temperature and density using parallel tempering. The parallel tempering technique is described in the appendix.

## IV Specific Heat

The specific heat is a thermodynamic quantity which undergoes a change signaling the glass transition. In experimental systems under constant pressure the specific heat exhibits a smooth step down as the temperature is lowered through the glass transition. In our simulations which are done at constant volume, the specific heat has a peak at the glass transition. It is a useful check of our calculation to see if the peak occurs at the same temperature (or density) as the drop in the linear generalized compressibility. There are two ways to compute the specific heat $`C_V`$ per particle at constant volume $`V`$. The first is by taking a derivative of the average energy $`E`$ per particle with respect to temperature: $`C_V=dE/dT`$. Since we study the system at discrete temperatures, we approximate the derivative by a finite difference:

$$C_V(T_n)=\frac{E(T_n)E(T_{n1})}{T_nT_{n1}}$$
(4)
where $`T_n>T_{n1}`$ for all integers $`n`$. The second way to calculate the specific heat is from the fluctuations:

$$C_V=Nk_B\beta ^2\left(E_P^2E_P^2\right)$$
(5)
where $`k_B`$ is Boltzmann’s constant, $`\beta `$ is the inverse temperature, and $`E_P`$ is the potential energy per particle. In our three dimensional simulations the kinetic energy per particle is given by $`3k_BT/2`$, so it is the fluctuations in the potential energy $`E_P`$ per particle which determine the temperature dependence of the specific heat. Thus

$$C_V=\frac{3}{2}k_B+Nk_B\beta ^2\left(E_P^2E_P^2\right)$$
(6)
In equilibrium these two ways of calculating the specific heat should agree. So we compare the results of calculating $`C_V`$ both ways as a check on our calculation and to make sure the system has equilibrated in all the basins that were visited in the energy landscape.

### IV.1 Specific Heat Versus Temperature

The specific heat at constant volume exhibits a peak at the glass transition as shown in Figure 3. The data in this figure is for 512 particles and was averaged over 6 runs with a measurement time of $`3\times 10^6`$ md steps. Notice that there is good agreement between calculating the specific heat by taking a derivative of the energy with respect to temperature (see eq. (4)) and by using fluctuations (see eq. (6)). This implies that the system has equilibrated within the basins that it visits in the energy landscape. We find similar agreement for other run times. At low temperatures the specific heat goes to 3$`k_B`$ as expected for classical oscillators while at high temperatures $`C_V`$ approaches $`3k_B/2`$ which corresponds to an ideal gas.

The peak in the specific heat occurs at $`T0.3`$ which corresponds to $`\mathrm{\Gamma }1.46`$. The temperature of the peak coincides with the mode coupling $`T_C=0.303`$ that we deduced from the intermediate scattering function data. Longer run times lead to a sharper peak in the specific heat as can be seen in Figure 4 which shows the specific heat for 512 particles for several different measuring times. The peaks would presumably be sharper if we had used a finer temperature scale. At high temperatures the agreement between the different times is very good. Perera and Harrowell Perera99  have found a specific heat peak in a two dimensional binary mixture of soft spheres. They argue that their peak is an equilibrium feature. However, in our case, at temperatures below the peak, the system has fallen out of equilibrium and has become trapped in a basin in the energy landscape. We shall see this later by examining the energy of the inherent structures (potential energy minima) as a function of temperature. Thus the fact that the peak in the specific heat occurs at or very close to the mode coupling $`T_C`$ is a result of the relaxation times (see Fig. 2) becoming comparable to and exceeding the simulation run times as the temperature drops below $`T_C`$. When this happens, the system falls out of equilibrium and undergoes a kinetic glass transition.

The specific heat $`C_P`$ of experimental systems at constant pressure exhibits a downward step at the glass transition during cooling and a peak at slightly higher temperatures upon heating DeBolt76 . As can be seen in Figure 5, in our warming up simulations, which are done at constant volume, the specific heat peak sharpens and moves toward higher temperatures compared to the cooling runs. This is consistent with what is seen in experiments. The hysteresis is consistent with the system falling out of equilibrium and getting stuck in a basin of the energy landscape.

As we mentioned in the introduction, some have suggested that the glass transition has an underlying second order phase transition Kirkpatrick89 ; Sethna91 ; Kivelson95 ; Ernst91 ; Dasgupta91 ; Menon95 . Unlike typical second order phase transitions, there is no experimental evidence that the specific heat diverges at the glass transition. This is consistent with our simulations. In simulations one looks for a divergence by examining whether the quantity increases systematically with system size. In Figure 6 we plot $`C_V`$ for systems with 64, 216, 512 and 1000 particles. As one can see, the specific heat does not exhibit any size dependence. However, we cannot rule out the possibility that a thermodynamic phase transition occurs at temperatures below where we fall out of equilibrium. Indeed theories which postulate a thermodynamic transition put the transition temperature well below the mode coupling $`T_C`$.

### IV.2 Specific Heat Versus Density

In Figure 7 we show the specific heat as a function of density. As the density increases, the specific heat rises to a peak at $`\rho _o^{\mathrm{peak}}=0.8`$. This corresponds to $`\mathrm{\Gamma }=1.44`$ which is in good agreement with the $`\mathrm{\Gamma }`$ value of 1.46 that we found for the specific heat peak when we varied the temperature. Going to higher densities corresponds to going to lower temperatures. At densities higher than 0.8, the system falls out of equilibrium.

## V Generalized Compressibilities

As we mentioned in the introduction, the generalized compressibilities are thermodynamic probes that are a function of the microscopic structure of the system. They are solely a function of the positions of the particles and do not depend on their histories. So one could take snapshots of the configurations of the particles at different instances, scramble the order of the snapshots, and still be able to calculate the generalized compressibilities. Averaging over these snapshots corresponds to ensemble averaging. In this sense the generalized compressibilities are thermodynamic quantities which can be calculated solely from the microstates of the system and do not depend on the system’s dynamics or kinetics.

We now derive expressions for the linear and nonlinear generalized compressibility. To probe the density fluctuations, we follow the approach of linear response theory and consider applying an external potential $`\frac{\mathrm{\Delta }P}{\rho _o}\varphi (\stackrel{}{r})`$ which couples to the local density $`\rho (\stackrel{}{r})=_{i=1}^N\delta (\stackrel{}{r}\stackrel{}{r}_i)`$ where $`\stackrel{}{r}_i`$ denotes the position of the $`i^{th}`$ particle. $`\rho _o`$ is the average density. $`\mathrm{\Delta }P`$ has units of pressure and sets the magnitude of the perturbation. $`\varphi (\stackrel{}{r})`$ is a dimensionless function of position that must be compatible with the periodic boundary conditions imposed on the system, i.e., it must be continuous across the boundaries, but is otherwise arbitrary. This adds to the Hamiltonian H of the system a term

$$U=\frac{\mathrm{\Delta }P}{\rho _o}_Vd^3r\varphi (\stackrel{}{r})\rho (\stackrel{}{r})=\frac{\mathrm{\Delta }P}{\rho _o}\underset{i}{}\varphi (\stackrel{}{r}_i)\frac{\mathrm{\Delta }P}{\rho _o}\rho _\varphi $$
(7)
where we have defined $`\rho _\varphi =_Vd^3r\varphi (\stackrel{}{r})\rho (\stackrel{}{r})=_i\varphi (\stackrel{}{r}_i)`$. $`\rho _\varphi `$ is the inner product of $`\varphi `$ and $`\rho (\stackrel{}{r})`$, and we can regard it as a projection of the density onto a basis function $`\varphi (r)`$, i.e., $`\rho _\varphi =<\rho |\varphi >`$. It weights the density fluctuations according to their spatial position. The application of the external potential will induce an average change $`\delta \rho _\varphi `$ in $`\rho _\varphi `$:

$$\delta \rho _\varphi =\rho _\varphi _U\rho _\varphi _{U=0}$$
(8)
where the thermal average $`\rho _\varphi _U`$ is given by

$$\rho _\varphi _U=\frac{1}{Z}\text{Tr}\left[e^{\beta (H+U)}\rho _\varphi \right]$$
(9)
The partition function $`Z=\text{Tr}e^{\beta (H+U)}`$ and $`\beta `$ is the inverse temperature. For small values of $`\mathrm{\Delta }P`$, this change can be calculated using perturbation theory Toda91 . Up to third order in $`\mathrm{\Delta }P`$, we find

$`\delta \rho _\varphi `$ $`=`$ $`{\displaystyle \frac{\beta \mathrm{\Delta }P}{\rho _o}}\rho _\varphi ^2_c+{\displaystyle \frac{\beta ^2\mathrm{\Delta }P^2}{2\rho _o^2}}\rho _\varphi ^3_c`$ (10)
$`{\displaystyle \frac{\beta ^3\mathrm{\Delta }P^3}{6\rho _o^3}}\rho _\varphi ^4_c,`$
where the cumulant averages are

$`\rho _\varphi ^2_c`$ $`=`$ $`\rho _\varphi ^2\rho _\varphi ^2`$ (11)
$`\rho _\varphi ^3_c`$ $`=`$ $`\rho _\varphi ^33\rho _\varphi \rho _\varphi ^2+2\rho _\varphi ^3`$ (12)
$`\rho _\varphi ^4_c`$ $`=`$ $`\rho _\varphi ^44\rho _\varphi \rho _\varphi ^33\rho _\varphi ^2^2+`$ (13)
$`12\rho _\varphi ^2\rho _\varphi ^26\rho _\varphi ^4`$
with the thermal average $`\rho _\varphi ^n=\rho _\varphi ^n_{U=0}`$. The third order cumulant, eq.(12), is zero in the liquid phase because for every configuration there exists an equivalent configuration with the opposite sign of $`\left(\rho _\varphi \rho _\varphi \right)`$ and so we will not consider this term any further. We can recast eq. (10) as a power series in the perturbation $`\mathrm{\Delta }P`$:

$$\frac{\delta \rho _\varphi }{N}=\frac{1}{6\rho _ok_BT}\chi _l\mathrm{\Delta }P+\frac{1}{6(\rho _ok_BT)^3}\chi _{nl}(\mathrm{\Delta }P)^3$$
(14)
where

$$\chi _l=\frac{6}{N}(\rho _\varphi )^2_c\chi _{nl}=\frac{1}{N}(\rho _\varphi )^4_c.$$
(15)
In this paper we will focus our attention on the linear ($`\chi _l`$) and nonlinear ($`\chi _{nl}`$) dimensionless generalized compressibilities defined by the above expressions.

We now discuss the choice of the function $`\varphi `$. We consider applying the potential along the direction $`\mu `$ of one of the coordinate axes so that $`\varphi (\stackrel{}{r})=\varphi (r^\mu )`$. A natural candidate for $`\varphi (r^\mu )`$ is $`\mathrm{cos}(k_\mu r^\mu )`$ (no implied sum over repeated indices) with $`k_\mu =2\pi n/L`$, where $`n=1,2,\mathrm{}`$ In this case, $`\rho _\varphi `$ is the $`k^{th}`$ mode of the cosine transform of the density. We will also consider the simpler function $`\varphi (r^\mu )=|r^\mu |/L`$. The absolute value corresponds to the case where all the particles feel a force along the $`\mu `$th direction pointing towards the origin. It gives results very similar to $`\varphi (r^\mu )=\mathrm{cos}(k_\mu r^\mu )`$ for small $`k`$ at a fraction of the computational cost. (No sum over repeated indices.) So our results in this paper correspond to two cases:

$$\rho _\varphi =\underset{i}{}|r_i^\mu |/L$$
(16)
which is rather like a center of mass, and

$$\rho _\varphi =\underset{i}{}\mathrm{cos}(k_\mu r_i^\mu )$$
(17)
Since the system is isotropic, we compute the compressibilities for each direction and then average over the direction $`\mu `$.

In most of our calculations we work in the canonical ensemble where we fix the volume $`V`$, the number $`N`$ of particles and the density $`\rho _o`$. However, it is straightforward to generalize our results to the grand canonical ensemble where the number of particles is not fixed. We simply replace the thermal average defined in equation (9) by

$$\rho _\varphi _U=\frac{1}{𝒵}\underset{N}{}e^{\mu N}\text{Tr}\left[e^{\beta (H_N+U_N)}\rho _\varphi \right]$$
(18)
where $`\mu `$ is the chemical potential, $`H_N`$ is the Hamiltonian with $`N`$ particles, $`U_N`$ is given by eq. (7) for a system with $`N`$ particles, and $`𝒵`$ is the grand canonical partition function given by

$$𝒵=\underset{N}{}e^{\mu N}\text{Tr}\left[e^{\beta (H_N+U_N)}\right]$$
(19)
The generalized compressibilities can be defined using equations (11) through (15) with the thermal averages $`\rho _\varphi ^n=\rho _\varphi ^n_{U=0}`$ defined in the grand canonical ensemble.

### V.1 Results for Linear Generalized Compressibility

We now turn to our results for the binary glass forming liquid.

#### V.1.1 $`\chi _l`$ from Absolute Value of Positions versus Temperature

We will first discuss the linear generalized compressibility calculated from the absolute values of the particle positions using eqs. (15) and (16).

Figure 8 shows the linear generalized compressibility as a function of temperature for different run times. The compressibility at high temperatures is independent of $`T`$. In the vicinity of the glass transition $`\chi _l`$ drops. Notice that as the measuring time $`t_M`$ increases (and hence as the cooling rate decreases), the temperature of the drop decreases and becomes more abrupt. The measuring time can be thought of as the number of snapshots at a single temperature that we use to calculate the compressibility. The linear compressibility is proportional to the width of the distribution of $`\rho _\varphi `$, so the drop in $`\chi _l`$ corresponds to the sudden narrowing of the distribution $`P(\rho _\varphi )`$. If we regard $`\rho _\varphi `$ as a generalized center of mass, then the drop in $`\chi _l`$ signals the sudden arrest in the fluctuations of the generalized center of mass. In other words, at the glass transition the motion of the particles is largely frozen and hence, the generalized center of mass does not move around much. This is consistent with recent observations of the colloidal glass transition in which the size of the clusters of “fast” particles drops dramatically at the glass transition Weeks00 .

Notice that at longer measuring times, the temperature $`T_{\mathrm{drop}}`$ at which the generalized linear compressibility drops is roughly at the mode coupling temperature $`T_C=0.303`$. Let us define $`T_{\mathrm{drop}}`$ as the temperature at which $`\chi _l`$ has dropped halfway down. For $`10^6`$ md steps, $`T_{\mathrm{drop}}0.33`$; for $`3\times 10^6`$ md steps, $`T_{\mathrm{drop}}0.30`$; and for $`10^7`$ md steps, $`T_{\mathrm{drop}}0.27`$. Thus we are able to stay in equilibrium down to the mode coupling temperature for our longer runs. This is what we would expect when we compare these run times, which are longer than 1 million time steps, to the $`\alpha `$ relaxation time $`\tau `$ which is about 1 million time steps at $`T=0.29`$ which is just below $`T_C`$. Thus the fact that the drop in the linear generalized compressibility occurs at or very close to the mode coupling $`T_C`$ is a result of the relaxation times (see Fig. 2) becoming comparable to and exceeding the simulation run times as the temperature drops below $`T_C`$. When this happens, the system falls out of equilibrium and undergoes a kinetic glass transition.

The behavior exhibited by $`\chi _l`$ can be quantified using a scaling ansatz: $`\chi _l(t_M,T)=g(\mu =t_M/\tau (T))`$, where the characteristic time has the Volgel–Fulcher form $`\tau (T)=\mathrm{exp}(A/(TT_o))`$. The inset of Figure 8 shows that the data collapse onto a single curve with $`A=0.75`$, $`T_o=0.15`$. (The data could not be fitted using $`\tau (T)=A(TT_o)^\gamma `$ as suggested by simple mode coupling theories Gotze92 .) Notice that $`T_o`$ lies below the Vogel–Fulcher temperature $`T_{VF}=0.21`$ and the MCT critical temperature $`T_C=0.303`$ deduced by fitting the temperature dependence of the relaxation times. This scaling suggests that $`\chi _l`$ becomes a step function for infinite $`t_M`$ and that the drop in the compressibility would become a discontinuity at infinitely long times. This abrupt drop is consistent with a sudden arrest of the motion of the particles in the liquid which is the kinetic view of the glass transition. The abrupt drop also appears to be in agreement with Mezard and Parisi’s proposal that the glass transition is a first order phase transition with a jump in the specific heat Mezard99 ; Parisi97b . Indeed the temperature at which $`\chi _l`$ drops agrees with the temperature of the peak in the specific heat shown in Figure 4. The specific heat provides an independent check of the glass transition temperature. However, the drop in $`\chi _l`$ and the specific heat peak are due to the system falling out of equilibrium, and therefore we cannot really tell if there is an underlying true thermodynamic transition.

One way to see that the system is falling out of equilibrium is to plot the inherent structure energy per particle Stillinger82 ; Stillinger84 ; Sastry98 . An inherent structure is a particular system configuration whose energy corresponds to the minimum of a basin in the energy landscape. The energy landscape is a 3N dimensional surface defined by the potential energy of the system which is a function of the particles’ coordinates. During each run we sampled the configurations found at each temperature. Each configuration lies somewhere in a basin and we used the method of conjugate gradients Press92  to find the inherent structure energy of that basin. The result is shown in Figure 9 where we plot the average inherent structure energies versus the temperature of the configuration that was originally saved. At high temperatures the inherent structure energy $`e_{IS}`$ per particle is flat as a function of temperature. As the system is cooled, $`e_{IS}`$ decreases rather steeply Sastry98 . The inherent structure energy flattens off at low temperatures where the system has fallen out of equilibrium and has become stuck in one basin. For each measuring time the temperature below which the generalized linear compressibility drops corresponds to the temperature below which the inherent structure energy flattens off at low temperatures. Thus the temperature of the drop in $`\chi _l`$ and the peak in the specific heat corresponds to the temperature below which the system falls out of equilibrium and ceases to explore deeper basins of the energy landscape.

The behavior of the linear generalized compressibility seen in Figure 8 is similar to that seen in measurements of the real part of the frequency dependent dielectric function $`\epsilon ^{}(\omega )`$ Menon95 . In that case as the frequency decreased, the temperature of the peak in $`\epsilon ^{}(\omega )`$ decreased and the drop in $`\epsilon ^{}(\omega )`$ below the peak became more abrupt. By extrapolating their data to $`\omega =0`$, Menon and Nagel Menon95  argued that $`\epsilon ^{}(\omega =0)`$ should diverge at the glass transition, signaling a second order phase transition. We have looked for evidence of this divergence by examining samples of different sizes to see if the linear generalized compressibility increased systematically with system size. As shown in Figure 10 we find no size dependence and no indication of a diverging linear generalized compressibility. However, as we mentioned earlier, this does not preclude the possibility that a thermodynamic phase transition occurs below $`T_{drop}`$. It simply means that if there is a growing correlation length, it is smaller than the size of our system at $`T>T_{drop}`$. Another possible reason for the absence of size dependence may be that if there is an underlying thermodynamic phase transition, then its order parameter may not couple to the local density $`\rho (\stackrel{}{r})`$.

So far we have shown the results of cooling the system. In order to look for hysteretic behavior we have done runs in which we heat a system of 512 particles by starting at our lowest temperature $`T=0.1`$ with a configuration obtained by cooling the system. We then increased the temperature in steps of $`\mathrm{\Delta }T=0.05`$. As before we equilibrate at each temperature for $`10^4`$ time steps and then measure quantities for an additional $`10^6`$ time steps. Our results are shown in Figure 11. Notice the slight hysteresis with the rise in $`\chi _l`$ upon warming being at a slightly higher temperature than the drop in $`\chi _l`$ upon cooling. This hysteresis is consistent with the kinetic arrest of motion and with the hysteresis found for the specific heat in Figure 5.

#### V.1.2 $`\chi _l`$ from Cosine of Positions versus Temperature

We now consider calculating the linear generalized compressibility from the cosine transform of the density using eqs. (15) and (17). So if we apply a cosine potential along, say the $`\mu =x`$ direction, then

$$\rho _\varphi =\underset{i}{}\mathrm{cos}(k_xx_i)$$
(20)
where the wave vector $`k_x=2\pi n/L`$ with $`n=1,2,\mathrm{}`$. The wave vectors are compatible with the periodic boundary conditions of our simulations. Since the system is isotropic, we average $`\chi _l`$ over the $`x`$, $`y`$, and $`z`$ directions. The resulting linear generalized compressibility is qualitatively similar in its temperature dependence to the linear compressibility calculated using the absolute values of the particles’ positions (eq. (16)). Figure 12 shows the linear generalized compressibility versus temperature for various values of the wave vector. The data is for a binary mixture of 512 particles with a measuring time of $`10^6`$ md steps and averaged over 10 runs.

Just as for the absolute value case, we find that as we increase the measuring time, the drop in the linear generalized compressibility calculated using cosine becomes sharper at the glass transition. This is shown in Figure 13 which shows $`\chi _l`$ as a function of temperature for measuring times of $`2\times 10^5`$, $`10^6`$, and $`3\times 10^6`$ md steps with $`k=2\pi /L`$, i.e., $`n=1`$.

Figure 14 shows the linear generalized susceptibility versus the wave vector $`k`$ in units of $`2\pi /L`$ for various temperatures. Note that the dependence is nonmonotonic.

#### V.1.3 $`\chi _l`$ from Absolute Value of Positions versus Density

In Figure 15 we plot the generalized linear compressibility versus density calculated from the absolute value of the positions of the particles. In Figure 16 we plot the generalized linear compressibility versus density calculated from the cosine of the particles’ positions. In both cases we see that $`\chi _l`$ drops with increasing density. The drop becomes more abrupt as the measuring time increases. This drop is similar to what we see when we cool the system at fixed density. The density ($`\rho 0.8`$, $`\mathrm{\Gamma }1.44`$) at which the drop occurs agrees with the density at which there is a peak in the specific heat as shown in Fig. 7.

#### V.1.4 $`\chi _l`$ from Absolute Value of Positions versus Temperature in a Slab Geometry

So far we have considered systems with a fixed number of particles, but as we mentioned earlier in this section, we can generalize our results to the grand canonical ensemble where the number $`N`$ of particles can vary. We have examined the generalized linear compressibility $`\chi _l`$ calculated from the absolute value of the particle positions using equations (15) and (16) for a slab of our system. In other words we have divided a system of $`8^3`$ particles into 8 slabs of equal thickness perpendicular to the $`x`$ axis. The number of particles in any given slab is not fixed. However, in eq. (15) we set the average number $`N`$ of particles in each slab equal to the total number of particles in the system divided by the number of layers, i.e., $`N=8^3/8=64`$. Such a slab geometry mimicks experiments on colloidal suspensions of binary mixtures in which the focal plane of the camera can essentially see only one monolayer of polystyrene balls comment:cherry . In figure 17 we show the generalized linear compressibility for a slab for two different measuring times. Again we see that the drop is sharper as the measurement time becomes longer. Thus allowing for fluctuations in the number of particles does not change the qualitative behavior of $`\chi _l`$ at the glass transition. Comparing Figures 8 and 17, we see that the temperature of the drop for the slab and the bulk agree. We also notice that the drop is sharper for the bulk where presumably the greater number of particles results in better statistical averaging.

#### V.1.5 $`\chi _l`$ from Absolute Value of Positions versus Density in a Slab Geometry

Since experiments on colloidal suspensions usually vary the density rather than the temperature, we have done simulations where we set the temperature $`T=1`$ and vary the density. Again we divide our system of $`N=8^3=512`$ particles into 8 slabs and measure $`\chi _l`$ in one of those slabs. The results are shown in Figure 18.

As one can see from the figure, $`\chi _l`$ drops as the density is increased and the drop becomes more abrupt as the measuring time lengthens. Comparing Figures 15, 16 and 18, we see that the drop occurs at the same density ($`\rho 0.8`$) as in the bulk.

### V.2 Nonlinear Generalized Compressibility

We now turn to the case of the nonlinear generalized compressibility $`\chi _{nl}`$ given by eq. (15). We are motivated by the case of spin glasses where the nonlinear magnetic compressibility diverges at the spin glass transition while the linear compressibility only has a cusp Bhatt88 ; Levy86 . There have been a few studies of nonlinear response functions in structural glasses Dasgupta91 ; Wu91 , but these have not found any divergences. Our results are consistent with this conclusion. In particular we find that the nonlinear generalized compressibility is zero above and below the glass transition temperature, though it does show a glitch at the glass transition temperature. There is no systematic increase with system size, indicating the absence of a divergence at temperatures above the glass transition. This does not rule out a divergence below the glass transition temperature where our system has fallen out of equilibrium. It also does not rule out a thermodynamic transition that does not couple to the local density. Because $`\chi _{nl}`$ is sensitive to the tails of the distribution of $`\rho _\varphi `$, one must be careful to obtain a good ensemble average in the liquid above the glass transition temperature. We have done this by doing 16 or 32 runs, each involving 200,000 time steps, with different initial conditions, stringing them together as though they were one long run and then taking the appropriate averages. In some sense this approach mixes molecular dynamics and Monte Carlo; the simulation follows the equations of motion for a given amount of time and then “jumps” to another configuration which again evolves according to the molecular dynamics equations until the next jump. We call this approach “global averaging.” It produces a better ensemble average of $`<\rho _\varphi ^2>^2`$ which enters into $`\chi _{nl}`$ in eq. (13). The resulting $`\chi _{nl}`$ is shown in Figure 19 which was calculated from the absolute values of the particles’ positions using eqs. (15) and (16).

$`\chi _{nl}`$ also took longer to equilibrate than $`\chi _l`$. By plotting $`\chi _{nl}`$ versus run time, we found that one had to run at least $`10^6`$ time steps at $`T=1`$ before $`\chi _{nl}`$ appeared to saturate (see Figure 20).

## VI Ordinary Isothermal Compressibility

The ordinary isothermal compressibility $`\kappa _T`$ can be related to the fluctuations in the number, volume or density of the system. In most of our calculations we fix the volume, number of particles and density, so that there are no such fluctuations and the system has $`\kappa _T=0`$. However in the grand canonical ensemble $`\kappa _T`$ is given by

$$\kappa _T=\frac{1}{\rho _ok_BT}\frac{N^2N^2}{N}$$
(21)
where $`\rho _o=<N>/V`$. To relate $`\kappa _T`$ to the linear generalized compressibility $`\chi _l`$, we choose a uniform potential $`\varphi (\stackrel{}{r})=1`$. Then $`\rho _\varphi =d^3r\rho (\stackrel{}{r})=N`$ and

$$\kappa _T=\frac{1}{6\rho _ok_BT}\chi _l$$
(22)
In principle one can also obtain $`\kappa _T`$ from the $`k0`$ limit of the static structure factor $`S(\stackrel{}{k})=(1/N)\rho _\stackrel{}{k}\rho _\stackrel{}{k}`$ where $`\rho _\stackrel{}{k}=_{i=1}^N\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i)`$ is the Fourier transform of the local density $`\rho (\stackrel{}{r})`$. The limit of $`S(\stackrel{}{k})`$ for $`k0`$ in an isotropic and homogeneous system is Hansen

$`S(0)`$ $`=`$ $`1+\rho _o{\displaystyle (g(r)1)d^3r}`$ (23)
$`=`$ $`\rho _ok_BT\kappa _T`$
where $`g(r)`$ is the radial distribution function. Note that in a system with fixed volume and particle number, the normalization of $`g(r)`$ leads to $`S(k0)=0`$. This is consistent with the fact that such a system has $`\kappa _T=0`$. Eq. (23) yields a nonzero value for $`\kappa _T`$ in a system which has fluctuations in volume, particle number or density. Even in such a compressible system taking the $`k0`$ limit of $`S(k)`$ can suffer from finite size effects Bernu85  because the farthest apart that any two particles can be along any given coordinate axis is $`L/2`$ when there are periodic boundary conditions. So at wave vectors $`k`$ with components smaller than 4$`\pi /L`$, $`S(k)`$ can have spurious results. (For example in our simulations we found that this manifests itself as a slight upturn in $`S(k)`$ at small $`k`$.) It is possible to extrapolate to distances larger than $`L/2`$ using various approaches Bernu85 . We chose not to use this approach to calculate $`\kappa _T`$ since we work in a system with fixed $`N`$ and $`V`$. We should note however that simulations Horbach99  working in the $`NVT`$ ensemble with fixed $`N`$, $`V`$, and $`T`$, have successfully used eq. (23) to find the isothermal compressibility. We can resolve this with the fact that $`S(k0)=0`$ in the $`NVT`$ ensemble by noting that for values of $`k>4\pi /L`$, $`S(k)`$ should give the same value in the $`NVT`$ ensemble as in the grand canonical ensemble. So if $`L`$ is large enough, fitting $`S(k)`$ to the small $`k`$ form $`S(k)=S(k0)+Ak^2`$, where $`A`$ is a constant, should yield the correct value of $`\kappa _T`$ as long as $`k>4\pi /L`$.

Rather than using eq. (23), we calculated $`\kappa _T`$ by monitoring a small subvolume inside of our system and keeping track of the fluctuations in the number of particles in the subvolume. Let us define a dimensionless ordinary isothermal compressibility $`K_T`$ by

$$\kappa _T=\frac{1}{\rho k_BT}K_T$$
(24)
where

$$K_T=\frac{N^2N^2}{N}$$
(25)
In order to calculate $`K_T`$, we have monitored a subvolume that had on average 25% of the total number of particles. Essentially we drew an imaginary boundary in the middle of our system that enclosed 25% of the total volume and counted the number of particles in this subvolume as a function of time. By monitoring the fluctuations in the number of particles in this subvolume, we could calculate $`K_T`$. The results for a subvolume which had on average 128 particles out of a total of 512 particles are shown in Figure 21 where $`K_T`$ is plotted versus temperature. We see that it has the same basic shape as the linear generalized compressibililty with a drop at the same temperature as $`\chi _l`$. As with $`\chi _l`$, the drop becomes sharper with increasing measuring time.

While $`K_T`$ shows behavior similar to the linear generalized compressibility $`\chi _l`$ as a function of temperature for a given size, it is unlike $`\chi _l`$ in that it suffers from finite size effects. We have demonstrated this by making measurements on systems with a total of 64, 216, 512, and 1000 particles. The measurements were made by counting the number of particles in a subvolume that was 25% of the total volume. Such a small subvolume has a large surface to volume ratio which produces large finite size effects. To understand this, we note that in such a small subvolume, a significant number of the particles are very close to the boundary of the subvolume. Fluctuations in the positions of these particles moves them in and out of the subvolume, producing large fluctuations in the number of particles in the subvolume. The smaller the system, the bigger this effect is. This produces large finite size effects in $`K_T`$ even at high temperatures where the system easily equilibrates. This can be seen in Figure 22. One can see that $`K_T`$ decreases with increasing system size at high temperatures above the drop in the compressibility. One of the advantages of the linear generalized compressibility is the absence of these finite size effects. In fact the linear generalized compressibility shows no size dependence at temperatures above the observed glass transition (see Figure 10).

## VII Diffusion Constant

The diffusion of the particles reflects the kinetics of the system and becomes very small below the glass transition temperature. We calculate the diffusion constant $`D`$ using the Einstein relation

$$D=\underset{t\mathrm{}}{lim}\frac{1}{6Nt}\underset{i=1}{\overset{N}{}}(𝐫_i(t)𝐫_i(0))^2$$
(26)
where $`𝐫_i`$ are true displacements of the $`i`$th particle. Since we are using periodic boundary conditions, if the particle has crossed the box several times, then this must be included in $`𝐫_i`$.

### VII.1 Diffusion Constant Versus Temperature

As the system is cooled through the glass transition, the diffusion constant calculated using equation (26) becomes very small. This is shown in Figure 23 where the diffusion constant for 512 particles is plotted on a logarithmic scale. The diffusion constant varies smoothly over the entire temperature range. The curves corresponding to different cooling rates begin to separate as the system falls out of equilibrium at the glass transition temperature where the specific heat peaks and where the linear generalized compressibility drops abruptly. Figure 23 also shows the diffusion constant for binary mixtures of several different sizes. Notice that there is no apparent size effect.

### VII.2 Diffusion Constant Versus Density

In Figure 24 we show the diffusion constant as a function of density. We see that the diffusion decreases smoothly as the density $`\rho `$ increases. The curves corresponding to different measurement times begin to separate at the glass transition density where the specific heat peaks and where the linear generalized compressibility drops abruptly.

## VIII Relation between $`\chi _l`$ and $`S(k)`$

When the system has translational invariance, we can relate the linear generalized compressibility $`\chi _l`$ to the static structure factor $`S(k)`$ which is measured in experiments such as neutron scattering. $`S(k)`$ is also used as an input for mode coupling theories Gotze92  and it is generally assumed that $`S(k)`$ does not show any essential variations near the glass transition as the temperature or density is varied. As we shall see, our calculation agrees with this.

By definition, the static structure factor $`S(\stackrel{}{k})=(1/N)\rho _\stackrel{}{k}\rho _\stackrel{}{k}`$. To relate $`S(k)`$ to $`\chi _l`$, we first note that $`S(k)`$ is the Fourier transform of the static density–density autocorrelation function $`G(\stackrel{}{r})`$, i.e., $`S(\stackrel{}{k})=\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r})G(\stackrel{}{r})d^3r`$, where $`G(\stackrel{}{r})=(1/N)\rho (\stackrel{}{r}^{}+\stackrel{}{r})\rho (\stackrel{}{r}^{})d^3r^{}`$. In the supercooled liquid above the glass transition, the system has translational invariance, and we can write

$$\rho (\stackrel{}{r})\rho (\stackrel{}{r}^{})=\frac{1}{V}g\left(\stackrel{}{r}\stackrel{}{r}^{}\right)$$
(27)
where $`V`$ is volume and $`g`$ is a function of the difference $`\left(\stackrel{}{r}\stackrel{}{r}^{}\right)`$. In this case $`G(\stackrel{}{r})=g(\stackrel{}{r})/N`$. In $`\chi _l`$ we meet

$$\rho _\varphi ^2=d^3rd^3r^{}\varphi (\stackrel{}{r})\varphi (\stackrel{}{r}^{})\rho (\stackrel{}{r})\rho (\stackrel{}{r}^{})$$
(28)
If there is translational invariance,

$`\rho _\varphi ^2`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d^3rd^3r^{}\varphi (\stackrel{}{r})\varphi (\stackrel{}{r}^{})g\left(\stackrel{}{r}\stackrel{}{r}^{}\right)}`$ (29)
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \frac{d^3k}{(2\pi )^3}\varphi (\stackrel{}{k})\varphi (\stackrel{}{k})g(\stackrel{}{k})}`$
$`=`$ $`\rho _o{\displaystyle \frac{d^3k}{(2\pi )^3}|\varphi (\stackrel{}{k})|^2S(\stackrel{}{k})}`$
where $`\varphi (\stackrel{}{k})`$ is the Fourier transform of $`\varphi (\stackrel{}{r})`$. Converting the integral $`Vd^3k/(2\pi )^3`$ to a sum $`_\stackrel{}{k}`$ and using eq. (15), we obtain

$$\chi _l=\frac{6}{N}\left[\frac{N}{V^2}\underset{\stackrel{}{k}}{}\left(|\varphi (\stackrel{}{k})|^2S(\stackrel{}{k})\right)\rho _\varphi ^2\right]$$
(30)
As an example, let us choose $`\varphi (\stackrel{}{r})=\mathrm{cos}(k_xx)`$ with the proviso that $`k_x0`$. (When $`\stackrel{}{k}=0`$, $`\chi _l=0`$ since the potential is uniform and there are no fluctuations allowed with fixed $`N`$ and $`V`$.) With this choice of $`\varphi (\stackrel{}{r})`$, $`\rho _\varphi =_i\mathrm{cos}(k_xx_i)`$. Then one can show by explicitly calculating $`\varphi (\stackrel{}{k})`$ and by using eqs. (29) and (30) that

$$\chi _l(\stackrel{}{k})=3S(\stackrel{}{k})\frac{6}{N}[\mathrm{Re}\rho _\stackrel{}{k}]^2$$
(31)
where translational invariance allows us to write

$`S(\stackrel{}{k})`$ $`=`$ $`(2/N)\rho _\varphi ^2`$ (32)
$`=`$ $`(2/N){\displaystyle \underset{ij}{}}\left[\mathrm{cos}(k_xx_i)\mathrm{cos}(k_xx_j)\right]`$
Note that the value of $`k`$ that we use to probe the glass transition is typically of order $`k_L2\pi /L`$ which is much smaller than the value of $`k_{\mathrm{peak}}`$ at which the first peak in $`S(k)`$ appears. $`k_{\mathrm{peak}}2\pi /\sigma _Ak_L`$ where $`\sigma _A`$ is the diameter of type A spheres in the binary mixture.

We have numerically calculated $`[S(k)]_{\mathrm{run}}`$ for our binary mixture using equation (32) with $`\varphi (\stackrel{}{r})=\mathrm{cos}(k_\mu r_\mu )`$. (No sum over repeated indices. $`\mu =x`$, $`y`$, or $`z`$.) $`[\mathrm{}]_{\mathrm{run}}`$ is an average over all the runs and over $`\stackrel{}{k}`$ being parallel to $`x`$, $`y`$, and $`z`$. The result is shown in Figure 25 and one can see that $`[S(k)]_{\mathrm{run}}`$ does not vary much through the glass transition which is consistent with what is assumed in mode coupling theory. The figure also shows $`[\mathrm{Re}<\rho _k>]_{\mathrm{run}}`$ where $`<\mathrm{}>`$ is a thermal average over a single run. If $`[S(k)]_{\mathrm{run}}`$ and $`[\mathrm{Re}<\rho _k>]_{\mathrm{run}}`$ do not vary much through the glass transition, how can the difference between the two terms in (31) decrease and produce a drop in $`\chi _l`$? To answer this, note that there are two inequivalent ways in which one can calculate $`\chi _l`$. So far we have calculated $`\chi _l`$ for each run and then averaged over the different runs. This approach is what we used in Figs. 811 and results in a sharp drop in the linear generalized compressibility at the glass transition. Let us call this a run–by–run average for which we can write:

$$\chi _l=\frac{6}{N}\left[\rho _\varphi ^2\rho _\varphi ^2\right]_{\mathrm{run}}$$
(33)
The drop in $`\chi _l`$ comes about because the width of the distribution of $`\mathrm{Re}(\rho _\varphi )`$ becomes much smaller below the transition. At low temperatures structural arrest hinders the exploration of phase space and reduces the fluctuations in Re$`(\rho _\varphi )`$.

The other way to calculate the generalized linear compressibility is with global averaging in which we string together a series of separate runs, treat it as one giant run, and then do the averaging required to calculate the generalized linear compressibility $`\chi _l^{\mathrm{global}}`$.

$$\chi _l^{\mathrm{global}}=\frac{6}{N}\left\{\left[\rho _\varphi ^2\right]_{\mathrm{run}}\left(\left[\rho _\varphi \right]_{\mathrm{run}}\right)^2\right\}$$
(34)
Careful inspection of eqs. (33) and (34) reveals that the difference is in whether we square the thermal average $`\rho _\varphi `$ and then average over runs to obtain $`\left[\rho _\varphi ^2\right]_{\mathrm{run}}`$ in $`\chi _l`$, or average $`\rho _\varphi `$ over runs and then square it to obtain $`\left(\left[\rho _\varphi \right]_{\mathrm{run}}\right)^2`$ in $`\chi _l^{\mathrm{global}}`$. If we now return to Figure 25 and take the difference of $`[S(k)]_{\mathrm{run}}`$ and $`\left([\mathrm{Re}<\rho _k>]_{\mathrm{run}}\right)^2`$, we obtain the global average:

$$\chi _l^{\mathrm{global}}=3[S(k)]_{\mathrm{run}}\frac{6}{N}\left([\mathrm{Re}<\rho _k>]_{\mathrm{run}}\right)^2$$
(35)
The result of both types of averaging is shown in Figures 25 and 26. Notice that $`\chi _l^{\mathrm{global}}`$ does not exhibit a drop with decreasing temperature while $`\chi _l`$ does. To understand why there is no drop in $`\chi _l^{\mathrm{global}}`$, note that by combining several different runs, very different configurations are sampled which produces much larger fluctuations in the generalized center of mass at low temperatures compared to $`\chi _l`$. As a result $`\chi _l^{\mathrm{global}}`$, which is a measure of the size of these fluctuations, does not have an abrupt drop.

## IX Summary

To summarize, we have introduced a new quantity called the generalized compressibility which depends solely on the positions of the particles and not on their histories. The generalized compressibility can easily be calculated in the canonical (e.g., $`NVT`$) and grand canonical ensembles. In particular it is well defined in a system which has particle number and volume fixed. In addition it does not suffer from the finite size effects often encountered in calculating the ordinary compressibility. The linear generalized compressibility drops abruptly at the observed glass transition due to the kinetic arrest of motion. This makes it an good quantity to calculate or measure in order to find the observed glass transition as a function of density or temperature. The generalized compressibility can be experimentally measured in several ways. It can be directly measured in colloidal experiments which monitor the positions of the particles. Measurements of the width of the distribution of $`\rho _\stackrel{}{k}`$, the spatial Fourier transform of the density, would also yield the linear generalized compressibility.

We thank Sharon Glotzer, Walter Kob, Andrea Liu and Francesco Sciortino for helpful discussions. This work was supported in part by DOE grant DE-FG03-00ER45843 as well as by CULAR funds provided by the University of California for the conduct of discretionary research by Los Alamos National Laboratory.

## X Appendix: Parallel Tempering

In calculating the intermediate scattering function at a given temperature, we initialized the run using a configuration at that temperature generated by parallel tempering. In this appendix we briefly describe the parallel tempering method.

We implement parallel tempering (PT) Hukushima96 ; Hukushima98 ; Kob01 ; Michele01  by choosing the temperatures at which we wish to have measurements made. We then run molecular dynamics simulations in parallel at these temperatures using a temperature constraint algorithm Rapaport95  to keep the temperature of each simulation constant. At 100 time step intervals we attempt to switch the configurations of two neighboring temperatures using a Metropolis test which ensures that the energies of the configurations sampled at any given temperature have a Boltzmann distribution. Let $`\beta _1`$ and $`\beta _2`$ be two neigboring inverse temperatures, and let $`U_1`$ and $`U_2`$ be the corresponding potential energies of the configurations at these temperatures at a time step just before the possible swap. If $`\mathrm{\Delta }=(\beta _1\beta _2)(U_2U_1)`$, then the switch is accepted with probability unity if $`\mathrm{\Delta }0`$ and with probability $`\mathrm{exp}(\mathrm{\Delta })`$ if $`\mathrm{\Delta }>0`$. The temperatures are chosen so that the acceptance ratio is between 30% and 75%. At the temperatures in the vicinity of the mode coupling $`T_C`$, the acceptance ratio was typically above 75% for L=6 and above 60% for L=8. After a swap is accepted, the velocities of the particles in each configuration are rescaled to suit their new temperature. Each configuration is then evolved using molecular dynamics for another 100 time steps. Switching configurations allows a given simulation to do a random walk in temperature space in which it visits both low temperatures and high temperatures. This helps to prevent a simulation from becoming trapped in a valley of the energy landscape at low temperatures. Typically we equilibrate for 2 million time steps and then do measurements for an additional 4 million time steps.

Present address: Internap, Seattle, WA 98101.


text:
# Circuit Theory for Full Counting Statistics in Multi-Terminal Circuits

## Abstract

We propose a theory that treats the current, noise, and, generally, the full current statistics of electron transfer in a mesoscopic system in a unified, simple and efficient way. The theory appears to be a circuit theory of $`2\times 2`$ matrices associated with Keldysh Green functions. We illustrate the theory by considering the big fluctuations of currents in various three-terminal circuits.

The field of quantum noise in mesoscopic systems has been exploded during the last decade, most achievements being summarized in a recent review article. BlanterReview  Measurement of fractional charge in Quantum Hall regime Frac\_Charge , noise measurements in atomic-size junctions Cron  and superconductors Kozhe  are milestones of the field and demonstrate the importance of quantum noise as a unique tool to study electron correlations and entanglements of different kinds. A very important step has been made in Levitov  where an elegant theory of full counting statistics (FCS) has been presented. This theory encompasses not only noise, but all higher momenta of the charge transfer.

Starting from pioneering work of BüttikerBuettiker , a special attention has been paid to noise and statistics of electron transfer in multi-terminal circuits. The correlations of currents flowing to different terminals reveal Fermi statistics of electrons. These cross-correlations have been recently observed. cross-corr  Although the noise correlations for several relevant layouts have been understood BlanterReview , the evaluation of FCS still encountered difficulties. For instance, an attempt to build up FCS with ”minimal correlation approach” BlanterSukhor  has lead to contradictions .BlanterSchomerus  This is unfortunate, since higher-order current correlations supply information about higher-order electron correlations and multi-particle interference. This information is of fundamental importance and can be hardly obtained by any other means.

In this letter, we present a calculational scheme that allows for easy evaluation of FCS in a multi-terminal mesoscopic systems. It is of a great intellectual enjoyment that this scheme is a simple and universal one. In fact, it is hardly more complicated than a conventional circuit theory of electric transport and is based on slight extension of Kirchoff rules to $`2\times 2`$ matrix structures.

We start by introducing current operators $`\widehat{I}_i`$, each being associated with the current to a certain terminal $`i`$. Extending the method of RefYuli  we introduce a Keldysh-type Green function defined by

$$\left(i\frac{}{t}\widehat{H}+\frac{1}{2}\overline{\tau }_3\underset{i}{}\chi _i(t)\widehat{I}_i\right)\stackrel{ˇ}{G}(t,t^{})=\delta (11^{})$$
(1)
Here we follow notations of a comprehensive review Rammer , $`\chi _i`$ are time-dependent parameters, $`\overline{\tau }_3`$ is a $`2\times 2`$ matrix in Keldysh space, $`\widehat{H}`$ is the one-particle Hamiltonian that incorporates all information about the system layout, including boundaries, defects and all kinds of elastic scattering. We use ”hat”, ”bar” and ”check” to denote operators in coordinate space, matrices in Keldysh space and operators in direct product of these spaces respectively. The Eq. 1 defines the Green function unambiguously provided boundary conditions are satisfied: $`\stackrel{ˇ}{G}(t,t^{})\overline{G}(x,x^{};t,t^{})`$ approaches the common equilibrium Keldysh Green functions $`\stackrel{ˇ}{G}_i^{(0)}(tt^{})`$ provided $`x,x^{}`$ are sufficiently far in the terminal $`i`$. These $`\stackrel{ˇ}{G}_i^{(0)}(tt^{})`$ incorporate information about the state of the terminals: their voltages $`V_i`$ and temperatures $`T_i`$.

One can easily see by traditional diagrammatic methods Rammer  that the expansion of $`\stackrel{ˇ}{G}`$ in $`\chi _i(t)`$ generates all possible diagrams for higher order correlators of $`\widehat{I}_i(t)`$ and thereby incorporates all the information about statistics of charge transfer. If we limit our attention to low-frequency limit of current correlations, we can keep time-independent $`\chi _i`$. In this case, the Green functions are functions of time difference only and the Eq. 1 separates in energy representation. It is convenient to introduce the following $`\chi _i`$-dependent action defined as a sum of closed diagrams:

$$\frac{S}{\chi _i}=it_0\frac{d\epsilon }{2\pi }\mathrm{Tr}\left(\overline{\tau }_3\widehat{I}_i\stackrel{ˇ}{G}(\epsilon )\right)$$
(2)
This allows us to express the probability for $`N_i`$ electrons to be transferred to the terminal $`i`$ during time interval $`t_0`$

$$P(\{N_i\})=_\pi ^\pi \underset{i}{}\frac{d\chi _i}{2\pi }e^{S(\{\chi _i\})i_iN_i\chi _i}.$$
(3)
(Higher-order) derivatives of $`S`$ with respect to $`\chi _i`$ give (higher-order) moments of $`P(\{N_i\})`$. First derivatives yield average currents to terminals, second derivatives correspond to the noises and noise correlations.

Using special properties of current operators, $`\chi `$-dependent terms in Eq. 1 can be gauged away.RefYuli ; Wolfgang  The $`\chi `$ dependence of $`\stackrel{ˇ}{G}`$ is thereby transferred to the boundary conditions: the gauged Green function far in each terminal shall approach $`\stackrel{ˇ}{G}_i(ϵ)`$ defined as

$$\stackrel{ˇ}{G}_i(ϵ)=\mathrm{exp}(i\chi _i\overline{\tau }_3/2)\stackrel{ˇ}{G}_i^{(0)}(ϵ)\mathrm{exp}(i\chi _i\overline{\tau }_3/2)$$
(4)
In the present form, the Eq. (1) with relations (4,2) solves the problem of determination of the FCS for any arbitrary system layout: one just has to find exact quantum-mechanical solution of a Green function problem. This is hardly constructive, and we proceed further by deriving a simplified semiclassical approach. First, we note that even in its exact quantummechanical form the Eq. 1 possesses an important property. We consider the quantity defined similar to standard definition of current density, $`\overline{j}^\alpha (x,ϵ)lim_{xx^{}}(^\alpha ^\alpha )\overline{G}(x,x^{};ϵ)/m`$. By virtue of Eq. 1 this quantity conserves so that

$$\frac{\overline{j}^\alpha (x,ϵ)}{x^\alpha }=0$$
(5)
This looks like the conservation of particles at a given energy. However, this relation contains more information since it is a conservation law for a $`2\times 2`$ matrix current.

Next, we construct a theory which makes use of this conservation law. We concentrate on the Green function in coinciding points, $`\overline{G}(x,ϵ)i\overline{G}(x,x^{};ϵ)/\pi \nu `$. So defined Green function has been introduced in several semiclassical theories. Rammer ; Larkin ; Noneq  It satisfies the normalization condition $`\overline{G}^2=\overline{1}`$. We relate the ”current density” $`\overline{j}`$ to gradients and/or changes of $`\overline{G}(x)`$, very much like the electric current density is related to the voltage in circuit theory of electric conductance. Following the approach of the circuit theory, we separate a mesoscopic layout into elements: nodes and connectors, so that the $`\overline{G}(x)`$ is constant across the nodes and drops across the connectors. One may associate a graph with each circuit, so that its lines $`(i,j)`$ would denote the connectors, and internal and external vertices correspond to the nodes and terminals, respectively. (See Fig. 1) This separation of actual layout is rather heuristic, similar to separation of an electric conductor of a complicated geometry onto nodes and circuit theory elements. The bigger the number and the finer the mesh of the nodes and connectors, the better the circuit theory approximates the actual layout. The nodes are similar to the terminals, the difference is that $`\overline{G}`$ is fixed in the terminals and yet to be determined in the nodes. The $`\overline{G}`$ in nodes are determined from Kirchoff rules reflecting the conservation law (5): sum of the matrix currents from the node over all connectors should equal zero at each energy. For this, we should be able to express the matrix current via each connector as a function of two matrices $`\overline{G}_{i,j}`$ at its ends.

The connector $`(i,j)`$ can be quite generally characterized by a set of transmission eigenvalues $`T_n^{(ij)}`$.General ; Noneq  The problem to solve is to express matrix current via the connector in terms of $`\overline{G}_{i(j)}`$. This problem shall be addressed by a more microscopic approach and was solved in Noneq  for Keldysh-Nambu matrix structure of $`\stackrel{ˇ}{G}`$. It is a good news that the derivation made in Noneq  does not depend on concrete matrix structure and can be used for the present problem without any modification yielding

$$\overline{I}_{ij}=\frac{1}{2\pi }\underset{n}{}𝑑E\frac{T_n^{(ij)}[\overline{G}_i,\overline{G}_j]}{4+T_n^{(ij)}\left(\{\overline{G}_i,\overline{G}_j\}2\right)}.$$
(6)
Each connector $`(i,j)`$ in the layout contributes to the total $`\chi _i`$-dependent action (2). The corresponding $`S_{ij}`$ contribution reads:Wolfgang

$$S_{ij}(\chi )=\frac{t_0}{2\pi }\underset{n}{}𝑑E\text{Tr}\mathrm{ln}\left[1+\frac{1}{4}T_n^{(ij)}\left(\{\overline{G}_i,\overline{G}_j\}2\right)\right].$$
(7)
Now we are ready to present a set of circuit theory rules that enables us to evaluate the FCS for an arbitrary mesoscopic layout. i. The layout is separated onto terminals, nodes, and connectors. ii. The $`\overline{G}_j`$ in each terminal $`j`$ is fixed by relation (4) thus incorporating information about voltage, temperature and counting field $`\chi `$ in each node. iii. For each node $`k`$, the matrix current conservation yields a Kirchoff equation $`_i\overline{I}_{ik}=0`$, where the summation is going over all connectors $`(i,k)`$ attached to node $`k`$, and $`\overline{I}_{ik}`$ are expressed with (6) in terms of $`\overline{G}_{i(k)}`$. iv. The solution of resulting equations with condition $`\overline{G}_k^2=1`$ fixes $`\overline{G}_k`$ in each node. v. The total action $`S(\chi )`$ is obtained by summing up the contributions $`S_{ij}(\{\chi _i\})`$ of individual connectors, those are given by (7): $`S(\{\chi _i\})=_{(i,j)}S_{ij}(\{\chi _i\})`$ vi. The statistics of electron transfer is obtained from the relation (3).

It is time to discuss the limits of applicability of the whole scheme. By virtue of semiclassical approach, the mesoscopic fluctuations coming from interference of electrons penetrating different connectors are disregarded. So that, we assume that conductivities of all connectors are much bigger than conductance quantum $`e^2/\pi \mathrm{}`$. The same condition provides the absence of Coulomb blockade effects in the system. Besides, we have disregarded the possible processes of inelastic relaxation in the system. The latter can be eventually taken into account but it would considerably complicate the scheme. The point is that the inelastic scattering would mix up the $`\overline{G}`$ at different energies, so that one can not solve the circuit theory equations separately at each energy.

As an illustration of the presented scheme, we will consider in the rest of the paper the FCS of the 3-terminal chaotic quantum dot. The system is sketched in the inset of Fig. 2. The heuristic circuit, associated with this mesoscopic system is shown by dashed lines. It includes only 3 external vertices, corresponding to terminals, 3 arbitrary connectors, associated with the contacts, and the node $`\{4\}`$, representing the quantum dot itself. This separation is valid provided the cavity is in the quantum chaotic regime. (See Aleiner  for definition). This regime corresponds to full isotropization of the Green function $`\stackrel{ˇ}{G}(x,x^{},ϵ)`$ within the dot, so that $`\overline{G}_4(ϵ)`$ can be regarded as a constant at a given energy.

Since the normalization $`\overline{G}_k^2=1`$ holds for each vertex, we use the parametrization $`\overline{G}_k=𝐠_k𝝉`$, $`𝐠_k𝐠_k=1`$. Here $`𝐠_k`$ is a 3-D vector, and $`𝝉=(\overline{\tau }_1,\overline{\tau }_2,\overline{\tau }_3)`$. In the absence of counting fields the Green functions in the terminals are given by a zero condition $`\overline{G}_k^{(0)}=\left(\begin{array}{cc}12f_k& 2f_k\\ 2(1f_k)& 2f_k1\end{array}\right),`$ where Fermi distribution function $`f_k(E)=\{\mathrm{exp}[(EeV_k)/T_k]+1\}^1`$ accounts for the bias voltages $`V_k`$ and the temperatures $`T_k`$ in the terminals. The $`\chi _i`$-dependence of $`\overline{G}_k(\chi )`$ is then given by Eq. 4.

We see that $`\overline{G}_4(\chi )=𝐠_4𝝉`$ is in fact the only function to be found. For that, we proceed by applying the current conservation law, $`_{k=1}^3\overline{I}_{k,4}=0`$, inside the dot. We present the currents $`\overline{I}_{k,4}`$ given by (6) in the form $`\overline{I}_{k,4}=\frac{1}{2}Z_k(𝐠_k𝐠_4)[\overline{G}_k,\overline{G}_4]`$, the scalar function $`Z_k(x)`$ incorporating the information about transmission eigenvalues in each connector $`k`$: $`Z_k(x)\underset{n}{}T_n^{(k,4)}/[2+T_n^{(k,4)}(x1)].`$ It can be evaluated for any particular distribution $`\rho (T)`$ of transmission eigenvalues in the given connector and completely defines its scattering properties. For a example, if we denote $`R_0=\pi \mathrm{}/e^2`$, then $`R_k^1=2R_0^1Z_k(1)`$ is an inverse resistance of the connector. One can also express the Fano factor $`F=T(1T)/T`$ as $`F=12(d/dx)\mathrm{log}Z(x)|_{x=1}`$. With the use of $`Z_k(x)`$ the conservation law can be efficiently rewritten as $`[\underset{k=1}{\overset{3}{}}p^k\overline{G}_k,\overline{G}_4]=0`$, where $`p^k=Z_k\left(𝐠_k𝐠_4\right)`$. The latter enables one to look for the vector $`𝐠_4`$ in the form $`𝐠_4=M^1\underset{k=1}{\overset{3}{}}p^k𝐠_k,`$ with $`M(\chi )`$ being an unknown normalization constant. Using the condition $`𝐠_4𝐠_4=1`$ we obtain the set of equations

$$p^i=Z_i\left(M^1\underset{j=1}{\overset{3}{}}g_{ij}(\chi )p^j\right),M^2=\underset{i,j=1}{\overset{3}{}}g_{ij}(\chi )p^ip^j$$
(8)
where

$`g_{ij}(\chi )`$ $`=`$ $`𝐠_i(\chi )𝐠_j(\chi )=(12f_i)(12f_j)`$
$`+`$ $`2e^{i(\chi _i\chi _j)}f_i(1f_j)+2e^{i(\chi _i\chi _j)}f_j(1f_i)`$
The fixed point $`\{p_{}^i(\chi ),M_{}(\chi )\}`$ of the mapping (8) finally yields $`\overline{G}_4`$ in question.

The total action can be found by applying the rule (v) at page 1 and reads

$$S(\chi )=\frac{t_0}{\pi }\underset{i=1}{\overset{3}{}}𝑑ϵS_i\left(M_{}^1(\chi )\underset{j=1}{\overset{3}{}}g_{ij}(\chi )p_{}^i(\chi )\right)$$
(9)
The partial contributions $`S_i(x)`$ in the above equation should be determined from the relation $`\frac{}{x}S_i(x)=Z_i(x)`$, $`S_i(1)=0`$.

We specifically consider three particular types of connectors: tunnel(T), ballistic(B) and diffusive(D). The corresponding contributions to action are: $`S_T(x)=\frac{1}{2}(R_0/R)(x1),`$ $`S_B(x)=(R_0/R)\mathrm{log}[(1+x)/2]`$, $`S_D(x)=\frac{1}{4}(R_0/R)\mathrm{log}^2(x+\sqrt{x^21})`$ RefYuli , $`R`$ being a resistance of the connector. For tunnel connector $`T_n1`$ for all $`n`$. For ballistic connector $`N`$ channels are opened ($`T_n1`$ for $`nN`$), and the other are closed. In the diffusive connector the transmission eigenvalues are distributed according to universal law $`\rho (T)=R_0/2RT\sqrt{1T}`$.

Analytical results for FCS (9) are plausible only for the system with tunnel connectors. To assess general situation we solve the Eqs. (8) for given $`\chi _i`$ numerically. To find the probability distribution, we evaluated the integral (3) in the saddle point approximation, assuming $`\chi _i`$ to be complex numbers. Saddle point approximation is always valid in the low frequency limit we consider, since in this case both action $`S`$ and number of transmitted particles $`N_i=I_it_0/e1`$. Due to the current conservation law only two of three counting fields $`\chi _i`$ are independent, and one can set $`\chi _3=0`$. The relevant saddle point of the function $`\mathrm{\Omega }(\chi )=S(\chi )+i\chi _1I_1t_0/e+i\chi _2I_2t_0/e`$ always corresponds to purely imaginary numbers $`\{\chi _1^{},\chi _2^{}\}`$ in the upper half plane. The probability reads $`P(I_1,I_2)\mathrm{exp}[\mathrm{\Omega }(\chi ^{})]`$. Evidently, $`\mathrm{\Omega }(\chi ^{})`$ is the Legendre transform of the action, and it can be regarded as implicit function on $`I(\chi ^{})`$.

In the following we assume the short noise regime $`eVkT`$ when the thermal fluctuations can be disregarded. The energy integration in (9) becomes trivial, since $`f_i(ϵ)=0`$ or 1, and it is sufficient to consider only the case $`V_1=V_2=0`$, $`V_3=V`$. Any other possible setup can be reduced to the number of previous ones by appropriate subdividing a relevant energy strip. The results of these calculations are shown in Fig.2 and 3. We see that the maximum of probability occurs at $`I_1=I_2=V/3R`$, $`I_3=2V/3R`$ that simply reflects the usual Kirchoff rules. The current distribution $`P(I_1,I_2)`$ for a ballistic system is bounded. It is due to the fact, that $`Z_B(x)`$ contains the finite number of open channels, contrary to the tunnel or diffusive type configurations, where it is not the case. From (8) and (9) we can also find a zero noise and noise correlations matrix $`\stackrel{~}{S}_{ij}=eR^1VF_{ij}`$, $`F_{11}=F_{22}=(4+3F)/27`$, $`F_{33}=(4+6F)/27`$, $`F_{12}=2/27`$, $`F_{13}=F_{23}=(2+3F)/27`$, where $`F`$ is a Fano factor. Since $`F_B=0`$, $`F_D=1/3`$ and $`F_T=1`$ one concludes, that for a fixed average currents through connectors the Gaussian’s currents fluctuations will increase in the sequence ballistic$``$diffusive$``$tunnel. Fig.2 and 3 show, that the similar behavior is also traced in the regime of the large current fluctuations. The essential point here is that the cross-correlations always persist regardless the concrete construction of the connectors. For the case of multilead chaotic cavities the results for shot noise in our theory coincides with those, obtained by means of random matrix theory Langen , and with use of ”minimal correlation approach” BlanterSukhor .

In conclusion, we present a constructive theory for counting statistics for electron transfer in mesoscopic systems. With this theory, one can easily make theoretical predictions for all FCS, thereby facilitating experimental activities in this direction. Up to now, only the noise has been measured. In our opinion, the measurements of FCS can be easily performed with threshold detectors that produce a signal provided the current in a certain terminal exceeds the threshold value. If the threshold value exceeds the average current, the detector will be switched by relatively improbable fluctuations of the current described by FCS, its signal being proportional to the probability $`P(I_1,I_2)`$ of these fluctuations.

We appreciate useful discussions with G. E. W. Bauer, C. W. J. Beenakker, L. P. Kouwenhoven and L. S. Levitov. This work is a part of the research program of the ”Stichting voor Fundamenteel Onderzoek der Materie” (FOM), and I acknowledge the financial support from the ”Nederlandse Organisatie voor Wetenschappelijk Onderzoek” (NWO).