Dynamics of interacting Brownian particles: a diagrammatic formulation

# Dynamics of interacting Brownian particles: a diagrammatic formulation

## Abstract

We present a diagrammatic formulation of a theory for the time dependence of density fluctuations in equilibrium systems of interacting Brownian particles. To facilitate derivation of the diagrammatic expansion we introduce a basis that consists of orthogonalized many-particle density fluctuations. We obtain an exact hierarchy of equations of motion for time-dependent correlations of orthogonalized density fluctuations. To simplify this hierarchy we neglect contributions to the vertices from higher-order cluster expansion terms. An iterative solution of the resulting equations can be represented by diagrams with three and four-leg vertices. We analyze the structure of the diagrammatic series for the time-dependent density correlation function and obtain a diagrammatic interpretation of reducible and irreducible memory functions. The one-loop self-consistent approximation for the latter function coincides with mode-coupling approximation for Brownian systems that was derived previously using a projection operator approach.

## I Introduction

There has been a lot of interest in recent years in the dynamics of interacting Brownian particles (1); (2). The reason for this interest is twofold. First, experiments have provided a wealth of information about the motion of individual colloidal particles (3). A system of interacting Brownian particles is the simplest model of a colloidal suspension. Second, interacting Brownian particles constitute the simplest model system on which one can test techniques and approximations of non-equilibrium statistical mechanics. It is a simpler model system than a simple fluid for both fundamental reasons (irreversibility is built in) and technical reasons (particles have fewer degrees of freedom due to the overdamped character of their motion).

In this paper we present a diagrammatic approach to the description of time-dependent density fluctuations in equilibrium systems of interacting Brownian particles (here an equilibrium system means a stable or a metastable (e.g. supercooled) equilibrium system). The original inspiration for this work was a series of three papers (4); (5); (6) by Hans Andersen in which a general framework of a diagrammatic approach to the dynamics of fluctuations in equilibrium simple fluids was presented. An important feature of Andersen’s approach was adoption of a specific set of basis functions, termed Boley (7) basis. As lucidly explained in Ref. (5), one of the advantages of this set of basis functions is an enormous simplification of the initial condition for the whole hierarchy of equations for time-dependent correlation functions (8). Additional motivation for our work comes from renewed interest (9); (10); (11); (12) in developing field theories for systems of strongly interacting particles and in using these theories to generate approximate self-consistent approaches to the dynamics of these systems. Such field theories usually lead to diagrammatic series for so-called response and correlation functions. Our work might be the first step in a reverse procedure: constructing a field theory from a diagrammatic approach. Also, our approach provides a very simple derivation of the mode-coupling theory (13) that has been extensively used to describe colloidal systems. This theory has been previously derived using a projection operator approach (14). More recent, field-theoretical derivations have either been found unsatisfactory (15) or are quite involved (10); (11). Finally, new techniques have recently been developed for strongly correlated many-body quantum systems that allow one to numerically integrate (16) and approximately analyze (17) whole diagrammatic series. It is hoped that these methods could be adopted to classical many-body systems and, in particular, that they could be used to evaluate diagrammatic series presented in this paper.

Our diagrammatic approach to the dynamics of equilibrium systems of interacting Brownian particles is similar to that developed by Andersen (4); (5); (6) for simple fluids. In spite of the fact that a Brownian system is simpler than a simple fluid, in the present problem it is advantageous to introduce two different sets of basis functions. As a consequence, a general structure that leads to the emergence of so-called irreducible memory function appears naturally in the diagrammatic expansion. Our approach uses one important aspect of Andersen’s work: in Ref. (4) the existence of a basis of orthogonalized many-particle phase-space density fluctuations was established. We use a consequence of this result: we assume the existence of a basis consisting of orthogonalized many-particle density fluctuations in the Fourier space. We also assume the existence of a second, closely-related orthogonalized basis of many-particle self-density fluctuations in the Fourier space. The latter basis was used previously in the description of self-diffusion in Brownian systems (18); (19). Our main, formal result is a hierarchy of equations for time-dependent correlations of the orthogonalized many-particle densities. An important feature of this hierarchy is that all the interactions are renormalized: they are expressed in terms of equilibrium correlation functions. To simplify the structure of the hierarchy, we neglect the contributions to the terms describing inter-particle interactions (i.e. vertices) coming from higher-order cluster expansion terms. An iterative solution of the simplified hierarchy can be interpreted in terms of diagrams. After some simplifications we obtain an expansion in terms of diagrams consisting of lines corresponding to free diffusion and three-, and four-leg vertices. We analyze the structure of the diagrammatic expression for the density correlation function and show that so-called irreducible memory function appears in a very natural way. Finally, we present a diagrammatic derivation of the standard, Götze-like (13); (14) mode-coupling approximation.

The paper is organized as follows. In Sec. II, we introduce two sets of basis functions. In Sec. III, we derive exact, formal equations of motion for time-dependent correlations of orthogonalized many-particle densities and in Sec. IV, we simplify these equations by neglecting contributions to the vertices from higher order cluster expansion terms. Sec. V is devoted to the derivation of diagrammatic representation: first, the approximate equations of motion are re-written as integral equations; then, the iterative solution of the latter equations is interpreted in terms of labeled diagrams; finally, a series expansion in terms of labeled diagrams is rewritten in terms of unlabeled diagrams. In Sec. VI, we analyze the series expansion and present diagrammatic expressions for so-called memory and irreducible memory functions. In Sec. VII, we show that a self-consistent one-loop approximation for the irreducible memory function is equivalent to the mode-coupling approximation. We close in Sec. VIII with a discussion of our results, a comparison with other approaches, and an outline of future research.

## Ii Basis functions: orthogonalized many-particle densities

We consider a system of interacting Brownian particles in volume . The average density is . The brackets indicate a canonical ensemble average at temperature . In Secs. II and III, we consider a large but finite system and in Sec. IV, we take the thermodynamic limit,

We start by introducing a set of Fourier transforms of many-particle densities,

 Nk(k1,...,kk)=N∑i1≠...≠ik=1e−ik1⋅ri1−...−ikk⋅rik. (1)

Here and , denote positions of the particles. For simplicity, we will henceforth use term many-particle densities for the Fourier transforms of these densities. Also, we will sometimes use abbreviated notation. Hence may be written as or even as . Also, sum over wavevectors, may be written as . It should be noted that densities are symmetric functions of their arguments.

Following Andersen (4), we introduce orthogonalized many-particle densities using the language of a Hilbert space. The densities are mapped onto vectors,

 Nk(k1,...,kk)↔|Nk(k1,...,kk)⟩, (2)

and the scalar product is defined as

 ⟨Nk|Nl⟩=⟨NkN∗l⟩, (3)

where the asterisk denotes complex conjugation.

To define a set of vectors corresponding to the orthogonalized many-particle densities we start from the 0-particle density,

 |n0⟩↔n0≡N. (4)

Next, we introduce a projection operator onto a subspace spanned by , and define as the part of that is orthogonal to ,

 |n1⟩≡(1−P0)|N1⟩. (5)

Having introduced we can define a projection operator onto a subspace spanned by it. This allows us to define as the part of that is orthogonal to and ,

 |n2⟩≡(1−P0−P1)|N2⟩. (6)

Higher-order orthogonalized many-particle densities can be introduced by continuing this recursive procedure. The orthogonalized densities are symmetric functions of their arguments. The set of the orthogonalized densities constitutes the Boley (7) basis for the present problem.

It should be emphasized that the orthogonalization procedure described above implicitly assumes the existence of the projection operators (4). The simplest, trivial example is that of . We can write as

 P1=∑1,2|n1(1)⟩K1(1;2)⟨n1(2)| (7)

where is the inverse of the ,

 ∑k3⟨n1(k1)n∗1(k3)⟩K1(k3,k2)=δk1,k2. (8)

One notes immediately that , where is the static structure factor, and thus .

In general, we can formally write

 Pk = 1(k!)2∑1,...,k,1′,...,k′|nk(1,..,k)⟩Kk(1,..,k;1′,...,k′) (9) ×⟨nk(1′,...,k′)∣∣.

Here is the inverse of ,

 Fk(1,...,k;1′,...,k′) = ⟨nl(1,...,k)|nk(1′,...,k′)⟩ (10) ≡ ⟨nk(1,...,k)n∗k(1′,...,k′)⟩,
 1k!∑1′′,...,k′′⟨nk(1,...,k)n∗k(1′′,...,k′′)⟩ (11) ×Kl(1′′,...,k′′;1′,...,k′)=Il(1,...,k;1′,...,k′).

In Eq. (11) is an identity defined as

 Il(1,...,k;1′,...,k′)=∑℘(1′,...,k′)δ1,1′...δk,k′, (12)

where denotes a permutation of the arguments , and the sum is over distinct permutations.

The question of the existence of functions is related to the question of the existence of similar functions that was discussed and answered affirmatively in Sec. 3 of Ref. (4) (a careful reader will have by now noticed that we partially follow notation introduced in that paper). The only, minor difference is that the functions considered in this work are Fourier transforms of the many-particle densities in position space whereas the functions considered in Ref. (4) are many-particle densities in phase-space.

It will become clear in the next section that in addition to the set of densities , it is advantageous to introduce another set of orthogonalized densities. This set of densities was implicitly used in investigations of self-diffusion in Brownian systems (18); (19).

 Ns1(k1)=e−ik1⋅r1. (13)

depends on the particle number 1; note that there is nothing special about selecting this particular particle and any other particle can be used in its place.

Next, we define analogous many-particle self-densities,

 Nsk(k1,...,kk)=N∑i1≠...≠ik−1=2e−ik1⋅r1−ik2⋅ri1−...−ikk⋅rik−1, (14)

and associated vectors in the Hilbert space,

 Nsk(k1,...,kk)↔∣∣Nsk(k1,...,kk)⟩, (15)

It should be noted that self-densities are symmetric functions of their last arguments.

Finally, we perform a recursive orthogonalization. To make this procedure similar to that used for many-particle densities we start with the 0-particle self-density,

 ∣∣ns0⟩↔ns0≡1, (16)

and then we define the 1-particle self-density,

 |ns1⟩≡(1−Ps0)|Ns1⟩, (17)

where is a projection operator on a subspace spanned by . Next, we introduce a projection operator onto a subspace spanned by and we define as the part of that is orthogonal to and ,

 ∣∣ns2⟩≡(1−Ps0−Ps1)|Ns2⟩. (18)

Again, higher-order orthogonalized self-densities can be introduced by continuing this procedure.

As before, the orthogonalization procedure relies upon the existence of projection operators . Formally we can write them as

 Psk = 1((k−1)!)2∑1,...,k,1′,...,k′∣∣nsk(1,..,k)⟩ (19) ×Ksk(1,..,k;1′,...,k′)⟨nsk(1′,...,k′)∣∣

Here is the inverse of the ,

 Fsk(1,...,k;1′,...,k′) = ⟨nsk(1,...,k)|nsk(1′,...,k′)⟩ (20) ≡ ⟨nsk(1,...,k)ns∗k(1′,...,k′)⟩,
 Missing or unrecognized delimiter for \left (21) ×Ksk(1′′,...,k′′;1′,...,k′)=Isk(1,...,k;1′,...,k′).

In Eq. (21) is an identity defined as

 Isk(1,...,k;1′,...,k′)=∑℘(2′,...,k′)δ1,1′...δk,k′ (22)

where denotes a permutation of the arguments , and the sum in Eq. (22) is over distinct permutations.

The question of the existence of projection operators is equivalent to that of the existence of functions . Here, we assume here that these functions exists and we leave the proof of this fact for a future study (such a proof probably can be done following the analysis presented in the Appendix B of Ref. (4)).

It should be noted that the bases and are not independent. For example,

 ⟨ns1(k1)|n1(k2)⟩=S(k1)δk1,k2. (23)

However, it is easy to see that

 ⟨nsk|nl⟩=0fork

## Iii Exact, formal equations of motion

We start with a formal expression for the time-dependent correlation function of a -particle density at time and an -particle density at time ,

 ⟨Nkexp(Ωt)N∗l⟩. (25)

Here denotes the Smoluchowski operator,

 Ω=D0N∑i=1∇i⋅(∇i−βFi), (26)

where is the diffusion coefficient of an isolated Brownian particle, denotes a partial derivative with respect to ,

 ∇i=∂∂ri, (27)

with being the Boltzmann constant, and denotes a force acting on particle ,

 Fi=∑j≠iFij=−∑j≠i∇iV(rij). (28)

with being the inter-particle potential. Finally, in expression (25) the equilibrium distribution stands to the right of the quantity being averaged and the Smoluchowski operator, and all other operators act on it as well as on everything else (unless parentheses indicate otherwise).

The orthogonalized many-particle densities are linear combinations of densities and thus we can easily define the following time-dependent correlation functions,

 ⟨nkexp(Ωt)n∗l⟩. (29)

As emphasized in Ref. (5), the advantage of dealing with time-dependent correlation functions (29) is that the initial condition is diagonal, i.e.

 ⟨nkn∗l⟩=0ifk≠l. (30)

Another advantage of using functions (29) is that in equations of motion for bare interactions (i.e. forces , ) are automatically renormalized by equilibrium correlation functions.

To derive a hierarchy of equations of motion for correlation functions (29) we follow Andersen (5) and ascribe the time-dependence to vectors . Explicitly, is defined as the vector associated with

 nk(k1,...,kk;t)≡exp(Ω†t)nk(k1,...,kk) (31)

where denotes the adjoint Smoluchowski operator,

 Missing or unrecognized delimiter for \right (32)

It should be emphasized that the adjoint operator acts only on the densities.

We decompose the time derivative of into a linear combination of ,

 ∂∂t|nk(1,..,k;t)⟩≡∣∣Ω†nk(1,...,k;t)⟩ (33) =∞∑l=01l!∑1,...,lQkl(1,...,k;1,...,l)|nl(1,...,l;t)⟩.

The formulas for the coefficients can be obtained in a number of ways (see, e.g., Ref. (5)). The result is

 Qkl(1,...,k;1,...,l) = 1l!∑1′,...,l′⟨nk(1,...,k)Ωn∗l(1′,...,l′)⟩ (34) ×Kl(1′,...,l′;1,...,l).

Next, we analyze matrix elements of the Smoluchowski operator, . Since all the particles are the same and the equilibrium distribution is symmetric with respect to the particle exchange, we can re-write matrix element in the following way

 ⟨nk(k1,...,kk)Ωn∗l(q1,...,ql)⟩= (35) −D0N⟨(∇1nk(k1,...,kk))⋅(∇1n∗l(q1,...,ql))⟩,

where, as emphasized by the parentheses, derivatives act only on the densities.

It is clear that is a linear combination of with . This allows us to insert projection operators into expression (35) for matrix element ,

 ⟨nk(k1,...,kk)Ωn∗l(q1,...,ql)⟩=−D0N (36) ×min(k,l)∑m=0⟨(∇1nk(k1,...,kk))Psm⋅(∇1n∗l(q1,...,ql))⟩.

Finally, if then, unless or ,

 ⟨(∇1nk(k1,...,kk))ns∗m(q1,...,qm)⟩=0. (37)

Eq. (37) follows from integrating by parts and then using the fact that is orthogonal to all for . As a consequence of Eq. (37), the only non-vanishing matrix elements of the Smoluchowski operator are , and :

 ⟨nk(k1,...,kk)Ωn∗k+1(q1,...,qk+1)⟩= (38) iD0Nk∑i=1ki⋅⟨nsk(ki,k2,...,ki−1,ki+1,...,kk)(∇1n∗k+1(q1,...,qk+1))⟩,
 ⟨nk(k1,...,kk+1)Ωn∗k(q1,...,qk)⟩= (39) −iD0Nk∑i=1⟨(∇1nk+1(k1,...,kk+1))ns∗k(qi,q2,...,qi−1,qi+1,...,qk)⟩⋅qi,
 ⟨nk(k1,...,kk)Ωn∗k(q1,...,qk)⟩= (40) −D0Nk∑i,j=1ki⋅qj⟨nsk(ki,k2,...,ki−1,ki+1,...,kk)ns∗k(qj,q2,...,qj−1,qj+1,...,qk)⟩ −D0N⟨(∇1nk(k1,...,kk))Psk−1(∇1n∗k(q1,...,qk))⟩.

One should note that the diagonal matrix element consists of two different parts. This decomposition, which appears here in a very natural way, will lead to the emergence of an irreducible memory function.

To derive the formulas for coefficients , we contract the expressions for matrix elements with functions . It is obvious that the only non-vanishing coefficients are , and .

We are now in a position to write down a hierarchy of equations of motion for vectors , . This hierarchy could be a starting point for a theory for time-dependent many-particle density correlations. In this paper we are only concerned with the time-dependent single-particle density correlation function, . Thus, rather than presenting the most general hierarchy, we write down an equation of motion for ,

 ∂∂t⟨n1(1;t)n∗1(1′)⟩=∑1′′Q11(1;1′′)⟨n1(1′′;t)n∗1(1′)⟩ (41) +12!∑1′′,2′′Q12(1;1′′,2′′)⟨n2(1′′,2′′;t)n∗1(1′)⟩,

and a hierarchy of equations of motion for functions that couple to , i.e. time-dependent many-particle correlations , ,

 ∂∂t⟨nk(1,...,k;t)n∗1(1′)⟩= (42) 1(k−1)!∑1′′,...,(k−1)′′Qkk−1(1,...,k;1′′,...,(k−1)′′)⟨nk−1(1′′,...,(k−1)′′;t)n∗1(1′)⟩ + 1k!∑1′′,...,k′′Qkk(1,...,k;1′′,...,k′′)⟨nk(1′′,...,k′′;t)n∗1(1′)⟩ + 1(k+1)!∑1′′,...,(k+1)′′Qkk+1(1,...,k;1′′,...,(k+1)′′)⟨nk+1(1′′,...,(k+1)′′;t)n∗1(1′)⟩.

The hierarchy (4142) is the main formal result of this paper. One could now follow Andersen (5) and use Eqs. (4142) as a starting point for a formally exact diagrammatic approach. Here, we follow a different route: first we approximate vertices and then we formulate a diagrammatic approach.

Before introducing approximations, let us comment on general structure of Eqs. (4142). First, a given correlation function couples, via equations of motion, to (except for ) and . Second, the initial condition for this hierarchy is very simple,

 ⟨nk(t=0)n∗1⟩=0ifk>1. (43)

Thus, in a hierarchy of integral equations that is equivalent to Eqs. (4142), and in an iterative solution of this hierarchy, there are no terms related to correlations except for . Third, it can easily be shown that the vertices can be expressed in terms of equilibrium correlation functions. Thus, bare interactions present in a hierarchy of equations of motion for correlation functions have been renormalized (5). In particular, within a simple approximation discussed in the next section, the bare force is replaced by a derivative of a direct correlation function.

## Iv Approximate equations of motion: lowest order cluster expansion terms

Vertices that enter into the exact, formal equations of motion (4142) can be expressed in terms of equilibrium correlation functions. In general, exact expressions for higher order vertices include higher order correlation functions, i.e. correlation functions beyond the pair correlation function . Such higher order correlation functions are not readily available and are usually approximated and/or neglected once formal expressions for time-dependent functions of interest have been derived. Here, we follow an alternative route: we approximate vertices before deriving a diagrammatic expansion.

A complete cluster expansion of vertices can be performed following Sec. II and Appendix A of Ref. (6). We only give expressions for the lowest order terms in the complete cluster expansion. To get these terms it is sufficient to retain only the lowest order terms in the cluster expansions of the matrix elements (38-40) and of functions . The analysis is straightforward albeit the intermediate formulas are rather long. We need the lowest order cluster expansion terms for (and its complex conjugate), and . Including only the lowest order cluster expansion terms, the first quantity is given by the following expression

 ⟨nsk(ki,k2,...,kk)(∇1n∗k+1(q1,...,qk+1))⟩= (44) k+1∑l>j=1∑℘(q1,...,qk+1[qj,ql])⟨ns1(ki)(∇1n∗2(qj,ql))⟩ ×⟨k∏m=1m≠in1(km)k+1∏n=1j≠n≠ln∗1(qn)⟩fac.

Here the notation means remove and from the preceding list and thus denotes a permutation of wavevectors , , and . Finally, in Eq. (44) the following shorthand notation is used,

 ⟨n1(k1)...n1(kk)n∗1(q1)...n∗1(qk)⟩fac≡ (45) ⟨n1(k1)n∗1(q1)⟩...⟨n1(kk)n∗1(qk)⟩.

The second quantity, , is given by

 ⟨nsk(ki,k2,...,ki−1,ki+1,...,kk) (46) ×ns∗k(qj,q2,...,qj−1,qj+1,...,qk)⟩= ∑℘(q1,...,qk[qj])⟨ns1(ki)ns∗1(qj)⟩ ×⟨k∏m=1m≠in1(km)k∏n=1n≠jn∗1(qn)⟩fac,

where the notation means remove from the preceding list. Finally, including only the lowest order cluster expansion terms, has the following simple form

 Kk(k1,...,kk;q1,...,qk)= (47) ∑℘(q1,...,qk)K1(k1;q1)...K1(kk;qk).

We substitute expressions (46-47) into the formulas for vertices , , and and, after some calculations, we obtain

 Qkk+1(k1,...,kk;q1,...,qk+1) = iD0Nk∑i=1k+1∑l>j=1∑q′j,q′lki⋅⟨ns1(ki)(∇1n∗2(q′j,q′l))⟩ (48) Extra open brace or missing close brace
 Qk+1k(k1,...,kk+1;q1,...,qk) = −iD0Nk+1∑j>i=1k∑l=1∑q′l⟨(∇1n2(ki,kj))ns∗1(q′l)⟩⋅q′l (49) ×K1(q′l;ql)Ik−1(k1,...,kk+1[ki,kj]|q1,...,qk[ql]),
 Qkk(k1,...,kk;q1,...,qk)= (50) −D0Nk∑i,j=1∑q′jki⋅q′j⟨ns1(ki)ns∗1(q′j)⟩K1(q′j;qj)Ik−1(k1,...,kk[ki]|q1,...,qk[qj]) −D0Nk∑j>i=1k∑m>l=1∑℘(q1,...,qk[ql,qm])∑q′′l,q′′m∑q′1⟨(∇1n2(ki,kj))ns∗1(q′1)⟩ ×⟨ns1(q′1)(∇1n∗2(q′′l,q′′m))⟩K1(q′′l;ql)K1(q′′m;qm)Ik−2(k1,...,kk[ki,kj]|q1,...,qk[ql,qm]).

The right-hand-sides of expressions (48-50) involve two-particle correlation function (more precisely, its Fourier transform, i.e. the static structure factor) and function . The exact expression for the latter function involves a three-particle correlation function. As is customary (13); (14), we use the convolution approximation for the three-particle contribution to