text:
# Condensation in the zero range process: stationary and dynamical properties

## 1 Introduction

Low dimensional stochastic particle systems far from equilibrium have a much richer structure than their equilibrium counterparts. In particular, even in one spatial dimension there is the possibility of a phase transition. On the other hand, we do not have available general criteria which would allow us to determine, for example, the phase diagram. A recent attempt in this direction is a proposal by Kafri et al. , who study phase separation in one dimension. Roughly speaking, they map the domain wall dynamics of clusters to a zero range process, for which analytical tools are available. Under suitable conditions, the clusters tend to grow. This coarsening is an intriguing phenomenon already on the level of the zero range process itself, and we will investigate it in much greater detail than available so far.

The zero range process is a stochastic particle system on the lattice $`^d`$ where the jump rate $`g(k)`$ of a given particle depends only on the occupation number $`k`$ at its current position. This model was originally introduced as a simple example of an interacting Markov process . Various properties have been established, among them the existence of the dynamics under very general conditions, classification of invariant measures, and hydrodynamic limits . If $`g(k)`$ is decreasing in $`k`$, then this induces an effective attraction between particles, as first noted in . As a result, there is a critical background density and excess particles condense on a non-extensive fraction of the volume.

Such a clustering phenomenon can be studied on two distinct levels. Firstly, the phenomenon is present already in the steady state. Based on some results for large deviations of independent, identically distributed random variables  and under general assumptions on $`g(k)`$, we will prove that for a typical steady state configuration there is a background phase at some critical density. Any additional mass is concentrated at a single, randomly located site. For the model introduced in , we will analyze the statistical properties of the background phase in detail. Secondly, there is the dynamics of clustering with the steady state appearing only in the long time limit. In a very recent preprint Godrèche  addresses this problem. He assumes a uniform initial density and investigates numerically how the probability distribution of the number of particles at some given site evolves in time, with particular attention to the macroscopic component of that distribution. In view of our static result we study the dynamics of Evans’ model , pursuing a somewhat different approach: During the initial nucleation process random sites are selected, at which a macroscopic number of particles accumulates. We investigate the effective dynamics of the number of particles at such cluster sites, in particular how the smaller occupation numbers become extinguished to the benefit of the larger ones.

Our contribution should be understood as a case study in the context of phase transitions in one-dimensional systems far from equilibrium, which has been a topic of major interest in the past decade (see  and references therein). Of particular interest is the occurrence of phase separation in systems with two conservation laws (see e.g.  and ) whose macroscopic behavior has only recently been examined to some extent . Given the correspondence to the zero-range process our results also provide new information on the stability of domain walls (shocks) which separate macroscopic regions of different phases in two-component systems. Domain wall stability already proved to be a key ingredient in the theory of boundary-induced phase transitions in systems with one conservation law  and thus may shed light on boundary-induced spontaneous symmetry breaking  in two-component systems.

To give a brief outline: in the following Section we discuss the stationary measures for the zero range process and in Section 3 the two main results on the equivalence of ensembles and the structure of the condensed phase are established. In Sections 4 and 5 we study the condensation transition for the model introduced in , first on the level of the stationary measure, and secondly through the dynamics of the clustering.

## 2 Zero range process and its infinite volume stationary measures

For notational simplicity we restrict ourselves to a description in one space dimension, but our results on the equivalence of ensembles hold for arbitrary dimension. Rather than defining the zero range process directly on an infinite lattice (cf. ), we first consider a finite system, compute the (unique) stationary measure and analyze it in the limit of infinite system size.

Let us consider a zero range process on the one-dimensional lattice $`\mathrm{\Lambda }_L=\{1,\mathrm{},L\}`$ of $`L`$ sites with periodic boundary conditions. Let $`\eta _x`$ be the number of particles on site $`x\mathrm{\Lambda }_L`$, where $`=\{0,1,2,\mathrm{}\}`$. The state space is given by $`\mathrm{\Omega }_L=^{\mathrm{\Lambda }_L}`$ and we denote a particle configuration by $`𝜼=(\eta _x)_{x\mathrm{\Lambda }_L}\mathrm{\Omega }_L`$. At a given site $`x\mathrm{\Lambda }_L`$, the number $`\eta _x`$ of particles decreases by one after an exponential waiting time with rate $`g(\eta _x)`$ and the leaving particle jumps to site $`x+y`$ with probability $`p(y)`$. The jump probabilities $`p:[0,1]`$ are normalized, $`_yp(y)=1`$, $`p(0)=0`$, assumed to be of finite range, $`p(y)=0`$ for $`|y|>R`$, and irreducible, $`p(1)>0`$. For the dynamics to be well defined in the limit $`L\mathrm{}`$ and to be nondegenerate the rate function $`g:[0,\mathrm{})`$ has to satisfy

$`\underset{k}{sup}|g(k+1)g(k)|<\mathrm{},g(k)>g(0)=0\text{for all }k>0.`$ (2.1)
The generator of the zero range process is then given by

$`(f)(𝜼)={\displaystyle \underset{x=1}{\overset{L}{}}}{\displaystyle \underset{y=R}{\overset{R}{}}}g(\eta _x)p(y)\left(f(𝜼^{x,x+y})f(𝜼)\right),`$ (2.2)
regarded as a linear operator on $`C(\mathrm{\Omega }_L,)`$. We used the shorthand $`\eta _z^{x,x+y}=\eta _z\delta (z,x)+\delta (z,((x+y1)modL)+1)`$ for all $`z\mathrm{\Lambda }_L`$, where $`\delta (.,.)`$ denotes the Kronecker delta function. The process conserves the number of particles $`\mathrm{\Sigma }_L(𝜼)=_{x=1}^L\eta _x`$, thereby dividing the configuration space into the finite, invariant subsets $`\mathrm{\Omega }_{L,N}=\{𝜼\mathrm{\Omega }_L|\mathrm{\Sigma }_L(𝜼)=N\}`$ with $`N`$. Restricted to $`C(\mathrm{\Omega }_{L,N},)`$ with $`L,N`$ fixed, $``$ is a finite dimensional matrix and the process is well defined. However for $`L=\mathrm{}`$ this is true only for “reasonable” initial conditions and under the assumption (2.1), see .

The following results for the stationary measures are well known and taken from . The zero range process (2.2) on $`\mathrm{\Omega }_{N,L}`$ is an irreducible Markov jump process with the unique stationary measure

$`\mu ^{N,L}(𝜼)={\displaystyle \frac{1}{Z(N,L)}}{\displaystyle \underset{x=1}{\overset{L}{}}}W(\eta _x)\delta (\mathrm{\Sigma }_L(𝜼),N).`$ (2.3)
The weight $`W`$ is given by

$`W(k):={\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{1}{g(i)}}`$ (2.4)
and the normalizing partition function is

$`Z(N,L)={\displaystyle \underset{𝜼\mathrm{\Omega }_L}{}}{\displaystyle \underset{x=1}{\overset{L}{}}}W(\eta _x)\delta (\mathrm{\Sigma }_L(𝜼),N).`$ (2.5)
Clearly (2.3) resembles a canonical ensemble in statistical mechanics. Therefore, in the limit of large system size $`L,N\mathrm{}`$ with fixed particle density $`\rho =N/L`$, (2.3) is expected to be equivalent to a grand canonical product measure, which is defined through

$`\nu _\varphi ^L(𝜼)={\displaystyle \underset{x=1}{\overset{L}{}}}\nu _\varphi (\eta _x)\text{with }\nu _\varphi (k)={\displaystyle \frac{1}{Z(\varphi )}}W(k)\varphi ^k,`$ (2.6)
and where the fugacity $`\varphi 0`$ is adjusted to fix the average density.

Let $`\varphi _c`$ be the radius of convergence of the grand canonical (one site) partition function

$`Z(\varphi )={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}W(k)\varphi ^k.`$ (2.7)
The measure (2.6) is well defined for fugacities $`\varphi [0,\varphi _c)`$ and its average particle density $`\rho (\varphi )`$ as a function of $`\varphi `$ is given by

$`\rho (\varphi )={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}k\nu _\varphi (k)=\varphi {\displaystyle \frac{\mathrm{log}Z(\varphi )}{\varphi }}.`$ (2.8)
The range of $`\rho `$ is the interval $`[0,\rho _c)`$, with $`\rho (0)=0`$ and $`\rho _c=lim_{\varphi \varphi _c}\rho (\varphi )`$ the critical density. $`\varphi \rho (\varphi )`$ is strictly increasing and we denote the inverse function on $`[0,\rho _c)`$ by $`\varphi (\rho )`$. If $`\varphi _c=\mathrm{}`$, then $`\rho _c=\mathrm{}`$ (see , Lemma 2.3.3), whereas for $`\varphi _c<\mathrm{}`$, both $`\rho _c=\mathrm{}`$ and $`\rho _c<\mathrm{}`$ are possible. In the second case $`Z(\varphi _c)<\mathrm{}`$ (see , Lemma 2.3.3) and $`\nu _{\varphi _c}`$ is a well defined probability measure with $`\eta _x_{\nu _{\varphi _c}}=\rho _c`$. Thus we set

$`\varphi (\rho )=\{\begin{array}{cc}\text{inverse of }\rho (\varphi )\hfill & ,\text{for }\rho <\rho _c\hfill \\ \varphi _c\hfill & ,\text{for }\rho \rho _c\hfill \end{array},`$ (2.11)
where $`\rho _c`$ can be either finite or infinite, cf. Figure 1 in Section 4. The reason for this particular convention will become clear in the next section. In this way we may regard the measure $`\nu _\varphi `$ also as a function of $`\rho `$ through $`\nu _{\varphi (\rho )}`$.

The link between canonical and grand canonical measures is given through the pointwise limit of the $`n`$-point marginal,

$`\underset{L\mathrm{}}{lim}\mu _n^{[\rho L],L}(𝐤)={\displaystyle \underset{i=1}{\overset{n}{}}}\nu _{\varphi (\rho )}(k_i).`$ (2.12)
Here $`[a]`$ denotes the integer part of $`a`$ and for every $`n`$, $`𝐱=(x_1,\mathrm{},x_n)\mathrm{\Lambda }_L^n`$ with $`x_ix_j`$ for $`ij`$, $`𝐤=(k_1,\mathrm{},k_n)^n`$ the $`n`$-point marginal is defined as $`\mu _{n,𝐱}^{N,L}(𝐤):=\mu ^{N,L}(\{\eta _{x_1}=k_1,\mathrm{},\eta _{x_n}=k_n\})`$. Since the measure $`\mu ^{N,L}`$ is permutation invariant (see (2.3)), the marginals do not depend on the sites $`𝐱`$ individually, but only on their number $`n`$. For this reason only $`n`$ is specified in our notation.

A rigorous result on this equivalence is available in , Appendix 2, but it does not cover the supercritical case $`\rho _c<\mathrm{}`$ and $`N/L\rho _c`$. It will be discussed in Section 3, where particular attention is given to the statistical properties of the excess density. If $`\rho _c<\mathrm{}`$, then $`\nu _{\varphi _c}`$ is well defined and by necessity decays subexponentially, which will be explained in (3.14). Examples with power-law decay are given in ,  (Example 2.3.4). In  Evans studies the dependence on the power-law exponent in more detail and introduces a generic model where $`\rho _c`$ can be either finite or infinite, depending on a system parameter. Its stationary and dynamical properties will be studied in detail in Sections 4 and 5.

The jump probabilities $`p`$ do not influence the stationary measures, but they play an important role for the relaxation dynamics, see Section 5. The zero range process is reversible if and only if $`p`$ is symmetric. The stationary current is given by $`j=m(p)g(\eta _x)`$, where $`m(p)=_{y=R}^Ryp(y)`$ is the first moment of $`p`$, which for non-symmetric jump probabilities is generically non-zero. The average jump rate $`g`$ of the two ensembles is given by

$`g(\eta _x)_{\mu ^{N,L}}={\displaystyle \frac{Z(N1,L)}{Z(N,L)}}`$ $`\text{(canonical)},`$
$`g(\eta _x)_{\nu _\varphi ^L}={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}g(k)\nu _\varphi (k)=\varphi `$ $`\text{(grand-canonical)}.`$ (2.13)
Thus for $`m(p)>0`$ the grand canonical current $`j(\rho )=m(p)\varphi (\rho )`$ is monotone increasing in $`\rho `$, approaching its maximum value $`m(p)\varphi _c`$ as $`\rho \rho _c`$, and correspondingly for $`m(p)<0`$.

## 3 Equivalence of supercritical measures

We consider the zero range process with $`\rho _c<\mathrm{}`$ and supercritical canonical measures with $`\rho =N/L>\rho _c`$. The heuristic picture, developed in  and , is that most sites of the system are distributed according to $`\nu _{\varphi _c}`$ with mean occupation number $`\rho _c`$. For large $`L`$, the $`(\rho \rho _c)L`$ excess particles presumably condense on a few sites. If so, locally one will observe the grand-canonical ensemble with $`\varphi =\varphi _c`$. This picture is made precise in

Theorem 1. (Equivalence of ensembles)
Let $`\rho \varphi (\rho )`$ be defined as in (2.8), (2.11). Then for every $`\rho [0,\mathrm{})`$ the $`n`$-th marginal has the pointwise limit

$`\underset{L\mathrm{}}{lim}\mu _n^{[\rho L],L}(𝐤)={\displaystyle \underset{i=1}{\overset{n}{}}}\nu _{\varphi (\rho )}(k_i).`$ (3.1)
The canonical partition functions converge as

$`\underset{L\mathrm{}}{lim}{\displaystyle \frac{1}{L}}\mathrm{log}Z([\rho L],L)=\mathrm{log}Z(\varphi (\rho ))\rho \mathrm{log}\varphi (\rho ).`$ (3.2)
Proof. It is convenient to characterize the distance between the canonical and grand-canonical measures through the relative entropy $`S`$. For two arbitrary probability measures $`\mu ,\nu `$ on a countable set $`\mathrm{\Omega }`$ it is defined as

$`S(\mu |\nu )=\{\begin{array}{cc}\underset{\omega \mathrm{\Omega }}{}\mu (\omega )\mathrm{log}\frac{\mu (\omega )}{\nu (\omega )}\hfill & ,\text{if }\mu \nu \hfill \\ \mathrm{}\hfill & ,\text{otherwise}\hfill \end{array}`$ (3.5)
where $`\mu \nu `$ means that $`\mu `$ is absolutely continuous w.r.t. $`\nu `$. It has the properties (see e.g. )

$`S(\mu |\nu )0,S(\mu |\nu )=0\mu (\omega )=\nu (\omega )\text{ for all }\omega \mathrm{\Omega }.`$ (3.6)
To prove (3.1) it is therefore enough to establish that

$`\underset{L\mathrm{}}{lim}S\left(\mu _n^{[\rho L],L}|\nu _{\varphi (\rho )}^n\right)=0`$ (3.7)
for every $`\rho [0,\mathrm{})`$.

We recall that $`\mathrm{\Sigma }_L(𝜼)=_{x=1}^L\eta _x`$ denotes the number of particles. By definition (2.3) one has for a fixed $`𝜼\mathrm{\Omega }_L`$

$`\mu ^{N,L}(𝜼)={\displaystyle \frac{W^L(𝜼)\delta (\mathrm{\Sigma }_L(𝜼),N)}{Z(N,L)}},`$ (3.8)
where $`W^L`$ denotes the product measure $`W^L(𝜼)=_{x=1}^LW(\eta _x)`$. Thus, using (3.5), it follows that for every $`\varphi [0,\varphi _c]`$

$`S\left(\mu ^{N,L}|\nu _\varphi ^L\right)={\displaystyle \underset{𝜼\mathrm{\Omega }_{L,N}}{}}\mu ^{N,L}(𝜼)\mathrm{log}{\displaystyle \frac{W^L(𝜼)}{\nu _\varphi ^L(𝜼)Z(N,L)}},`$ (3.9)
since $`\mu ^{N,L}`$ is absolutely continuous w.r.t. $`\nu _\varphi ^L`$. From (2.5) and (2.6) we conclude that

$`\nu _\varphi ^L(𝜼)Z(\varphi )^L`$ $`=`$ $`W^L(𝜼)\varphi ^N\text{for }𝜼\mathrm{\Omega }_{L,N}\text{and}`$
$`\nu _\varphi ^L\left(\{\mathrm{\Sigma }_L=N\}\right)Z(\varphi )^L`$ $`=`$ $`Z(N,L)\varphi ^N.`$ (3.10)
To simplify notation we use the shorthand $`\nu _\varphi ^L(A)=_{𝜼A}\nu _\varphi ^L(𝜼)`$ for a subset $`A`$ of the configuration space and $`\{𝜼|\mathrm{\Sigma }_L(𝜼)=N\}=\{\mathrm{\Sigma }_L=N\}`$. Inserting in (3.9) this yields

$`S\left(\mu ^{N,L}|\nu _\varphi ^L\right)=\mathrm{log}\nu _\varphi ^L\left(\{\mathrm{\Sigma }_L=N\}\right).`$ (3.11)
At this point we use the subadditivity of $`S`$, namely if two measures $`\mu `$, $`\nu `$ have marginals $`\mu _i`$, $`\nu _i`$, $`i=1,2`$, then

$`S(\mu |\nu )S(\mu _1|\nu _1)+S(\mu _2|\nu _2).`$ (3.12)
Therefore for every $`n\{1,\mathrm{},L\}`$ and $`\varphi [0,\varphi _c]`$

$`S\left(\mu _n^{N,L}|\nu _\varphi ^n\right){\displaystyle \frac{1}{[L/n]}}\mathrm{log}\nu _\varphi ^L\left(\{\mathrm{\Sigma }_L=N\}\right).`$ (3.13)
The key point is to maximize $`\nu _\varphi ^L\left(\{\mathrm{\Sigma }_L=N\}\right)`$ by appropriately adjusting $`\varphi =\varphi (\rho )`$ as defined in (2.11). In the subcritical case, $`\rho =N/L<\rho _c`$, we have $`\varphi (\rho )<\varphi _c`$ and $`\nu _{\varphi (\rho )}`$ has exponential moments (see (2.6)). Then the variance $`\sigma ^2`$ of $`\nu _\varphi `$ is finite and the limit distribution of $`(\mathrm{\Sigma }_L(\rho L))/(\sigma \sqrt{L})`$ is given by the normal distribution $`𝒩(0,1)`$ (cf. (4.11)). By the local limit theorem (see e.g. ) we get in this case $`\nu _{\varphi (\rho )}^L\left(\{\mathrm{\Sigma }_L=[\rho L]\}\right)1/\sqrt{L}`$ for large $`L`$. For $`\rho =\rho _c`$ the decay of $`\nu _{\varphi _c}`$ is subexponential, since $`\varphi =\varphi _c`$ is the radius of convergence of the partition function $`Z(\varphi )`$ in (2.7). Thus

$`\underset{k\mathrm{}}{lim}{\displaystyle \frac{1}{k}}\mathrm{log}\nu _{\varphi _c}(k)=0`$ (3.14)
and the second moment of $`\nu _{\varphi _c}`$ could be infinite, leading to a non-normal limit distribution (cf. (4.12)). Since the first moment of $`\nu _{\varphi _c}`$ equals $`\rho _c<\mathrm{}`$, by the local limit theorem for non-normal limit distributions (see also ) we get the lower bound $`\nu _{\varphi _c}^L\left(\{\mathrm{\Sigma }_L=[\rho _cL]\}\right)1/L`$. The supercritical case $`\rho >\rho _c`$, where $`\varphi (\rho )=\varphi _c`$, can be reduced to the critical one via

$`\nu _{\varphi _c}^L\left(\left\{\mathrm{\Sigma }_L=[\rho L]\right\}\right)`$ (3.15)
$`\nu _{\varphi _c}^L\left(\left\{\eta _L=[\rho L][\rho _c(L1)],\mathrm{\Sigma }_{L1}=[\rho _c(L1)]\right\}\right)=`$
$`=\nu _{\varphi _c}\left([\rho L][\rho _c(L1)]\right)\nu _{\varphi _c}^{L1}\left(\left\{\mathrm{\Sigma }_{L1}=[\rho _c(L1)]\right\}\right).`$
Both terms decay subexponentially, the first one using (3.14) and the second one as in the critical case. Thus in all cases we have a subexponential lower bound on $`\nu _{\varphi (\rho )}^L\left(\{\mathrm{\Sigma }_L=[\rho L]\}\right)`$ and the limit (3.7) follows for all $`\rho [0,\mathrm{})`$ from (3.13),

$`\underset{L\mathrm{}}{lim}S\left(\mu _n^{[\rho L],L}|\nu _{\varphi (\rho )}^n\right)\underset{L\mathrm{}}{lim}{\displaystyle \frac{1}{[L/n]}}\mathrm{log}\nu _{\varphi (\rho )}^L\left(\{\mathrm{\Sigma }_L=[\rho L]\}\right)=0.`$ (3.16)
To establish (3.2) we use the second line of (3) and immediately get

$`{\displaystyle \frac{1}{L}}\mathrm{log}Z(N,L){\displaystyle \frac{1}{L}}\mathrm{log}\nu _\varphi ^L\left(\{\mathrm{\Sigma }_L=N\}\right)=\mathrm{log}Z(\varphi ){\displaystyle \frac{N}{L}}\mathrm{log}\varphi ,`$ (3.17)
for all $`N,L`$. With $`\varphi =\varphi (N/L)`$ we can use the above estimates, so that the second term on the left vanishes in the limit $`L\mathrm{}`$ and the assertion (3.2) follows. $`\mathrm{}`$
Theorem 1 ensures us that the volume fraction of the condensed phase vanishes in the limit $`L\mathrm{}`$. In principle it could still contain an infinite number of sites, and the question remains, how many condensed sites there are in a typical configuration. The answer depends on the large-$`k`$ behavior of the critical distribution $`\nu _{\varphi _c}(k)`$. From (3.14) we know already that it decays subexponentially. We will show that for a large class of such subexponential distributions the excess particles condense on a single, randomly located site.

Let $`𝒞_2`$ be the set of distributions $`\nu _{\varphi _c}`$ which have finite second moment and for which the integrated tail $`\nu _{\varphi _c}\left(\{\eta _1k\}\right)`$ is heavier than $`\mathrm{exp}[k^\alpha ]`$ for some $`\alpha (0,1/2)`$. $`𝒞_p`$ denotes the set of distributions with power-law tail $`\nu _{\varphi _c}(k)k^b`$, $`b(2,3]`$, for which the second moment diverges. For a detailed description and results on subexponential distributions we refer the reader to .

Theorem 2. Let $`\nu _{\varphi _c}`$ be in the class $`𝒞_2𝒞_p`$ as defined above, with first moment $`\rho _c<\mathrm{}`$. Then for the sequence of the corresponding canonical measures $`\mu ^{[\rho L],L}`$ with $`\rho >\rho _c`$ one has

$`\underset{L\mathrm{}}{lim}\mu ^{[\rho L],L}\left(\left\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\right\}\right)=1.`$ (3.18)
The probability that there is a site which contains at least $`(\rho \rho _c)L`$ particles converges to one in the thermodynamic limit. But the occupation number cannot be substantially larger than $`(\rho \rho _c)L`$, since almost all sites are distributed according to $`\nu _{\varphi _c}`$, as proved in Theorem 1. Therefore all excess particles are condensed on a single site in the limit $`L\mathrm{}`$. The proof of Theorem 2 uses large deviation results on the asymptotic behavior of $`\nu _{\varphi _c}^L\left(\{\mathrm{\Sigma }_L\rho L\}\right)`$ for $`L\mathrm{}`$, which we summarize in the following Lemma for our purpose.

Lemma. Let $`\eta _1,\eta _2,\mathrm{}`$ be i.i.d. random variables with mean zero and probability distribution $`P𝒞_2𝒞_p`$. Then with $`\mathrm{\Sigma }_L=_{x=1}^L\eta _x`$ one has for any $`\rho >0`$

$`\underset{L\mathrm{}}{lim}\left|{\displaystyle \frac{P^L\left(\{\mathrm{\Sigma }_L\rho L\}\right)}{LP\left(\{\eta _1\rho L\}\right)}}1\right|=0.`$ (3.19)
Proof of the Lemma. In the case $`P𝒞_p`$ see  (Chapter 1, Corollary 1.1.1 to 1.1.3) and for distributions in $`𝒞_2`$ see .

The interpretation of the Lemma is that under a distribution in $`𝒞_2𝒞_p`$ the rare event $`\{\mathrm{\Sigma }_L\rho L\}`$ in the limit $`L\mathrm{}`$ is realized by the deviation of a single (randomly positioned) site with probability one. Thus for $`\nu _{\varphi _c}𝒞_2𝒞_p`$ this is also a typical configuration under the canonical measure $`\mu ^{[\rho L],L}`$, since for $`\rho >\rho _c`$ the latter is basically given by the grand canonical critical measure $`\nu _{\varphi _c}`$ under the condition $`\{\mathrm{\Sigma }_L=[\rho L]\}`$, cf. (3.21). This argument will be made precise in the following.

Proof of Theorem 2. Obviously one has

$`\mu ^{[\rho L],L}\left(\left\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\right\}\right)1`$ (3.20)
for all $`L`$. To find a lower bound that converges to $`1`$ for $`L\mathrm{}`$ we apply the results of the Lemma. We note that for all $`𝜼\mathrm{\Omega }_{L,N}`$ it is

$`\mu ^{[\rho L],L}(\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\})={\displaystyle \frac{\nu _{\varphi _c}^L(\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\})}{\nu _{\varphi _c}^L\left(\{\mathrm{\Sigma }_L=[\rho L]\}\right)}}`$ (3.21)
cf. (3.8), and for $`L\mathrm{}`$

$`\nu _{\varphi _c}^L(\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\})=L\nu _{\varphi _c}(\{\eta _1[(\rho \rho _c)L]\})(1+o(1)).`$ (3.22)
Shifting the expectation value to zero via $`\eta _x^{}=\eta _x\rho _c`$, $`\mathrm{\Sigma }_L^{}:=_{x=1}^L\eta _x^{}`$ and using $`\nu _{\varphi _c}^L\left(\{\mathrm{\Sigma }_L=[\rho L]\}\right)\nu _{\varphi _c}^L\left(\{\mathrm{\Sigma }_L[\rho L]\}\right)`$ we get from (3.21) and (3.22)

$`\mu ^{[\rho L],L}\left(\left\{\underset{1xL}{\mathrm{max}}\eta _x[(\rho \rho _c)L]\right\}\right)`$ (3.23)
$`{\displaystyle \frac{L\nu _{\varphi _c}\left(\left\{\eta _1^{}[(\rho \rho _c)L]\right\}\right)}{\nu _{\varphi _c}^L\left(\left\{\mathrm{\Sigma }_L^{}[(\rho \rho _c)L]\right\}\right)}}\left(1+o(1)\right).`$
By (3.19) the righthand side of (3.23) converges to $`1`$ in the limit $`L\mathrm{}`$ and (3.18) is shown. $`\mathrm{}`$

## 4 Stationary properties near criticality

In  Evans studies the zero range process with rates

$`g_b(k)=\theta (k)(1+b/k),`$ (4.1)
where $`\theta (0)=0`$ and $`\theta (k)=1`$ for $`k>0`$. He observes that for $`b>2`$ the critical density $`\rho _c<\mathrm{}`$. Our goal here is to study the properties of the invariant measures for any $`b>0`$. The stationary weight for the zero range process with rates $`g_b`$ is given by

$`W(k)={\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{1}{1+b/i}}={\displaystyle \frac{k!}{(1+b)_k}}={\displaystyle \frac{k!\mathrm{\Gamma }(1+b)}{\mathrm{\Gamma }(1+b+k)}},`$ (4.2)
where $`(a)_k=_{i=0}^{k1}(a+i)`$ denotes the Pochhammer symbol, $`a`$, $`k`$. The grand canonical partition function of (2.7) is

$`Z(\varphi )=_2F_1(1,1;1+b;\varphi ):={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)_k(1)_k}{(1+b)_k}}{\displaystyle \frac{\varphi ^k}{k!}}.`$ (4.3)
Its radius of convergence is $`\varphi _c=1`$, and $`{}_{2}{}^{}F_{1}^{}`$ denotes the hypergeometric function , which has the expansion

$`{}_{2}{}^{}F_{1}^{}(k,k;k+b;\varphi )={\displaystyle \frac{\mathrm{\Gamma }(k+b)\mathrm{\Gamma }(kb)}{\mathrm{\Gamma }(k)^2}}(1\varphi )^{bk}\left[1+𝒪(1\varphi )\right]+`$ (4.4)
$`{\displaystyle \frac{\mathrm{\Gamma }(k+b)\mathrm{\Gamma }(bk)}{\mathrm{\Gamma }(b)^2}}\left[1+{\displaystyle \frac{k^2}{1+kb}}(1\varphi )+𝒪(1\varphi )^2\right].`$
The particle density (2.8) is given by

$`\rho (\varphi )={\displaystyle \frac{\varphi _2F_1(2,2;2+b;\varphi )}{(1+b)_2F_1(1,1;1+b;\varphi )}}\text{for }\varphi <1.`$ (4.5)
In the following we analyze the grand-canonical single site measure $`\nu _\varphi `$ of (2.6) in the limit $`\varphi 1`$, i.e. near the critical density $`\rho _c`$. For $`\rho _c<\mathrm{}`$ the limit $`\nu _1`$ is well defined, as discussed in Section 2, and according to Theorem 1 it is the distribution of the non-condensed phase for supercritical systems with $`N/L=\rho >\rho _c`$. As long as $`\varphi <1`$ the distribution $`\nu _\varphi `$ has exponential moments. For $`\varphi =1`$ the exponential tail of $`\nu _\varphi `$ disappears and the tail becomes proportional to the weight $`W(k)`$ as defined in (4.2). Using Stirling’s formula, the behavior of the weight for large $`k`$ is given by the power-law $`W(k)\mathrm{\Gamma }(1+b)k^b`$. These distributions have moments up to order $`b1`$. Thus different scenarios are encountered as $`b`$ is varied.
The case $`0<b1`$:
For $`b<1`$ the leading order in the asymptotic expansion for $`Z`$ and $`\rho `$ is given by

$`Z(\varphi )`$ $``$ $`\mathrm{\Gamma }(1+b)\mathrm{\Gamma }(1b)(1\varphi )^{b1}\mathrm{},`$
$`\rho (\varphi )`$ $``$ $`{\displaystyle \frac{\varphi }{(1+b)^2(1b)}}(1\varphi )^1\rho _c=\mathrm{},`$ (4.6)
as $`\varphi 1`$. For every density, the stationary distribution in the limit $`L\mathrm{}`$ is given by the grand-canonical measure $`\nu _\varphi `$ (see Figure 1(a)). The probability to have a fixed number of particles on a given site vanishes with increasing density, as is shown for the example of an empty site in Figure 1(b). Thus in the limit there is an infinite number of particles on every site with probability one, as it should be for homogeneous systems with $`\rho \mathrm{}`$.

For $`b=1`$ this picture does not change qualitatively, except for the logarithmic corrections

$`Z(\varphi )`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}(1\varphi )}{\varphi }}\mathrm{},`$
$`\rho (\varphi )`$ $`=`$ $`{\displaystyle \frac{\varphi }{(\varphi 1)\mathrm{log}(1\varphi )}}1\rho _c=\mathrm{},`$ (4.7)
as $`\varphi 1`$.

The case $`1<b2`$:
For $`1<b<2`$ the leading order terms change, and

$`Z(\varphi )`$ $``$ $`{\displaystyle \frac{b}{b1}}+\mathrm{\Gamma }(1+b)\mathrm{\Gamma }(1b)(1\varphi )^{b1}Z(1)={\displaystyle \frac{b}{b1}},`$
$`\rho (\varphi )`$ $``$ $`\varphi (b1)\mathrm{\Gamma }(b)\mathrm{\Gamma }(2b)(1\varphi )^{b2}\rho _c=\mathrm{},`$ (4.8)
as $`\varphi 1`$. As before, for $`b=2`$ the first order terms have logarithmic corrections but the qualitative behavior does not change. In particular, $`\rho _c=\mathrm{}`$ and the stationary distribution is described by the grand-canonical ensemble for every density $`\rho `$ (see Figure 1(a)).

However, somewhat surprisingly, the character of this distribution for large $`\rho `$ differs from the case $`b1`$. Since $`Z(1)<\mathrm{}`$, $`\nu _1`$ is well defined and there is a non-zero probability to have a fixed number of particles at a given site,

$`\nu _1(0)`$ $`=`$ $`{\displaystyle \frac{1}{Z(1)}}={\displaystyle \frac{b1}{b}},`$
$`\nu _1(k)`$ $`=`$ $`{\displaystyle \frac{W(k)}{Z(1)}}\mathrm{\Gamma }(b)(b1)k^b\text{for large }k.`$ (4.9)
For example the probability of an empty site, given by $`\nu _{\varphi (\rho )}(0)=1/Z(\varphi (\rho ))`$, decreases monotonically with $`\varphi `$, i.e. with increasing density $`\rho `$. In contrast to the case $`b1`$, it does not vanish in the limit $`\rho \mathrm{}`$, however, it reaches the non-zero value $`\nu _1(0)=(b1)/b`$ (see Figure 1(b)). So no matter how large the density, the fraction of empty sites in a typical configuration is always greater than $`(b1)/b`$.

Distributions with power-law tails are well studied (see e.g.  and references therein). A typical configuration for this stationary distribution, i.e. a set of $`L`$ i.i.d. random variables $`\eta _x`$ drawn from $`\nu _1`$, is known to have a hierarchical structure. The $`n`$-th largest value of the set $`\{\eta _1,\mathrm{},\eta _L\}`$ scales as $`(\mathrm{\Gamma }(b1)L/n)^{1/(b1)}`$, which holds for every $`b>1`$. In our particular case $`1<b2`$, this means that the particle number $`\mathrm{\Sigma }_L`$ also scales as $`L^{1/(b1)}`$ and thus grows faster than the number of summands $`L`$. Therefore the particle density $`\mathrm{\Sigma }_L/L`$ diverges as $`L^{(2b)/(b1)}`$ and the highest occupied site contains a nonzero fraction of the particles in the system. This hierarchical structure of typical configurations can be understood as a precursor for the condensation phenomenon to be discussed in the next part.
The case $`b>2`$:
In this case $`b`$ is large enough so that besides the normalization also the first moment of the grand canonical distribution converges in the limit $`\varphi 1`$:

$`Z(\varphi )`$ $``$ $`{\displaystyle \frac{b}{b1}}{\displaystyle \frac{b}{(b1)(b2)}}(1\varphi )Z(1)={\displaystyle \frac{b}{b1}},`$
$`\rho (\varphi )`$ $``$ $`{\displaystyle \frac{1}{b2}}+\varphi (b1)\mathrm{\Gamma }(b)\mathrm{\Gamma }(2b)(1\varphi )^{b2}\rho _c={\displaystyle \frac{1}{b2}},`$ (4.10)
as $`\varphi 1`$. We note that for $`b>3`$ also the second moment $`\sigma ^2`$ of the distribution $`\nu _1`$ exists and the number of particles satisfies the usual central limit theorem

$`\underset{L\mathrm{}}{lim}\nu _1^L\left(\xi _1{\displaystyle \frac{\mathrm{\Sigma }_L\rho _cL}{\sigma \sqrt{L}}}\xi _2\right)={\displaystyle _{\xi _1}^{\xi _2}}G(\xi )𝑑\xi ,`$ (4.11)
where $`G`$ denotes the Gaussian probability density with zero mean and unit variance. The density for $`b>3`$ is of order $`\rho (\varphi )=1/(b2)+𝒪(1\varphi )`$ and its first derivative is finite at $`\varphi =1`$ and given by $`\rho ^{}(1)=(b1)^2/((b3)^2(b2)^2)`$ (see Figure 1(a)).

As explained already in Section 3, the most occupied site contains of order $`L^{1/(b1)}`$ particles, and for $`b<3`$ this fluctuation is larger than $`\sqrt{L}`$. Therefore the scaling limit leads to a self-similar distribution, which is given by the completely asymmetric Lévy distribution $`L_{(b1)}`$ (for details see  or ),

$`\underset{L\mathrm{}}{lim}\nu _1^L\left(\xi _1{\displaystyle \frac{\mathrm{\Sigma }_L\rho _cL}{\left(b\mathrm{\Gamma }(b1)L\right)^{1/(b1)}}}\xi _2\right)={\displaystyle _{\xi _1}^{\xi _2}}L_{(b1)}(\xi )𝑑\xi .`$ (4.12)
With (4) we have $`\rho ^{}(1)=\mathrm{}`$ for $`b<3`$, leading to a differentiable function $`\varphi (\rho )`$, as shown in Figure 1(a).

## 5 Dynamics of the condensation

The stationary distribution investigated so far carries no information on the kinetics of the condensation. A natural set-up is to start with particles uniformly distributed at the supercritical density $`\rho >\rho _c`$. In the beginning the excess particles condense at a few random sites. Such a site containing many excess particles is called a cluster site. Thus there are several clusters which are essentially immobile as will be discussed below. On the remaining sites, called bulk sites, the distribution relaxes to $`\nu _1`$. With increasing time the larger clusters will gain particles at the expense of the smaller ones, causing some of the clusters to disappear. Eventually only a single cluster containing all excess particles survives, which is typical for the stationary distribution, as has been discussed already in Section 3.

In the following we will study the kinetics of condensation and its dependence on system parameters in detail. Compared to Section 4 the jump probabilities play an important role for the dynamical properties of the system. We will focus on nearest-neighbor jumps which are either totally asymmetric, i.e. particles only jump to the right with $`p_a(y)=\delta (1,y)`$, or symmetric, i.e. $`p_s(y)=\left(\delta (1,y)+\delta (1,y)\right)/2`$. In contrast to previous sections we are not able to rigorously prove our statements, but use heuristic considerations which are corroborated through comparison with simulation data.

### 5.1 Cluster formation

To distinguish cluster and bulk sites, we define a site $`x\{1,\mathrm{},L\}`$ to be a cluster if it contains a macroscopic fraction of the excess particles $`\eta _x>\alpha (\rho \rho _c)L`$. The prefactor $`\alpha (0,1]`$ is rather arbitrary but for simulations it is important that the clusters are well separated from the bulk fluctuations, which are of order $`(\mathrm{\Gamma }(b1)L)^{1/(b1)}`$. Since fluctuations grow only sublinearly with $`L`$ this separation is clearly guaranteed in the limit $`L\mathrm{}`$, and for finite systems it holds for sufficiently large $`L`$ depending on $`\alpha `$ and $`\rho `$. In our simulations we choose $`\alpha =1/40`$, requiring system sizes of about 200 sites minimum for the values of $`\rho `$ and $`b`$ considered. Let $`n(t)`$ be the number of cluster sites at time $`t`$ and $`m_i(t)`$, $`i=1,\mathrm{},n(t)`$ be the size of the $`i`$-th largest cluster, i.e. $`m_1(t)\mathrm{}m_{n(t)}(t)`$. These quantities depend also on the system parameters $`b`$, $`L`$, and $`\rho `$.

By definition a typical cluster has a size of order $`\alpha (\rho \rho _c)L`$ and so there are of the order of $`1/\alpha `$ cluster sites. The time scale for the formation of such clusters is very roughly estimated as follows: $`𝒪(L)`$ particles have to move a distance of order $`L`$ to form the cluster. So in the asymmetric case the time scale for cluster formation is $`𝒪(L^2)`$. The dependence on $`\rho `$ and $`b`$ is not so obvious, since the bulk has not yet relaxed to $`\nu _1`$ and the speed of particles still changes. In the symmetric case the time scale is $`𝒪(L^3)`$, since the particles diffuse without a drift.

In Figure 2(a), in the totally asymmetric case, we plot the average number of particles in the condensed phase, normalized by the number of excess particles,

$`f(t)={\displaystyle \frac{\underset{i=1}{\overset{n(t)}{}}m_i(t)}{(\rho \rho _c)L}}.`$ (5.1)
The time axis is scaled proportional to $`\tau _a=(\rho \rho _c)^2L^2/b`$ (cf. Section 5.2), since this choice gives the best data collapse when $`\rho `$ and $`b`$ are varied. On this time scale most excess particles become trapped in a cluster and the bulk relaxes to $`\nu _1`$. In Figure 2(b) the average number of clusters is plotted as a function of time. The number of clusters grows for a short time and then starts decreasing again.

### 5.2 Coarsening

Once the bulk has relaxed to $`\nu _1`$, each bulk site loses particles at the average rate $`g_b_{\nu _1}=1`$. In the asymmetric case this results in a particle current $`j=\varphi =1`$ (see (2) and Figure 1). On top of that, excess particles are exchanged between clusters. The bulk can be seen as a homogeneous medium where the excess particles move, and the cluster sites as boundaries where they enter and exit. A cluster of size $`m>0`$ loses excess particles with rate $`g(m)1=b/m`$ and gains particles from neighboring clusters. Since this rate decreases with increasing cluster size, smaller clusters lose particles to the larger ones.

To quantitatively describe this coarsening process we study the normalized mean cluster size $`\overline{m}(t)=f(t)/n(t)`$ as a function of time for large system sizes $`L`$. The ensemble average of this quantity (denoted by $`\mathrm{}`$) is expected to grow according to a scaling law

$`\overline{m}(t)t^\beta ,`$ (5.2)
with a scaling exponent $`\beta `$ . To estimate this exponent we notice that the time scale for the coarsening process is determined by two factors: Firstly, the rate at which excess particles leave a cluster of size $`m`$ and enter the bulk. Secondly, the typical time of such a particle to reach the neighboring cluster. The mobility of excess particles in the bulk is characterized by the average exit rate from an occupied site,

$`g(k)|k>0_{\nu _1}={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}g(k){\displaystyle \frac{\nu _1(k)}{\nu _1(k0)}}={\displaystyle \frac{g(k)_{\nu _1}}{1\nu _1(0)}}=(1{\displaystyle \frac{b1}{b}})^1=b,`$ (5.3)
as follows from (2) and (4). The excess particles perform a random walk in the bulk, which is either biased or unbiased, depending on the first moment of the jump probabilities. The time it takes to reach the next cluster at a distance of order $`𝒪(L/n)=𝒪(m/(\rho \rho _c))`$ is given by the mean first passage time for a random walk . In the following we distinguish between asymmetric and symmetric jump probabilities, to estimate the time scale for a particle to leave a cluster and the one for reaching the neighboring cluster.

In the totally asymmetric case excess particles leave a cluster of size $`m`$ with rate $`b/m`$ and move towards the right neighboring cluster without returning, so the time to lose one particle scales like $`𝒪(m/b)`$. Since the mean first passage time of a biased random walk is proportional to the distance, this particle spends a time $`𝒪(m/(b(\rho \rho _c)))`$ in the bulk, because $`b`$ is the speed of the particle. Thus the typical times for exiting a cluster and entering the next neighbor are of the same order in $`L`$ and there are $`𝒪(1/(\rho \rho _c))`$ excess particles in the bulk. The coarsening time scale, determined by the typical time for a cluster to lose all $`m`$ particles, is thus proportional to

$`t_a(m):={\displaystyle \frac{m^2}{b}},\text{ i.e.}\beta _a=1/2.`$ (5.4)
So for the asymmetric case the mean cluster size is predicted to grow like $`\overline{m}(t)(bt)^{1/2}`$.

In the case of symmetric jump probabilities $`p_s`$, excess particles perform an unbiased random walk in the bulk with diffusion constant $`b`$. Thus the mean first passage time to reach the next cluster is proportional to the square of the distance, i.e. $`𝒪(m^2/(b(\rho \rho _c)^2))`$. In contrast to the asymmetric case, it is very likely that particles return to the cluster they left. The probability that they do not return but reach the neighboring cluster, which is the relevant event for coarsening, is inverse proportional to the diffusion distance, i.e. $`𝒪((\rho \rho _c)/m)`$. So the typical time of a particle to leave a cluster is $`𝒪(m^2/(b(\rho \rho _c)))`$ and as before there are $`𝒪(1/(\rho \rho _c))`$ excess particles in the bulk. The coarsening time scale is therefore proportional to

$`t_s(m):=t_a(m){\displaystyle \frac{m}{(\rho \rho _c)}}={\displaystyle \frac{m^3}{b(\rho \rho _c)}},\text{ i.e.}\beta _s=1/3.`$ (5.5)
Thus the mean cluster size is predicted to grow like $`\overline{m}(t)((\rho \rho _c)bt)^{1/3}`$.

In general, the time scale for the coarsening regime is determined by the largest clusters with size $`𝒪((\rho \rho _c)L)`$ and is thus of order

$`\tau _a`$ $`:=`$ $`t_a\left((\rho \rho _c)L\right)=(\rho \rho _c)^2L^2/b`$
$`\tau _s`$ $`:=`$ $`t_s\left((\rho \rho _c)L\right)=(\rho \rho _c)^2L^3/b`$ (5.6)
for asymmetric resp. symmetric jump probabilities. The growth exponents for the two cases are confirmed by simulations and shown in Figure 3(a) for the totally asymmetric and 3(b) for the symmetric jump probabilities. Using the time scale $`\tau _a`$ resp. $`\tau _s`$ and normalizing $`\overline{m}`$ by the number of excess particles $`(\rho \rho _c)L`$ the data for different system sizes collapse. The measured growth exponents from these data are

$`\beta _a=0.514\pm 0.005\text{and}\beta _s=0.334\pm 0.004,`$ (5.7)
which agree with the above predictions. Independently from us, these exponents have been obtained in  by numerical simulations.

Note that the clusters coarsen on the same time scale as they nucleate. However, looking at the number of time steps in Figures 2 and 3, the time it takes to nucleate clusters is by a factor of 10 shorter than the coarsening regime.

### 5.3 Saturation

Eventually all clusters except for two will have disappeared and finite size effects become dominant. The scaling law (5.2) is no longer valid in this regime, since the mean cluster size $`\overline{m}`$ saturates towards its limiting value. The two clusters exchange particles until one of them vanishes and the system has reached its stationary state where all excess particles are concentrated at a single cluster site. In Figure 4 we plot the average size of the three largest clusters $`m_i(t)`$, $`i=1,2,3`$ normalized by $`(\rho \rho _c)L`$. Note that the coarsening regime ends at latest when the third largest cluster has disappeared, and thus takes only about a tenth of the total equilibration time. In the following, we focus on the totally asymmetric jump probabilities, but the symmetric choice would lead to an effective evolution equation of the same form.

Let $`M=m_1+m_2`$ be the total number of particles at the two largest cluster sites. On the time scale $`t^{}=t/((\rho \rho _c)L/b)`$ the two clusters exchange single particles with effective rates $`(\rho \rho _c)L/m_i`$, $`i=1,2`$ (see discussion in Section 5.2). The fluctuations of $`M`$ on this time scale are only $`𝒪(1)`$. Thus it is $`M=(\rho \rho _c)L+𝒪(1)`$, since the bulk is relaxed to $`\nu _1`$ and all other clusters have disappeared. Let $`q(m,t^{})`$ be the probability of having $`m=0,\mathrm{},M`$ particles at one cluster site and $`Mm`$ at the other one. The dynamics is then governed by the effective master equation

$`{\displaystyle \frac{}{t^{}}}q(m,t^{})=q(m,t^{})\left[{\displaystyle \frac{\theta (m)}{m/M}}+{\displaystyle \frac{\theta (Mm)}{1m/M}}\right]+`$ (5.8)
$`q(m1,t^{}){\displaystyle \frac{\theta (m)}{1(m1)/M}}+q(m+1,t^{}){\displaystyle \frac{\theta (Mm)}{(m+1)/M}}.`$
The exchange rates on the right hand side only depend on the rescaled variable $`m/M`$ and not on the system parameters $`\rho `$, $`L`$ and $`b`$. For large $`M`$, $`m/M`$ varies on the time scale $`t^{}/Mt/((\rho \rho _c)^2L^2/b)`$, confirming that $`\tau _a`$ (5.2) is the appropriate time scale for the two-cluster situation. Therefore the plots in Figure 4 are independent of the system parameters. However, in the following discussion we stick to the time scale $`t^{}`$ and the discrete variable $`m`$.

For any initial condition the solution of (5.8) tends to the inverse binomial distribution $`q^{}(m)1/\left(\genfrac{}{}{0pt}{}{M}{m}\right)`$. It is symmetric around $`m=M/2`$ and for small $`m`$ we have $`q^{}(m)=𝒪(M^m)`$. Thus the two extreme occupation numbers $`m=0`$ resp. $`m=M`$ are the most probable ones and in the limit $`L\mathrm{}`$ both have probability $`1/2`$, consistent with the results of Section 3. For $`m=\alpha M`$, $`\alpha (0,1)`$ it is $`q^{}(m)=𝒪\left(\sqrt{M}\left(\alpha ^\alpha (1\alpha )^{1\alpha }\right)^M\right)`$ using Stirling’s formula. Thus, in the stationary state, the typical time for a macroscopic fluctuation of the cluster size diverges exponentially with the system size $`L`$.

To study the relaxation dynamics, we write (5.8) in the canonical form, using the discrete derivative $`_mf(m):=f(m+1)f(m)`$,

$`{\displaystyle \frac{}{t^{}}}q(m,t^{})`$ $`=`$ $`_m\left(a(m)q(m,t^{})\right)+_m^2\left(d(m)q(m,t^{})\right)`$
$`a(m)`$ $`=`$ $`{\displaystyle \frac{1}{(1m/M)}}{\displaystyle \frac{M}{m}}={\displaystyle \frac{2m/M1}{m/M(1m/M)}}`$
$`d(m)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{(1m/M)}}+{\displaystyle \frac{M}{m}}\right)={\displaystyle \frac{1}{2m/M(1m/M)}}.`$ (5.9)
For ease of notation, we take $`0<m<M`$ and ignore the boundary terms. Note that (5.3) is symmetric around $`m=M/2`$. It describes diffusive motion in a double well potential with drift $`a(m)`$ and diffusion coefficient $`d(m)`$, with the slightly unusual feature that the minima of the potential are located close to the boundaries at $`m=1`$ and $`m=M1`$.

The master equation must be supplied with a suitable initial condition $`q(m,0)`$, which, since resulting from a complex coarsening process, is not readily available. A crude estimate can be found by noting that $`q(m,0)`$ is roughly proportional to the lifetime of the occupation number $`m`$. It is determined by the inverse exit rate for $`m`$ taken from Equation (5.8), and thus we expect $`q(m,0)6m/M(1m/M)`$. This is a symmetric single hump distribution with mean $`M/2`$ and standard deviation $`M/(2\sqrt{5})0.22M`$. Comparing with the simulation data at the time when the third largest cluster has just disappeared, $`m_3(t)/M<0.01`$, we indeed find a single hump distribution with mean $`M/2`$ and standard deviation $`0.166M`$. When solving Equation (5.3) with this initial distribution, the expectation of the larger cluster size, given by $`1/2+|m/M1/2|_{q(m,t)}`$, is indistinguishable from $`m_1(t)`$ in Figure 4. Note that, except for the time scales, Figure 4(a) and 4(b) are almost identical, confirming that the effective master equation for the symmetric case is of the same form as (5.8).

## Acknowledgments

S.G. acknowledges the support of the Graduiertenkolleg “Mathematik im Bereich ihrer Wechselwirkung mit der Physik” and the support of theDAAD/CAPES program “PROBRAL”. He is also grateful for a fruitful research visit at the “Forschungszentrum Jülich”, where this work was initiated.


format: text

text:
# Using bond-length dependent transferable force constants to predict vibrational entropies in Au-Cu, Au-Pd, and Cu-Pd alloys

## I Introduction

In the last decade, a clear prescription has emerged to obtain phase diagrams of materials from first-principles. When two compounds mix, some amount of site disorder occurs. First-principles alloy theory has been mostly preoccupied with defining appropriate energy models for such partially disordered systems, and obtaining the parameters for them from increasingly accurate Density Function Theory (DFT) methods. In particular, the cluster expansion approach 1 has been highly successful, as it allows one to parameterize the energy of systems with disorder on a fixed framework of sites. Many binary phase diagrams have been calculated in this way 2-6.

Most of these first-principles phase diagram calculations do not include vibrational entropy effects. The assumption that vibrational entropy differences between phases are small, and can therefore be neglected, has been recently examined experimentally 7-11 and theoretically 12-18. In Al2Cu 18 vibrational entropy contributions were shown to be essential in reversing the stability of $`\theta `$ and $`\theta `$’ phases at finite temperatures. In the Al-Sc system, vibrational entropy contributions were shown to increase the solubility limits 27-fold 16,17.

For alloys with partial site disorder, vibrational entropy in first-principles phase diagrams is computationally intensive to determine. In conventional first-principles alloy theory, a cluster expansion has to be fit to the ground state energies of a large number of A-B ordered states (for a binary alloy). Vibrational effects can be formally included by fitting to vibrational free energies, rather than ground state energies 19. This requires that one determines the phonon spectrum for many configurations. Unfortunately, this is computationally expensive.

One idea that has been proposed to alleviate this problem is the use of force constants that can be transferred between different environments. The existence of transferable force constants would greatly reduce the number of calculations required to determine vibrational entropies from first-principles: currently, almost all of the work in calculating vibrational entropies is involved with the numerical determination of force constants for every environment in a large number of structures.

The idea of using transferable force constants to reduce computational expense is tremendously appealing. Despite this, until recently, there was little clear evidence supporting the validity of their use. For instance, experimental data in Fe3Al shows that force constants can have a clear configurational dependence 20. Strictly speaking, it is not possible to rigorously define configuration-independent force constant matrices as the form of the force constant matrix depends on the symmetry of the configuration 21. Nevertheless, configuration-independent transferable force constants have been used in computations to study segregation in Ni3Al 22 and elastic properties in Li-Al 23. In Li-Al, these force constants predicted elastic constants well, but reproduced vibrational entropies poorly. Transferable force constants have also been explored in oxides 24 and semiconductors 25. In these materials, transferability is examined in the short-ranged force constants, which are obtained after subtracting out an analytical long-ranged dipole-dipole contribution. None of these studies examined the use of transferable force constants to calculate vibrational entropy. This is a challenging issue, as vibrational entropy differences between structures are a fraction of the total vibrational entropy and hence, high relative accuracy is required.

Previous attempts at defining transferable force constant matrices 23 defined force constant matrices for a given pair type. Recently, it was shown that high-accuracy in predicting the vibrational entropy could be achieved by using bond length-dependent transferable force constants 12,26. This simple bond length-dependent force constant model was shown to work well in predicting vibrational entropies in Ni3Al and Pd3V 12,13. The force constants were found to be transferable between different FCC structures at a single composition; their transferability between configurations with different atomic ratios or different chemical systems was untested.

In this paper, we evaluate the accuracy of the length-dependent transferable force constants approach in calculating vibrational entropies in the Au-Cu, Au-Pd, and Cu-Pd systems, and also investigate the transferability of such relations between different chemical systems. Our results show that in the Au-Cu, Au-Pd, and Cu-Pd systems, transferable force constants can be defined that depend on length and pair type. We show that these force constants can be exchanged between different structures and chemical systems, while still retaining good accuracy for calculating the vibrational entropy of intermetallic structures.

## II Methodology

A supercell method was used to calculate force constants 12,26-29. This method consists of creating a series of supercells with well-chosen displacements, calculating the resulting forces on atoms, then using a least-squares fit to obtain the force constants. All total energy and force calculations were performed using the Vienna ab initio Simulation Package (VASP) 30,31, which implements Blöchl’s projector augmented wave (PAW) approach 32,33 within the local density approximation (LDA). The calculations used energy cutoffs of 300-400 eV and $``$150-200 unique $`\stackrel{}{k}`$-points per 4 atoms. This resulted in energies converged to $``$5 meV, and forces converged to 1%. Special care was taken to make sure forces on the initial structure vanished before supercell configurations were generated for the force constant calculation.

The superstructures and their perturbations were determined using an efficient algorithm 26, which generates supercells with perturbed atoms and their corresponding set of linear equations. Third-order anharmonic terms, which are present when positive and negative displacements are not symmetrically equivalent, were eliminated by subtracting forces taken from calculations of perturbations of equal magnitude but opposite sign 26. Vibrational entropy calculations were converged to within 0.015 $`k_B`$: this value is obtained by calculating the configuration-dependent part of the high-temperature limit of the vibrational entropy per atom, $`S_{vib}=k_B\underset{0}{\overset{\mathrm{}}{}}g(\omega )\mathrm{ln}\omega d\omega `$, for increasing force constant range until $`S_{vib}`$changed by less than 0.015 $`k_B`$. Typically, this requires $``$4-5 neighbors for L12 structures and 2 neighbors for SQS-8 structures.

The above procedure was applied to calculate vibrational entropies of L12-ordered and disordered phases in the Au3Cu, Cu3Au, Au3Pd, Pd3Au, Cu3Pd, Pd3Au, and Cu3Pt systems. In addition, we calculated vibrational entropies of L12-ordered Ag3Au and Au3Ag. The disordered structure was approximated by an 8 atom special quasirandom structure (SQS) 34. These structures have been shown to give the best possible approximation to a disordered structure, within a given number of sites. In Ni3Al, Morgan et. al. tested the convergence of vibrational properties with SQS size, and found that a SQS with 8 atoms is a good approximation of the disordered state 35. The SQS-8 used in this study has been used successfully in first-principles vibrational entropy calculations of disordered Ni3Al and Pd3V 12,13.

To define transferable force constant matrices, some approximations are necessary. The form of a force constant matrix, and hence the number of non-zero force constants, depends on the symmetry of structure. In this work, we use a stretching-bending force constant model. In this model, the coordinate system of each force constant matrix is transformed, so that the $`z`$-axis is aligned along the segment joining the two atoms in question. Two further approximations are necessary to obtain transferable force constants. First, bending terms are averaged, so that they are orientation independent. Second, off-diagonal terms are constrained to be zero. The resulting force constant matrix has only two independent terms - a stretching term $`s`$, and a bending term $`b`$ \- and the form of the matrix can be written

$$\left(\begin{array}{cccccccccccccccccccc}b& 0& 0& & & & & & & & & & & & & & & & & \\ 0& b& 0& & & & & & & & & & & & & & & & & \\ 0& 0& s& & & & & & & & & & & & & & & & & \end{array}\right).$$
(1)
The stretching term ($`s)`$ can be used as a qualitative measure of the bond strength. However, it is not possible to keep only the stretching terms ($`s)`$ in the force constant matrix: doing so can result in large errors (0.2 $`k_B)`$ in the vibrational entropy 12,26. We examined two types of errors associated with using stretching-bending force constant matrices: errors introduced in the force constants and errors introduced in the vibrational entropy. For all force constant matrices examined, the rms error associated with constraining the bending terms to be equal was 0.052 eV/A2 per bending term. The rms error associated with setting off-diagonal terms to zero was 0.009 eV/A2 per off-diagonal term. Stretching terms contain no errors associated with the stretching-bending force constant model: they are directly obtained from the coordinate transformation, with no further approximations.

The errors introduced into the vibrational entropy by using the stretching-bending force constant model are shown in Table 1. For 14 out of the 16 structures tested, $`S_{vib}`$ calculated using full-force constant matrices and $`S_{vib}`$ calculated using the stretching-bending force constant model agreed to within 0.01 $`k_B`$. For the L12 Au3Cu structure, the stretching-bending force constant approximation failed when using nearest-neighbor force constants; for the SQS-8 Au3Cu structure, the stretching-bending force constant approximation failed when using both nearest-neighbor and all-neighbor force constants (where “all-neighbor” is defined as the force constant range at which vibrational entropies were converged to within 0.015 $`k_B)`$. In these structures, using stretching-bending force constant matrices resulted in a dynamically unstable structure. Tests of each approximation revealed that the negative phonon modes in the L12 Au3Cu structure arise because the off-diagonal terms to set to zero; negative phonon modes in the SQS-8 Au3Cu arise due to both the neglect of the off-diagonal terms and the constraint on the bending terms. In previous work, the errors in vibrational entropy by using the simplified stretching-bending force constant model have been tested for five structures in the Ni3Al and Pd3V system 12,13. In those systems, the errors on vibrational entropy for all structures were less than 0.01 $`k_B`$/atom.

## III Results and Discussion

Calculated vibrational entropies for all structures examined in the Ag-Ag, Au-Cu, Au-Pd, Cu-Pt, and Cu-Pd system are listed in Table 1. For the Cu-Au system, our results can be compared with previous experimental 10,36 and theoretical 14 work. We calculated $`S_{vib}\left(Cu\right)=`$ 4.87 $`k_B`$ and $`S_{vib}\left(Au\right)=`$2.85 $`k_B`$, which yields for the formation entropies$`S_{vib}^{form}(L1_2Cu_3Au)=`$ 0.10 $`k_B`$ and $`S_{vib}^{form}(L1_2Au_3Cu)=`$ 0.11 $`k_B`$. Our calculated $`S_{vib}^{form}(L1_2Cu_3Au)`$agrees well with previous experimental work, which found $`S_{vib}^{form}(L1_2Cu_3Au)=`$ 0.07 $`\pm `$ 0.03 $`k_B^{36}`$. Our calculated formation entropies also agree well with previous theoretical work, which found $`S_{vib}^{form}(L1_2Cu_3Au)=`$ 0.10 $`k_B`$ and $`S_{vib}^{form}(L1_2Au_3Cu)=`$ 0.14 $`k_B^{14}`$. The entropy change upon disordering in Cu3Au $`\left(\mathrm{\Delta }S_{vib}^{orderdisorder}\left(Cu_3Au\right)\right)`$ is calculated to be 0.07 $`\pm `$ 0.045 $`k_B`$. This value is lower than the experimental value of 0.14 $`\pm `$ 0.05 $`k_B^{10}`$. However, this calculation is in good agreement with previous theoretical work, in which $`\mathrm{\Delta }S_{vib}^{orderdisorder}\left(Cu{}_{3}{}^{}Au\right)=`$ 0.08 $`k_B`$ was obtained 14. For Au3Cu, we calculate $`\mathrm{\Delta }S_{vib}^{orderdisorder}\left(Au{}_{3}{}^{}Cu\right)=0.01\pm 0.03k_B`$, which is lower than the previously calculated $`\mathrm{\Delta }S_{vib}^{orderdisorder}\left(Au_3Cu\right)=`$ 0.05 $`k_B`$ 14.

Calculated full force constant matrices were transformed to the stretching-bending force constant model, with the results plotted in Figures 1-6. The different bond lengths correspond to equilibrium bond lengths in various structures. For all pair types, the force constant stiffness for both stretching and bending terms decreases with increasing bond length. In particular, the stiffness of first nearest-neighbor force-constant matrices shows a strong dependence on bond length. This dependence of bond stiffness on bond length shows that one can not rely too much upon simple bond-counting arguments 26 when predicting or explaining vibrational entropy differences between phases because the bond counting effect does not take into account bond-length changes upon disordering. Figures 1-6 show that small changes in bond length upon disordering can have a large effect on force constants and vibrational entropy. The relationship between atomic relaxations and vibrational entropy has been previously noted 12,16-18.

For a given bond type, stretching and bending force constants for all systems and structures are found to lie on a single curve. Thus, the dependence of bond stiffness on ordering, composition, or chemical system can be explained almost entirely in terms of bond length changes. Remarkably, other effects on force constant stiffness, such as changes in charge density associated with different configurations and chemical systems, are small. This may be because our study is limited to noble-metal intermetallics.

For each pair type, we constructed a relationship between force constant stiffness and bond length. Our goal is to parameterize the stiffness of the stretching and bending force constants as a function of bond length, then use “fitted” force constants from these functions to predict vibrational entropy. Stretching terms were fit to a second-order polynomial; bending terms were fit to a linear function. The fit was restricted to first-nearest neighbors force constants. The rms fitting error over all stretching terms was 0.102 eV/A2; the rms fitting error over all bending terms was 0.048 eV/A2.

Although attempts to use longer-ranged fitted force constants were unsuccessful, the errors introduced by using only nearest-neighbor force constants was small. In the systems studied, vibrational entropy converged quickly with respect to neighbors: the difference between $`S_{vib}`$ calculated using only first nearest-neighbors and $`S_{vib}`$ calculated using all neighbors was typically around $``$0.00-0.02 $`k_B`$,with a maximum difference of 0.03 $`k_B`$ in L12 Pd3Cu. Previous studies on intermetallics and group-IV semiconductors have also shown that vibrational entropy can converge quickly with respect to neighbors 12,13,37.

The results of using fitted force constants to predict vibrational entropies are shown in Table 1. For all structures, the vibrational entropies obtained using fitted force constants agree well with the vibrational entropies obtained using directly calculated force constants. The rms error between $`S_{vib}^{calc,allneigh,fullfc}`$, and $`S_{vib}^{fit,1nn,sb}`$ for all structures investigated was 0.032 $`k_{B,}`$with a maximum error of 0.058 $`k_B`$ in L12 Pd3Au. The errors are $``$1-2% of the absolute vibrational entropy. These calculations show the predictive power of using fitted stretching-bending force constants to predict vibrational entropies.

We also tested the accuracy of using fitted force constants, in representing the effect of homogeneous volume changes. In metals, pairwise expansion can usually not capture the effect of such electron density changes of the energy. The force constants and vibrational entropy of L12 and SQS-8 structures for Pd3Cu, Pd3Au, Cu3Pd, Cu3Au, and Au3Pd were recalculated at a volume $``$2% larger than the equilibrium volume. For the L12 Pd3Cu structure, the vibrational entropies were also calculated at a 4% and 6% volume increase were also calculated. These direct results were then compared with predicted results obtained from using force constants derived from equilibrium volume calculations. None of the force constants at these increased volumes were included in the fit. The results of using these fitted force constants to predict vibrational entropy are shown in Table 2. For all structures at 2% increased volume, the rms error between $`S_{vib}^{calc,allneigh,fullfc}`$ and $`S_{vib}^{fit,1nn,sb}`$ was 0.033 $`k_B`$, with a maximum error of 0.053 $`k_B`$ in L12 Pd3Au. Thus, even when the force constant data is not included in the fit, the length dependent force constant function gives force constants that accurately predict vibrational entropies.

Often, one is interested in vibrational entropy differences, rather than absolute vibrational entropies. We used vibrational entropy data from Tables 1-2 to obtain two types of vibrational entropy differences that are of interest: the vibrational entropy difference between ordered and disordered phases, and vibrational entropy difference between a structure at two volumes separated by $``$2%. The latter is related to the thermal expansion by the expression $`\alpha _L=\frac{\mathrm{\Delta }S}{3B\mathrm{\Delta }V}`$, where B is the bulk modulus. This data is shown in Table 3. In the systems studied, we estimated the effect of vibrations on calculated phase boundaries by using the following procedure: The change in transition temperature when vibrations are included in the phase diagram calculation is given by $`T_{config+vib}^{\alpha \beta }=T_{configonly}^{\alpha \beta }\left(1+\frac{\mathrm{\Delta }S_{vib}^{\alpha \beta }}{\mathrm{\Delta }S_{config}^{\alpha \beta }}\right)^{14}`$. Thus, $`\frac{\mathrm{\Delta }S_{vib}^{\alpha \beta }}{\mathrm{\Delta }S_{config}^{\alpha \beta }}`$ determines the effect of lattice vibrations on calculated phase diagrams. The configurational entropy per atom for a binary solid solution depends on the state of short-range order. It has a maximum value for an ideal (random) solution of $`S_{config}^\alpha =k\left[c\mathrm{ln}c+\left(1c\right)\mathrm{ln}(1c)\right]`$. Hence, the maximum value of $`\mathrm{\Delta }S_{config}^{\alpha \beta }`$ in the systems studied is 0.562 $`k_B`$/atom: this value occurs when a fully ordered state ($`S_{config}^\alpha =0)`$ at $`c`$ = 0.25 transforms to a fully random solid solution. In most cases, the actual value of $`\mathrm{\Delta }S_{config}^{\alpha \beta }`$ will be smaller than this, owing to short-range order in the disordered state and also some disorder in the low temperature phase. Nevertheless, our approximate estimate shows that phase boundaries in the systems studied would change by $``$5-13% if vibrations were included.

The data in Table 3 can be used to examine the ability of transferable force constants to predict vibrational entropy differences. The rms error between calculated and fit differences $`\left(\mathrm{\Delta }S_{vib}^{calc,allneigh,fullfc}\mathrm{\Delta }S_{vib}^{fit,\mathrm{\hspace{0.17em}1}nn,sb}\right)`$ was 0.030 $`k_{B,}`$ with a maximum error of 0.063 $`k_B`$ in Cu3Au. For vibrational entropy differences between structures at volumes 2% apart, the rms error between calculated and fit differences was 0.016 $`k_B`$ with a maximum error of 0.030 $`k_B`$ in L12 Au3Pd. Thus, the transferable force constants used here are able to predict vibrational entropy differences reasonably well.

In examining the data, one can make a few observations with respect to the accuracy of using transferable force constants. First, errors in absolute entropies from using transferable force constants (Tables 1-2) were typically $``$0.00-0.05 $`k_B`$for all structures. This was true for structures with both large and small absolute vibrational entropies. Second, in all cases, the transferable force constants were able to correctly predict the relative hierarchy of entropy differences. Thus, large entropy differences are predicted to be large; small entropy differences are predicted to be small. A small entropy difference was never predicted to be large, or vice versa. Third, when using transferable force constants, errors in entropy differences between the same structure at different volumes tended to cancel; errors in entropy differences between different structures (at the same composition) did not tend to cancel as much (Table 3). Thus, it is likely that there are small structure-dependent contributions to the force constant stiffness that are not captured by our length-dependence force constant model. Fourth, using transferable force constants gave small absolute errors on vibrational entropy differences ($``$0.00-0.07 $`k_B)`$. For small vibrational entropy differences, this sometimes resulted in $``$15-20% errors, whereas the percentage errors in larger vibrational entropy differences were smaller. Lastly, this method can greatly reduce the computational cost of calculating vibrational properties of intermetallics.

## IV Conclusion

The vibrational entropies for 16 structures in the Ag-Au, Au-Cu, Au-Pd, Cu-Pd, and Cu-Pt systems have been calculated. A simplified model was used with only stretching and bending terms, making the force constant matrices independent of symmetry. With these approximations, the form of the force constant matrices was independent of crystal symmetry. We found that the variation of force constants with ordering, composition, or chemical system can be explained almost entirely through changes in bond length. This method represents a promising way to include vibrational effects in phase diagram calculations at a moderate computational cost.

## Acknowledgements

This work was supported by Department of Energy, Office of Basic Energy Sciences under Contract No. DE-FG02-96ER45571. We gratefully acknowledge computing resources provided by NPACI through the Texas Advanced Computing Center. We would also like to thank postdoctoral researcher Dane Morgan for insightful observations.

1J. M. Sanchez, F. Ducastelle, and D. Gratias, Physica A, 128a, 334 (1984).

2F. Ducastelle, Order and phase stability in alloys, Elsevier Science Pub. Co., Amsterdam, 1991.

3D. de Fontaine, Solid State Physics, 47, 33 (1994).

4G. Ceder, M. Asta, and D. de Fontaine, Physica C, 177, 106 (1991).

5A. Zunger, in NATO ASI on Statics and Dynamics of Alloy Phase Transformations,edited by P. E. Turchi and A. Gonis Plenum Press, New York, 1994, Vol. p.361.

6M. Asta, R. McCormack, and D. de Fontaine, Phys. Rev. B, 48, 748 (1993).

7B. Fultz, L. Anthony, L. J. Nagel, et al., Phys. Rev. B, 52, 3315 (1995).

8L. Anthony, J. K. Okamoto, and B. Fultz, Phys. Rev. Lett., 70, 1128 (1993).

9L. Anthony, L. J. Nagel, J. K. Okamoto, et al., Phys. Rev. Lett., 73, 3034 (1994).

10L. J. Nagel, L. Anthony, and B. Fultz, Philos. Mag. Lett., 72, 421 (1995).

11M. E. Manley and B. Fultz, Philos. Mag. B, 80, 1167 (2000).

12A. van de Walle and G. Ceder, Phys. Rev. B, 61, 5972 (2000).

13A. van de Walle, G. Ceder, and U. V. Waghmare, Phys. Rev. Lett., 80, 4911 (1998).

14V. Ozolins, C. Wolverton, and A. Zunger, Phys. Rev. B, 58, R5897 (1998).

15P. D. Tepesch, A. F. Kohan, G. D. Garbulsky, et al., J. Am. Ceram. Soc., 79, 2033 (1996).

16V. Ozolins and M. Asta, Phys. Rev. Lett., 86, 448 (2001).

17M. Asta and V. Ozolins, Phys. Rev. B, 64, art. no. 094104 (2001).

18C. Wolverton and V. Ozolins, Phys. Rev. Lett., 86, 5518 (2001).

19G. D. Garbulsky and G. Ceder, Phys. Rev. B, 49, 6327 (1994).

20I. M. Robertson, J. Phys.: Condens. Matter , 3, 8181 (1991).

21M. Sluiter, M. Weinart, and Y. Kawazoe, Europhys. Lett., 43, 183 (1998).

22M. Sluiter and Y. Kawazoe, Philos. Mag. A, 78, 1353 (1998).

23M. H. F. Sluiter, M. Weinert, and Y. Kawazoe, Phys. Rev. B, 59, 4100 (1999).

24P. Ghosez, E. Cockayne, U. V. Waghmare, et al., Phys. Rev. B, 60, 836 (1999).

25P. Giannozzi, S. de Gironcoli, P. Pavone, et al., Phys. Rev. B, 43, 7231 (1991).

26A. van de Walle and G. Ceder, Rev. Modern Phys., 74, 11 (2002).

27S. Wei and M. Y. Chou, Phys. Rev. Lett., 69, 2799 (1992).

28S. Wei and M. Y. Chou, Phys. Rev. B, 50, 2221 (1994).

29K. Kunc and R. M. Martin, Phys. Rev. Lett., 48, 406 (1982).

30G. Kresse and J. Furthmuller, Comp. Mater. Sci., 6, 15 (1996).

31G. Kresse and J. Furthmuller, Phys. Rev. B, 54, 11169 (1996).

32G. Kresse and D. Joubert, Phys. Rev. B, 59, 1758 (1999).

33P. E. Blochl, Phys. Rev. B, 50, 17953 (1994).

34A. Zunger, S. H. Wei, L. G. Ferreira, et al., Phys. Rev. Lett., 65, 353 (1990).

35D. Morgan, J. D. Althoff, and D. de Fontaine, J. Phase Equilib., 19, 559 (1998).

36P. D. Bogdanoff and B. Fultz, Philos. Mag. B, 79, 753 (1999).

37G. Garbulsky, Ph.D. Thesis, http://africa.mit.edu/papers.


text:
# Two-loop self-energy correction in high-𝐙 hydrogen-like ions

( Feb. 28, 2003)

## Abstract

A complete evaluation of the two-loop self-energy diagrams to all orders in $`Z\alpha `$ is presented for the ground state of H-like ions with $`Z40`$.

The calculation of the two-loop self-energy correction to the Lamb shift is at present one of the most challenging problems in bound-state QED. Until very recently, this project has been addressed to mainly within the $`Z\alpha `$-expansion approach. In it, the two-loop self-energy contribution is represented as an expansion over $`Z\alpha `$ and $`\mathrm{ln}(Z\alpha )`$,

$`F(Z\alpha )`$ $`=`$ $`B_{40}+(Z\alpha )B_{50}+(Z\alpha )^2[L^3B_{63}`$ (1)
$`+L^2B_{62}+LB_{61}+B_{60}]+\mathrm{},`$
where $`F(Z\alpha )=\mathrm{\Delta }E/[m(\alpha /\pi )^2(Z\alpha )^4]`$, $`L=\mathrm{ln}(Z\alpha )^2`$. Whereas the lowest-order term $`B_{40}`$ has been known for a long time, calculations of higher-order contributions have not been accomplished until recently. The correction $`B_{50}`$ was found to be surprisingly large Pachucki94 ; Eides95 , $`B_{50}=24.27`$. This result has significantly changed the theoretical prediction for the Lamb shift in hydrogen and resolved the disagreement with the experimental value existing at that time. The leading logarithmic term $`B_{63}`$ was derived first in Karshenboim93  and later confirmed in Manohar00 ; Pachucki01 . (It should be noted that the first evaluation Karshenboim93 , while yielding the right result, is not completely correct, as will be discussed below.) The two remaining logarithmic corrections $`B_{62}`$ and $`B_{61}`$ have also been elaborated lately by Pachucki Pachucki01 . Again, as in order $`(Z\alpha )^5`$, the result obtained turned out to be surprisingly large. The numerical value of $`B_{61}`$ is $`50.3`$, which reverses the sign of the overall logarithmic contribution for hydrogen. This indicates that the convergence of the $`Z\alpha `$ expansion for the two-loop self-energy correction is remarkably slow, and a conclusion has been drawn in Pachucki01  that a numerical evaluation with Dirac-Coulomb propagators is desirable even for hydrogen.

The calculation of the two-loop self-energy diagrams (Fig. 1) without an expansion in $`Z\alpha `$ started with the irreducible contribution of the diagram (a) (known also as the loop-after-loop (LAL) correction), which is by far the simplest part of the total set. Such an evaluation was first accomplished in Mitrushenkov95  for high-$`Z`$ ions, and later in Mallampalli98  for all ions, including hydrogen. The latter investigation demonstrated a rather peculiar behaviour of the LAL correction in the low-$`Z`$ region. It was shown that for hydrogen its actual value was of a different sign and magnitude than the value based on first two terms of the $`Z\alpha `$ expansion. In addition, a different result was found in Mallampalli98  for the leading logarithmic contribution $`B_{63}`$ as compared to the analytical evaluation Karshenboim93 . (We note that in the latter work the $`B_{63}`$ term was evaluated for the whole set of two-loop self-energy diagrams. However, it was argued that it originated from the LAL contribution only.) As a result, a question was raised in Mallampalli98  about the possibility that the $`Z\alpha `$ expansion for the two-loop self-energy could be inadequate even for hydrogen. This speculation attracted attention and several investigations followed. The subsequent numerical calculation Goidenko99  claimed to be compatible with the analytical result. However, the third numerical evaluation by one of us Yerokhin00lal  confirmed the first result Mallampalli98 . At the same time, the total value of the $`B_{63}`$ contribution was confirmed independently by several groups, e.g., in Manohar00 . To throw light on this intricate situation we performed Yerokhin01lal  an analytic calculation of the $`B_{63}`$ term separately for the LAL correction and found agreement both with our previous numerical result and with that of Mallampalli98 . Our conclusion was that the LAL correction provided an additional cubed logarithmic contribution that had been omitted in the original analytical calculation Karshenboim93 . However, this additional term vanishes when the whole set of two-loop self-energy diagrams is taken into account. Recently, analogous additional terms were reported for the leading logarithmic contribution for $`P`$ states Jentschura02 .

The evaluation of the remaining contributions in Fig. 1 is by far more difficult. These contributions are: the reducible part of the diagram (a), the overlapping diagram (b), and the nested diagram (c). The first attempt to evaluate them to all orders in $`Z\alpha `$ was made by Mallampalli and Sapirstein Mallampalli98a . In that work, the contribution of interest was rearranged in 3 parts, referred to by the authors as the ”$`M`$”, ”$`P`$”, and ”$`F`$” terms. (We will discuss this separation in more detail below.) Mallampalli and Sapirstein calculated only the $`M`$ and $`F`$ terms, while the $`P`$ term was left out since a new numerical technique had to be developed for its computation. In addition, since the numerical procedure turned out to be very time consuming, the actual calculation of the $`M`$ term was carried out only for two ions, uranium and bismuth. Subsequently, in the investigation by two of us Yerokhin01sese  we accomplished the computation of the remaining $`P`$ term for $`Z=83`$, 90, and 92, which formally completed the calculation of the two-loop self-energy. However, as we will see, the rearrangement of the whole correction into the $`M`$, $`P`$, and $`F`$ terms is artificial since all the three are divergent. A proper treatment should include these terms simultaneously. In addition, more than two points (in $`Z`$) are needed in order to analyze the $`Z`$ dependence of the correction and to compare it with the known terms of the $`Z\alpha `$ expansion. All these issues are addressed to in the present investigation.

Let us now turn to the evaluation of the two-loop self-energy diagrams. The first problem to be solved is the separation of ultraviolet (UV) divergences. The standard method of renormalization in QED is developed for diagrams involving only free-electron propagators, treating them in momentum space. Therefore, our strategy is to subtract similar diagrams with electron propagators containing zero or one interaction with the nuclear Coulomb field in order to make the corresponding point-by-point difference UV finite. The subtracted diagrams can be then evaluated in momentum space or in the mixed momentum-coordinate representation. For the first-order self-energy, this approach was first implemented in Snyderman91Blundell91 . The situation is much more difficult in the case of the two-loop self-energy. Here, for the first time, we encounter the overlapping UV divergences. For example, the diagram in Fig. 1(b) can be considered as consisting of two overlapping vertex subdiagrams, each of which is UV divergent. The presence of the overlapping divergences makes the structure of subtraction terms much more elaborate than that in the first order. Moreover, some of these terms contain both bound-electron propagators and UV-divergent subdiagrams. Such situation had never been encountered before, and a new numerical technique had to be developed for the evaluation of these subtraction terms.

Following Mallampalli98a , we rearrange the contribution of the diagrams in Fig. 1 in 4 parts: the LAL part, the $`M`$, $`P`$, and $`F`$ terms. The LAL correction is defined by the irreducible part of the diagram (a). Since its evaluation is relatively easy and has been performed by several groups, we do not discuss it here. The $`M`$ term is diagrammatically represented by Fig. 2. It consists of 3 parts originating from the nested diagram, the overlapping diagram, and the reducible part of the diagram (a). The subtractions in the $`M`$ term are chosen so that each of these 3 parts is separately UV finite. Next, we should account for the subtracted terms. Those that contain only free-electron propagators can be treated in momentum space using the standard Feynman-parametrization technique. For those that involve bound-electron propagators, we introduce additional subtractions that remove the overlapping UV divergences. This is graphically represented by Fig. 3. The corresponding contribution is referred to as the $`P`$ term. It consists of 3 parts, each containing only single UV-divergent subgraphs. Taking the first part as an example, we see that the difference shown in the figure is UV divergent only due to the inner self-energy loop, while the divergence due to the outer self-energy loop is canceled. Finally, we collect all terms we have subtracted and denote them as the $`F`$ term depicted in Fig. 4. It consists of Feynman diagrams that contain free-electron propagators only.

We should also mention the infrared (IR) reference-state divergences that are present in the $`M`$ and $`P`$ terms. These singularities can occur in bound-state QED calculations when energies of the intermediate states in the spectral decomposition of electron propagators coincide with the valence-state energy. An analysis given, e.g., in Mallampalli98a  shows that the IR-divergent terms cancel each other in the sum of the $`M`$ and the $`P`$ term. To sum up our discussion of divergences, we separately write divergent contributions to the terms under consideration:

$`\mathrm{\Delta }E_M`$ $`=`$ $`\mathrm{\Delta }E_M^f\mathrm{\Delta }E_{\mathrm{IR}},`$ (2)
$`\mathrm{\Delta }E_P`$ $`=`$ $`\mathrm{\Delta }E_P^f+\mathrm{\Delta }E_{\mathrm{IR}}+L^{(1)}\mathrm{\Delta }E_{\mathrm{SE}}^{(2+)},`$ (3)
$`\mathrm{\Delta }E_F`$ $`=`$ $`\mathrm{\Delta }E_F^f+B^{(1)}\mathrm{\Delta }E_{\mathrm{SE}}^{(2+)},`$ (4)
where the index $`f`$ labels finite contributions, $`\mathrm{\Delta }E_{\mathrm{IR}}`$ is the IR-divergent contribution, $`L^{(1)}`$ and $`B^{(1)}`$ are the one-loop renormalization constants that fulfill the Ward identity $`L^{(1)}=B^{(1)}`$, and $`\mathrm{\Delta }E_{\mathrm{SE}}^{(2+)}`$ is the many-potential part of the one-loop self-energy correction.

We now turn to the numerical evaluation of these terms. It was carried out in the Feynman gauge. The $`P`$ term was evaluated along the lines described in detail in our previous investigation Yerokhin01sese . The calculation of the $`F`$ and $`M`$ terms was performed in a way, in many respects similar to that of Mallampalli98a . The details of the calculation will be published elsewhere. Here, we focus on major novel features of our evaluation as compared to Mallampalli98a . The first point is a different treatment of the reference-state IR divergences. In Mallampalli98a , they were regulated by altering the valence energy $`\epsilon _a`$ to $`\stackrel{~}{\epsilon }_a=\epsilon _a(1\delta )`$. The actual calculations were performed keeping a finite regulator $`\delta `$, and the limit $`\delta 0`$ was evaluated numerically. According to our experience, that approach, while being technically easy to handle, does not allow to control the accuracy of the computation effectively. In our approach, we introduce some subtractions in order to make the terms under consideration finite, separating IR divergences in the form $`\mathrm{\Delta }E_{\mathrm{IR}}`$. The divergent contributions cancel each other explicitly in the sum, and we can perform the whole computation without introducing any actual IR regulator. However, in order to allow the term-by-term comparison with the previous evaluation Mallampalli98a , we performed our calculations with the regulator $`\delta `$ as well.

The second new feature of our approach is a different procedure employed for the double summation over the partial waves in the $`M`$ term. In Mallampalli98a , the photon angular momenta $`l_1`$ and $`l_2`$ were chosen as independent expansion parameters. We found it technically more convenient to employ for this purpose the absolute values of the relativistic angular parameter $`\kappa `$ of two electron propagators, $`|\kappa |=j+1/2`$. Thus, we turn the nested and overlapping contributions to the $`M`$ term into tables of values $`X_{|\kappa _1|,|\kappa _2|}`$. Next, we perform a resummation of the table: $`X_{|\kappa _1|,|\kappa _2|}Y_{ij}`$, where $`i=\left||\kappa _1||\kappa _2|\right|`$, $`j=|\kappa _1|+|\kappa _2|`$. Finally we sum up the table: first, we fix $`i`$ and extrapolate the sum over $`j`$ to infinity, and then sum over $`i`$ and estimate the tail of the expansion.

Now we discuss the computer time necessary for the evaluation of the $`M`$ term. In the previous evaluation by Mallampalli and Sapirstein, a total time of 7323 h was required for a given value of $`Z`$. In our numerical approach, the typical time of the evaluation of one element $`X_{\kappa _1\kappa _2}`$ is about 1 h for the IBM PWR3 processor with 350 MHz, both for the nested and the overlapping diagram. The typical number of elements for a given $`Z`$ was 440 for the nested diagram and 320 for the overlapping diagram. This shows that the time consumption in our numerical procedure is smaller than that of Mallampalli98a , although it is still very large.

In Table I we present the numerical results for finite parts of the $`M`$, $`P`$, and $`F`$ terms. The table shows that our numerical values for the LAL and $`F`$ terms agree very well with the ones from Mallampalli98 ; Mallampalli98a  but there is a significant deviation for the $`M`$ term. More specifically, our calculation for $`Z=92`$ yields $``$2.137, 4.679(2), and $``$3.837(2) for the reducible, nested, and overlapping contributions to the $`M`$ term, respectively. These results should be compared correspondingly with $``$2.137, 4.669(5), and $``$4.387(5) from Mallampalli98a . We see that the leading source of discrepancy is the overlapping diagram. Taking into account the complexity of the computation, it is difficult to suppose what the reason for this disagreement could be.

As in the case of the one-loop self-energy, the evaluation becomes problematic very fast as $`Z`$ decreases. It is due in part to the fact that some individual contributions exhibit a nearly $`Z`$-independent behaviour, while the total correction scales as $`(Z\alpha )^4`$. Numerical problems restricted our calculation to the region $`Z40`$. In Fig. 5, we compare our non-perturbative results with the known terms of the $`Z\alpha `$ expansion. Although we can not as yet say anything conclusive about the higher-order terms, the figure suggests that the results obtained by two different methods could be compatible.

Summing up, we have evaluated all contributions to the two-loop self-energy correction for H-like ions with $`Z40`$. As this correction has been the major source of the uncertainty of theoretical values for the ground-state Lamb shift in these systems, our evaluation improves their accuracy by an order of magnitude Yerokhin01sese . While the experimental precision for H-like uranium is not presently sufficient to probe the new contribution, this should become possible when the experiment currently planned at GSI Stoehlker00  is completed. The question of utmost importance is to extend the present evaluation to low-$`Z`$ ions, where higher-order terms could enter at the level of experimental interest even at $`Z=`$1 Pachucki01 , as well as to excited states. For the $`2p_{1/2}`$-$`2s`$ transition in Li-like high-$`Z`$ ions the two-loop self-energy presently defines the uncertainty of the theoretical prediction Yerokhin01sese  and can be probed by comparing with experimental data available.

This study was supported in part by RFBR (Grant No. 01-02-17248), by Ministry of Education (Grant No. E02-3.1-49), and by the program ”Russian Universities” (Grant No. UR.01.01.072). V.Y. acknowledges the support from the Ministère de l’Education Nationale et de la Recherche, the foundation ”Dynasty”, and International Center for Fundamental Physics. The computation was performed on the CINES and IDRIS national computer centers.