text:
# Quantum Dot in the Kondo Regime coupled to p-wave superconductors

## I introduction

The Kondo effect is one of the simplest manifestations of many-body physics exhibiting strong correlationshewson . Originally, it has been introduced in order to explain the occurrence of a shallow minimum of the resistivity $`\rho (T)`$ of bulk metals doped with a small concentration of magnetic impurities. It has been confirmed that at low temperatures (below the Kondo temperature $`T_K`$) scattering from magnetic impurities is strong, leading to an enhanced resistivity. It is now well established that at this strong coupling regime the theory of perturbation fails due to logarithmic divergences of higher (than second) order diagrams.

The Kondo physics also plays an important role in electron transport through quantum dots, where instead of a magnetic impurity one now encounters localized electrons glazman88 ; hershfield91 ; ournca . Inspection of the dot-electron Green function reveals an enhancement of the local density of states at the Fermi energy (the Abrikosov-Suhl resonance). Thus, unlike the case of Kondo effect in bulk metals, its hallmark experimental manifestation in quantum dots is an enhancement of the conductance. An important milestone in this field has been recorded by recent observations of the Kondo effect in transport through quantum dotsgoldhaber ; reviews  and in a finite length carbon nanotubes lin . These experiments have paved the way for a new class of experimental investigations into the physics of strongly correlated electrons in general and the Kondo effect in particular. Together with more recent experiments unitary  demonstrate the feasibility of exploiting tunable physical parameters of a quantum dot system in order to yield important information on Kondo systems and other many-body related phenomena. Such studies, for example, encompass the entire crossover region between the Kondo limit, the mixed valence regime and the non-Kondo (weak-coupling) domain.

While the Kondo physics in a quantum dot attached on both its sides to normal ($`N`$) metallic leads received much recent attention, it has been realized that novel physical effects emerge if (one or both) electrodes attached through a tunneling barrier to a quantum dot in the Kondo regime ($`K`$) is a superconductor ($`S`$). The central electron-transport mechanism in such $`SKS`$ and $`SKN`$ junctions is that of Andreev reflections when two particles tunnel together coherently to form a Cooper pair in the superconductor. In fact, the physics of $`SKS`$ and $`SKN`$ junctions is determined by the interplay of these Andreev reflections with the Kondo-formed resonance in the spectral density of states of the dot electron BTK . This interplay between the Kondo effect and superconductivity in the sub-gap region is rather effective here. We then expect the the $`IV`$ characteristics in the sub-gap region to be most important. Andreev reflections play a dominant role both in and out of equilibrium. In the former case they are responsible for occurrence of direct Josephson current while in the latter case they are the cause of dissipative currents at sub-gap voltages as well as for suppressing the zero bias anomaly in $`NS`$ junctions. The effect becomes even more significant when the effective Kondo temperature $`T_K`$ exceeds the superconducting energy gap $`\mathrm{\Delta }`$. A natural candidate for such $`SKN`$ junction is that in which the role of the quantum dot or the Kondo impurity is played by a Carbon nano-tube $`(CNT)`$, where relatively high values of $`T_K`$ can be achieved. Indeed, Kondo effect with rather high Kondo temperature $`(T_K>1K`$) has been reported in $`NCNTN`$ junctions with a few conducting channels lin . Moreover, superconducting junctions with a weak link formed by $`CNT`$ have already been fabricated Kas ; Christ .

Let us illuminate the somewhat subtle distinction between $`swave`$ and $`pwave`$ $`SKN`$ junctions. In the case of an $`SKN`$ junction when one electrode is an $`swave`$ superconductor and the Kondo impurity is weakly coupled to the $`S`$ and $`N`$ electrodes, is determined by the competition between two important effects. The first one is the formation of a Kondo singlet which screens the bare impurity spin and drives the system toward the unitary limit at very low temperatures. The second one is the existence of the superconducting gap which implies a vanishingly small density of low energy electron states fazio ; raimondi ; ambe . These are precisely the electron states which are needed in order to screen the Kondo impurity. Note that if the Kondo temperature is larger than the superconducting gap ($`T_K>\mathrm{\Delta }`$) then electron states outside the gap can participate in the screening interaction. However, this condition cannot be easily realized. Consider, on the other hand, an $`SKN`$ junction in which the $`S`$ electrode consists of an unconventional superconductor which is oriented relative to the interface in such a way that the pair potential reverses its sign on the Fermi surface. In this case, mid-gap zero-energy states $`(ZES)`$ are formed which are localized near the surface of the unconventional superconductor. These states can now participate in screening the impurity spin through the Kondo effect and emergence of sub-gap current is expected. In the experimentally feasible setup of $`SCNTN`$ where just a few tunneling channels are present, charge is carried mainly by quasiparticles moving perpendicular to the interface. This restricts the possible values of the angle $`\theta `$ between the superconducting surface and the direction of the injected quasiparticles tanaka . Formation of $`ZES`$ is possible when $`\mathrm{\Delta }(\theta )=\mathrm{\Delta }(\pi \theta )`$. If the impurity is almost point-like, the relevant injection angle is of course $`\theta =0`$. Unconventional superconductors for which the symmetry of the pair potential is that of a triplet ($`pwave`$ superconductors) satisfy the above condition. These have recently been discovered by Maeno, et al maeno  in $`Sr_2RuO_4`$. For impurities of finite extent, one may also consider formation of $`ZES`$ in $`dwave`$ superconductors. Thus, the physics of $`SKN`$ junctions with an $`S`$ electrode whose order parameter has a non-trivial symmetry is affected by the formation of $`ZES`$ in the Kondo regime.

Beyond investigating the conductance dependence on the applied bias we propose and explore theoretically other novel experimental tools suitable for probing the Kondo regime. These include shot-noise measurements and the Josephson (direct) current. Like the conductance, we will calculate and analyze shot-noise power spectrum in $`SKN`$ junctions at very low temperature where one lead is either an $`swave`$ superconductor (mid-gap surface states are absent) or a $`pwave`$ superconductor (mid-gap surface states are present). We will also study the temperature dependence of the Josephson current in these junctions.

In order to study this strongly interacting non-equilibrium problem in the framework of workable approximation, we restrict ourselves to voltages that do not exceed the superconducting gap so that the Kondo temperature can satisfy the inequality $`T_K>max(\mathrm{\Delta },eV)`$. Under this condition we can apply the well known approach familiar in the physics of the Anderson impurity Hamiltonian bickers , namely the mean field slave boson approximation (MFSBA) (see below). It consists of a workable scheme for calculating the conductance and the zero frequency shot-noise power. (Here we also apply the MFSBA to calculate Josephson current in $`SKS`$ junctions.)

The combined effect of multiple Andreev reflections (MAR) and electron-electron interactions on the shot-noise in $`SKS`$ junctions was studied in Ref. us1, ; compagnano, ; arovas2, . It was shown that in the strong coupling Kondo regime the effective transmission is enhanced, and a ballistic-like channel opens up. Hence, an interplay between MAR and the Kondo resonance yields an excess current at zero bias and the $`IV`$ curve behaves similarly to that of noninteracting ballistic junctions. At very high values of $`T_K`$ and in the low voltage limit the current approaches the noninteracting value $`I_{AR}=4e\mathrm{\Delta }/h`$. Analogously, the shot-noise power is shown to display a pronounced maximum at $`V=0`$ and a decay as $`1/V`$ at small bias,familiar in the standard theory of noninteracting $`SNS`$ junctions averin .

In section II the Hamiltonian of $`SKN`$ and $`SKS`$ junctions is defined, and the bare action is introduced. After integrating out the fermion fields pertaining to the superconducting leads an effective action is obtained, depending solely on the dot variables. Calculations and presentations of conductance, shot-noise power and Josephson current are respectively detailed in sections III, IV and V.

## II Model Hamiltonian and effective action

The dynamics of systems like $`SKN`$ or $`SKS`$ junctions is governed by the Hamiltonian

$`H=H_L+H_R+H_d+H_t+H_c,`$ (1)
in which $`H_j`$, ($`j=L,R`$) are the Hamiltonians of the electrodes which depend on the electron field operators $`\psi _{a\sigma }(𝐫,t)`$ where $`𝐫=(x,y)`$ and $`\sigma =\pm `$ is the spin index.

$`𝑯_j={\displaystyle 𝑑r[\mathrm{\Psi }_{j\sigma }^{}(𝒓)\xi ()\mathrm{\Psi }_j\gamma \mathrm{\Psi }_j^{}(𝒓)\mathrm{\Psi }_j^{}\mathrm{\Psi }_j(𝒓)\mathrm{\Psi }_j(𝒓)]}.`$ (2)
Here $`\gamma `$ is the BCS coupling constant and $`\xi ()=^2/2m\mu `$ with $`\mu `$ being the chemical potential at temperature $`T`$. The Planck constant is set $`\mathrm{}=1`$ and whenever appropriate, the spin, space and time dependence of all the field operators will not be explicitly displayed. As in Refs. us1, ; us,  the dot is represented by a single level Anderson impurity with energy $`ϵ_0<0`$ and Hubbard repulsion parameter $`U`$. In the Kondo regime of interest here we set $`U\mathrm{}`$ and assume $`|ϵ_0|`$ to exceed any other energy scale except $`U`$. In this case it is convenient to express the dot and the tunneling Hamiltonians $`H_d`$ and $`H_t`$ via slave boson (operators $`b,b^{}`$) and slave fermion (operators $`c,c^{}`$) auxiliary fieldscoleman . Explicitly,

$`H_d=ϵ_0{\displaystyle \underset{\sigma }{}}c_\sigma ^{}c_\sigma ,`$
$`H_t={\displaystyle \underset{j\sigma }{}}𝒯_jc_\sigma ^{}b\psi _{j\sigma }(\mathrm{𝟎},t)+h.c.,`$ (3)
where $`𝒯_j`$ is the tunneling amplitude. Finally, the Hamiltonian of the system must also include a term which prevents double occupancy in the limit $`U\mathrm{}`$. This term reads,

$`H_c=\lambda ({\displaystyle \underset{\sigma }{}}c_\sigma ^{}c_\sigma +b^{}b1),`$ (4)
where $`\lambda `$ is a Lagrange multiplier.

Following Ref. us1  let us consider the dynamical “partition function”

$`Z{\displaystyle 𝒟[F]\mathrm{exp}(i𝒮)},`$ (5)
where the path integral is carried out over all fields $`[F]`$ and the action $`𝒮`$ is obtained by integrating the Lagrangian pertaining to the Hamiltonian (1) along the Keldysh contour. In performing the functional integrations the boson field operators are treated as $`c`$-numbers. As a result one arrives at an effective action expressed in terms of the Green functions of the leads.

$`S_{\mathrm{eff}}=i\mathrm{Tr}\mathrm{ln}\widehat{G}^1{\displaystyle 𝑑t[\widehat{\lambda }\sigma _z(\widehat{b}\widehat{b}1)]}.`$ (6)
Here $`\widehat{\lambda }=(\lambda _1,\lambda _2)`$, $`\widehat{b}=(b_1,b_2)`$ and $`\sigma _z`$ are diagonal matrices acting in Keldysh space. The inverse propagator $`\widehat{G}^1`$ depends on the Green functions of the electrodes.us, . Performing the standard basis rotation in Keldysh space one finds,

$`\widehat{G}^1(ϵ,ϵ^{})=\delta (ϵϵ^{})(ϵ\tau _z\stackrel{~}{ϵ}{\displaystyle \frac{\mathrm{\Gamma }b^2}{2}}\tau _z\widehat{g}_+(ϵ)\tau _z),`$ (7)
where $`\stackrel{~}{ϵ}=ϵ_0+\lambda `$ is the renormalized level position (in the Kondo limit one has $`\stackrel{~}{ϵ}0`$) and $`\mathrm{\Gamma }=(\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)/2𝒯_{,}^{}{}_{}{}^{2}`$ is the usual transparency parameter. The $`2\times 2`$ matrix representation (in Keldysh space) for $`g`$ is composed of diagonal elements $`\widehat{g}^{R/A}(ϵ)`$ and an upper off-diagonal element $`\widehat{g}^K(ϵ)=(\widehat{g}^R\widehat{g}^A)th(ϵ/2T)`$. Here and below we define

$`\widehat{g}_\pm =\gamma _L\widehat{g}_L\pm \gamma _R\widehat{g}_R,`$ (8)
with asymmetry parameters $`\gamma _j=\mathrm{\Gamma }_j/\mathrm{\Gamma }`$. The matrix $`\widehat{g}_R`$ has the standard structure with retarded and advanced superconductor Green functions which in the $`swave`$ case reads,

$`\widehat{g}^{R/A}(ϵ)=i{\displaystyle \frac{(ϵ\pm i0)+|\mathrm{\Delta }|\tau _x}{\sqrt{(ϵ\pm i0)^2|\mathrm{\Delta }|^2}}},`$ (9)
while for the $`pwave`$ case with incidence angle $`\alpha `$ it reads,

$`\widehat{g}^{R/A}(ϵ)=\widehat{g}_1^{R/A}+\widehat{g}_2^{R/A}\tau _x,`$
with (11)
$`\widehat{g}_1^{R/A}={\displaystyle \frac{i\sqrt{(ϵ\pm 0)^2|\mathrm{\Delta }|^2}\mathrm{cos}\alpha ϵ\mathrm{sin}\alpha }{ϵ\mathrm{cos}\alpha +i\sqrt{(ϵ\pm 0)^2|\mathrm{\Delta }|^2}\mathrm{sin}\alpha }}`$
$`\widehat{g}_2^{R/A}={\displaystyle \frac{i|\mathrm{\Delta }|}{ϵ\mathrm{cos}\alpha +i\sqrt{(ϵ\pm 0)^2|\mathrm{\Delta }|^2}\mathrm{sin}\alpha }}.`$
Note that the Pauli matrices $`\tau _{x,y,z}`$ act in Nambu space. For a one channel system one can take $`\alpha `$=0 and simplify the expression for the Green functions,

$`\widehat{g}_1^{R/A}={\displaystyle \frac{i\sqrt{(ϵ\pm 0)^2|\mathrm{\Delta }|^2}}{ϵ\pm 0}}`$
$`\widehat{g}_2^{R/A}={\displaystyle \frac{i|\mathrm{\Delta }|}{ϵ\pm 0}}.`$ (12)
In the case of $`NKS`$ junction, the left lead is represented by the Keldysh Green’s function of a normal metal $`\widehat{g}_L`$: $`\widehat{g}_L^{A/R}=\pm i`$ ; the Keldysh Green function $`\widehat{g}_L^K`$ has only diagonal matrix elements equal to $`2ith[(ϵ\pm eV)/2T]`$. Performing the variation of the effective action with respect to the fields $`b`$ and $`\lambda `$ a couple of self-consistency equations are obtained that determine these fields. In order to explicitly write down these self-consistency equations let us introduce the bare Kondo temperature $`T_K^0=Dexp[\pi ϵ_0|/(2\mathrm{\Gamma })]`$ and define a parameter $`X`$ by $`\mathrm{\Gamma }b^2=T_K^0X`$, where $`D`$ is the energy bandwidth. Then the MFSBA equations take the form

$`X={\displaystyle \frac{i\mathrm{\Gamma }}{2T_K^0}}\mathrm{Tr}\widehat{G}^K\tau _z,`$ (13)
$`\lambda ={\displaystyle \frac{i\mathrm{\Gamma }}{8}}\mathrm{Tr}[\widehat{G}^K\tau _z(\widehat{g}_+^R+\widehat{g}_+^A)+(\widehat{G}^R+\widehat{G}^A)\tau _z\widehat{g}_+^K]\tau _z,`$ (14)
where the trace includes energy integration as well. Eq. (13) effectively determines the Kondo temperature (through the parameter $`X`$), and reflects the constraint which prevents double occupancy in the limit $`U\mathrm{}`$. The second self-consistency equation (14) defines the renormalized energy level position $`\stackrel{~}{ϵ}`$. Let us briefly discuss the validity range of the present analysis. The MFSBA is known to encode the Kondo Fermi-liquid behavior at low temperatures. An important parameter here is the ratio between the Kondo temperature and the superconducting gap $`t_KT_K^0/\mathrm{\Delta }`$ fazio ; raimondi ; ambe . For $`t_K1`$ a Fermi liquid behavior is expected. Accordingly, in this regime Eq. (13) should have a nonzero solution $`X0`$ which corresponds to nonzero $`T_K`$. On the other hand, in the limit of large $`\mathrm{\Delta }`$ the only possible solution is the trivial one $`b=0`$ (and, hence, $`T_K=0`$). Quantitatively, the MFSBA is reliable only for sufficiently large values of $`t_K`$. We believe, however, that it can provide useful qualitative information also for moderate values of $`t_K`$ describing a crossover between the Kondo regime and the Coulomb blockade dominant domain us . It is worth noting here that, strictly speaking, the applied bias voltage $`V`$ attenuates the Kondo resonance and lowers $`T_K`$. Hence, for the reliability of the MFSBA in non-equilibrium situations, both $`\mathrm{\Delta }`$ and $`eV`$ should not exceed the Kondo temperature. Attention below is mainly focused on the sub-gap voltage regime $`eV\mathrm{\Delta }`$ in which case $`t_K`$ appears to be the only relevant parameter.

## III Conductance

The expression for the tunneling current obtained in Ref. us1  has a simple representation in terms of the dot Green function,

$`I=i{\displaystyle \frac{eXt_K}{8\mathrm{}}}\mathrm{Tr}[(\widehat{G}^R\tau _z\tau _z\widehat{G}^A)\widehat{g}_{}^K\widehat{G}^K\stackrel{~}{g}],`$ (15)
where for $`SKN`$ junctions we denote

$`\stackrel{~}{g}=\gamma _R(\widehat{g}^R\tau _z\tau _z\widehat{g}^A)2i\gamma _L\tau _z,`$
$`\stackrel{~}{g}^{R,A}=\gamma _R(\tau _z\widehat{g}^{R,A}+\widehat{g}^{R,A}\tau _z).`$ (16)
The last quantity will be used later on as it enters the expression for the noise power spectrum. Being combined with eqs. (13) and (14) the result (15) can be conveniently used for computing the transport current and the differential conductance of an $`SKN`$ junction in the Kondo regime for different values of $`t_K`$.

For sufficiently large $`t_K`$ we anticipate a strong Kondo resonance and the $`GV`$ curves for both $`swave`$ and $`pwave`$ superconductors is expected to resemble that of purely ballistic junctions (for $`\gamma _R=\gamma _L=1`$ ) without interactionBTK . Indeed, in the limit of large $`t_K1`$ which corresponds to the unitary, pure ballistic case, our expression for the current (15) reduces to that derived in Ref. BTK . This agreement is further supported by our numerical calculations carried out for $`t_K=100`$ (see Fig.1). Calculations of the current and conductance were also performed for $`t_K=5`$,$`3`$, and $`2`$. The results are displayed in figure 1. For voltages in the sub-gap region, the conductance of an S(p-wave)KN junction is distinct from that of an S(s-wave)KN junction. Such difference shows that the ZES supports the formation of a Kondo singlet for the lower values $`t_K=5`$, $`3`$ and $`2`$ and effectively turn the junction to be more transparent, approaching the BTK limit BTK .

Fig. 1 The conductance $`G`$ (in units of $`e^2/h`$) versus the bias $`V`$ (in units of $`\mathrm{\Delta }/e`$) for an $`swave`$ SKN (dash curves) and $`SpwaveKN`$ (solid lines) junction at sub-gap voltages with $`\mathrm{\Gamma }/T_K^0=200`$. The parameter $`t_K`$ =2,3,5,100 (from down to top). The upper line corresponding to $`t_K`$ =100 coincides for s and p wave superconducting leads.

Let us briefly summarize our results for the current and conductance of SKN junctions for $`s`$ and $`p`$ wave superconductors. In the limit of large $`t_K`$ the $`GV`$ curve is practically independent of $`t_K`$ and resembles that of a ballistic junction, as indicated in the upper curve in Fig.1. For lower values of $`t_K`$, junctions with s-wave superconductor lead are driven away from the unitary limit and the underlying physics becomes much richer. It reflects the influence of both $`\mathrm{\Delta }`$ and $`V`$ on the Kondo resonance and on the actual value of the Kondo temperature. For $`t_K=5`$ the $`GV`$ curve noticeably deviates from that obtained in the non-interacting limit. For $`t_K=2`$ the competition between gap-related suppression of the Kondo effect and the effective transparency of the junction becomes essential, leading to further decrease of the conductance. However, it is interesting to note that for and S(s-wave)KN junction at $`t_K=2`$ the conductance displays a small peak at the gap edge. Its interpretation is that the Kondo correlations strongly compete with superconductivity and influence the quasiparticle correlations when the energy exceeds the gap. Such an effect takes place even when $`t_K<1`$ (see ambe ). On the other hand, the S(p-wave)KN junction is less influenced to variations of $`t_K`$. This may be explained by the fact that due to the presence of ZES, superconductivity plays a minor role in the formation of the Kondo resonance and BTK behavior persists for smaller values of $`t_K`$.

## IV shot-noise

The shot-noise spectrum is usually defined as the symmetrized current-current correlation function lesovik ; rev

$`K(t_1,t_2)`$ $`=`$ $`\mathrm{}[I(t_1)I(t_2)I^2],`$ (17)
where $`\mathrm{}`$ denotes quantum averaging with the Hamiltonian (1). The general expression for the zero frequency shot-noise power in junctions with one (or even two) superconducting leads (in MFSBA approximation) was obtained in our recent work us1 . Like in the case of SKS junctions, it is useful to write $`K=(K_1+K_2)e^2\mathrm{\Delta }/(8\mathrm{})`$ for which the expressions derived are,

$`K_1={\displaystyle \frac{Xt_K}{2}}\mathrm{Tr}\{(\widehat{g}_+^R\widehat{g}_+^A)(\widehat{G}^R\widehat{G}^A)\widehat{g}_+^K\widehat{G}^K\},`$ (18)
$`K_2={\displaystyle \frac{(Xt_K)^2}{8}}\mathrm{Tr}\{(\widehat{G}^K\stackrel{~}{g})^22\tau _z\stackrel{~}{g}\tau _z\widehat{G}^A\stackrel{~}{g}\widehat{G}^R`$
$`[2\stackrel{~}{g}\widehat{G}^R\tau _z\widehat{g}_{}^K\widehat{G}^K+(\stackrel{~}{g}^R\widehat{G}^R)^2(\widehat{G}^A\widehat{g}_{}^K\tau _z)^2+\mathrm{h}.\mathrm{c}.]\}.`$ (19)
Expressions (18) and (19) (supplemented by the self-consistency eqs. (13) and (14)) are then solved numerically for the same set of parameters $`\mathrm{\Gamma }/T_K^0=200`$, $`t_K=100`$, 5, 3 and $`2`$. The results for the shot-noise power spectrum $`K`$ versus the applied voltage $`V`$ are displayed in figure 2. These results are clearly correlated with those for the $`GV`$ curve and can be summarized as follows: In the limit $`t_K1`$ the characteristics of shot-noise power spectrum for both $`s`$ and $`p`$ wave superconductors are consistent with those obtained for purely ballistic junctions which exhibit strong suppressing of the shot-noise power in the sub-gap region. At lower $`t_K`$ the physics is distinct. For $`t_K=5`$ the noise spectrum for $`s`$ wave superconducting lead still shows features typical for a junction with relatively high transparency, while the results for $`t_K=2`$ are more similar to those for a low transparency junction. Such dependence is explicitly exposed in the plot of the Fano factor versus the applied voltage (see Fig.3). Though the Fano factor does not reach the maximum value of 2, it is strongly enhanced for the smaller value of $`t_K=2`$. For $`pwave`$ superconductor the shot-noise power (as function of voltage) reflects the same physics as in the conductance: ZES makes the Kondo resonance less vulnerable to the impact of superconductivity and the junction remains close to the unitary limit nearly in the whole range of values of $`t_K`$ considered here (see the solid curves on Fig.2 and Fig.3).

Fig. 2 The shot-noise power $`K`$ (in units of $`e^2\mathrm{\Delta }/h`$) as a function of $`V`$ (in units of $`\mathrm{\Delta }/e`$) for an $`SKN`$ junction. The parameters and notations are the same as in Fig. 1.($`t_K`$ grows from top to down. The curves with $`t_K`$=100 coincide for both s and p wave superconductors.

Fig. 3 The Fano factor $`K/2eI`$ as a function of $`V`$ (in units of $`\mathrm{\Delta }/e`$) for an $`SKN`$ junction. The parameters and notations are the same as in Fig. 1.

## V Josephson current

In this section we study an equilibrium property (Josephson effect) of an SKS junction in which both electrodes are $`pwave`$ superconductors (for comparison we also represent the results for s-wave superconducting leads). In the Kondo regime, the self-consistent equations for the $`pwave`$ case(13) and (14) now read

$`\stackrel{~}{ϵ}+{\displaystyle \frac{2\mathrm{\Gamma }}{\pi }}\mathrm{log}{\displaystyle \frac{\stackrel{~}{ϵ}}{T_K}}=`$
$`{\displaystyle \underset{\omega }{}}[2\mathrm{\Gamma }{\displaystyle \frac{(1+\alpha (\omega ))\sqrt{\omega ^2+\mathrm{\Delta }^2}+\beta (\omega )\mathrm{\Delta }^2/\omega \mathrm{cos}^2\frac{\delta }{2}}{(1+\alpha (\omega ))^2\omega ^2+\stackrel{~}{ϵ}^2+\beta (\omega )^2\mathrm{\Delta }^2\mathrm{cos}^2\frac{\delta }{2}}}`$
$`{\displaystyle \frac{2\mathrm{\Gamma }|\omega |}{\omega ^2+\stackrel{~}{ϵ}^2}}]`$ (20)
$`\stackrel{~}{\mathrm{\Gamma }}`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \left[\frac{2\stackrel{~}{ϵ}}{(1+\alpha (\omega ))^2\omega ^2+\stackrel{~}{ϵ}^2+\beta (\omega )^2\mathrm{\Delta }^2\mathrm{cos}^2\frac{\delta }{2}}\right]}`$ (21)
Here $`\alpha (\omega )=\stackrel{~}{\mathrm{\Gamma }}\sqrt{\omega ^2+\mathrm{\Delta }^2}/\omega ^2`$, $`\beta (\omega )=\stackrel{~}{\mathrm{\Gamma }}/\omega `$ and $`\delta `$ is the phase difference between two superconductors. Self consistent equations for the case of $`swave`$ superconducting leads were derived in Ref. arovas2

The Josephson current is obtained by variation of free energy with respect to the phase of the order parameter,

$$I=\frac{2e}{\mathrm{}}\frac{F}{\delta }.$$
(22)
Explicitly we then find,

$`I={\displaystyle \frac{e}{\mathrm{}}}{\displaystyle \underset{\omega }{}}{\displaystyle \frac{(\beta (\omega )\mathrm{\Delta })^2\mathrm{sin}\delta }{\omega ^2(1+\alpha (\omega ))^2+\stackrel{~}{ϵ}^2+(\beta (\omega )\mathrm{\Delta })^2\mathrm{cos}^2\frac{\delta }{2}}}`$ (23)
The self-consistency equation and the expression for the Josephson current can easily be extended to the case of an anisotropic coupling :$`\mathrm{\Gamma }_{L/R}=\mathrm{\Gamma }(1\pm p)`$. For this, one should replace $`\mathrm{cos}^2\frac{\delta }{2}(\mathrm{cos}^2\frac{\delta }{2}+p^2\mathrm{sin}^2\frac{\delta }{2})`$ and $`I\mathrm{sin}\delta I(1p^2)\mathrm{sin}\delta `$.

For a small anisotropy parameter $`p=0.1`$ and for values of $`t_K`$ =$`100,5,3,2`$, we calculate the Josephson current as a function of the phase difference $`\delta `$ at temperatures close to $`T=0`$. The plot is depicted in Fig.4.

Fig. 4 Josephson current versus phase difference $`\delta `$ at $`T0`$. Dash and solid curves correspond to s and p-wave superconductors, correspondingly.The parameters are the same as in Fig.1. The value of $`t_K`$ decreases from top downward.
In the $`swave`$ case (dashed lines in figure 4) the current at $`t_K`$ =$`100`$ qualitatively corresponds to the unitary limit glazman , while its amplitude is decreased at smaller values of $`t_K`$. Nevertheless, in this region of parameters, the Kondo effect strongly competes with superconductivity and therefore there are no traces of a $`\pi `$-junction. The mid-gap states for $`pwave`$ superconductors act as if they increase the effective normal region of the junction and therefore, the weight of the contribution to the Josephson current which is proportional to $`sin\delta `$ is reduced. When $`t_K`$ becomes smaller, the Josephson current for $`pwave`$ junction tends to be similar to that in an SIS junction.

Fig. 5 Temperature dependence of the maximum Josephson current for s-wave (dashed curves) and p-wave(solid curves) cases when $`t_K=`$100, 5, 3 and 2 (from the top to the bottom). The phase $`\delta `$ of each curve is chosen so as to give the maximum value of Josephson current at $`T=0`$ in Fig. 4. The rest parameters are the same as in Fig. 1.
By the same reasoning the temperature dependence (see Fig. 5) of the Josephson current is more similar to the usual SIS junction for the $`pwave`$ case, while for $`swave`$ superconductors the current behavior is closer to the temperature dependence of the Josephson Current in SNS junctions. To be more precise, Fig. 5 displays the temperature dependence of the maximum Josephson current $`I(T)`$ for s-wave (dashed lines) and p-wave (solid lines) cases. For p-wave junctions, $`I(T)`$ assumes relatively large value near $`T=T_C`$ while for s-wave junctions, it rapidly deceases as $`T`$ increases. This difference stems from the question of whether the mid-gap state appears or not. The situation is similar to the one encountered in Josephson current through in an SIS system tanaka&kashiwaya . The difference becomes more prominent when the Kondo effect is suppressed; $`I(T)`$ drops more rapidly for s-wave junctions than for p-wave ones as $`t_K`$ decreases.

In conclusion, we have analyzed an important physical problem involving strong correlations, the Kondo effect and superconductivity. These aspects can be combined in an $`SKN`$ junction consisting of an Anderson impurity (in the Kondo regime). We have developed a theoretical framework by which it is possible to investigate an interplay between Andreev reflections and the formation of Kondo resonance in the Kondo regime $`T<\mathrm{\Delta }<T_K`$. In this limit we calculated Josephson Current. We have also investigated non-equilibrium aspects and elucidated the nonlinear $`GV`$ characteristics. Moreover, we calculated the shot-noise power spectrum of $`SKN`$ junctions at voltages $`eV\mathrm{\Delta }`$. It is found that at sufficiently large $`t_K`$ the Kondo resonance effectively turns the junction behavior to be similar to that of highly transparent non-interacting weak links for both $`s`$ and $`pwave`$ superconductors. However, when the ratio of the Kondo temperature to the superconducting gap becomes smaller the behavior of these two types of junctions is quite different: the Kondo resonance persists much more effectively for junctions with $`pwave`$ leads than it does for $`swave`$ leads.

Acknowledgment We would like to thank Y. Tanaka for very helpful and stimulating discussions. This research is supported by DIP German Israel Cooperation project Quantum electronics in low dimensions by the Israeli Science Foundation grant Many-Body effects in non-linear tunneling and by the US-Israel BSF grant Dynamical instabilities in quantum dots. One of us (TA) is supported by JSPS Research Fellowships for Young Scientists.


format: text

text:
# Dynamics of a Massive Piston in an Ideal Gas: Oscillatory Motion and Approach to Equilibrium

## 1 Introduction

The time evolution of a gas filled container divided by a massive piston is an old problem which has attracted much attention recently from both a conceptual and computational point of view, \[Li, G, KBM\]. While we do not believe that there is any “paradox” associated with this problem, the basic mechanism of energy transport across the piston from the hot side to the cold side which are at equal pressure but different temperatures, was already described for an idealized version, by one of us in 1959 \[L1\], many intriguing questions remain. Not surprisingly the microscopic motion of the piston, and thus the time evolution of the system from an initial, far from equilibrium state, to the final true equilibrium state is not easy to compute analytically or even numerically for large systems \[KBM\]. It is therefore interesting to consider the greatly simplified case when the gas particles only interact via collisions with the massive piston \[L1, GF, GP\]. This serves, among other things, as a model for the approach to equilibrium in such a “weakly” interacting macroscopic system.

Here we carry out studies on this system under the “simplest” initial conditions. The positions and velocities of gas particles of unit mass, in an insulated box of size $`L^3`$, are picked, at $`t=0`$, from a Poisson process with some density function, i.e. they are assumed to be independent random variables with a given double peaked velocity distribution $`p(v)=p(v)`$ (see (2.4)) and a constant spatial density throughout the box. The piston, which acts as a single particle with area $`L^2`$ and mass $`M_LL^2`$, is then released at the position $`X(0)=L/2`$ with zero velocity $`V(0)=0`$.

It might be thought that the random fluctuations of the piston’s velocity due to the initial randomness in the positions and velocities of the gas particles will be relatively small when $`L`$ is large (the number of gas particles growing like $`L^3`$) and vanish in an appropriate hydrodynamic scaling limit (as $`L\mathrm{}`$) of time, space and piston mass. In fact, some rigorous results in this direction were obtained in \[LPS\] and \[CLS\] for a class of initial conditions more general than those considered here. It is proven in those papers that for those initial distributions the dynamics of the piston, in the hydrodynamic scaling limit, is governed by deterministic equations for as long as no particle collides with the piston more than twice. After that period, however, we could not control the random fluctuations in the particle configuration anymore. It is clear, though, that for an initial state which is spatially uniform and has the same velocity distribution on both sides of the piston, the deterministic macroscopic evolution would predict that the position of the piston remains stationary forever. On the other hand, when the initial velocity distribution of the particles is not Maxwellian, the system is not at thermal equilibrium, and we expect it to somehow evolve toward equilibrium for almost any initial configuration (with respect to the Liouville measure on the energy surface). This is indeed what we found for every initial distribution. The path to the Maxwellin taken by the system turned out however to be quite sensitive to the initial distribution.

Here we present the results of detailed numerical simulations for one particular initial density $`p(v)`$ that vanishes unless $`0.5|v|1`$. We found that the quiescent state becomes unstable after a few (5–10) recollisions of each gas particle with the piston. The system then develops large (on a macroscopic scale) oscillatory motions which last for a long time. They look like smooth harmonic oscillations with an (almost) constant amplitude and piston speed comparable to that of the gas particles. Eventually, though, the oscillations dampen and vanish, and the system approaches a stable equilibrium state with constant gas density and Maxwellian velocity distribution111The latter is of course expected on general grounds for, as pointed out by Boltzmann: “\[the Maxwell distribution\] is characterized by the fact that by far the largest number of possible velocity distributions have the characteristic properties of the Maxwell distribution, and compared to these there are only a relatively small number of possible distributions that deviate significantly form Maxwell’s. Whereas Zermelo says that the number of states that finally lead to the Maxwellian state is small compared to all possible states, I assert on the contrary that by far the largest number of possible states are “Maxwellian” and that the number that deviate from the Maxwellian state is vanishingly small”, \[B\], see also \[L2\]..

Our main conclusions from the numerical investigations, which we also “derive” heuristically, is that 1) while the dynamics are non chaotic in the technical sense of there being no positive Lyapunov exponents, there are for typical initial conditions222Typical here means of probability converging to one (presumably, exponentially fast) as $`N\mathrm{}`$ with respect to the invariant Liouville measure. There are clearly some exceptional initial configurations (e.g., those with a complete symmetry about $`x=L/2`$ with opposite velocities) when $`V(t)=0`$ at all times and so the speed of gas particles never changes. We expect (for the reasons given by Boltzman \[B\]) all such initial states to be unstable with respect to generic small perturbations. enough interactions to bring the system to equilibrium along “interesting” nontrivial pathways, and 2) that the time of onset of the instability grows with $`L`$ as $`L\mathrm{log}L`$. This means that even on the hydrodynamical scale $`\tau =t/L`$, the deterministic behavior in which nothing happens on the spatial scale, $`y=x/L`$, would remain valid when $`L\mathrm{}`$. In terms of $`\tau `$ however this onset time grows only like $`\mathrm{log}L`$ so that in “practice” this behavior extends to macroscopic systems. We believe that this instability is related to the fact that the deterministic solutions are unstable for the kind of $`p(v)`$ we consider here. For other initial $`p(v)`$, for which the deterministic solution is stable, we do not get such large oscillations but, as expected, we still get an approach to Maxwellian, see Section 7.

## 2 Description of the model

Consider a cubical container $`\mathrm{\Lambda }_L=[0,L]\times [0,L]\times [0,L]`$ filled with an ideal gas consisting of $`N`$ particles. The container is divided into two parts by a wall (piston) orthogonal to the $`x`$ axis. At time $`t=0`$ the wall is released and then it can move freely without friction inside $`\mathrm{\Lambda }_L`$ along the $`x`$-axis, under the action of elastic collisions with the gas particles, each of which has the same fixed mass $`m`$. Since the piston’s area is $`L^2`$, we assume its mass $`M_L`$ to be proportional to $`L^2`$ and given by $`M_L=bmL^2`$ with a fixed constant $`b>0`$ (we set $`m=1`$ and $`b=2`$ in our numerical simulations).

Since the components of the particle velocities perpendicular to the $`x`$-axis play no role in the dynamics of the piston, we may assume that each particle has only one component of velocity, $`v`$, directed along the $`x`$-axis. Hence, each gas particle can be specified by a pair $`(x_i,v_i)`$, where $`i=1,\mathrm{},N`$. We shall take the total number of particles $`N`$ and the total kinetic energy of the system proportional to $`L^3`$; and we are interested in the behavior of the system for $`L1`$.

The initial configuration of gas particles and their velocities $`\{(x_i,v_i)\}_{i=1}^N`$ is selected at random with statistics given by a (two-dimensional) Poisson process on the $`x,v`$ plane with a given density, $`p(x,v)=L^2p(v)`$. More precisely, for any domain $`D[0,L]\times \text{I}\text{R}^1`$ the number of particles in $`D`$, i.e. $`N_D=\mathrm{\#}\{i:(x_i,v_i)D\}`$, has a Poisson distribution with parameter

$$\lambda _D=L^2_Dp(v)𝑑x𝑑v$$

That is, for each $`k0`$,

$$P(N_D=k)=\frac{\lambda _D^k}{k!}e^{\lambda _D}$$

The numbers of particles in nonoverlapping domains are independent. The function $`p(v)=p(v)`$ takes values of order one, and the factor of $`L^2`$ is simply the cross-sectional area of the container $`\mathrm{\Lambda }_L`$. The piston is initially at rest at the midpoint, $`X(0)=L/2`$, $`V(0)=0`$ and the position and velocity at time $`t`$ are denoted by $`X`$, $`0XL`$, and $`V`$.

We note that the total number of gas particles $`N`$ in the container is random. So are the numbers of particles to the left and to the right of the piston, call them $`N_{}`$ and $`N_+`$, respectively ($`N=N_{}+N_+`$). The initial total kinetic energy of the system, $`\frac{1}{2}mv_i^2`$, is also a random variable. Once chosen, the values of $`N_{},N_+`$ and the total kinetic energy $`E`$ of gas plus piston are (presumably, the only) integrals of motion.

The gas particles and the piston move freely (with constant velocity) between elastic collisions of particles with the piston and the walls. When a particle collides with a wall at $`x=0`$ or $`x=L`$, its velocity simply reverses. If a particle with velocity $`v`$ hits the piston whose velocity is $`V`$, then their velocities after the collision, call them $`v^{}`$ and $`V^{}`$, respectively, are given by

$$V^{}=(1\epsilon )V+\epsilon v$$
(2.1)
$$v^{}=(1\epsilon )v+(2\epsilon )V$$
(2.2)
where

$$\epsilon =\frac{2m}{M+m}=\left(\frac{2}{bL^2+1}\right)$$
(2.3)
In order to avoid recollisions of gas particles with the piston, for at least some initial period of time, we impose a velocity cutoff333As explained in Section 7, the choice of $`p(v)`$ plays an important role in the later time evolution of the system.

$$p(v)0\mathrm{if}|v|v_{\mathrm{max}}\mathrm{or}|v|v_{\mathrm{min}}$$
(2.4)
with some $`0<v_{\mathrm{min}}<v_{\mathrm{max}}<\mathrm{}`$, cf. \[CLS, LPS\]. Under these conditions, there will be an initial time interval of length $`O(L)`$ during which each particle colliding with the piston will have to travel to the wall, bounce off it, and then travel back to the piston before it can hit it again. Therefore, during that interval, there will be no recollisions of any gas particle with the piston. We call such an interval of time the zero-recollision interval. Likewise, if the piston remains slow enough after recollisions occur, there will be another interval of time of order $`L`$ during which each gas particle experiences at most one recollision with the piston. We call it the one-recollision interval.

The hydrodynamic limit is now obtained by rescaling space, $`y=x/L`$, $`0<y<1`$, and time, $`\tau =t/L`$. In the new space-time coordinates $`y,\tau `$ the zero-recollision interval and the one-recollision interval will be of order one, and we will denote them by $`(0,\tau _0)`$ and $`(\tau _0,\tau _1)`$ respectively. As already mentioned, we prove in \[CLS\], for a more general class of initial densities $`p(x,v)`$ that during the interval $`(0,\tau _1)`$ the functions $`Y_L(\tau )`$ and $`W_L(\tau )`$

$$Y_L(\tau )=X(\tau L)/L,W_L(\tau )=dY_L/d\tau =V(\tau L)$$

converge uniformly in probability, as $`L\mathrm{}`$, to some deterministic functions $`\overline{Y}(\tau )`$ and $`\overline{W}(\tau )`$; $`\overline{W}`$ satisfies a certain algebraic equation and $`d\overline{Y}/d\tau =\overline{W}`$. In the case considered here, when the density function $`p(x,v)`$ is independent of $`x`$ and symmetric in $`v`$, i.e. $`p(v)=p(v)`$ for all $`x,v`$, then the hydrodynamic evolution is trivial: $`\overline{Y}(\tau )0.5`$ and $`\overline{W}(\tau )0`$ for all $`\tau >0`$. Hence, $`Y_L(\tau )0.5`$ and $`W_L(\tau )0`$ as $`L\mathrm{}`$, uniformly for all $`\tau (0,\tau _1)`$.

In this paper we investigate numerically what happens to $`Y_L(\tau )`$ and $`W_L(\tau )`$ beyond the interval $`(0,\tau _1)`$. More precisely, for how long does the stochastic trajectory $`\{Y_L(\tau ),W_L(\tau )\}`$ remain close to the deterministic one $`\{\overline{W}(\tau ),\overline{Y}(\tau )\}`$? And what does it look like in the long run?

## 3 Numerical results

In the computer simulations described here and in Sections 4-6, we set

$$p(x,v)=p(v)=\{\begin{array}{cc}1\hfill & \mathrm{if}0.5|v|1\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}$$
(3.1)
so $`v_{\mathrm{min}}=0.5`$ and $`v_{\mathrm{max}}=1`$. The $`x`$ and $`v`$ coordinates of all the particles are then independent random variables uniformly distributed in their ranges $`0<x<L`$ and $`v_{\mathrm{min}}|v|v_{\mathrm{max}}`$. Our computer program first selects $`N`$ according to the Poisson law with mean $`L^3`$, and then generates all $`(x_i,v_i)`$, $`1iN`$, independently according to their uniform distributions. We used the random number generator described in \[MN\]. The parameter $`L`$ changed in our simulations from $`L=30`$ to $`L=300`$. For $`L=300`$ the system contains $`L^3=27,000,000`$ particles.

Once the initial data is generated randomly, the program computes the dynamics by using the elastic collision rules (2.1), (2.2). All calculations were performed in double precision, with coordinates and velocities of all particles stored and computed individually. We note that memory requirements alone can be enormous – the program needs over 430Mb RAM in order to run the model with $`L=300`$. Also, the task of determining which particle is to collide with the piston next is nontrivial. To avoid a long search through the entire set of $`N`$ particles after each collision, we assembled a smaller group of particles located in a vicinity of the piston, where the search is performed, and updated this group periodically, as time goes on. The calculations were done in C++ and HPF (High Performance Fortran) on a Dell Power Edge machine with Dual 733MHz processors at the University of Alabama in Birmingham. Our code is available on the web page referred to in the next section.

Figure 1 presents a typical trajectory of the piston. Here $`L=100`$. The position and time are measured in hydrodynamic variables $`Y=X/L`$, $`0<Y<1`$, and $`\tau =t/L`$.

Initially, the piston barely moves about its equilibrium position: recall that the hydrodynamic trajectory of the piston is $`\overline{Y}(\tau )0.5`$ for all $`\tau >0`$ and that this holds exactly for $`\tau <2`$ as $`L\mathrm{}`$, \[CLS\]. Then, at times $`\tau `$ between $`3`$ and $`5`$, the random vibrations of the piston grow and become quite visible on the $`y`$-scale, but for a short while they look random. After that the piston starts travelling back and forth along the $`y`$ axis, making excursions farther and farther away from the equilibrium point $`y=0.5`$. Very soon, at $`\tau =\tau _{\mathrm{max}}8`$, the swinging motion of the piston reaches its maximum, $`(\mathrm{\Delta }Y)_{\mathrm{max}}=\mathrm{max}|Y(\tau )0.5|0.1`$. Then the oscillations of the piston dampen in size and seem to stabilize at an amplitude $`A0.04`$. At the same time the trajectory of the piston smoothes out and enters an oscillatory mode with a period $`\tau _{\mathrm{per}}1.63`$.

The velocity of the piston $`W(\tau )`$ follows a similar pattern. Its random fluctuations grow after $`\tau 2.5`$, then it starts swinging up and down, reaches the maximum value of $`W_{\mathrm{max}}=\mathrm{max}|W(\tau )|0.4`$ at time $`\tau 9.5`$. After that the oscillations of $`W(\tau )`$ dampen and seem to stabilize.

Note that the graph of the function $`Y(\tau )`$ looks much smoother than that of $`W(\tau )`$, as would be expected from the fact that $`Y(\tau )`$ is the integral of $`W(\tau )`$. Interestingly, both functions $`Y`$ and $`W`$ smooth out as time goes on.

This cycle of the gas motion in the container between the walls and the piston continues for a long time with the same period $`\tau _{\mathrm{per}}1.63`$, independent of $`L`$, but the amplitudes of both $`Y(\tau )`$ and $`W(\tau )`$ are slowly decreasing, see Fig. 3.

The oscillations of the piston with decaying amplitude can be described, in the interval $`30<\tau <1000`$, approximately by

$$Y_1(\tau )Ae^{\lambda (\tau 20)}\mathrm{sin}\omega (\tau \alpha )$$
(3.2)
with $`A=0.046`$ and some constant $`\lambda >0`$. Correspondingly, $`W_1(\tau )=dY_1/d\tau `$ in the same interval $`30<\tau <1000`$ is

$`W_1(\tau )`$ $``$ $`\lambda Y_1+Ae^{\lambda (\tau 20)}\omega \mathrm{cos}\omega (\tau \alpha )`$ (3.3)
$`=`$ $`Ae^{\lambda (\tau 20)}[\lambda \mathrm{sin}\omega (\tau \alpha )+\omega \mathrm{cos}\omega (\tau \alpha )]`$
$`=`$ $`A_1e^{\lambda (\tau 20)}\mathrm{sin}\omega (\tau \beta )`$
with $`A_1=A\sqrt{\omega ^2+\lambda ^2}`$ and some $`\beta `$ related to $`\alpha `$.

To check how well our formula (3.2) agrees with the data, we computed the amplitude $`A(\tau )`$ as a function of time $`\tau `$, by fitting a sine function $`Y_0(\tau )=A\mathrm{sin}\omega (\tau \alpha )`$ “locally”, on the interval $`(\tau 5,\tau +5)`$ for each $`\tau `$. Fig. 4 shows $`A(\tau )`$ on the logarithmic scale, which looks almost linear on the interval $`30<\tau <800`$. (After that, $`Y(\tau )`$ becomes quite unstable, with already small amplitude $`A(\tau )`$ decreasing abruptly, possibly due to the interference from random fluctuations, so we left that part out.)

We used the least squares fit to estimate $`\lambda =0.00264`$ for the run shown on Figs. 1-4. Since $`\lambda `$ is small, the oscillations indeed die out very slowly. The “half-life” time (the time it takes to reduce the amplitude by a factor of two) is $`\tau _{1/2}=\lambda ^1\mathrm{ln}2263`$. Note that over the time interval $`30<\tau <800`$ we only observed the reduction of the amplitude by a factor of about 10, and the periodic oscillations were still visible on the plot at $`\tau >900`$. Also, $`\lambda `$ and hence $`\tau _{1/2}`$ depend on the system size $`L`$, see below.

## 4 Dependence of the Time Evolution on System Size

Many of the characteristics of the piston trajectory described above ($`(\mathrm{\Delta }Y)_{\mathrm{max}}`$, $`W_{\mathrm{max}}`$, $`A`$, $`\tau _{\mathrm{per}}`$) appear to be independent of $`L`$. The following table444The numerical data here are given primarily for demonstrating the typical scale of oscillations. They are not meant to be estimates of physical parameters, so we do not provide error bars. For $`L200`$, where more than one experimental trajectory was generated, we have estimated that the accuracy of our numerical values in Table 1 is at least within 10%. presents computed values of all these parameters for different $`L`$’s.

| $`L`$ | $`(\mathrm{\Delta }Y)_{\mathrm{max}}`$ | $`W_{\mathrm{max}}`$ | $`A`$ | $`\tau _{\mathrm{per}}`$ |
| --- | --- | --- | --- | --- |
| 30 | 0.114 | 0.40 | 0.042 | 1.70 |
| 50 | 0.121 | 0.39 | 0.045 | 1.65 |
| 80 | 0.105 | 0.37 | 0.042 | 1.65 |
| 100 | 0.100 | 0.34 | 0.041 | 1.64 |
| 120 | 0.122 | 0.39 | 0.041 | 1.62 |
| 150 | 0.111 | 0.37 | 0.045 | 1.62 |
| 200 | 0.102 | 0.37 | 0.037 | 1.62 |
| 250 | 0.100 | 0.41 | 0.042 | 1.64 |
| 300 | 0.122 | 0.42 | 0.045 | 1.65 |

Table 1. Principal characteristics of the piston dynamics.

In each case, we averaged over several experimentally generated trajectories (for $`L=250`$ and 300, just one trajectory was used).

There are however other quantities such as $`\tau _{\mathrm{max}}`$, $`\tau _{1/2}`$, and the related $`\lambda `$, which depend in a systematic way on $`L`$. In particular, we estimated numerically that $`\tau _{1/2}L^{1.3}`$, hence $`\lambda L^{1.3}`$, see Fig. 5. We will denote $`\lambda =\lambda _L`$ and discuss it further in Section 6. We also noticed that in some runs with larger $`L`$’s (such as $`L=150`$ and $`L=200`$) the exponent $`\lambda `$ changes with time, it is higher when $`\tau <100`$ and lower when $`\tau >100`$.

But most importantly, the time of the largest oscillations $`\tau _{\mathrm{max}}`$ and the related time of the onset of the instability $`\tau _c`$, see below, seem to slowly grow with $`L`$, very likely as $`\mathrm{log}L`$. To understand this fact, we looked into the mechanism of the build-up of random fluctuations of the piston position and velocity, which eventually result in their large nearly harmonic oscillations. To this end we plotted the histogram of the (empirical) density of gas particles in the $`y,v`$ plane at various times $`0<\tau <30`$, see samples in Fig. 6. The initial density (at time zero) is almost uniform over the domain $`0<x<L`$ and $`v_{\mathrm{min}}|v|v_{\mathrm{max}}`$ (variations in the initial configuration always exist, because it is generated randomly). Then, for $`0<\tau <1`$, the piston experiences random collisions with particles and acquires a speed of order $`M_L^{1/2}=O(1/L)`$, see \[L1, Ho, DGL\]. These small fluctuations of the piston velocity result in changes of the velocities of the particles which leave the piston after collision. Thus the outgoing particles on the right hand side of the piston have velocities in the interval $`(v_{\mathrm{min}}+2W(\tau ),v_{\mathrm{max}}+2W(\tau ))`$ while those on the left hand side of the piston have velocities in the interval $`(v_{\mathrm{min}}+2W(\tau ),v_{\mathrm{max}}+2W(\tau ))`$. Hence, the region in the $`y,v`$ plane where the density of the particles is positive is no longer a rectangle with straight sides, now its boundaries are curves whose shape nearly repeats the graph of the randomly evolving piston velocity $`W(\tau )`$. While the variations of $`O(1/L)`$ of these boundary curves may seem small, it is crucial that on opposite sides of the piston they go in opposite directions. Indeed, when $`W(\tau )>0`$, then the outgoing particles on the right hand side accelerate and those on the left hand side slow down. When $`W(\tau )<0`$ the opposite happens.

Next, the particles that have collided with the piston travel to the wall and come back to the piston. Now their densities are less regular than they were initially – the regions in the $`x,v`$ plane where the density is positive, are curvilinear domains. When they hit the piston, they shake it back and forth more forcefully than before, because the velocities of the incoming particles on the opposite sides of the piston are now negatively correlated. When particles on the right hand side are fast, those on the left hand side are slow, and vice versa. This produces a resonance-type effect destabilizing the piston dramatically and the velocity of the piston $`W(\tau )`$ experiences larger fluctuations than before. The velocities of the newly outgoing particles will again go up and down in opposite direction, on a greater scale than before.

As time goes on, the above phenomenon repeats over and over, with larger and larger fluctuations of the gas and piston velocities, until the distribution of gas particles completely breaks down. At times $`\tau 10`$, two large clusters of particles are formed, one on each side of the piston. When one cluster bombards the piston, the other moves away from it and hits the wall, then they exchange their roles. The clusters have sizes of about 0.3–0.5 in the $`y`$ direction and the particle velocities range from about 0.2 to just over 1. The average velocity is about 0.5–0.6 and so the clusters hammer the piston periodically with period 1.6–2.0, which is close to the experimentally determined period of piston oscillations, see above.

Fig. 6 shows six snapshots of the empirical density of gas particles taken at different times. At $`\tau =0`$ the gas fills (almost uniformly) two rectangles $`\{(y,v):0.5<|v|<1,0<y<1\}`$. At $`\tau =2.3`$ one can see some ripples on the boundaries of these rectangles. At time $`\tau =4.2`$ the irregularities grow and at $`\tau =5.9`$ the rectangular shape is broken down. Two large clusters of particles are formed, both appear in the upper half-plane $`v>0`$, i.e. at that time both clusters move to the right (one toward the piston, the other away from it). Later the density undergoes strange formations ($`\tau =7.4`$) but eventually smoothes out and enters a slow process of convergence to Maxwellian ($`\tau =18.6`$) described below. Note a narrow white line around $`v=0`$, meaning the total lack of very slow particles at time $`\tau =18.6`$.

A longer sequence of snapshots at times $`0<\tau <30`$ is posted on the web page

$$\mathrm{www}.\mathrm{math}.\mathrm{uab}.\mathrm{edu}/\mathrm{chernov}/\mathrm{piston}/\mathrm{pictures}/\mathrm{piston}.\mathrm{html}$$

It gives a spectacular view of the entire system evolution.

The above analysis may suggest that the fluctuations of the piston velocity roughly increase by a constant factor during each time interval of length one. Indeed, initial random fluctuations $`W_aO(1/L)`$ result in additional changes of velocities of outgoing particles by $`2W_a`$. When those particles come back to the piston (in time $`\mathrm{\Delta }\tau 1`$), they kick its velocity to the level of $`2W_a`$. Then the newly outgoing particles acquire an additional velocity $`4W_a`$, etc. Over each time interval of length one the fluctuations double in size. This is an obvious oversimplification of the real dynamics, but it leads to a reasonable conjecture

$$W_a(\tau )\frac{CR^\tau }{L}$$
(4.4)
where $`W_a(\tau )`$ are typical fluctuations of the piston velocity at time $`\tau `$ and $`C,R>0`$ are constants.

We tested the above formula numerically as follows. Let $`W_c>0`$ be some preset critical value of the piston speed and $`\tau _c=inf\{\tau >0:|W(\tau )|W_c\}`$ the (random) time when $`W_c`$ is first reached. This time plays the role of the “onset” of large fluctuations of the piston velocity. One would expect, based on (4.4) that

$$\tau _c\mathrm{ln}(W_cL/C)/\mathrm{ln}R$$
(4.5)
i.e. $`\tau _c`$ grows as $`\mathrm{ln}L`$ when $`L`$ increases.

We found $`\tau _c`$ experimentally for $`W_c=0.1`$ and $`W_c=0.15`$ and checked that (4.5) agreed well with the data, see Fig. 7. By the least squares fit we estimated $`C=0.45`$ and $`R=1.6`$.

## 5 Approach to equilibrium

We examined the convergence of the velocity distribution of gas particles to a Maxwellian. At any given time $`\tau >0`$, let

$$F_\tau (u)=\mathrm{\#}\{i:v_i<u\}/N$$

be the empirical (cumulative) distribution function of particle velocities. At equilibrium, it should be close to the normal distribution function $`\mathrm{\Phi }(x)`$. As a measure of their closeness, we used the supremum of the difference

$$D_\tau =\underset{\mathrm{}<u<\mathrm{}}{sup}|F_\tau (u)\mathrm{\Phi }(u)|$$

Initially, $`D_00.245`$ for our choice of $`p(v)`$. One can expect that $`D_\tau `$ converges to 0 as $`\tau `$ grows when $`N`$ is large. In fact, if the velocities $`v_i`$ were independent and had Maxwellian distribution (which it would be in true statistical equilibrium), then $`D_\tau `$ would be of order $`O(1/\sqrt{N})`$, and the product $`D_\tau \sqrt{N}`$ would have a certain limit distribution, see the theory of the Kolmogorov-Smirnov statistical test \[Lu\]. In particular, it is known that for a Maxwellian the probability $`P(D_\tau \sqrt{N}>1)0.2`$. Based on this, we define the time of convergence to equilibrium by

$$\tau _{\mathrm{eq}}=inf\{\tau >0:D_\tau \sqrt{N}<1\}$$
(5.1)
Here the constant 1 as a critical value is chosen arbitrarily. We estimated $`\tau _{\mathrm{eq}}`$ for various $`L`$’s and found that $`\tau _{\mathrm{eq}}aL^b`$ with some constants $`a,b>0`$. By a least squares fit we found $`a=0.18`$ and $`b=2.47`$, see Fig. 8. (Note that the accuracy of our experimental data seems to increase with $`L`$, as the points are getting closer to each other and to the fitting line for larger $`L`$’s on Fig. 8.)

The plot of the product $`S=D_\tau \sqrt{N}`$ versus $`\tau `$ is given on Fig. 9 (for a particular run with $`L=40`$). It shows that, after an initial sharp drop for $`0<\tau <20`$, the statistic $`S`$ decreases exponentially in $`\tau `$. Another commonly used statistic to measure closeness to Maxwellian is

$$S^{}=3\frac{M_4}{M_2^2}$$

where $`M_2`$ and $`M_4`$ are the second and the fourth sample moments of the empirical velocity distribution, respectively. Fig. 9 shows that $`S^{}`$ converges to zero in a similar manner (for the same run with $`L=40`$).

Theoretical Considerations

As we noted in Section 2, the total energy $`E`$ and the numbers of particles in the left and right compartments ($`N_{}`$ and $`N_+`$) are integrals of motion. With these quantities fixed, the model can be reduced to a billiard system in a high-dimensional polyhedron by standard techniques, as we show next.

Let $`\{x_i\}`$, $`i=1,\mathrm{},N_+`$, denote the $`x`$-coordinates of the particles to the right of the piston, and $`\{x_i\}`$, $`i=1,\mathrm{},N_{}`$, those to the left of it (ordered arbitrarily). Put $`x_0=X\sqrt{M}`$, where $`X`$ is the coordinate of the piston and $`M`$ is its mass. Then the configuration space of the system (in the coordinates $`x_i`$, $`N_{}iN_+`$) is a polyhedron $`Q\text{I}\text{R}^{N+1}`$ (recall that $`N=N_{}+N_+`$) defined by inequalities

$$0x_N_{},\mathrm{},x_1x_0/\sqrt{M}x_1,\mathrm{},x_{N_+}L$$

It is known that the dynamics of our mechanical system (“gas+piston”) corresponds to the billiard dynamics in $`Q`$, see \[CFS\]. That is, the configuration point $`𝐪Q`$ moves freely and experiences specular reflections at the boundary $`Q`$. The velocity vector

$$𝐩=\dot{𝐪}=\{v_N_{},\mathrm{},v_1,V\sqrt{M},v_1,\mathrm{},v_{N_+}\}$$

has constant length, since $`𝐩^2=2E=`$const. Therefore, the phase space of the billiard system is $`=Q\times S_\rho ^N`$ where $`S_\rho ^N`$ is the $`N`$-dimensional sphere of radius $`\rho =\sqrt{2E}`$.

The billiard system has a natural equilibrium state given by the Liouville measure $`\mu `$ on $``$, which is the product of a uniform measure on the polyhedron $`Q`$ and a uniform (Lebesgue) measure on the sphere $`S_\rho ^N`$, i.e. $`d\mu =dqdp`$. The properties of billiard dynamics depend heavily on the curvature of the boundary $`Q`$. In our case $`Q`$ is a polyhedron, hence its boundary consists of flat sides with zero curvature. A prototype of such systems is billiard in a polygon. It is well known that (see, e.g., \[C\])

Fact. For billiards in polygons and polyhedra (and hence, for our mechanical model of a piston in the ideal gas) all Lyapunov exponents vanish, and so does the Kolmogorov-Sinai entropy.

Systems with zero Lyapunov exponents and zero entropy are not regarded as chaotic, but they still may be ergodic. In fact, billiards in generic polygons are ergodic \[KMS\]. Moreover, for many nonergodic polygons, the phase space is foliated by invariant subsurfaces on which the dynamics is ergodic.

Even though there are no similar results, to our knowledge, for billiards in high-dimensional polyhedra, one can expect that they, too, have similar properties. That is, they are generically ergodic or become ergodic after trivial reductions. In our case, the billiard in $`Q`$ is, perhaps, ergodic for typical values of $`M`$, or else the phase space is foliated by invariant submanifolds on which the dynamics is ergodic, and that those submanifolds fill $``$ pretty densely. In the latter case, one would hardly distinguish experimentally between such a nonergodic system and a truly ergodic one.

Hence, we can assume that our system is ergodic or very close to ergodic in the above sense. Then almost every trajectory eventually behaves according to the invariant measure $`\mu `$, independently of the initial state. In particular, for any initial gas density and velocity distribution (given by the function $`p(x,v)`$, see Section 2) the hydrodynamic regime for a finite $`L`$ is only valid on a finite interval of time – eventually the system will relax to equilibrium. We expect in fact that in terms of the “macroscopic” variables, say, the one particle distribution function, the system will relax to an effective equilibrium, as defined by (5.1) in terms of $`\tau _{\mathrm{eq}}`$, which is much smaller than the exponentially long time (in $`L`$) required for the ergodic theorem. So the real question is how does this time depend on $`L`$. According to our earlier discussion $`\tau _c\mathrm{log}L`$ and $`\tau _{\mathrm{eq}}L^{5/2}`$.

At equilibrium, the distribution of coordinates $`x_i`$ and velocities $`v_i`$ are determined by the Liouville measure $`\mu `$, which is uniform in the phase space. Physically interesting (and only observable) are its marginal measures, i.e. projections, on lower-dimensional subspaces. The marginal measures of the velocities are normal (Gaussian) for large $`N`$ (as $`N\mathrm{}`$).

In particular, each individual velocity $`v_i`$ converges in law to a Maxwellian (i.e., normal) distribution with zero mean and variance $`2E/N=`$const. The same holds for the “piston” component of the velocity, $`\dot{x}_0=V\sqrt{M}`$, hence the piston velocity $`V`$ will be normally distributed with zero mean and standard deviation const$`/\sqrt{M}=`$const$`/L`$, as $`L\mathrm{}`$. In our case $`V`$ has standard deviation $`\sqrt{7/24}/L0.5/L`$. This conclusion agrees well with our numerical data555It also allows us to estimate, in a peculiar way, the time of convergence to equilibrium, $`\tau _{\mathrm{eq}}`$, discussed earlier. Assume all the initial velocities are of order one, and note that the largest Maxwellian velocities are of order $`\sqrt{N}=L^{3/2}`$. Since each collision with the piston adds $`O(1/L)`$ to a particle’s velocity, it takes $`L^{5/2}`$ collisions to reach the maximum. An amazing agreement with the estimate $`\tau _{\mathrm{eq}}L^{2.47}`$ reported near Fig. 8..

The equilibrium distribution of the piston coordinate $`X`$ is also determined by the projection of the uniform measure $`dq`$ on $`Q`$ onto the $`x_0`$ axis. Before we do that, let us get rid of $`M`$ in the definition of both $`Q`$ and $`x_0`$. A simple change of variable $`X=x_0/\sqrt{M}`$ allows us to redefine $`Q`$ by

$$0x_N_{},\mathrm{},x_1Xx_1,\mathrm{},x_{N_+}L$$

Furthermore, rescaling $`Y=X/L`$ and $`y_i=x_i/L`$ gives a new, simpler, definition of $`Q`$:

$$0y_N_{},\mathrm{},y_1Yy_1,\mathrm{},y_{N_+}1$$

“Integrating away” the variables $`y_i`$ yields the following equilibrium density for $`Y`$:

$$f(Y)=cY^N_{}(1Y)^{N_+}$$

for $`0<Y<1`$, where $`c`$ is the normalizing factor that can be computed explicitly. Put $`z=(Y0.5)\sqrt{8K}`$, then the density of $`z`$ is given asymptotically by $`f(z)\frac{1}{\sqrt{2\pi }}e^{z^2/2}`$. Hence, $`Y`$ is asymptotically gaussian with mean $`0.5`$ and variance $`(4N)^1=(4L^3)^1`$. Therefore, in equilibrium

$$|Y0.5|\frac{1}{2L\sqrt{L}}\frac{1}{2\sqrt{N}}$$

and the probability of observing larger fluctuations is exponentially small in $`N`$. Note that fact is independent of the piston mass, as it has to be for an equilibrium (classical) system. This is consistent with our observation of the piston coordinate reaching a positive constant value, $`|Y0.5|0.1`$, during a short time interval $`0<\tau <20`$, since our initial conditions were selected according to the density function $`p(x,v)`$ which has probability exponentially small in $`N`$.

## 6 Remarks

1. We checked the accuracy of our computer program in various ways. A simple one was based on the effect of round-offs on the energy of the system: the total energy was found to be practically constant over the whole period $`0<\tau <1000`$. A more sensitive test of the accuracy of the computations consists in using the time-reversal symmetry of the dynamics. Suppose at some time $`\overline{\tau }>0`$ the velocities of all gas particles and the piston are reversed ($`v_iv_i`$). Then the system is supposed to trace back its past trajectory and arrive to its initial state with all velocities reversed at time $`2\overline{\tau }`$. We verified this property numerically for various $`\overline{\tau }<50`$ and always found that the system did repeat its past trajectory in the sense that the graphs of $`Y(\tau )`$ and $`W(\tau )`$ over the interval $`(\overline{\tau },2\overline{\tau })`$ looked like perfect mirror images of the corresponding graphs over the interval $`(0,\overline{\tau })`$. Therefore, one can assume that round off errors remain negligibly small during times $`\tau <50`$. However, this accuracy test failed for larger times, $`\overline{\tau }>100`$, indicating that the system then “loses memory” of its initial state due to round-offs. This does not seem to affect the overall picture, though. As yet another test, we carried out some computations with single rather than double precision, and found that the changes were little. Still, some estimates requiring long runs, such as the equilibrium time $`\tau _{\mathrm{eq}}`$ on Fig. 8, might not be very accurate.

2. We tested the dependence of the piston oscillations on the $`b`$ factor involved in the formula $`M_L=bmL^2`$. Remember that we set $`b=2`$ in our main experiments. When we changed it to $`b=20`$ (this made the piston 10 times heavier), then the oscillations started slightly later and their amplitude was slightly larger, but otherwise the picture was very much the same. When we changed $`b`$ to 0.2 (and this made the piston 10 times lighter), then the oscillations started at about the same time as for $`b=2`$, but they dampened somewhat faster.

We also tried to change the piston mass $`M_L`$ even more drastically. When we set it to $`L^3`$ (insteas of $`L^2`$), it appeared that oscillations started only after a very long initial quiescent period. But we did not examine this fact in detail, see \[GPL\]. On the contrary, when we set the piston mass $`M_L`$ to $`L`$ (instead of $`L^2`$), large oscillations started very soon, but very quickly dampened and disappeared.

3. We note that Eq. (3.3) describes the time evolution of a damped harmonic oscillator. Accepting (3.3) over some time range, say, $`\tau [30,800]`$ for $`L=100`$, we can then look at the “inverse problem” of finding the effective spring constant and damping coefficient.

Using the original variables, $`t`$ and $`X_L(t)`$, we write

$$M_L\frac{d^2X_L}{dt^2}+K_L(X_LL/2)+\eta _L\frac{dX_L}{dt}=0$$

which has the solution $`X_LL/2e^{\alpha t}`$ with

$$\alpha =\frac{\eta _L}{2M_L}\pm i\sqrt{\frac{K_L}{M_L}}\left[1\frac{\eta _L^2}{4K_LM_L}\right]^{1/2}$$

This yields, to the lowest order in $`\eta _L^2/4K_LM_L`$, remembering that $`M_L=2L^2`$ and that $`\omega =2\pi /\tau _{\mathrm{per}}`$, with $`\tau _{\mathrm{per}}=t_{\mathrm{per}}/L1.63`$, $`K_L8\pi ^2/\tau _{\mathrm{per}}^2`$, i.e. the effective “restoring force” $`K_L`$ is independent of $`L`$. On the other hand, the damping coefficient is $`\eta _L=4L\lambda _L`$, were $`\lambda _L`$ was found experimentally to decrease as $`\lambda _LL^\gamma `$ with $`\gamma =1.3`$, see Section 4. Hence $`\eta _LL^{1\gamma }=O(L^{0.3})`$, i.e. the damping gets weaker as $`L\mathrm{}`$.

## 7 Discussion

We have presented here numerical results concerning the time evolution of a system with many degrees of freedom (up to $`27\times 10^6`$) one of which, the position of the piston, plays a very special role. Starting with the particular initial particle distribution given by (3.1) we found two striking features of the evolution: (1) the velocity distribution of the particles approaches a Maxwellian, i.e. the system goes toward thermal equilibrium and (2) the time evolution of the piston followed closely, after some initial period, that of a damped harmonic oscillator over an extended time interval, with initial oscillations as large as 1/10 of the system size.

As already noted, we expect from general considerations \[B\] that (1) should be true for any initial density $`p(x,v)`$ provided the number of particles $`N`$ ($`L^3`$) is large enough and one waits long enough. But what about (2)? Clearly, if we choose for $`p(|v|)`$ a Maxwellian, we expect only thermal fluctuations of order $`O(L^{3/2})`$ in the position of the piston. This is indeed what we found numerically. For other, non-maxwellian, initial densities the question turns out to be far from trivial. We are currently working with more general initial densities and will report results in a separate paper \[CCLP\]. Below we outline our program and mention some preliminary findings.

For simplicity, we assume that $`p(x,v)=p(|v|)`$ and $`X(0)=L/2`$, $`V(0)=0`$. In this case there is no apriori bias for the piston to move at all, and as already noted, the hydrodynamical (deterministic) equations, see \[CLS\], predict that, in the limit $`L\mathrm{}`$, the system would remain frozen in the initial state: $`Y(\tau )Y(0)=0.5`$ and $`W(\tau )W(0)=0`$ for all $`\tau >0`$ (in the variables $`y`$ and $`\tau `$, see Section 2). The density $`p(y,v,t)`$ will also remain constant in time. But what about the particle system with a large but finite $`L`$? How will the piston and the gas behave while the particle velocities make their way to a Maxwellian?

To answer this question we note that since the initial configuration of particles is generated randomly from a Poisson process with the density $`p(y,v,0)=p(|v|)`$, the actual (empirical) density of the particles, such as the one shown on Fig. 6 at $`\tau =0`$, does not exactly coincide with $`p(y,v,0)`$. Random fluctuations of the empirical density are typically of order $`O(1/L)`$. Hence, the actual initial distribution of particles can be thought of as a small perturbation of the function $`p(|v|)`$ and can be written as $`p(|v|)+\epsilon p_1(y,v)`$ with $`\epsilon =1/L`$ and some (random) function $`p_1(y,v)`$ of order one.

Now, we conjecture that for large $`L`$ the evolution of the particle system closely follows the solutions of the hydrodynamical equations derived in \[CLS, LPS\] with a perturbed initial density $`p(|v|)+\epsilon p_1(y,v)`$, rather than the stationary solution corresponding to the unperturbed density $`p(|v|)`$. In particular, if the latter solution is unstable, then small perturbations grow exponentially in time (in the units of $`\tau =t/L`$, of course), hence the system can be destabilized in $`\tau \mathrm{log}\epsilon =\mathrm{log}L`$ units of time. This would be in agreement with our analysis and numerical estimates in the end of Section 4. Hence, the instability of the hydrodynamical equations becomes a key issue.

When we were finishing the present paper, we received a message from E. Caglioti and E. Presutti who (a) proved that the hydrodynamical equations are stable when $`p(|v|)`$ is monotonically nonincreasing in $`|v|`$, i.e. $`p^{}(|v|)0`$, and (b) suggested that they might be unstable for our class on non-monotone $`p(|v|)`$. We checked the suggestion (b) for our particular density (3.1) and found that it was indeed correct; we proved that small perturbations grow exponentially in $`\tau `$. Furthermore, when starting with a perturbed initial density with $`\epsilon =0.01`$, we found numerically that the corresponding solution of the hydrodynamical equations resembled very well the evolution of the particle system described here, including large nearly harmonic oscillations of the piston during the interval $`10<\tau <30`$. Further work in this direction is currently under way \[CCLP\].

Conversely, when we simulated a particle dynamics with a nonincreasing initial density $`p(|v|)`$ the oscillations essentially disappeared666We tried a uniform “flat” function given by $`p(|v|)=1`$ for $`|v|v_{\mathrm{max}}`$ and a triangular one $`p(|v|)=1|v|/v_{\mathrm{max}}`$.. On the other hand, the particle velocity distribution still approached a Maxwellian, albeit at a somewhat slower pace.

Acknowledgement. We thank F. Bonetto, R. Dorfman, Ch. Gruber, J. Piasecki, N. Simanyi, Ya. Sinai, and V. Yakhot for many useful discussions. We thank H. van den Bedem and C. Lesort for technical assistance in numerical experiments. N. Chernov was partially supported by NSF grant DMS-0098788. J. Lebowitz was partially supported by NSF grant DMR-9813268 and by Air Force grant F49620-01-0154. This work was completed when the authors stayed at the Institute for Advanced Study with partial support by NSF grant DMS-9729992.


text:
# Fluctuation-dissipation relation for the Ising-Glauber model with arbitrary exchange couplings

## 1 Introduction

The knowledge about out-of-equilibrium processes is far from being as advanced as for systems at equilibrium. In particular, the fluctuation-dissipation theorem (FDT) which holds at equilibrium is known to be violated out-of-equilibrium. This theorem states that at equilibrium the response function $`R_{\mathrm{eq}}(t,s)`$ at time $`t`$ to an infinitesimal field branched to the system at time $`s`$ is related to the time-derivative of the two-time autocorrelation function $`C_{\mathrm{eq}}(t,s)`$:

$$R_{\mathrm{eq}}(t,s)=\beta \frac{C_{\mathrm{eq}}(t,s)}{s}.$$
(1)
In the Ising case, the response function reads $`R(t,s)=\frac{\delta \sigma _i(t)}{\delta h_i(s)}`$ and the correlation function $`C(t,s)=\sigma _i(t)\sigma _i(s)`$. Based on a mean-field study of spin-glasses, Cugliandolo et al  have conjectured that for asymptotically large times the FDT has to be corrected by a multiplicative factor $`X(t,s)`$ which moreover depends on time only through the correlation function:

$$R(t,s)=\beta X(C(t,s))\frac{C(t,s)}{s}.$$
(2)
Exact results have been obtained for the ferromagnetic Ising chain  that confirm this conjecture. The violation ratio $`X(t,s)`$ has been also numerically computed for many systems: 2d and 3d-Ising ferromagnets , 3d Edwards-Anderson model , 3d and 4d Gaussian Ising spin glasses , 2d Ising ferromagnet with dipolar interactions $`\mathrm{}`$ The response function has been exactly calculated in the case of the Ising chain  but the result is not written in terms of thermodynamic averages making difficult any conjecture on the way the FDT is violated.

We present an analytic study of the Glauber dynamics of the Ising model with arbitrary exchange couplings. The response function to an infinitesimal magnetic field is first calculated at equilibrium and the compatibility with the fluctuation-dissipation theorem is verified. In the last section, the response function is derived out-of-equilibrium and the expression is compared to the equilibrium equivalent.

## 2 Model, dynamics and general response function

### 2.1 The model and its dynamics

We consider a classical Ising model whose degrees of freedom are scalar variables $`\sigma _i=\pm 1`$ located on the nodes of a $`d`$-dimensional lattice. Let us denote $`\mathrm{}(\{\sigma \},t)`$ the probability to observe the system in the state $`\{\sigma \}`$ at time $`t`$. The dynamics of the system is supposed to be governed by the master equation:

$$\frac{\mathrm{}(\{\sigma \},t)}{t}=\underset{\{\sigma ^{}\}}{}\left[\mathrm{}(\{\sigma ^{}\},t)W(\{\sigma ^{}\}\{\sigma \})\mathrm{}(\{\sigma \},t)W(\{\sigma \}\{\sigma ^{}\})\right]$$
(3)
where $`W(\{\sigma \}\{\sigma ^{}\})`$ is the transition rate from the state $`\{\sigma \}`$ to the state $`\{\sigma ^{}\}`$ per unit of time. The normation is written $`_{\{\sigma ^{}\}}W(\{\sigma \}\{\sigma ^{}\})=1`$. The transition rates are defined by the stationarity condition

$$\underset{\{\sigma ^{}\}}{}\left[\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})W(\{\sigma ^{}\}\{\sigma \})\mathrm{}_{\mathrm{eq}}(\{\sigma \})W(\{\sigma \}\{\sigma ^{}\})\right]=0$$
(4)
where $`\mathrm{}_{\mathrm{eq}}(\{\sigma \})`$ is the equilibrium probability distribution which reads for the Ising model with general exchange-couplings:

$$\mathrm{}_{\mathrm{eq}}(\{\sigma \})=\frac{1}{𝒵}e^{\beta (\{\sigma \})}=\frac{1}{𝒵}e^{\beta _{k,l<k}J_{kl}\sigma _k\sigma _l}.$$
(5)
Equation (4) is satisfied when the detailed balance holds:

$$\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})W(\{\sigma ^{}\}\{\sigma \})=\mathrm{}_{\mathrm{eq}}(\{\sigma \})W(\{\sigma \}\{\sigma ^{}\}).$$
(6)
This last unnecessary but sufficient condition is fulfilled by the heat-bath single-spin flip dynamics defined by the following transition rates:

$`W(\{\sigma \}\{\sigma ^{}\})`$ $`=\left[{\displaystyle \underset{lk}{}}\delta _{\sigma _l,\sigma _l^{}}\right]{\displaystyle \frac{\underset{\sigma }{}\delta _{\sigma _k^{},\sigma }\mathrm{}_{\mathrm{eq}}(\{\sigma _l\}_{lk},\sigma _k=\sigma )}{_\sigma \mathrm{}_{\mathrm{eq}}(\{\sigma _l\}_{lk},\sigma _k=\sigma )}}`$ (7)
$`=\left[{\displaystyle \underset{lk}{}}\delta _{\sigma _l,\sigma _l^{}}\right]{\displaystyle \frac{e^{\beta _{lk}J_{kl}\sigma _l^{}\sigma _k^{}}}{_{\sigma =\pm 1}e^{\beta _{lk}J_{kl}\sigma _l^{}\sigma }}}.`$
Only the single-spin flip $`\sigma _k\sigma _k^{}`$ is allowed. The product of Kronecker deltas ensures that all other spins are not modified during the transition. The spin $`\sigma _k`$ takes after the transition the new value $`\sigma _k^{}`$ chosen according to the equilibrium probability distribution $`\mathrm{}_{\mathrm{eq}}(\{\sigma \})`$. All spins are successively updated during the evolution of the system. These transition rates reduce to Glauber’s ones  in the case of the Ising chain.

### 2.2 The response function

A magnetic field $`h_i`$ is coupled to the spin $`\sigma _i`$ between the times $`s`$ and $`s+\mathrm{\Delta }s`$ where $`\mathrm{\Delta }s`$ is supposed to be small. Between these instants, the transition rates are changed to $`W_h(\{\sigma \}\{\sigma ^{}\})`$ satisfying equation (7) with an additional term $`\beta h_i\sigma _i`$ in the Hamiltonian of the equilibrium probability distribution (5):

$$W_h(\{\sigma \}\{\sigma ^{}\})=\left[\underset{lk}{}\delta _{\sigma _l,\sigma _l^{}}\right]\frac{e^{\beta \left[_{lk}J_{kl}\sigma _l^{}\sigma _k^{}h_i\sigma _k^{}\delta _{k,i}\right]}}{_{\sigma =\pm 1}e^{\beta \left[_{lk}J_{kl}\sigma _l^{}\sigma h_i\sigma \delta _{k,i}\right]}}.$$
(8)
The transition rates are identical to the case $`h_i=0`$ apart from the single-spin flip involving the spin $`\sigma _i`$. Using the Bayes relation and the master equation (3), the average of the spin $`\sigma _i`$ at time $`t`$ can be expanded to lowest order in $`\mathrm{\Delta }s`$ under the following form:

$`\sigma _i(t)`$ $`={\displaystyle \underset{\{\sigma \}}{}}\sigma _i\mathrm{}(\{\sigma \},t)`$ (9)
$`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s+\mathrm{\Delta }s)\mathrm{}(\{\sigma ^{}\},s+\mathrm{\Delta }s)`$
$`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s+\mathrm{\Delta }s)[(1\mathrm{\Delta }s)\mathrm{}(\{\sigma ^{}\},s).`$
$`.+\mathrm{\Delta }s{\displaystyle \underset{\{\sigma ^{\prime \prime }\}}{}}\mathrm{}(\{\sigma ^{\prime \prime }\},s)W_h(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})].`$
$`W_h`$ being the only quantity depending on the magnetic field in equation (9), the first term disappears after a derivative with respect to the magnetic field. In the last term, the single-spin flip involving the spin $`\sigma _i`$, i.e. the one coupled to the magnetic field, is the only one that gives a non-vanishing contribution. The integrated response function between the time $`s`$ and $`s+\mathrm{\Delta }s`$ is $`\left[\frac{\sigma _i(t)}{h_i}\right]_{h_i0}`$ and the limit $`\mathrm{\Delta }s0`$ leads to the response function:

$$\underset{\mathrm{\Delta }s0}{lim}\frac{1}{\mathrm{\Delta }s}\left[\frac{\sigma _i(t)}{h_i}\right]_{h_i0}=\underset{\mathrm{\Delta }s0}{lim}\frac{1}{\mathrm{\Delta }s}_s^{s+\mathrm{\Delta }s}R_{ii}(t,u)𝑑u=R_{ii}(t,s).$$
(10)
According to the definition of the response function (10), the calculation can be limited to the lowest order in $`\mathrm{\Delta }s`$. As a consequence, the conditional probability in equation (9) can be replaced by $`\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)`$.

## 3 Equilibrium fluctuation-dissipation relation

### 3.1 The response function

If the system is in thermal equilibrium at time $`s`$, the probability distribution $`\mathrm{}(\{\sigma \},s)`$ has to be replaced in equation (9) by $`\mathrm{}_{\mathrm{eq}}(\{\sigma \})`$. Then the sum over $`\{\sigma ^{\prime \prime }\}`$ can be calculated:

$`{\displaystyle \underset{\{\sigma ^{\prime \prime }\}}{}}`$ $`\mathrm{}_{\mathrm{eq}}(\{\sigma ^{\prime \prime }\})W_h(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})`$ (11)
$`={\displaystyle \underset{\{\sigma ^{\prime \prime }\}}{}}{\displaystyle \frac{1}{𝒵}}e^{\beta _{k,l<k}J_{kl}\sigma _k^{\prime \prime }\sigma _l^{\prime \prime }}\times \left[{\displaystyle \underset{ji}{}}\delta _{\sigma _j^{\prime \prime },\sigma _j^{}}\right]{\displaystyle \frac{e^{\beta \left[_{ji}J_{ij}\sigma _i^{}\sigma _j^{}h_i\sigma _i^{}\right]}}{_{\sigma _i^{(3)}=\pm 1}e^{\beta \left[_{ji}J_{ij}\sigma _i^{(3)}\sigma _j^{}h_i\sigma _i^{(3)}\right]}}}`$
$`=\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})e^{\beta h_i\sigma _i^{}}{\displaystyle \frac{\mathrm{cosh}\left(\beta \underset{ji}{}J_{ij}\sigma _j^{}\right)}{\mathrm{cosh}\left(\beta \left[_{ji}J_{ij}\sigma _j^{}h_i\right]\right)}}.`$
The derivative with respect to $`h_i`$ followed by the limit $`h_i0`$ gives $`\beta \mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})\left[\sigma _i^{}\mathrm{tanh}\left(\beta _{ji}J_{ij}\sigma _j^{}\right)\right]`$. When replaced in equation (9), one obtains the following integrated response function

$$\left[\frac{\sigma _i(t)}{h_i}\right]_{h_i0}=\beta \mathrm{\Delta }s\left[\sigma _i(t)\sigma _i(s)\sigma _i(t)\mathrm{tanh}\left(\beta \underset{ji}{}J_{ij}\sigma _j(s)\right)\right]$$
(12)
and according to equation (10) the response function reads

$$R_{ii}(t,s)=\beta \left[\sigma _i(t)\sigma _i(s)\sigma _i(t)\mathrm{tanh}\left(\beta \underset{ji}{}J_{ij}\sigma _j(s)\right)\right].$$
(13)
### 3.2 Equivalence with the usual expression

The equilibrium fluctuation-dissipation theorem states that the response function is related to the time-derivative of the correlation function. We will now show that the expression (13) is compatible with this statement. At equilibrium, the correlation function is invariant under a time-translation. As a consequence, one has the relation

$$\frac{C_{ii}(t,s)}{t}=\frac{C_{ii}(t,s)}{s}.$$
(14)
Moreover, it can be shown that the conditional probability $`\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)`$ satisfies the master equation (3) for the time-variable $`t`$ and a similar equation for $`s`$:

$$\left(1\frac{}{s}\right)\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)=\underset{\{\sigma ^{\prime \prime }\}}{}\mathrm{}(\{\sigma \},t|\{\sigma ^{\prime \prime }\},s)W(\{\sigma ^{}\}\{\sigma ^{\prime \prime }\}).$$
(15)
The time-derivative of the correlation function is then

$`{\displaystyle \frac{C_{ii}(t,s)}{t}}`$ $`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\sigma _i^{}{\displaystyle \frac{\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)}{t}}\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})`$
$`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\sigma _i^{}{\displaystyle \frac{\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)}{s}}\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})`$
$`=C_{ii}(t,s)+{\displaystyle \underset{\{\sigma \},\{\sigma ^{}\},\{\sigma ^{\prime \prime }\}}{}}\sigma _i\sigma _i^{}\mathrm{}(\{\sigma \},t|\{\sigma ^{\prime \prime }\},s)W(\{\sigma ^{}\}\{\sigma ^{\prime \prime }\})\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\}).`$
The sum over $`\{\sigma ^{\prime \prime }\}`$ can be performed on the lines of equation (11) for the Ising-Glauber model. One obtains

$`{\displaystyle \underset{\{\sigma ^{}\}}{}}\sigma _i^{}W(\{\sigma ^{}\}\{\sigma ^{\prime \prime }\})\mathrm{}_{\mathrm{eq}}(\{\sigma ^{}\})`$ $`={\displaystyle \frac{1}{𝒵}}{\displaystyle \underset{\{\sigma ^{}\}}{}}\sigma _i^{}{\displaystyle \frac{\left[\underset{ji}{}\delta _{\sigma _j^{},\sigma _j^{\prime \prime }}\right]e^{\beta _{ji}J_{ij}\sigma _i^{\prime \prime }\sigma _j^{\prime \prime }}}{_{\sigma _i^{(3)}=\pm 1}e^{\beta _{ji}J_{ij}\sigma _i^{(3)}\sigma _j^{\prime \prime }}}}e^{\beta _{k,l<k}J_{kl}\sigma _k^{}\sigma _l^{}}`$ (17)
$`=\mathrm{}_{\mathrm{eq}}(\{\sigma ^{\prime \prime }\})\mathrm{tanh}\left(\beta {\displaystyle \underset{ji}{}}J_{ij}\sigma _j^{\prime \prime }\right).`$
Replacing in (3.2) and identifying with (13), one gets the expression of the equilibrium fluctuation-dissipation theorem for the Ising-Glauber model:

$$R_{ii}(t,s)=\beta \frac{C_{ii}(t,s)}{s}=\beta \left[\sigma _i(t)\sigma _i(s)\sigma _i(t)\mathrm{tanh}\left(\beta \underset{ji}{}J_{ij}\sigma _j(s)\right)\right].$$
(18)
## 4 General fluctuation-dissipation relation

Calculations have been made up to now within the hypothesis that the system is at equilibrium at time $`s`$. We will now derive a relation valid even if the system is far from equilibrium. We start with the expansion (9) of the average $`\sigma _i(t)`$. $`W_h`$ being the only quantity depending on the magnetic field, the integrated response function between the time $`s`$ and $`s+\mathrm{\Delta }s`$ is to lowest order in $`\mathrm{\Delta }s`$

$$\left[\frac{\sigma _i(t)}{h_i}\right]_{h_i0}=\mathrm{\Delta }s\underset{\genfrac{}{}{0pt}{}{\{\sigma \},\{\sigma ^{}\},}{\{\sigma ^{\prime \prime }\}}}{}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)\left[\frac{W_h}{h_i}(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\right]_{h_i0}\mathrm{}(\{\sigma ^{\prime \prime }\},s).$$
(19)
The derivative of the transition rate $`W_h`$ defined by equation (8) is easily taken and reads

$$\left[\frac{W_h}{h_i}(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\right]_{h_i0}=\beta W(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\left[\sigma _i^{}\mathrm{tanh}\left(\beta \underset{ji}{}J_{ij}\sigma _j^{}\right)\right].$$
(20)
Combining (19) and (20), the integrated response function is rewritten as

$`\left[{\displaystyle \frac{\sigma _i(t)}{h_i}}\right]_{h_i0}=`$ $`\mathrm{\Delta }s\beta {\displaystyle \underset{\genfrac{}{}{0pt}{}{\{\sigma \},\{\sigma ^{}\},}{\{\sigma ^{\prime \prime }\}}}{}}\sigma _i\sigma _i^{}\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)W(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\mathrm{}(\{\sigma ^{\prime \prime }\},s)`$ (21)
$`\mathrm{\Delta }s\beta {\displaystyle \underset{\genfrac{}{}{0pt}{}{\{\sigma \},\{\sigma ^{}\},}{\{\sigma ^{\prime \prime }\}}}{}}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)W(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})`$
$`\times \mathrm{tanh}\left(\beta {\displaystyle \underset{ji}{}}J_{ij}\sigma _j^{}\right)\mathrm{}(\{\sigma ^{\prime \prime }\},s).`$
In the first term of (21), the sum over $`\{\sigma ^{\prime \prime }\}`$ can be performed using the master equation (3) and the expression reorganised using the fact that after the limit $`h_i0`$, the conditional probability becomes invariant under a time-translation and depends only on the difference $`ts`$. One obtains

$`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\{\sigma \},\{\sigma ^{}\},}{\{\sigma ^{\prime \prime }\}}}{}}`$ $`\sigma _i\sigma _i^{}\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)W(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\mathrm{}(\{\sigma ^{\prime \prime }\},s)`$ (22)
$`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\sigma _i^{}\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)\left(1+{\displaystyle \frac{}{s}}\right)\mathrm{}(\{\sigma ^{}\},s)`$
$`=C_{ii}(t,s)+{\displaystyle \frac{C_{ii}(t,s)}{s}}{\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\sigma _i^{}{\displaystyle \frac{\mathrm{}}{s}}(\{\sigma \},t|\{\sigma ^{}\},s)\mathrm{}(\{\sigma ^{}\},s)`$
$`=C_{ii}(t,s)+{\displaystyle \frac{C_{ii}(t,s)}{s}}+{\displaystyle \frac{C_{ii}(t,s)}{t}}.`$
Similarly, in the second term of equation (21) the sum over $`\{\sigma ^{\prime \prime }\}`$ can be performed using the master equation (3) and leads to

$`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\{\sigma \},\{\sigma ^{}\},}{\{\sigma ^{\prime \prime }\}}}{}}`$ $`\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)\mathrm{tanh}\left(\beta {\displaystyle \underset{ji}{}}J_{ij}\sigma _j^{}\right)W(\{\sigma ^{\prime \prime }\}\{\sigma ^{}\})\mathrm{}(\{\sigma ^{\prime \prime }\},s)`$ (23)
$`={\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)\mathrm{tanh}\left(\beta {\displaystyle \underset{ji}{}}J_{ij}\sigma _j^{}\right)\left(1+{\displaystyle \frac{}{s}}\right)\mathrm{}(\{\sigma ^{}\},s)`$
$`=\left(1+{\displaystyle \frac{}{s}}+{\displaystyle \frac{}{t}}\right){\displaystyle \underset{\{\sigma \},\{\sigma ^{}\}}{}}\sigma _i\mathrm{}(\{\sigma \},t|\{\sigma ^{}\},s)\mathrm{tanh}\left(\beta {\displaystyle \underset{ji}{}}J_{ij}\sigma _j^{}\right)\mathrm{}(\{\sigma ^{}\},s).`$
Putting together equations (22) and (23) into (21), one can then evaluate the response function (10):

$$R_{ii}(t,s)=\beta \left(1+\frac{}{s}+\frac{}{t}\right)\left[C_{ii}(t,s)\sigma _i(t)\mathrm{tanh}\left(\beta \underset{ji}{}J_{ij}\sigma _j(s)\right)\right].$$
(24)
The hyperbolic tangent is the equilibrium value of a single spin $`\sigma _i`$ in the Weiss field created by all other spins at time $`s`$. Moreover, the hyperbolic tangent being an odd function of its argument and since $`\sigma _i^2=1`$, one can rewrite (24) in the following form:

$$R_{ii}(t,s)=\beta \left(1+\frac{}{s}+\frac{}{t}\right)\left[C_{ii}(t,s)\sigma _i(t)\sigma _i(s)\mathrm{tanh}\left(\beta ϵ_i(s)\right)\right].$$
(25)
where $`ϵ_i(s)=_{ji}J_{ij}\sigma _i(s)\sigma _j(s)`$ is the energy at site $`i`$.

## 5 Conclusion

The expression (24) (or (25)) of the response function makes apparent the way the FDT is violated. At equilibrium, all averages are invariant under a time-translation, i.e. depend only on $`ts`$ so the derivatives over $`s`$ and $`t`$ cancel and at equilibrium the average in equation (24) is equal to the time-derivative of the correlation function as shown in section 3.2. The expression (24) is very general. Since the response function is expressed as thermodynamic averages, it does not depend on the initial conditions. Moreover, it does not depend on the dimension of space or whether or not, the lattice is regular. It applies to any configuration of exchange couplings so for a disordered system, the average response function is just the average of expression (24) over all couplings configurations. Equations (22) and (23) suggest that the presence of the term $`1+\frac{}{s}+\frac{}{t}`$ in the expression of the response function may be a general feature of the response function out-of-equilibrium. Exact results for other models are highly desirable. The last term of (24) whose equilibrium limit is the time-derivative of the correlation function may also be evaluated perturbatively for certain models.

This work may also help to improve numerical simulations. Up to now, the FDT has been numerically checked by calculating the time-integrated response function. One can note that equation (19) can be used to calculate the response function directly during a Monte Carlo simulation.

The laboratoire de Physique des Matériaux is Unité Mixte de Recherche CNRS number 7556. The author gratefully acknowledges Dragi Karevski, Loïc Turban, Bertrand Berche and Alan Picone for critical reading of the manuscript.

## References