text:
# 1 Background: on dynamical urn models

## 1 Background: on dynamical urn models

Dynamical urn models are simplified models of physical reality, which have always played an important role in the elucidation of conceptual problems of statistical mechanics and probability theory.

The ancestor and prototype of this class of models is the Ehrenfest urn model . It was devised by P. and T. Ehrenfest, in their attempt to critically review Boltzmann’s $`H`$theorem. Consider $`N`$ balls, labeled from 1 to $`N`$, which are distributed in two urns (or boxes). At each time step a ball is chosen at random (i.e., an integer between 1 and $`N`$ is chosen at random), and moved from the box in which it is, to the other box. Let $`N_1`$ (respectively, $`N_2`$) be the numbers of balls in box number 1 (respectively, number 2) ($`N_1+N_2=N`$). If the process is repeated indefinitely, for any initial condition the system will relax to equilibrium, characterized by a binomial distribution of balls in, say, box number 1:

$$f_{k,\mathrm{eq}}=𝒫(N_1=k)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)\frac{1}{2^N}.$$
(1.1)
This result is both intuitively translucent, and easy to derive (see below). The partition function of the system is equal to (for $`2`$ boxes and $`N`$ balls)

$$Z(2,N)=\frac{2^N}{N!}.$$
(1.2)
There are $`2^N`$ possibilities of distributing the balls amongst the 2 boxes, however all $`N!`$ labelings of the balls are equivalent.

Finding the distribution of balls after $`n`$ steps, $`f_k(n)=𝒫(N_1(n)=k)`$, requires more effort. Kohlrausch and Schrödinger  found the master equation for $`f_k(n)`$, which they interpreted as the probability distribution of the position at time $`n`$ of a random walker (played here by $`N_1(n)`$). The full solution of this master equation was given later on by Kac, Siegert and Hess .

The Ehrenfest model is at the origin of a whole class of dynamical urn models, which we name the *Ehrenfest class*. They generalize the original Ehrenfest model in two ways, which we detail successively.

The first generalization consists in considering $`M`$ boxes instead of two. Then, at equilibrium, the joint distribution of the occupation numbers $`N_1`$, $`N_2`$, $`\mathrm{}`$, $`N_M`$, with

$$\underset{i=1}{\overset{M}{}}N_i=N,$$

is the multinomial distribution

$$𝒫(N_1=k_1,\mathrm{},N_M=k_M)=\left(\genfrac{}{}{0pt}{}{N}{k_1\mathrm{}k_M}\right)\frac{1}{M^N},$$

as a simple reasoning shows. The marginal distribution of $`N_1`$ is obtained by summation upon the other variables, and reads

$$f_{k,\mathrm{eq}}=𝒫(N_1=k)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)\frac{1}{M^N},$$

which is a simple generalization of (1.1). In the thermodynamic limit $`N\mathrm{}`$, $`M\mathrm{}`$, with fixed density $`\rho =N/M`$, this yields a Poisson law:

$$f_{k,\mathrm{eq}}=\mathrm{e}^\rho \frac{\rho ^k}{k!}.$$

The partition function of the system is equal to (for $`M`$ boxes and $`N`$ balls)

$$Z(M,N)=\frac{M^N}{N!},$$
(1.3)
as a simple generalization of the reasoning leading to (1.2) shows.

The dynamics of the Ehrenfest model and of its generalization to $`M`$ boxes takes place at infinite temperature because there is no constraint on the move of the drawn ball. The second direction of generalization consists in defining these models at finite temperature, by introducing energy. We assume that the energy is a sum of contributions of independent boxes:

$$E(N_1,\mathrm{},N_M)=\underset{i=1}{\overset{M}{}}E(N_i),$$

and choose a rule obeying detailed balance for the move of the drawn ball. For instance, for the Metropolis rule, the move is allowed with probability $`\mathrm{min}(1,\mathrm{exp}(\beta \mathrm{\Delta }E))`$. Heat-bath dynamics is another possible choice (see below).

The fundamental ingredients for the definition of the models belonging to the Ehrenfest class are therefore

the statistics: a ball is chosen at random, and put in a box chosen at random,

the choice of the energy of a box $`E(N_i)`$, and of a dynamical rule (Metropolis or heat-bath),

the geometry: for instance, boxes may be ordered on a line, or on the contrary be all connected. For short, we designate the latter geometry as the mean-field case.

The backgammon model  is a representative of the Ehrenfest class, corresponding to the choice (where $`\delta `$ is Kronecker symbol)

$$E(N_i)=\delta (N_i,0).$$
(1.4)
This model has been extensively studied, mainly in its mean-field formulation .

The dynamical (and equilibrium) properties of the Ehrenfest class depend crucially on the choice of statistics described in (i), which we will briefly refer to as the *ball-box statistics*. However, other choices are possible, which define new classes of dynamical urn models. The class which we will refer to for short as the *Monkey class*, because it corresponds to the image of a monkey playing at exchanging balls between boxes, is defined by:

the statistics: a box is chosen at random, from which any ball is drawn, and put in another box, chosen at random (*box-box statistics*),

the choice of energy and dynamical rule (as above),

the geometry (as above).

A first example of a model belonging to this class corresponds to taking definition (1.4) for the energy . This model, referred to as model B in , possesses non-trivial dynamical properties .

A second example, inspired from quantum gravity, corresponds to taking

$$E(N_i)=\mathrm{ln}(N_i+1).$$
(1.5)
It presents interesting properties both at equilibrium  and in nonequilibrium situations . In contrast with the backgammon model, or with model B, it possesses a phase transition between a fluid phase and a condensed phase at finite temperature . We will refer to this model as the *zeta urn model*, for reasons which will appear clear in the sequel. The present work is entirely devoted to the study of the nonequilibrium behavior of the zeta urn model in the mean-field geometry.

Before specializing to this model, let us first present, in parallel, some formalism which applies to the two classes of models defined above, in order to underline the fundamental role played by the choice of statistics for both the equilibrium and nonequilibrium properties of the models. Note that the equilibrium properties of the dynamical urn models defined above are independent of the geometry, because boxes are independent.

For the Ehrenfest class, the partition function reads

$$Z(M,N)=\underset{N_1}{}\mathrm{}\underset{N_M}{}\frac{p_{N_1}}{N_1!}\mathrm{}\frac{p_{N_M}}{N_M!}\delta (\underset{i}{}N_i,N),$$
(1.6)
where

$$p_{N_i}=\mathrm{e}^{\beta E(N_i)}$$

is the unnormalized Boltzmann weight attached to box number $`i`$. For the Monkey class we have

$$Z(M,N)=\underset{N_1}{}\mathrm{}\underset{N_M}{}p_{N_1}\mathrm{}p_{N_M}\delta (\underset{i}{}N_i,N).$$
(1.7)
Using the integral representation $`2\mathrm{i}\pi \delta (m,n)=dzz^{mn1}`$, we obtain

$$Z(M,N)=\frac{\mathrm{d}z}{2\mathrm{i}\pi z^{N+1}}\left[P(z)\right]^M,$$
(1.8)
where

$`P(z)={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{p_k}{k!}}z^k`$ $`(\mathrm{Ehrenfest}),`$ (1.9)
$`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}p_kz^k`$ $`(\mathrm{Monkey}).`$ (1.10)
The equilibrium properties of the models are therefore entirely encoded in the temperature-dependent generating series $`P(z)`$. The presence or absence of the factorial term $`k!`$ has direct implication on the analytic structure of the series, and by consequence on the possible existence of a phase transition at finite temperature.

The equilibrium probability distribution of the occupation number $`N_1`$ reads, for the Ehrenfest class,

$`f_{k,\mathrm{eq}}`$ $`=`$ $`𝒫(N_1=k)=\delta (N_1,k)`$ (1.11)
$`=`$ $`{\displaystyle \frac{1}{Z(M,N)}}{\displaystyle \underset{N_1}{}}\mathrm{}{\displaystyle \underset{N_M}{}}\delta (N_1,k){\displaystyle \frac{p_{N_1}}{N_1!}}\mathrm{}{\displaystyle \frac{p_{N_M}}{N_M!}}\delta ({\displaystyle \underset{i}{}}N_i,N)`$
$`=`$ $`{\displaystyle \frac{p_k}{k!}}{\displaystyle \frac{Z(M1,Nk)}{Z(M,N)}}.`$
For the Monkey class one obtains

$$f_{k,\mathrm{eq}}=p_k\frac{Z(M1,Nk)}{Z(M,N)}.$$
(1.12)
At infinite temperature, equation (1.8), together with (1.9) and (1.10), respectively lead to (1.3) and to

$$Z(M,N)=\frac{(M+N1)!}{(M1)!N!}.$$
(1.13)
In the thermodynamic limit at fixed density $`\rho `$, the free energy per box is defined as

$$\beta F=\underset{M\mathrm{}}{lim}\frac{1}{M}\mathrm{ln}Z(M,N),NM\rho .$$

At infinite temperature, equations (1.3) and (1.13) yield

$`\underset{\beta 0}{lim}\beta F=\rho \mathrm{ln}\rho \rho `$ $`(\mathrm{Ehrenfest}),`$
$`=\rho \mathrm{ln}\rho (\rho +1)\mathrm{ln}(\rho +1)`$ $`(\mathrm{Monkey}).`$
At finite temperature, the free energy can be obtained by evaluating the contour integral in (1.8) by the saddle-point method. The saddle-point value $`z_s`$ of $`z`$ is a function of temperature and density through the equation

$$\frac{z_sP^{}(z_s)}{P(z_s)}=\rho ,$$
(1.14)
and the free energy per box reads

$$\beta F=\rho \mathrm{ln}z_s\mathrm{ln}P(z_s).$$
(1.15)
Similarly, we obtain the following expressions for the equilibrium occupation probabilities in the thermodynamic limit

$`f_{k,\mathrm{eq}}={\displaystyle \frac{p_k}{k!}}{\displaystyle \frac{z_s^k}{P(z_s)}}`$ $`(\mathrm{Ehrenfest}),`$ (1.16)
$`=p_k{\displaystyle \frac{z_s^k}{P(z_s)}}`$ $`(\mathrm{Monkey}).`$ (1.17)
This formalism will be illustrated below on the zeta urn model, studied in this work.

For both classes of models, the temporal evolution of the occupation probability

$$f_k(t)=𝒫(N_1(t)=k)$$

is given by the master equation

$$\frac{\mathrm{d}f_k(t)}{\mathrm{d}t}=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\left(\pi _{k+1,\mathrm{}}+\pi _{\mathrm{},k1}\pi _{k,\mathrm{}}\pi _{\mathrm{},k}\right),$$
(1.18)
where $`\pi _{k,\mathrm{}}`$ denotes the contribution of a move from the departure box, containing $`k`$ balls, to the arrival box, containing $`\mathrm{}`$ balls. Restricting our study to the mean-field case, we have (for $`k,\mathrm{}0`$)

$`\pi _{k,\mathrm{}}=kf_kf_{\mathrm{}}W_{k,\mathrm{}}(1\delta _{k,0})`$ $`(\mathrm{Ehrenfest}),`$
$`=f_kf_{\mathrm{}}W_{k,\mathrm{}}(1\delta _{k,0})`$ $`(\mathrm{Monkey}),`$ (1.19)
where the term $`1\delta _{k,0}`$ accounts for the fact that the departure box cannot be empty, and where the acceptance rate $`W_{k,\mathrm{}}`$ depends on the dynamics chosen. With the Metropolis rule we have

$$W_{k,\mathrm{}}=\mathrm{min}(1,\frac{p_{k1}p_{\mathrm{}}}{p_kp_{\mathrm{}+1}}),$$

while the heat-bath rule leads to

$$W_{k,\mathrm{}}=\frac{p_{\mathrm{}+1}}{p_{\mathrm{}}}\left(\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}f_{\mathrm{}}\frac{p_{\mathrm{}+1}}{p_{\mathrm{}}}\right)^1,$$
(1.20)
which only depends on the label $`\mathrm{}`$ of the arrival box.

In all cases, equation (1.18) can be seen as the master equation of a random walk for $`N_1`$, i.e., over the positive integers $`k=0,1,\mathrm{}`$,

$`{\displaystyle \frac{\mathrm{d}f_k(t)}{\mathrm{d}t}}`$ $`=`$ $`\mu _{k+1}f_{k+1}+\lambda _{k1}f_{k1}\left(\mu _k+\lambda _k\right)f_k(k1),`$ (1.21)
$`{\displaystyle \frac{\mathrm{d}f_0(t)}{\mathrm{d}t}}`$ $`=`$ $`\mu _1f_1\lambda _0f_0,`$ (1.22)
generalizing the result of Kohlrausch and Schrödinger for the Ehrenfest model. In these equations,

$$\lambda _k=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{\pi _{\mathrm{},k}}{f_k},\mu _k=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{\pi _{k,\mathrm{}}}{f_k}$$

are, respectively, the hopping rate to the right, corresponding to $`N_1=kN_1=k+1`$, and to the left, corresponding to $`N_1=kN_1=k1`$. The equation for $`f_0`$ is special because one cannot select an empty box as a departure box, nor can $`N_1`$ be negative, hence $`\lambda _1=\mu _0=0`$. In other words a partially absorbing barrier is present at site $`k=0`$. The random walk is locally biased, to the right or to the left, according to whether its velocity $`\lambda _k\mu _k`$ is positive or negative, respectively. It is easy to verify that (1.21) and (1.22) preserve the sum rules

$`{\displaystyle \underset{k}{}}f_k(t)`$ $`=`$ $`1,`$ (1.23)
$`{\displaystyle \underset{k}{}}kf_k(t)`$ $`=`$ $`N_1(t)=\rho ,`$ (1.24)
expressing respectively the conservation of probability and of the number of particles.

The equilibrium occupation probabilities (1.16), (1.17) are recovered as the unique stationary state ($`\mathrm{d}f_k/\mathrm{d}t=0`$) of the master equations (1.21), (1.22). In contrast, except at infinite temperature, where $`W_{k,\mathrm{}}=1`$, so that the rates $`\lambda _k`$ and $`\mu _k`$ simplify, the master equations cannot be solved explicitly. The difficulty comes from the fact that the rates $`\lambda _k`$ and $`\mu _k`$ are functions of the $`f_k`$, hence the master equations are non-linear. However, as will be illustrated by the present work, the long-time behavior of these equations is amenable to analytic computations.

To summarize, for the two classes of dynamical urn models described above, finding the properties of their equilibrium states is in general easy, but attaining dynamical properties is much more difficult to achieve. However if one restricts the study to the fluctuating number of balls in a given box, denoted by $`N_1(t)`$ all throughout this paper, it is, in some cases, possible to predict the long-time behavior of its probability distribution, $`f_k(t)`$, and also of its two-time correlation and response functions. This has been done in previous studies for the backgammon model .

In the present work we focus our interest on the mean-field dynamical urn model defined by the choice of energy (1.5), and box-box statistics –the zeta urn model. We pursue the investigation of its nonequilibrium properties, initiated in . The next section is devoted to a more complete presentation of the model and to an outline of this paper.

## 2 The zeta urn model

At equilibrium, the zeta urn model is defined by its partition function (1.7), where the Boltzmann weight

$$p_{N_i}=\mathrm{e}^{\beta E(N_i)}=(N_i+1)^\beta $$
(2.1)
corresponds to the choice of energy (1.5. The equilibrium phase diagram of the model  easily follows from the analysis of the previous section (see (1.8)–(1.17)).

At low enough temperature $`(\beta >2)`$, there is a finite critical density:

$$\rho _c=\frac{P^{}(1)}{P(1)}=\frac{\zeta (\beta 1)\zeta (\beta )}{\zeta (\beta )},$$

where $`\zeta `$ denotes Riemann’s zeta function.

In the fluid phase $`(\rho <\rho _c)`$, the equilibrium distribution (1.17) decays exponentially, since $`z_s<1`$.

At the critical density $`(\rho =\rho _c)`$, we have $`z_s=1`$, hence

$$f_{k,\mathrm{eq}}=\frac{p_k}{P(1)}=\frac{(k+1)^\beta }{\zeta (\beta )},$$
(2.2)
which is known as the zeta distribution. In the regular part of the critical line ($`\beta >3`$), the mean squared population is finite, and equal to

$$N_i^2=\underset{k=0}{\overset{\mathrm{}}{}}k^2f_{k,\mathrm{eq}}=\mu _c=\frac{\zeta (\beta 2)2\zeta (\beta 1)+\zeta (\beta )}{\zeta (\beta )},$$

while it is infinite in the strong-fluctuation case ($`2<\beta <3`$). Throughout the following, we shall restrict the study to the regular part of the critical line.

In the condensed phase $`(\rho >\rho _c)`$, a macroscopic condensate of particles appears. Indeed, equation (2.2) still applies to all the boxes but one, in which an extensive number of particles, of order $`NM\rho _c=M(\rho \rho _c)`$, is condensed.

The dynamical definition of the model was given in the previous section. For heat-bath dynamics, equation (1.18), together with (1.19) and (1.20), leads to

$`{\displaystyle \frac{\mathrm{d}f_k(t)}{\mathrm{d}t}}=f_{k+1}(t)+\sigma (t)r_{k1}f_{k1}(t)\left(1+\sigma (t)r_k\right)f_k(t)(k1),`$
$`{\displaystyle \frac{\mathrm{d}f_0(t)}{\mathrm{d}t}}=f_1(t)\sigma (t)r_0f_0(t),`$ (2.3)
with

$$r_k=\frac{p_{k+1}}{p_k}=\frac{f_{k+1,\mathrm{eq}}}{f_{k,\mathrm{eq}}}=\left(\frac{k+1}{k+2}\right)^\beta $$

and

$$\sigma (t)=\frac{1f_0(t)}{{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}r_kf_k(t)}.$$

We assume that at time $`t=0`$ the system is quenched from its infinite-temperature equilibrium state to a finite temperature $`T=1/\beta `$. Hence, by (1.14) and (1.17), the initial occupation probabilities read

$$f_k(0)=\frac{\rho ^k}{(1+\rho )^{k+1}}.$$
(2.4)
In the fluid phase ($`\rho <\rho _c`$), the equilibrium distribution $`f_{k,\mathrm{eq}}`$ (1.17) is a stationary solution of (2.3), corresponding to $`\sigma _{\mathrm{eq}}=z_s`$. The convergence of $`f_k(t)`$ towards $`f_{k,\mathrm{eq}}`$ is characterized by a finite relaxation time, depending on $`\beta `$ and $`\rho `$.

The long-time behavior of the distribution $`f_k(t)`$, both at criticality ($`\rho =\rho _c`$) and in the condensed phase ($`\rho >\rho _c`$), is the subject of the next section. In section 4 we establish the dynamical equations obeyed by the two-time correlation and response functions of the population of a given box. It is indeed well-known that nonequilibrium properties are more fully revealed by two-time observables . Section 5 is devoted to the analysis of the equilibrium properties of these functions at criticality. Section 6 is devoted to the analysis of their nonequilibrium properties at criticality, including the violation of the fluctuation-dissipation theorem. In section 7 we briefly investigate the nonequilibrium properties of the model in the condensed phase.

## 3 Long-time behavior of occupation probabilities

### 3.1 At criticality ($`\rho =\rho _c`$)

We investigate how, starting from the disordered initial condition (2.4), the occupation probabilities $`f_k(t)`$ converge toward their critical equilibrium values $`f_{k,\mathrm{eq}}`$, given by (2.2), which are the stationary solutions of equations (2.3) corresponding to $`\sigma _{\mathrm{eq}}=1`$. In analogy with the analysis done in , we anticipate that $`\sigma (t)`$ converges to this value as a power law:

$$\sigma (t)1+At^\omega ,$$

and consider two regimes:

Regime I: $`k`$ fixed and $`t1`$

This is the “short-distance” regime, considering $`k`$ as the position of a fictitious random walker, as in section 1. It is therefore analogous to the Porod regime for phase-ordering systems .

Setting

$$f_k(t)f_{k,\mathrm{eq}}\left(1+v_kt^\omega \right),$$
(3.1)
equation (2.3) yields

$$v_k=v_0+Ak,$$
(3.2)
where $`v_0`$ and $`A`$ are determined below.

Regime II: $`k`$ and $`t`$ simultaneously large (scaling regime)

In this regime, we look for a similarity solution to (2.3) of the form

$$f_k(t)f_{k,\mathrm{eq}}F(u),u=kt^{1/2}.$$
(3.3)
The structure of the master equations (2.3) indeed dictates that the scaling variable is the combination $`kt^{1/2}`$. Starting from a random initial condition, for a large but finite time $`t`$, and for $`k`$ much smaller than an ordering size of order $`t^{1/2}`$, the system looks critical, i.e., the distribution $`f_k(t)`$ has essentially converged toward the equilibrium distribution $`f_{k,\mathrm{eq}}`$. This implies $`F(0)=1`$. To the contrary, for $`kt^{1/2}`$, the system still looks disordered, i.e., the $`f_k(t)`$ fall off very fast. Hence $`F(\mathrm{})=0`$. It will indeed be shown below that $`f_k(t)\mathrm{exp}(k^2/(4t))`$, as the scaling function falls of very fast for $`u1`$: $`F(u)\mathrm{exp}(u^2/4)`$. Note the close analogy between the present situation and critical coarsening for ferromagnetic spin systems , where the scaling variable is $`\left|𝐫\right|t^{1/z}`$, and where the role of $`f_k(t)`$ is played by the equal-time correlation function $`C(\left|𝐫\right|,t)`$.

In order to determine the exponent $`\omega `$, we use the sum rules (1.23) and (1.24), which yield respectively

$`t^\omega (v_0+A\rho _c)=t^{(\beta 1)/2}I_1,I_1={\displaystyle \frac{1}{\zeta (\beta )}}{\displaystyle _0^{\mathrm{}}}u^\beta \left(1F(u)\right)du,`$ (3.4)
$`t^\omega (v_0\rho _c+A\mu _c)=t^{(\beta 2)/2}I_2,I_2={\displaystyle \frac{1}{\zeta (\beta )}}{\displaystyle _0^{\mathrm{}}}u^{1\beta }\left(1F(u)\right)du.`$ (3.5)
These equations are compatible only if the right-hand side of equation (3.4) is subleading. In the case of a regular critical point ($`\beta >3`$), we thus obtain

$$\omega =\frac{\beta 2}{2}>\frac{1}{2}$$

and

$$A=\frac{I_2}{\mu _c\rho _c^2},v_0=\frac{\rho _cI_2}{\mu _c\rho _c^2}.$$

Inserting the form (3.3) into (2.3), and using the fact that $`\omega >1/2`$, leads to the differential equation

$$𝒟F(u)=0,$$
(3.6)
where $`𝒟`$ is the linear differential operator

$$𝒟=\frac{\mathrm{d}^2}{\mathrm{d}u^2}+\left(\frac{u}{2}+\frac{\beta }{u}\right)\frac{\mathrm{d}}{\mathrm{d}u}.$$
(3.7)
The solution of (3.6) is

$$F(u)=\frac{2^\beta }{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)}_u^{\mathrm{}}y^\beta \mathrm{e}^{y^2/4}dy.$$
(3.8)
We present in Appendix A an alternative way of solving equation (3.6), using the Mellin transformation. In particular (A.2) yields the explicit expression

$$I_2=\frac{M_{1F}(\beta 2)}{\zeta (\beta )}=\frac{\pi ^{1/2}\mathrm{\hspace{0.17em}2}^{1\beta }}{(\beta 2)\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\zeta (\beta )}.$$

As an illustration of the above, let us determine how the variance of the population of box number 1, $`VarN_1(t)=N_1(t)^2N_1(t)^2=N_1(t)^2\rho _c^2`$, converges at long times to its equilibrium value $`\mu _c\rho _c^2`$. Only the scaling regime matters for the long-time behavior of this quantity, and of all the quantities at criticality to be considered hereafter. Using (3.3), we obtain

$$VarN_1(t)(\mu _c\rho _c^2)\mathrm{\Delta }t^{(\beta 3)/2},$$

with, by (A.2),

$$\mathrm{\Delta }=\frac{M_{1F}(\beta 3)}{\zeta (\beta )}=\frac{2^{3\beta }}{(\beta 3)\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\zeta (\beta )}.$$

### 3.2 In the condensed phase ($`\rho >\rho _c`$)

We set

$$\sigma (t)1+At^{1/2},$$

and consider the same two regimes as at criticality (see  for more details).

Regime I: $`k`$ fixed and $`t1`$ (short-distance regime)

Equations (3.1) and (3.2) still hold, but now with $`\omega =1/2`$, for any $`\beta >2`$.

Regime II: $`k`$ and $`t`$ simultaneously large (scaling regime)

Looking for a similarity solution of equations (2.3), of the form

$$f_k(t)\frac{F(u)}{t},u=kt^{1/2},$$
(3.9)
we obtain for the scaling function $`F(u)`$ the linear differential equation

$$\frac{\mathrm{d}^2F}{\mathrm{d}u^2}+\left(\frac{u}{2}A+\frac{\beta }{u}\right)\frac{\mathrm{d}F}{\mathrm{d}u}+\left(1\frac{\beta }{u^2}\right)F=0.$$
(3.10)
Following , we notice that the amplitude $`A`$ is determined by the fact that equation (3.10) has an acceptable solution $`F(u)`$, vanishing as $`u0`$ and $`u\mathrm{}`$. The normalization of the solution $`F(u)`$ is determined by the sum rule (1.24), which yields

$$_0^{\mathrm{}}uF(u)du=\rho \rho _c.$$
(3.11)
The parameter $`v_0`$ entering Regime I is determined by the sum rule (1.23), leading to

$$v_0=A\rho _c_0^{\mathrm{}}F(u)du.$$

At variance with equation (3.6), the differential equation (3.10) cannot be solved in closed form. It can be recast in Schwarzian form, without first derivative, by setting

$$F(u)=uY(u)H(u),$$

with

$$Y(u)=u^{\beta /21}\mathrm{e}^{Au/2u^2/8}.$$
(3.12)
We thus obtain for $`H(u)`$ a differential equation of the form

$$H=0,$$
(3.13)
with

$$=\frac{\mathrm{d}^2}{\mathrm{d}u^2}+W(u),$$
(3.14)
and where the potential $`W(u)`$ reads

$$W(u)=\frac{u^2}{16}\frac{Au}{4}+\frac{\beta 3+A^2}{4}\frac{A\beta }{2u}+\frac{\beta (\beta +2)}{4u^2}.$$
(3.15)
Equation (3.13) is a biconfluent Heun equation . We will therefore refer to $``$ as the Heun operator, denoting its discrete eigenvalues by $`E_n`$ and the corresponding eigenfunctions bu $`H_n(u)`$. Equation (3.13) implies that the ground-state eigenvalue reads $`E_0=0`$, while the associated eigenfunction $`H_0(u)`$ is simply related to the scaling function $`F(u)`$:

$$F(u)=cuY(u)H_0(u).$$
(3.16)
By (3.11) we have

$$c=\frac{\rho \rho _c}{{\displaystyle _0^{\mathrm{}^{}}}u^2Y(u)H_0(u)du}.$$
(3.17)
The spectrum of the operator $``$, and related quantities such as the scaling function $`F(u)`$, can be further investigated in the limiting regimes of high and low temperature.

High temperature $`(\beta 0)`$

The analysis of this limiting situation will be helpful in section 7, although it is of no direct physical relevance, since the condensed phase only exists at low enough temperature $`(\beta >2)`$.

For $`\beta =0`$ and $`A=0`$, the potential (3.15) becomes $`W(u)=u^2/163/4`$, so that $``$ is the Hamiltonian of a harmonic oscillator on the half-line $`u0`$, up to a scale, with Dirichlet boundary condition at $`u=0`$. Its spectrum is $`E_n=n`$ $`(n=0,1,\mathrm{})`$, and the (unnormalized) ground-state eigenfunction reads

$$H_0(u)=u\mathrm{e}^{u^2/8},$$
(3.18)
hence

$$F(u)=\frac{\rho \rho _c}{2\pi ^{1/2}}u\mathrm{e}^{u^2/4}.$$
(3.19)
Perturbation theory can then be used to determine the small-$`\beta `$ behavior of the amplitude $`A`$. Expressing that the lowest eigenvalue reads $`E_0=0`$ with no correction, we obtain

$$_0^{\mathrm{}}\left(Au+\beta +\frac{2\beta }{u^2}\right)H_0^2(u)du+\mathrm{}=0,$$

hence

$$A\frac{\pi ^{1/2}}{2}\beta .$$
(3.20)
Low temperature $`(\beta \mathrm{})`$

In this other limiting situation, the operator $``$ simplifies as follows. If we rescale $`A`$ and $`u`$ according to $`Au\beta ^{1/2}`$, all terms in the expression (3.15) for the potential $`W(u)`$ scale as $`\beta `$. In other words, $`\beta `$ plays the role of $`1/\mathrm{}`$ in Quantum Mechanics, and $`\beta \mathrm{}`$ is the semi-classical regime. Expressing that the minimum of the potential is at zero yields the first estimates $`Au(2\beta )^{1/2}`$.

A more refined analysis consists in expanding the potential $`W(u)`$ around its minimum. Setting

$$A=(2\beta )^{1/2}a,u=(2\beta )^{1/2}+\beta ^{1/8}v,$$
(3.21)
the operator $``$ becomes

$$\frac{a^22}{4}+\beta ^{1/4}\left(\frac{\mathrm{d}^2}{\mathrm{d}v^2}+\frac{a}{2^{5/2}}v^2\right).$$
(3.22)
The expression in the parentheses is again proportional to the Hamiltonian of a harmonic oscillator. Expressing that the lowest eigenvalue of the right-hand-side of (3.22) is $`E_0=0`$ determines $`a=2^{1/2}2^{1/2}\beta ^{1/4}`$, hence

$$A(2\beta )^{1/2}2^{1/2}+2^{1/2}\beta ^{1/4}.$$

The spectrum of $``$ reads $`E_n\beta ^{1/4}n`$ ($`n=0,1,\mathrm{}`$), and the (unnormalized) ground-state eigenfunction reads $`H_0(u)\mathrm{e}^{v^2/4}`$, i.e.,

$$H_0(u)\mathrm{exp}\left(\frac{(u(2\beta )^{1/2})^2}{4\beta ^{1/4}}\right).$$
(3.23)
Both this expression and the function $`Y(u)`$ become singular in the $`\beta \mathrm{}`$ limit, so that a direct analysis of equation (3.10) is needed in order to derive the behavior of the scaling function $`F(u)`$ at low temperature. The above analysis suggests to look for a solution depending on the scaling variable $`y=(2\beta )^{1/2}u`$. The term involving the second-order derivative $`\mathrm{d}^2F/\mathrm{d}u^2`$ is then negligible. The simplified form of equation (3.10), namely

$$y(y1)^2\frac{\mathrm{d}F}{\mathrm{d}y}+(2y^21)F=0,y=(2\beta )^{1/2}u,$$

admits the normalized solution

$$F(u)=\frac{\rho \rho _c}{2E_1(1)\beta }\frac{y}{(1y)^3}\mathrm{exp}\left(\frac{1}{1y}\right)(0<y<1),$$

where $`E_1`$ is the first exponential integral $`(E_1(1)=0.219383934)`$. The scaling function $`F(u)`$ is therefore nonzero in the $`\beta \mathrm{}`$ limit only for $`y<1`$, i.e., $`u<(2\beta )^{1/2}`$. The upper bound $`y=1`$, i.e., $`u=(2\beta )^{1/2}`$, coincides with the point where the eigenfunction $`H_0(u)`$, as given by equation (3.23), peaks at low temperature.

## 4 Two-time observables: dynamical equations

In this section we will successively establish dynamical equations for the two-time correlation function $`C(t,s)`$, for its derivative $`C(t,s)/s`$, and for the response function $`R(t,s)`$. (See  for similar techniques.)

We consider the (connected) two-time correlation function $`C(t,s)`$ between the population of box number 1 at times $`s`$ (waiting time) and $`t`$ (observation time), with $`0st`$:

$$C(t,s)=N_1(t)N_1(s)N_1(t)N_1(s)=N_1(t)N_1(s)\rho ^2.$$

This definition can be recast as

$$C(t,s)=\underset{k=1}{\overset{\mathrm{}}{}}k\gamma _k(t,s)\rho ^2,$$

with

$$\gamma _k(t,s)=\underset{j=1}{\overset{\mathrm{}}{}}jf_j(s)𝒫\{N_1(t)=kN_1(s)=j\}.$$

The evolution of the $`\gamma _k(t,s)`$ with respect to $`t`$ is given, for $`ts`$, by a master equation similar to (2.3):

$`{\displaystyle \frac{\gamma _k(t,s)}{t}}=\gamma _{k+1}(t,s)+\sigma (t)r_{k1}\gamma _{k1}(t,s)\left(1+\sigma (t)r_k\right)\gamma _k(t,s)(k1),`$
$`{\displaystyle \frac{\gamma _0(t,s)}{t}}=\gamma _1(t,s)\sigma (t)r_0\gamma _0(t,s).`$ (4.1)
These equations preserve the sum rule $`_k\gamma _k(t,s)=\rho `$. At $`t=s`$, the initial conditions are $`\gamma _k(s,s)=kf_k(s)`$, implying

$$C(s,s)=\underset{k=1}{\overset{\mathrm{}}{}}k^2f_k(s)\rho ^2,$$

which is the variance of the population of box number 1 at time $`s`$, as it should.

In the discussion of the fluctuation-dissipation theorem, we will need expressions of the time derivative $`C(t,s)/s`$. We have

$$\frac{C(t,s)}{s}=\underset{k=1}{\overset{\mathrm{}}{}}k\phi _k(t,s),$$

with

$$\phi _k(t,s)=\frac{\gamma _k(t,s)}{s}.$$

The evolution of the $`\phi _k(t,s)`$ with respect to $`t`$ is again given, for $`ts`$, by equations similar to (2.3):

$`{\displaystyle \frac{\phi _k(t,s)}{t}}=\phi _{k+1}(t,s)+\sigma (t)r_{k1}\phi _{k1}(t,s)\left(1+\sigma (t)r_k\right)\phi _k(t,s)(k1),`$
$`{\displaystyle \frac{\phi _0(t,s)}{t}}=\phi _1(t,s)\sigma (t)r_0\phi _0(t,s),`$ (4.2)
with initial conditions at $`t=s`$

$`\phi _k(s,s)=f_{k+1}(s)+\sigma (s)r_{k1}f_{k1}(s)(k1),`$
$`\phi _0(s,s)=f_1(s),`$
so that

$$\left(\frac{C(t,s)}{s}\right)_{t=s}=\sigma (s)\underset{k=1}{\overset{\mathrm{}}{}}kr_kf_k(s)+2(1f_0(s))\rho .$$

The two-time response function $`R(t,s)`$ is a measure of the change in the mean population of box number 1 at time $`t`$, induced by an infinitesimal modulation of the conjugate variable, i.e., the local chemical potential acting on the same box, at the earlier time $`s`$.

In the presence of an arbitrary local, time-dependent chemical potential $`\mu (t)`$, the energy of box number 1 at time $`t`$ reads

$$E\left(N_1(t)\right)=\mathrm{ln}(N_1(t)+1)\mu (t)N_1(t).$$

The occupation probabilities of this box now depend on $`\mu (t)`$: we denote them by $`f_k^\mu (t)`$. In the thermodynamic limit, i.e., to leading order as $`M\mathrm{}`$, the occupation probabilities of all the other boxes $`(i=2,\mathrm{},M)`$ are still given by the $`f_k(t)`$.

The response function reads

$$R(t,s)=\left(\frac{\delta N_1(t)}{\delta \mu (s)}\right)_{\mu =0}=\underset{k=1}{\overset{\mathrm{}}{}}kh_k(t,s),$$

with

$$h_k(t,s)=\left(\frac{\delta f_k^\mu (t)}{\delta \mu (s)}\right)_{\mu =0}.$$

The modified occupation probabilities $`f_k^\mu (t)`$ obey the dynamical equations

$`{\displaystyle \frac{\mathrm{d}f_k^\mu (t)}{\mathrm{d}t}}=f_{k+1}^\mu (t)+\sigma (t)\mathrm{e}^{\beta \mu (t)}r_{k1}f_{k1}^\mu (t)\left(1+\sigma (t)\mathrm{e}^{\beta \mu (t)}r_k\right)f_k^\mu (t)(k1),`$
$`{\displaystyle \frac{\mathrm{d}f_0^\mu (t)}{\mathrm{d}t}}=f_1^\mu (t)\sigma (t)\mathrm{e}^{\beta \mu (t)}r_0f_0^\mu (t).`$ (4.3)
The initial values $`f_k^\mu (0)=f_k(0)`$ and the parameter $`\sigma (t)`$ are unchanged. The dynamical equations (4.3) preserve the sum rule $`_kf_k^\mu (t)=1`$.

Equations (4.3) imply that the $`h_k(t,s)`$ obey, for $`t>s`$,

$`{\displaystyle \frac{h_k(t,s)}{t}}=h_{k+1}(t,s)+\sigma (t)r_{k1}h_{k1}(t,s)\left(1+\sigma (t)r_k\right)h_k(t,s)(k1),`$
$`{\displaystyle \frac{h_0(t,s)}{t}}=h_1(t,s)\sigma (t)r_0h_0(t,s),`$ (4.4)
with initial conditions at $`t=s`$

$`h_k(s,s)=\beta \sigma (s)\left(r_{k1}f_{k1}(s)r_kf_k(s)\right)(k1),`$
$`h_0(s,s)=\beta \sigma (s)r_0f_0(s),`$
so that

$$R(s,s)=\beta \left(1f_0(s)\right).$$

The behavior of the two-time observables will now be successively investigated in the next three sections, first at criticality ($`\rho =\rho _c`$), both at equilibrium and in the nonequilibrium regime, and then in the condensed phase ($`\rho >\rho _c`$).

## 5 Equilibrium critical dynamics

When the waiting time $`s`$ becomes very large, keeping the difference $`\tau =ts0`$ fixed, two-time quantities reach their equilibrium values, which only depend on $`\tau `$, both in the fluid phase ($`\rho <\rho _c`$) and along the critical line ($`\rho =\rho _c`$). This section is devoted to the latter case.

The equilibrium correlation function $`C_{\mathrm{eq}}(\tau )`$ reads

$$C_{\mathrm{eq}}(\tau )=\underset{k=1}{\overset{\mathrm{}}{}}k\gamma _{k,\mathrm{eq}}(\tau )\rho _c^2,$$

where the equilibrium values $`\gamma _{k,\mathrm{eq}}(\tau )`$ of the $`\gamma _k(t,s)`$ obey

$`{\displaystyle \frac{\mathrm{d}\gamma _{k,\mathrm{eq}}(\tau )}{\mathrm{d}t}}=\gamma _{k+1,\mathrm{eq}}(\tau )+r_{k1}\gamma _{k1,\mathrm{eq}}(\tau )\left(1+r_k\right)\gamma _{k,\mathrm{eq}}(\tau )(k1),`$
$`{\displaystyle \frac{\mathrm{d}\gamma _{0,\mathrm{eq}}(\tau )}{\mathrm{d}t}}=\gamma _{1,\mathrm{eq}}(\tau )r_0\gamma _{0,\mathrm{eq}}(\tau ),`$ (5.1)
with initial conditions $`\gamma _{k,\mathrm{eq}}(0)=kf_{k,\mathrm{eq}}`$. This implies consistently that $`C_{\mathrm{eq}}(0)=\mu _c\rho _c^2`$ is the equilibrium variance of the population of a generic box.

As the time difference $`\tau `$ gets large, we have $`\gamma _{k,\mathrm{eq}}(\tau )\rho _cf_{k,\mathrm{eq}}`$, so that $`C_{\mathrm{eq}}(\tau )0`$. The decay of $`C_{\mathrm{eq}}(\tau )`$ for large $`\tau `$ can be investigated along the lines of section 3.1. The two regimes I and II are again to be considered separately, although results concerning the former will not be needed explicitly. In the scaling regime, we look for a similarity solution to equations (5.1) of the form

$$\gamma _{k,\mathrm{eq}}(\tau )f_{k,\mathrm{eq}}\tau ^{1/2}G_{\mathrm{eq}}(u),u=k\tau ^{1/2},$$
(5.2)
for which equations (5.1) yield the differential equation

$$\left(𝒟+\frac{1}{2}\right)G_{\mathrm{eq}}(u)=0,$$
(5.3)
with boundary conditions $`G_{\mathrm{eq}}(0)=0`$ and $`G_{\mathrm{eq}}(u)u`$ as $`u\mathrm{}`$. The differential operator $`𝒟`$ was defined in equation (3.7).

Equation (5.3) can be solved by the method of Appendix A. With the definition (A.1), we have the functional equation

$$\frac{M_{G_{\mathrm{eq}}}(z+2)}{M_{G_{\mathrm{eq}}}(z)}=\frac{z1}{2(z+2)(\beta 1z)},$$

whose suitably normalized solution reads

$$M_{G_{\mathrm{eq}}}(z)=\frac{\pi ^{1/2}\mathrm{\Gamma }\left(\frac{z1}{2}\right)\mathrm{\Gamma }\left(\frac{\beta +1z}{2}\right)}{2^{z+1}\mathrm{\Gamma }\left(\frac{z+2}{2}\right)\mathrm{\Gamma }\left(\frac{\beta }{2}\right)}(1<Rez<\beta +1).$$

The decay of the equilibrium correlation function at large $`\tau `$ is again dominated by the scaling regime, and the scaling form (5.2) yields

$$C_{\mathrm{eq}}(\tau )A_{\mathrm{eq}}\tau ^{(\beta 3)/2},$$
(5.4)
with

$$A_{\mathrm{eq}}=\frac{M_{G_{\mathrm{eq}}}(\beta 2)}{\zeta (\beta )}=\frac{\pi \mathrm{\Gamma }\left(\frac{\beta 3}{2}\right)}{2^\beta \mathrm{\Gamma }\left(\frac{\beta }{2}\right)^2\zeta (\beta )}.$$

As expected, the fluctuation-dissipation theorem holds at equilibrium. Indeed the equilibrium values $`\phi _{k,\mathrm{eq}}(\tau )`$ and $`h_{k,\mathrm{eq}}(\tau )`$ of $`\phi _k(t,s)`$ and $`h_k(t,s)`$ obey the same equations, identical to equations (5.1), with initial conditions:

$`h_{k,\mathrm{eq}}(0)=\beta \phi _{k,\mathrm{eq}}(0)=\beta \left(f_{k,\mathrm{eq}}f_{k+1,\mathrm{eq}}\right)(k1),`$
$`h_{0,\mathrm{eq}}(0)=\beta \phi _{0,\mathrm{eq}}(0)=\beta f_{1,\mathrm{eq}},`$
hence the identities

$$h_{k,\mathrm{eq}}(\tau )=\beta \phi _{k,\mathrm{eq}}(\tau )$$
(5.5)
and

$$R_{\mathrm{eq}}(\tau )=\beta \frac{\mathrm{d}C_{\mathrm{eq}}(\tau )}{\mathrm{d}\tau }.$$

The last formula is the fluctuation-dissipation theorem in its usual differential form.

## 6 Nonequilibrium critical dynamics

We now turn to the nonequilibrium critical behavior of the two-time correlation function $`C(t,s)`$, response function $`R(t,s)`$, and fluctuation-dissipation ratio

$$X(t,s)=\frac{R(t,s)}{\beta {\displaystyle \frac{C(t,s)}{s}}},$$
(6.1)
in the scaling regime where both time variables $`s`$ and $`t`$ are large and comparable. Hereafter $`x`$ will denote the dimensionless time ratio

$$x=\frac{t}{s}1.$$

We first analyze the scaling behavior of the correlation function $`C(t,s)`$. Looking for a two-variable scaling solution to equations (4.1) of the form

$$\gamma _k(t,s)f_{k,\mathrm{eq}}t^{1/2}G(u,x),u=kt^{1/2},x=t/s,$$
(6.2)
we obtain the partial differential equation

$$\left(x\frac{}{x}+𝒟+\frac{1}{2}\right)G(u,x)=0$$
(6.3)
for the scaling function $`G(u,x)`$, with initial condition

$$G(u,1)=uF(u),$$
(6.4)
and boundary condition $`G(\mathrm{},x)=0`$ for all $`x1`$. The function $`F(u)`$ is known from (3.8), and the operator $`𝒟`$ is given by (3.7).

As shown in Appendix B, equation (6.3) can be solved explicitly by the method of separation of variables. Equation (6.2) implies the scaling law

$$C(t,s)s^{(\beta 3)/2}\mathrm{\Phi }(x),$$
(6.5)
where

$$\mathrm{\Phi }(x)=\frac{x^{(\beta 3)/2}}{\zeta (\beta )}_0^{\mathrm{}}u^{1\beta }G(u,x)du.$$

Using (B.4), we obtain

$$\mathrm{\Phi }(x)=x^{\beta /2}\underset{n=0}{\overset{\mathrm{}}{}}A_nx^n,$$
(6.6)
with

$$A_n=\frac{a_nM_{G_n}(\beta 2)}{\zeta (\beta )}.$$

More explicitly, equations (B.3) and (B.5) imply that the leading coefficient of the expansion (6.6) reads

$$A_0=\frac{\pi ^{1/2}\mathrm{\hspace{0.17em}2}^{3\beta }\mathrm{\Gamma }\left(\frac{\beta +4}{2}\right)}{3(\beta +1)\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)^2\zeta (\beta )},$$
(6.7)
while the coefficient ratios are rational functions of $`\beta `$:

$$\frac{A_n}{A_0}=\frac{3}{2^n}\underset{j=0}{\overset{n1}{}}(2j+\beta )\underset{k=0}{\overset{n}{}}\frac{()^k}{(2k+3)k!(nk)!}\underset{\mathrm{}=0}{\overset{k1}{}}\frac{2\mathrm{}+\beta +4}{2\mathrm{}+\beta +3},$$
(6.8)
i.e.,

$$\frac{A_1}{A_0}=\frac{\beta (2\beta +3)}{10(\beta +3)},\frac{A_2}{A_0}=\frac{\beta (\beta +2)(8\beta ^2+52\beta +45)}{280(\beta +3)(\beta +5)},\text{etc.}$$

The expansion (6.6) is convergent over the whole physical domain $`(x>1)`$. For $`x1`$, i.e., $`\tau =tss`$, the equilibrium result (5.4) is recovered as

$$\mathrm{\Phi }(x)A_{\mathrm{eq}}(x1)^{(\beta 3)/2}.$$
(6.9)
These properties can be checked by noticing that the expression (B.1) for the Laguerre polynomials simplifies to

$$L_n^{(\beta +1)/2}(u^2/4)\left(\frac{2n^{1/2}}{u}\right)^{(\beta +1)/2}J_{(\beta +1)/2}\left(un^{1/2}\right),$$

where $`J`$ is the Bessel function, in the scaling regime where the order $`n`$ is large and $`u`$ is small. The subsequent integrals can be estimated for $`n`$ large, yielding

$$A_n\frac{\pi }{2^\beta \mathrm{\Gamma }\left(\frac{\beta }{2}\right)^2\zeta (\beta )}n^{(\beta 5)/2}.$$

This asymptotic expression ensures the convergence of the series (6.6) for all $`x>1`$, and establishes (6.9), including the prefactor.

We now turn to the derivative $`C(t,s)/s`$ of the correlation function. Looking for a two-variable scaling solution to equations (4.2) in the scaling regime, of the form

$$\phi _k(t,s)f_{k,\mathrm{eq}}t^{1/2}G^{(1)}(u,x),$$
(6.10)
we get the partial differential equation

$$\left(x\frac{}{x}+𝒟\frac{1}{2}\right)G^{(1)}(u,x)=0,$$
(6.11)
with initial condition

$$G^{(1)}(u,1)=\beta \frac{F(u)}{u}2F^{}(u).$$
(6.12)
Equation (6.11) can again be solved by the method of Appendix B. We thus obtain

$$\frac{C(t,s)}{s}s^{(\beta 1)/2}\mathrm{\Phi }^{(1)}(x),$$
(6.13)
where

$$\mathrm{\Phi }^{(1)}(x)=\frac{x^{(\beta 1)/2}}{\zeta (\beta )}_0^{\mathrm{}}u^{1\beta }G^{(1)}(u,x)du.$$
(6.14)
This expression will not be made more explicit, as equation (6.5) yields more directly

$$\mathrm{\Phi }^{(1)}(x)=\frac{3\beta }{2}\mathrm{\Phi }(x)x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x},$$
(6.15)
i.e.,

$$\mathrm{\Phi }^{(1)}(x)=x^{\beta /2}\underset{n=0}{\overset{\mathrm{}}{}}A_n^{(1)}x^n,$$
(6.16)
with

$$A_n^{(1)}=\left(n+\frac{3}{2}\right)A_n.$$
(6.17)
The scaling behavior of the response function $`R(t,s)`$ can be determined by the same approach. Looking for a two-variable scaling solution to equations (4.4) in the scaling regime, of the form

$$h_k(t,s)f_{k,\mathrm{eq}}t^{1/2}G^{(2)}(u,x),$$

equations (4.4) yield the partial differential equation

$$\left(x\frac{}{x}+𝒟\frac{1}{2}\right)G^{(2)}(u,x)=0,$$
(6.18)
with initial condition

$$G^{(2)}(u,1)=\beta ^2\frac{F(u)}{u}\beta F^{}(u).$$
(6.19)
We thus obtain

$$R(t,s)s^{(\beta 1)/2}\mathrm{\Phi }^{(2)}(x),$$
(6.20)
where

$$\mathrm{\Phi }^{(2)}(x)=\frac{x^{(\beta 1)/2}}{\zeta (\beta )}_0^{\mathrm{}}u^{1\beta }G^{(2)}(u,x)du,$$
(6.21)
i.e., explicitly

$$\mathrm{\Phi }^{(2)}(x)=x^{\beta /2}\underset{n=0}{\overset{\mathrm{}}{}}A_n^{(2)}x^n.$$
(6.22)
The leading coefficient of this expansion reads

$$A_0^{(2)}=\frac{\beta \pi ^{1/2}\mathrm{\hspace{0.17em}2}^{1\beta }\mathrm{\Gamma }\left(\frac{\beta +2}{2}\right)}{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)^2\zeta (\beta )}=\frac{3\beta (\beta +1)}{2(\beta +2)}A_0,$$

while the coefficient ratios are again rational functions of $`\beta `$:

$$\frac{A_n^{(2)}}{A_0^{(2)}}=\frac{1}{(\beta +1)2^n}\underset{j=0}{\overset{n1}{}}(2j+\beta )\underset{k=0}{\overset{n}{}}\frac{()^k(2k+\beta +1)}{(2k+1)k!(nk)!}\underset{\mathrm{}=0}{\overset{k1}{}}\frac{2\mathrm{}+\beta +2}{2\mathrm{}+\beta +3},$$

i.e.,

$$\frac{A_1^{(2)}}{A_0^{(2)}}=\frac{\beta (2\beta +1)}{6(\beta +1)},\frac{A_2^{(2)}}{A_0^{(2)}}=\frac{\beta (\beta +2)(8\beta ^2+28\beta +9)}{120(\beta +1)(\beta +3)},\text{etc.}$$

The scaling results (6.13) and (6.20) imply that the fluctuation-dissipation ratio $`X(t,s)`$, defined in equation (6.1), only depends on the time ratio $`x`$ in the nonequilibrium scaling regime. We have indeed

$$X(t,s)𝒳(x)=\frac{\mathrm{\Phi }^{(2)}(x)}{\beta \mathrm{\Phi }^{(1)}(x)}.$$
(6.23)
The scaling function $`𝒳(x)`$ turns out to be universal, whereas $`\mathrm{\Phi }(x)`$, $`\mathrm{\Phi }^{(1)}(x)`$, and $`\mathrm{\Phi }^{(2)}(x)`$, which respectively enter equations (6.5), (6.13), (6.20), are only universal up to a scale fixing. (A definition of universal quantities in the present context will be recalled in the Discussion.)

Equations (6.16) and (6.22) yield the expression

$$𝒳(x)=\frac{{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}A_n^{(2)}x^n}{\beta {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}A_n^{(1)}x^n}=\underset{n=0}{\overset{\mathrm{}}{}}\xi _nx^n,$$
(6.24)
which is again convergent for $`x>1`$. The limit fluctuation-dissipation ratio, $`X_{\mathrm{}}=𝒳(\mathrm{})=\xi _0`$, takes the simple value

$$X_{\mathrm{}}=\frac{\beta +1}{\beta +2}.$$
(6.25)
Equation (6.24) also yields

$$\xi _1=\frac{\beta ^2}{3(\beta +2)(\beta +3)},\xi _2=\frac{\beta ^2(2\beta ^3+2\beta ^2+78\beta +117)}{45(\beta +2)(\beta +3)^2(\beta +5)},\text{etc.}$$

The behavior of the fluctuation-dissipation ratio $`𝒳(x)`$ close to equilibrium (i.e., for $`x1`$) is also of interest. Equations (6.14) and (6.21) for $`x=1`$, together with (6.12) and (6.19), imply

$$\underset{x1}{lim}\left(\beta \mathrm{\Phi }^{(1)}(x)\mathrm{\Phi }^{(2)}(x)\right)=\frac{\beta }{\zeta (\beta )}M_F^{}(\beta 2)=\frac{2^{1\beta }\beta }{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\zeta (\beta )}.$$

This expression, together with equations (6.9) and (6.15), yields

$$𝒳(x)1\frac{2\mathrm{\Gamma }\left(\frac{\beta }{2}\right)^2}{\pi \mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\mathrm{\Gamma }\left(\frac{\beta 1}{2}\right)}(x1)^{(\beta 1)/2},$$

confirming that the fluctuation-dissipation theorem is restored for $`x1`$, and explicitly giving the leading violation of that theorem in the scaling regime, which is proportional to $`(\tau /s)^{(\beta 1)/2}`$.

## 7 Nonequilibrium dynamics in the condensed phase

We finally turn to the long-time behavior of the two-time quantities $`C(t,s)`$, $`R(t,s)`$, and $`X(t,s)`$ in the condensed phase ($`\rho >\rho _c`$).

In Regime II, we look for two-variable scaling solutions to equations (4.1), (4.2), and (4.4), inspired by (3.9), of the form

$$\gamma _k(t,s)\frac{k}{t}G(u,x),\phi _k(t,s)\frac{k}{t^2}G^{(1)}(u,x),h_k(t,s)\frac{k}{t^2}G^{(2)}(u,x),$$
(7.1)
where $`u=kt^{1/2}`$ and $`x=t/s`$.

The three scaling functions obey partial differential equations of the form (6.3), (6.11), or (6.18), where the differential operator $`𝒟`$ is replaced by

$$𝒟=\frac{\mathrm{d}^2}{\mathrm{d}u^2}+\left(\frac{u}{2}+A\frac{\beta +2}{u}\right)\frac{\mathrm{d}}{\mathrm{d}u}\frac{3}{2}+\frac{A}{u}.$$
(7.2)
This differential operator is related to the Heun operator $``$ (3.14) and to the function $`Y(u)`$ defined in equation (3.12) by the conjugation

$$Y(u)=𝒟Y(u).$$
(7.3)
Hence $`𝒟`$ and $``$ share the same eigenvalues $`E_n`$ ($`n=0,1,\mathrm{}`$), which are not known explicitly, as already mentioned, except $`E_0=0`$.

The two-time scaling functions defined in (7.1) are then determined in analogy with the critical case. We thus obtain

$$G(u,x)=Y(u)\underset{n=0}{\overset{\mathrm{}}{}}a_nx^{E_n1/2}H_n(u),G^{(1,2)}(u,x)=Y(u)\underset{n=0}{\overset{\mathrm{}}{}}a_n^{(1,2)}x^{E_n+1/2}H_n(u),$$
(7.4)
where the coefficients $`a_n`$, $`a_n^{(1,2)}`$ are determined by the initial conditions

$$G(u,1)=F(u),G^{(1)}(u,1)=\left(\frac{A}{u}\frac{\beta }{u^2}\right)F(u)\frac{2}{u}F^{}(u),G^{(2)}(u,1)=\frac{\beta }{u}F^{}(u).$$
(7.5)
Using the expansions (7.4), as well as the linear dependence of $`F(u)`$ in $`(\rho \rho _c)`$ (see equations (3.16), (3.17)), we obtain the following scaling forms for the correlation and response functions:

$$\begin{array}{cc}C(t,s)s^{1/2}(\rho \rho _c)\mathrm{\Phi }(x),& \mathrm{\Phi }(x)=\underset{n=0}{\overset{\mathrm{}}{}}A_nx^{E_n},\\ \frac{C(t,s)}{s}s^{1/2}(\rho \rho _c)\mathrm{\Phi }^{(1)}(x),& \mathrm{\Phi }^{(1)}(x)=\underset{n=0}{\overset{\mathrm{}}{}}A_n^{(1)}x^{E_n},\\ R(t,s)s^{1/2}(\rho \rho _c)\mathrm{\Phi }^{(2)}(x),& \mathrm{\Phi }^{(2)}(x)=\underset{n=0}{\overset{\mathrm{}}{}}A_n^{(2)}x^{E_n}.\end{array}$$
(7.6)
The scaling functions $`\mathrm{\Phi }(x)`$, $`\mathrm{\Phi }^{(1,2)}(x)`$ only depend on $`\beta `$. They take finite limit values, both at $`x=0`$ and at $`x=\mathrm{}`$. The fluctuation-dissipation ratio

$$X(t,s)𝒳(x)=\frac{\mathrm{\Phi }^{(2)}(x)}{\beta \mathrm{\Phi }^{(1)}(x)}$$

again only depends on the time ratio $`x`$ in the nonequilibrium scaling regime. The whole scaling function $`𝒳(x)`$ is universal, just as in the critical case. In particular, the limit fluctuation-dissipation ratio $`X_{\mathrm{}}=𝒳(\mathrm{})=A_0^{(2)}/A_0^{(1)}`$ in the condensed phase is universal, and only depends on $`\beta `$. Using equations (3.16) and (7.5), we can derive the following expression for $`X_{\mathrm{}}`$ in terms of the ground-state eigenfunction $`H_0(u)`$:

$$X_{\mathrm{}}=\frac{{\displaystyle _0^{\mathrm{}}}\left({\displaystyle \frac{u}{2}}A+{\displaystyle \frac{\beta }{u}}\right)H_0^2(u)du}{{\displaystyle _0^{\mathrm{}^{}}}uH_0^2(u)du}.$$
(7.7)
The above expressions cannot be made more explicit in general. More quantitative results can be derived at high temperature and at low temperature, using results of section 3.2.

Let us consider the limit fluctuation-dissipation ratio $`X_{\mathrm{}}`$, as given by equation (7.7). In the high-temperature case, by (3.18) and (3.20), we obtain

$$X_{\mathrm{}}=\frac{1}{2}(\pi 2)\beta +\mathrm{}(\beta 0).$$
(7.8)
In the low-temperature case, using (3.21) and (3.23), we obtain

$$X_{\mathrm{}}=\beta ^{1/2}\frac{\beta ^{3/4}}{4}+\mathrm{}(\beta \mathrm{}).$$
(7.9)
It is instructive to compare the limit fluctuation-dissipation ratios corresponding to the critical point and to the condensed phase. The value of $`X_{\mathrm{}}`$ at criticality is given by the analytical expression (6.25). The value of $`X_{\mathrm{}}`$ in the condensed phase is obtained by a numerical evaluation of (7.7): the ground-state eigenfunction $`H_0(u)`$ is obtained by numerically solving the differential equation (3.10) or (3.13), and then used to evaluate the integrals entering the result (7.7). The values thus obtained smoothly interpolate between the limiting laws (7.8) and (7.9). In Figure 1, thick full lines show the physical values of the fluctuation-dissipation ratios, while thin dashed lines show their continuation to high temperatures. Both fluctuation-dissipation ratios start from $`1/2`$ at infinite temperature, and converge to the limit values $`X_{\mathrm{}}=0`$ and $`X_{\mathrm{}}=1`$ at zero temperature, respectively in the condensed phase and at criticality (see Discussion).

## 8 Discussion

At the onset of this paper we gave a comparative presentation of two classes of dynamical urn models, the Ehrenfest class and the Monkey class. All these models have simple static properties, because their Hamiltonian is a sum of contributions of independent boxes. They possess however interesting nonequilibrium dynamical properties, even in the mean-field geometry. The backgammon model  is a prototypical example of the Ehrenfest class, while model B of Reference , and the zeta urn model  which is the subject of the present work, are examples of the other class.

Let us come back to the role of statistics in the definition and properties of such models. The essential difference between the Ehrenfest class and the Monkey class indeed resides in matters related to a priori statistics. Statistics enters the dynamical definition of models: the proposed moves for the Ehrenfest class (respectively, the Monkey class) are chosen according to ball-box statistics (respectively, box-box statistics). Consistently, statistics also enters their equilibrium definition: the partition functions (1.2), (1.3) of the original Ehrenfest model, and of its generalization to $`M`$ urns, involve a factorial of the total number of particles in their denominators. Inverse factorials $`1/N_i!`$, taking into account equivalent labelings of particles within each box, are also involved in the evaluation of the partition function (1.6).

In statistical mechanics with Maxwell-Boltzmann statistics, the presence of inverse factorials has its origin in the indiscernibility of identical classical objects. These factorials are absent for the Monkey class, in which the populations $`N_i`$ are involved in flat sums, such as (1.7), just as occupation numbers are in quantum-mechanical statistical mechanics with Bose-Einstein statistics. There is of course nothing quantum-mechanical in the urn models considered here. It can be said for short, following Reference , that the equilibrium statistics of the Ehrenfest class is Maxwell-Boltzmann, while that of the Monkey class is Bose-Einstein. Table 1 illustrates this discussion.

The main focus of this work concerns the nonequilibrium properties of the zeta urn model, with emphasis on the aging behavior of the two-time correlation and response functions of the fluctuating population of a given box, and of the corresponding fluctuation-dissipation ratio. We considered successively the critical line ($`\rho =\rho _c`$) and the condensed phase ($`\rho >\rho _c`$). We summarize these two cases below.

The critical line in the temperature-density plane corresponds to a line of fixed points, parametrized by inverse temperature $`\beta `$ . In other words, critical exponents, both static and dynamical, and more generally universal quantities, depend continuously on temperature. In the present context, a universal quantity ought to be independent:

– of the initial state (provided it is homogeneous and disordered, i.e., the $`f_k(0)`$ decay rapidly),

– of the specific form of the energy of each box (provided it diverges logarithmically, as $`E(N_i)\mathrm{ln}N_i`$, for large occupation numbers),

– and of details of the dynamics (such as Metropolis versus heat-bath).

Our results (6.5) and (6.20) for the two-time correlation and response function have the expected product form (see  and references therein), involving: a common non-universal prefactor, a negative power of the waiting time, related to the anomalous dimension of the observable, and a universal scaling function of the time ratio, or temporal aspect ratio, $`x=t/s`$. Accordingly, the fluctuation-dissipation ratio (6.23), $`X(t,s)𝒳(x)`$, only depends on $`x`$ in the nonequilibrium scaling regime. The limit value $`X_{\mathrm{}}=𝒳(\mathrm{})`$ has been recently emphasized  to be a new universal quantity, characteristic of nonequilibrium critical dynamics. For the zeta urn model, the result (6.25) applies to the regular part of the critical line ($`\beta >3`$), so that $`4/5<X_{\mathrm{}}<1`$, as shown by the upper curve in Figure 1. This range is unusual for a critical system. Indeed, statistical-mechanical models such as ferromagnets are observed to have $`0<X_{\mathrm{}}1/2`$ at their critical point. The upper bound $`X_{\mathrm{}}=1/2`$, corresponding to the mean-field situation , is also observed in a range of simpler models . It is worth noticing that the backgammon model also has a high fluctuation-dissipation ratio at its zero-temperature critical point, namely $`X(t,s)1C/(\mathrm{ln}s)^2`$ for $`st`$, where the amplitude $`C`$ depends both on the observable and on the dynamical rule .

The present analysis of condensation dynamics for $`\rho >\rho _c`$ extends and completes that begun in . We have investigated various quantities related to the occupation probabilities in the scaling regime, describing the growing condensate. An approximate analysis of the Heun operator (3.14) has allowed us to obtain asymptotic expressions for various quantities at high and low temperature. The fluctuation-dissipation ratio admits a non-trivial limit value $`X_{\mathrm{}}`$ throughout the condensed phase, shown by the lower curve in Figure 1. Expression (7.9) shows that $`X_{\mathrm{}}\beta ^{1/2}`$ slowly goes to zero at low temperature. This behavior is very different from that of conventional models. It is indeed currently accepted  that coarsening systems, such as ferromagnets quenched from a high-temperature initial state, have identically $`X_{\mathrm{}}=0`$ throughout their low-temperature phase, i.e., for any temperature below $`T_c`$. These unusual features of the zeta urn model are less of a surprise if one remembers that the condensation dynamics of the present model is basically different from a domain-growth or coarsening dynamics. In the latter case, phase separation takes place in a statistically homogeneous way, at least for an infinite system. To the contrary, in the present situation, condensation takes place in a very inhomogeneous fashion. The form of equation (3.9) indeed demonstrates that the condensate of particles is shared by an ever decreasing fraction of boxes, scaling as $`t^{1/2}`$, each of them having a population growing as $`t^{1/2}`$ (until size effects eventually become important for $`tM^2`$).

## Appendix A Solving equation (3.6) by Mellin transformation

The Mellin transformation provides an efficient alternative way of solving the differential equation (3.6).

We define the Mellin transform $`M_f(z)`$ of a function $`f(u)`$ by

$$M_f(z)=_0^{\mathrm{}}u^{z1}f(u)du,f(u)=\frac{\mathrm{d}z}{2\pi \mathrm{i}}u^zM_f(z).$$
(A.1)
Equation (3.6), together with the boundary condition $`F(0)=1`$, is equivalent to the following functional equation

$$\frac{M_{1F}(z+2)}{M_{1F}(z)}=\frac{z}{2(z+2)(\beta 1z)}$$

for the Mellin transform of $`1F(u)`$. The solution of this functional equation, with boundary condition $`F(\mathrm{})=0`$, reads

$$M_{1F}(z)=\frac{\mathrm{\Gamma }\left(\frac{\beta +1z}{2}\right)}{z\mathrm{\hspace{0.17em}2}^z\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)}(0<Rez<\beta +1).$$
(A.2)
Similarly, for $`Rez<0`$, the Mellin transform of $`F(u)`$,

$$M_F(z)=\frac{\mathrm{\Gamma }\left(\frac{\beta +1z}{2}\right)}{z\mathrm{\hspace{0.17em}2}^z\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)}(Rez<0),$$
(A.3)
is the formal opposite of equation (A.2).

We also give for further reference the expression of the Mellin transform of the derivative $`F^{}(u)`$:

$$M_F^{}(z)=\underset{1<Rez<\beta }{\underset{}{(z+1)M_{1F}(z+1)}}=\underset{Rez<1}{\underset{}{(z+1)M_F(z+1)}}=\frac{\mathrm{\Gamma }\left(\frac{\beta z}{2}\right)}{2^{z+1}\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)}(Rez<\beta ).$$
(A.4)
## Appendix B Solving equation (6.3) by separation of variables and spectral superposition

Partial differential equations such as (6.3), with given initial and boundary conditions, can be solved explicitly by the method of separation of variables and spectral superposition.

To do so, we need a basis of eigenfunctions of the differential operator $`𝒟`$ of equation (3.7). Inspired by the explicit form of the solution (3.8), we set

$$G(u)=u^{\beta +1}\mathrm{e}^{u^2/4}L(v),v=u^2/4.$$

The eigenvalue equation $`(𝒟E)G(u)=0`$ becomes

$$v\frac{\mathrm{d}^2L}{\mathrm{d}v^2}+\left(\frac{\beta +3}{2}v\right)\frac{\mathrm{d}L}{\mathrm{d}v}+(E1)L=0,$$

which is the differential equation obeyed by the Laguerre polynomials $`L_n^\alpha (v)`$ , with $`\alpha =(\beta +1)/2`$ and $`n=E1`$.

The eigenvalues of the operator $`𝒟`$ therefore read $`E_n=n+1`$, with $`n=0,1,\mathrm{}`$ The associated eigenfunctions,

$$G_n(u)=u^{\beta +1}\mathrm{e}^{u^2/4}L_n^{(\beta +1)/2}(u^2/4),$$

with

$$L_n^{(\beta +1)/2}(v)=\underset{k=0}{\overset{n}{}}\frac{\mathrm{\Gamma }\left(n+\frac{\beta +3}{2}\right)}{\mathrm{\Gamma }\left(k+\frac{\beta +3}{2}\right)}\frac{(v)^k}{k!(nk)!},$$
(B.1)
obey the orthogonality property

$$_0^{\mathrm{}}G_m(u)G_n(u)u^\beta \mathrm{e}^{u^2/4}du=N_n\delta _{m,n},N_n=\frac{2^{\beta +2}\mathrm{\Gamma }\left(n+\frac{\beta +3}{2}\right)}{n!}.$$
(B.2)
The Mellin transform $`M_{G_n}(z)`$ of these eigenfunctions can again be evaluated in closed form:

$$M_{G_n}(z)=\frac{2^{\beta z}\mathrm{\Gamma }\left(n+\frac{z+2}{2}\right)\mathrm{\Gamma }\left(\frac{\beta +1z}{2}\right)}{n!\mathrm{\Gamma }\left(\frac{z+2}{2}\right)}(Rez<\beta +1),$$
(B.3)
where the normalization has been fixed by the condition  $`L_n^{(\beta +1)/2}(0)=1`$.

The method of spectral superposition consists in looking for the solution $`G(u,x)`$ of equation (6.3) as a linear superposition of the form

$$G(u,x)=\underset{n=0}{\overset{\mathrm{}}{}}a_nx^{(n+3/2)}G_n(u).$$
(B.4)
The coefficients $`a_n`$ are determined by the initial condition (6.4). Using the orthogonality property (B.2), we obtain

$$a_n=\frac{1}{N_n}_0^{\mathrm{}}u^2F(u)L_n^{(\beta +1)/2}(u^2/4)du.$$

Then, using (B.1) and (B.3), we are left with

$$a_n=\frac{2^{1\beta }}{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)}n!\underset{k=0}{\overset{n}{}}\frac{()^k}{(2k+3)k!(nk)!}\frac{\mathrm{\Gamma }\left(k+\frac{\beta +4}{2}\right)}{\mathrm{\Gamma }\left(k+\frac{\beta +3}{2}\right)}.$$
(B.5)