Maximum Entropy Formulation of the Kirkwood Superposition Approximation

A. Singer^{∗}, Z. Schuss^{†}

April 19, 2004

Abstract

Using a variational formulation, we derive the Kirkwood superposition approximation for [[ COMMENT: the reader will not understand the phrase “entire space”]] systems at equilibrium. We define the entropy of particle triplets,[[ COMMENT: the reader will not understand the phrase “particle triplets”. Because this is an Abstract, you should use the full formal name if you possibly can, e.g., ‘triplet correlation function’ or whatever it is called in Berry, Rice, and Ross, or McQuarrie]] and show that the Kirkwood closure brings the entropy to its maximal value. [[ COMMENT: The following does not belong in an Abstract “The principle of maximum entropy is a well known principle in statistical mechanics and beyond.”]] Our approach leads to a different interpretation for the Kirkwood closure relation, usually explained by probabilistic consideration of dependency and independency of particles. [[COMMENT: the choice of the word ‘explain’ is exquisite, exactly right. Well done!]] find the generalization of the Kirkwood closure for finite volume systems with certain boundary conditions as needed to specify many nonequilibrium systems and most electrochemical or biophysical devices, e.g., a concentration cell, battery, or ionic channel. The generalized closure is slightly different [[ COMMENT: Grammatical English requires “different from” instead of “different than”]] the Kirkwood closure, and enables us to calculate the pair function near the boundary walls of the domain as well as defining characteristics of nonequilibrium quantities, such as flow and conductance. .

∗Department of Applied Mathematics, Tel-Aviv University, Ramat-Aviv, Tel-Aviv 69978, Israel, e-mail: amit sin@zahav.net.il, partially supported by research grants from the Israel Science Foundation, US-Israel Binational Science Foundation, and DARPA

†Department of Applied Mathematics, Tel-Aviv University, Ramat-Aviv, Tel-Aviv 69978, Israel, e-mail: schuss@post.tau.ac.il, partially supported by research grants from the Israel Science Foundation, US-Israel Binational Science Foundation, and DARPA

1 Introduction

The pair correlation function is one of the cornerstones in the theory of simple liquids [1, 3, 2, 4]. Many thermodynamic properties of the fluid can be derived from the pair function. I would cite Rowley, because it is so clear and concrete and computer oriented, even if it is unfashionable. There are mainly two approaches for finding the pair function. The first approach is based on the Ornstein-Zernike integral equation and a closure relation for the direct and indirect correlation functions. Many closure relations fall into this category, such as the Percus-Yevick approximation (PY), the hypernetted chain approximation (HNC) and the mean spherical approximation (MSA). The second approach relies on the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy, which relates the n-th correlation function with the n + 1-th correlation function, and assumes a closure relation that connects them. The Kirkwood superposition that approximates the triplet correlation function by the product of the three pair correlation functions, was the first suggested closure of this kind.

Both approaches and the underlying closure relations have their own advantages and disadvantages. Their success is often compared to either molecular dynamics (MD) or Monte-Carlo (MC) simulations. All closures succeed in the limit of (very) low density. However, when the density is increased they tend to fail, sooner or later. Choosing the “best” closure relation is an art in itself. Obviously, each choice of closure relation results in a different approximate solution for the pair function. The BBGKY equation (with or without the superposition approximation) has a great advantage over any

other approach. For hard spheres it predicts a point where \frac{dp}{dρ} = 0, and

hence it predicts solidification. Neither the PY or HNC theories can do this and so one expects the BBGKY approach to be superior at high densities, which is particularly important in applications to protein and channel biology REF Eisenberg, B., Proteins, Channels, and Crowded Ions. Biophysical Chemistry, 2003. 100: p. 507 - 517.

From a mathematical point of view, all approaches suffer from their lack of error estimates for the approximations being made. This issue is more serious than it may seem. Mathematically, there are no shortage of expansions that do not converge, or converge only after many many terms, or converge in some senses (uniformly, conditionally, asymptotically) and not others. If a few terms of such series are used, and error estimates are not available, it is easy to develop approximations with qualitative properties very different from the original function. Scientifically, mathematical issues can sometimes be sidestepped, if the social process of scientific investigation (‘guess and check’) converges to an agreed solution. If error estimates are not available, the scientific process is XXXX more difficult and may itself (socially) diverge. Our goal is to develop a mathematically defined closure approximation that has general meaning and allows error estimates.

In this paper we derive the Kirkwood superposition approximation from a variational (Euler-Lagrange) formulation. We define an entropy functional for the triplets function, and show that it produces the Kirkwood closure (for systems at equilibrium that fill all space and have ??? boundary conditions at infinity) [[COMMENT I think it important to say what kind of boundary conditions are used at infinity in the ‘fill all space’ case. Are there no boundary conditions, or certain kinds of boundary conditions, or any kind of boundary conditions?]] when it reaches its maximum value, under the constraint that its marginal is the pair function. We hope that this kind of formulation will enable us to find error estimates for this simple closure.

This maximum entropy formulation leads to a different closure for a finite volume system. We write the resulting closure as an integral equation. In the limit of entire space, the solution of this integral equation reduces to the Kirkwood superposition approximation. We describe an iterative procedure for solving the BBGKY equation with this generalized closure.

Using this generalized closure instead of the original Kirkwood closure produces different pair correlation functions. The generalized closure can be used to calculate the pair function in confined geometries with the kind of boundary conditions needed to define the devices of electrochemistry and transport biology. Specifically, we can deduce the resulting pair function near the domain walls (boundaries), and compare it to MC or MD simulations, and to the results of the original Kirkwood closure. Because these boundary conditions are spatially nonuniform, they usually imply flow, i.e., nonequilibrium conditions. The generalized closure thus also allows calculation of defining characteristics of nonequilibrium quantities, such as flow and conductance.

2 Maximum Entropy

Let Ω ⊂ ℝ^{3}. Suppose p_{2}(x_{1}, x_{2}) is a known symmetric probability distribution function (pdf)

∫_{Ω×Ω} p_{2}(x_{1}, x_{2}) dx_{1} dx_{2} = 1
p_{2}(x_{1}, x_{2}) = p_{2}(x_{2}, x_{1})
p_{2}(x_{1}, x_{2}) ≥ 0.

XXXX p_{2}(x_{1}, x_{2}) represents the joint pdf of finding two particles at locations x_{1} and x_{2}, as usually defined in the statistical mechanics of fluids.

Consider the following optimization problem. Find a pdf p_{3}(x_{1}, x_{2}, x_{3}) of the triplets, that satisfies the constraints

ϕ_{1}(p_{3}) = ∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{3} − p_{2}(x_{1}, x_{2}) = 0

ϕ_{2}(p_{3}) = ∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{2} − p_{2}(x_{1}, x_{3}) = 0

ϕ_{3}(p_{3}) = ∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{1} − p_{2}(x_{2}, x_{3}) = 0

which means that p_{2} is the marginal of p_{3}, by the definition of the marginal probability density, and brings the entropy functional

H(p_{3}) = − ∫_{Ω×Ω×Ω} p_{3}(x_{1}, x_{2}, x_{3}) log p_{3}(x_{1}, x_{2}, x_{3}) dx_{1} dx_{2} dx_{3} (1)

to its maximum. This is a variational formulation of the closure problem. To find its solution we use the Lagrange multipliers method. Let λ_{1}(x_{2}, x_{3}), λ_{2}(x_{1}, x_{3}), λ_{3}(x_{1}, x_{2}) be the multipliers of the constraints, and the functional

F (p_{3}, λ_{1}, λ_{2}, λ_{3}) = H + λ_{1}ϕ_{1} + λ_{2}ϕ_{2} + λ_{3}ϕ_{3}

= − ∫_{Ω×Ω×Ω} p_{3}(x_{1}, x_{2}, x_{3}) log p_{3}(x_{1}, x_{2}, x_{3}) dx_{1} dx_{2} dx_{3}

+λ_{1}(x_{2}, x_{3}) (∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{1} − p_{2}(x_{2}, x_{3}))

+λ_{2}(x_{1}, x_{3}) (∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{2} − p_{2}(x_{1}, x_{3}))

+λ_{3}(x_{1}, x_{2}) (∫_{Ω} p_{3}(x_{1}, x_{2}, x_{3}) dx_{3} − p_{2}(x_{1}, x_{2})) .

Since F is independent of the derivatives of p_{3}, the Euler-Lagrange equations read

∇_{x_{1}} [− log p_{3}(x_{1}, x_{2}, x_{3}) − 1 + λ_{1}(x_{2}, x_{3}) + λ_{2}(x_{1}, x_{3}) + λ_{3}(x_{1}, x_{2})] = 0,
∇_{x_{2}} [− log p_{3}(x_{1}, x_{2}, x_{3}) − 1 + λ_{1}(x_{2}, x_{3}) + λ_{2}(x_{1}, x_{3}) + λ_{3}(x_{1}, x_{2})] = 0,
∇_{x_{3}} [− log p_{3}(x_{1}, x_{2}, x_{3}) − 1 + λ_{1}(x_{2}, x_{3}) + λ_{2}(x_{1}, x_{3}) + λ_{3}(x_{1}, x_{2})] = 0.

Therefore

− log p_{3}(x_{1}, x_{2}, x_{3}) − 1 + λ_{1}(x_{2}, x_{3}) + λ_{2}(x_{1}, x_{3}) + λ_{3}(x_{1}, x_{2}) = C, (2)

where C is a constant, or equivalently

p_{3}(x_{1}, x_{2}, x_{3}) = γ_{1}(x_{2}, x_{3})γ_{2}(x_{1}, x_{3})γ_{3}(x_{1}, x_{2}),  (3)

where

γ_{1}(x_{2}, x_{3}) = e^{λ_{1}(x_{2},x_{3})−(C+1)/3}
γ_{2}(x_{1}, x_{3}) = e^{λ_{2}(x_{1},x_{3})−(C+1)/3}
γ_{3}(x_{1}, x_{2}) = e^{λ_{3}(x_{1},x_{2})−(C+1)/3}

Clearly, γ_{i} ≥ 0 (i = 1, 2, 3), and therefore p_{3} ≥ 0. Moreover, γ_{1} = γ_{2} = γ_{3} because p_{2} is symmetric. Let γ = γ_{1} = γ_{2} = γ_{3}. Then p_{3} takes the form

p_{3}(x_{1}, x_{2}, x_{3}) = γ(x_{1}, x_{2})γ(x_{2}, x_{3})γ(x_{1}, x_{3}).  (4)

We find γ from the constraint that p_{2} is the marginal of p_{3}

p_{2}(x_{1}, x_{2}) = γ(x_{1}, x_{2}) ∫_{Ω} γ(x_{1}, x_{3})γ(x_{2}, x_{3}) dx_{3}.  (5)

The symmetry of p_{2}, i.e. p_{2}(x_{1}, x_{2}) = p_{2}(x_{2}, x_{1}) follows that of γ, γ(x_{1}, x_{2}) = γ(x_{2}, x_{1}).

XXXX The pdf p_{2} is unknown. We know, however, that it satisfies the BBGKY equation

0 = fex(x_{1})p_{2}(x_{1}, x_{2}) + f(x_{2}, x_{1})p_{2}(x_{1}, x_{2}) − k_{B}T ∇_{x_{1}}p_{2}(x_{1}, x_{2})

+(N − 2) ∫_{Ω} f(x_{3}, x_{1})p_{3}(x_{1}, x_{2}, x_{3}) dx_{3},

where f(x_{2}, x_{1}) is the force XXXX [[COMMENT: ‘interacted’ is not English, ‘exerted’ is commonly used, but nothing is needed, and “less is more” in good writing, and thinking, as well, I suspect.]] on a particle located at x_{1}, by another particle located at x_{2}, and fex(x_{1}) is an external force field acting on a particle located at x_{1}. Substituting the maximum entropy closure (4) into the BBGKY equation (6), together with the integral equation (5), produces [[COMMENT: ‘results in’ is ungrammatical]] a system of integral equations for p_{2} and γ,

0 = fex(x_{1})p_{2}(x_{1}, x_{2}) + f(x_{2}, x_{1})p_{2}(x_{1}, x_{2}) − k_{B}T ∇_{x_{1}}p_{2}(x_{1}, x_{2})

+(N − 2) ∫_{Ω} f(x_{3}, x_{1})γ(x_{1}, x_{2})γ(x_{2}, x_{3})γ(x_{1}, x_{3}) dx_{3}

p_{2}(x_{1}, x_{2}) = γ(x_{1}, x_{2}) ∫_{Ω} γ(x_{1}, x_{3})γ(x_{2}, x_{3}) dx_{3}  (6)

Solving the system (6) will produce the pair function p_{2}.

3 Two Examples

In this section we discuss two examples in which the system (6) can be solved or simplified.

3.1 Non interacting particles in an external field

Non-interacting particles in an external field are described by f ≡ 0 and the system (6) is reduced to

0 = fex(x_{1})p_{2}(x_{1}, x_{2}) − k_{B}T ∇_{x_{1}}p_{2}(x_{1}, x_{2})  (7)

p_{2}(x_{1}, x_{2}) = γ(x_{1}, x_{2}) ∫_{Ω} γ(x_{1}, x_{3})γ(x_{2}, x_{3}) dx_{3}  (8)

Solving equation (7) yields

p_{2}(x_{1}, x_{2}) = e^{−Uex(x_{1})/k_{B}T} h(x_{2}),  (9)

where fex(x_{1}) = −∇_{x_{1}}Uex(x_{1}), and h(x_{2}) is an arbitrary function (the integration constant). Since p_{2}(x_{1}, x_{2}) = p_{2}(x_{2}, x_{1}), we have

p_{2}(x_{1}, x_{2}) = C^{−1}e^{−[Uex(x_{1})+Uex(x_{2})]/k_{B}T} ,  (10)

where C = [∫_{Ω} e^{−Uex(x)/k_{B}T} dx]^{2} is a normalization constant. As expected,

we find that x_{1}, x_{2} are independent random variables, p_{2}(x_{1}, x_{2}) = p_{1}(x_{1})p_{1}(x_{2}). In this case, the solution to (8) is given by

γ(x_{1}, x_{2}) = \sqrt{p_{2}(x_{1}, x_{2})} = \sqrt{p_{1}(x_{1})}\sqrt{p_{1}(x_{2})}.  (11)

Indeed,

γ(x_{1}, x_{2}) ∫_{Ω} γ(x_{1}, x_{3})γ(x_{2}, x_{3}) dx_{3} =

= \sqrt{p_{1}(x_{1})p_{1}(x_{2})} ∫_{Ω} \sqrt{p_{1}(x_{1})p_{1}(x_{3})}\sqrt{p_{1}(x_{2})p_{1}(x_{3})} dx_{3}

= p_{1}(x_{1})p_{1}(x_{2}) ∫_{Ω} p_{1}(x_{3}) dx_{3} = p_{1}(x_{1}, x_{2}).

3.2 Entire Space Systems

In the limit process Ω → ℝ^{3} all the pdfs tend to zero, so it is more convenient to work with densities. [[COMMENT: I know this will seem trivial to you, but the scientist–or mathematician unfamiliar with probability theory–will be confused here. He/she will think ‘pdf’ = ‘probability density function’ and will wonder how can a function which is already a density be converted again into another density? Please explain just enough so she/he is not confused. Anticipating confusion is a sign of good writing.]] First, let Ω ⊂ ℝ^{3} be a bounded domain. The previous example of non-interacting particles implies the definition

δ(x_{1}, x_{1}) = \frac{γ(x_{1}, x_{2})}{\sqrt{p_{1}(x_{1})p_{1}(x_{2})}},  (12)

so equation (8) reads

\frac{p_{2}(x_{1}, x_{2})}{p_{1}(x_{1})p_{1}(x_{2})} = δ(x_{1}, x_{2}) ∫_{Ω} p_{1}(x_{3})δ(x_{1}, x_{3})δ(x_{2}, x_{3}) dx_{3}.  (13)

We rewrite equation (13) as

\frac{p_{2}^{(Ω)}(x_{1}, x_{2})}{p_{1}^{(Ω)}(x_{1})p_{1}^{(Ω)}(x_{2})} =  (14)

δ^{(Ω)}(x_{1}, x_{2}) ∫_{Ω} p_{1}^{(Ω)}(x_{3})δ^{(Ω)}(x_{1}, x_{3})δ^{(Ω)}(x_{2}, x_{3}) dx_{3},

where p_{2}^{(Ω)}(x_{1}, x_{2}) = p_{2}(x_{1}, x_{2}), δ^{(Ω)}(x_{1}, x_{2}) = δ(x_{1}, x_{2}), p_{1}^{(Ω)}(x_{1}) = p_{1}(x_{1}), to emphasize their dependency on the specific domain Ω. Let

g_{2}(x_{1}, x_{2}) = lim_{Ω→ℝ^{3}} \frac{p_{2}^{(Ω)}(x_{1}, x_{2})}{p_{1}^{(Ω)}(x_{1})p_{1}^{(Ω)}(x_{2})}.  (15)

For example, if the two particles become independent when they are separated,

p_{2}^{(Ω)}(x_{1}, x_{2}) = p_{1}^{(Ω)}(x_{1})p_{1}^{(Ω)}(x_{2}) (1 + o(1)) , for |x_{1} − x_{2}| ≫ 1,  (16)

then lim_{|x_{2}|→∞} g_{2}(x_{1}, x_{2}) = 1.

Next, we show that lim_{Ω→ℝ^{3}} δ^{(Ω)}(x_{1}, x_{2}) = g_{2}(x_{1}, x_{2}), given the assumption (16). Indeed,

∫_{Ω} p_{1}(x_{3})\frac{p_{2}(x_{1}, x_{3})}{p_{1}(x_{1})p_{1}(x_{3})} \frac{p_{2}(x_{2}, x_{3})}{p_{1}(x_{2})p_{1}(x_{3})} dx_{3} = ∫_{Ω} p_{1}(x_{3}) (1 + o(1)) dx_{3}
= 1 + o(1).  (17)

Taking the limit Ω → ℝ^{3} [[COMMENT: introductory adverbial clauses and phrases require a comma in complex sentences and sentences involving mathematics are always in the grammatical class ‘complex’]] the o(1) term vanishes, and equation (14) follows

δ(x_{1}, x_{2}) = g_{2}(x_{1}, x_{2}),  (18)

as asserted.

We interpret equation (18) as the Kirkwood superposition approximation. Equations (4) and (12) implies that the triplet pdf satisfies

\frac{p_{3}^{(Ω)}(x_{1}, x_{2}, x_{3})}{p_{1}^{(Ω)}(x_{1})p_{1}^{(Ω)}(x_{2})p_{1}^{(Ω)}(x_{3})} = δ^{(Ω)}(x_{1}, x_{2})δ^{(Ω)}(x_{1}, x_{3})δ^{(Ω)}(x_{2}, x_{3}).  (19)

Taking the limit Ω → ℝ^{3}, using equation (18), we obtain

g_{3}(x_{1}, x_{2}, x_{3}) = lim_{Ω→ℝ^{3}} \frac{p_{3}^{(Ω)}(x_{1}, x_{2}, x_{3})}{p_{1}^{(Ω)}(x_{1})p_{1}^{(Ω)}(x_{2})p_{1}^{(Ω)}(x_{3})}
= g_{2}(x_{1}, x_{2})g_{2}(x_{1}, x_{3})g_{2}(x_{2}, x_{3}),  (20)

which is the Kirkwood closure relation for the triplet correlation function. [[COMMENT: This is not the version of Kirkwood recommended by Doug. For some reason, which he did not tell me clearly, he was sure that the Quadruplet correlation function was the one that should be closed. Doug said ‘If you use the quadruplet closure, you get the right form for the pair correlation function, and it will fit data very well and be an enormous step forward.’ Doug did not say what was wrong with closing at the triplet level.]] Usually, the motivation for using the Kirkwood superposition approximation is a probabilistic consideration of dependency and independency of particles (see also Section 5.[[COMMENT: If I were a reviewer, I would insist that you give a reference to the original Kirkwood paper, because your work is such a significant extension of the original.]]) Here we find another interpretation for the Kirkwood closure.

The Kirkwood closure is the (only) closure relation that brings the entropy of triplets of particles to its maximum value. The principle of maximum entropy is a well known principle in statistical mechanics, in testing statistical hypotheses [5, 6] and beyond I owe you the references and will give them to you Monday or Tuesday. In the following section we give further motivation for its use.

4 Minimum Helmholtz Free Energy

XXX Elementary textbooks in statistical mechanics mention that the Boltzmann distribution

p_{N}(x_{1}, …, x_{N}) = \frac{1}{Z_{N}} e^{−U(x_{1},...,x_{N})/k_{B}T} ,  (21)

brings the Helmholtz free energy

F (p) = U(p) − k_{B}T H(p)  (22)
= ∫_{Ω^{N}} U(x_{1}, …, x_{N})p(x_{1}, …, x_{N}) dx_{1} ⋯ dx_{N}
+k_{B}T ∫_{Ω^{N}} p(x_{1}, …, x_{N}) log p(x_{1}, …, x_{N}) dx_{1} ⋯ dx_{N},

to its minimum under the normalization constraint

∫_{Ω^{N}} p(x_{1}, …, x_{N}) dx_{1} ⋯ dx_{N} = 1.  (23)

For a pairwise additive potential together with an external field force

U(x_{1}, …, x_{N}) = ∑_{1≤i<j≤N} U(x_{i}, x_{j}) + ∑_{j=1}^{N} Uex(x_{j}),  (24)

the potential energy term U(p) of the Helmholtz free energy (22), takes the simple form

U(p) = ∫_{Ω^{N}} ( ∑_{1≤i<j≤N} U(x_{i}, x_{j}) + ∑_{j=1}^{N} Uex(x_{j})) p(x_{1}, …, x_{N}) dx_{1} ⋯ dx_{N}

= \frac{N(N − 1)}{2} ∫_{Ω^{2}} U(x_{1}, x_{2})p_{2}(x_{1}, x_{2}) dx_{1} dx_{2} + N ∫_{Ω} Uex(x_{1})p_{1}(x_{1}) dx_{1},

where

p_{2}(x_{1}, x_{2}) = ∫_{Ω^{N−2}} p_{N}(x_{1}, …, x_{N}) dx_{1} ⋯ dx_{N},
p_{1}(x_{1}) = ∫_{Ω} p_{2}(x_{1}, x_{2}) dx_{2},

are the marginal distributions. If the pdf p_{2}(x_{1}, x_{2}) is assumed to be known, as in Section 2, then the energy term of the Helmholtz free energy U(p) is also known. Therefore, minimizing the Helmholtz free energy, under the assumption that the pdf p_{2}(x_{1}, x_{2}) is known, is equivalent to maximizing the entropy, since U(p) is constant during the minimizing process. [[COM-MENT: Using the word ‘constant’ by itself is confusing or even frightening to a novice reader, and others as well, in a problem with so many variables. You need to say, nearly always, constant with respect to what. I may have the ‘what’ wrong here. If so, please correct my words, but you should not leave it out altogether, in my view.]]

5 Probabilistic Interpretation of the Kirkwood Closure

The Kirkwood superposition approximation (20) was the first closure relation to be suggested and tested in the theory of simple liquids. This fact might be explained by its simplicity and its intuitive origin. In this section we give a probabilistic interpretation of the Kirkwood superposition approximation, and find its generalization for closure relations of higher orders of the BBGKY hierarchy. The problem at level n (n ≥ 2) is to find an approximation for the n + 1-particle pdf in terms of the n-particle pdf. For example, the Kirkwood superposition approximation (20) closes the hierarchy at level n = 2. High orders closures are expected to be much more accurate, fitting better the experimental and simulated (MC or MD) data REFERENCES. However, the computational complexity increases drastically with n.

First, consider the case n = 2. We assume that particles become independent as they are separated (16), although we are certainly aware that long range forces like the electric field can produce strong correlations over all of even a very large domain. [[COMMENT: The last change is very very important, both because of what we will eventually do with these closure relations and because without this sentence the paper is dangerously misleading.]] In order to make the exact equality (16) into an approximation, we assume that there exists a distance d > 0 such that

p_{2}(x_{1}, x_{2}) = p_{1}(x_{1})p_{1}(x_{2}),  (25)

Figure 1: Four configurations of three particles. (a) no intersections (b) one intersection (c) two intersections (d) three intersections.

for |x_{1} − x_{2}| > d. Three interchangeable particles can be in four different configurations with respect to the distance d, depending on the number of intersections (see Figure 1). In all configurations but configuration (d), where all three particles intersect, there are at least two particles that do not intersect. Since the particles are interchangeable we may assume that |x_{1} − x_{3}| > d. Applying Bayes’ law we have

p_{3}(x_{1}, x_{2}, x_{3}) = p_{3}(x_{3}|x_{1}, x_{2})p_{2}(x_{1}, x_{2}).  (26)

By the independence assumption (25) we have

p_{3}(x_{3}|x_{1}, x_{2}) = p_{2}(x_{3}|x_{2}).  (27)

Therefore,

p_{3}(x_{1}, x_{2}, x_{3}) = p_{2}(x_{3}|x_{2})p_{2}(x_{1}, x_{2}) = \frac{p_{2}(x_{2}, x_{3})}{p_{1}(x_{2})}p_{2}(x_{1}, x_{2}). (28)

Multiplying by 1 = \frac{p_{2}(x_{1}, x_{3})}{p_{1}(x_{1})p_{1}(x_{3})} we obtain

p_{3}(x_{1}, x_{2}, x_{3}) = \frac{p_{2}(x_{1}, x_{2})p_{2}(x_{2}, x_{3})p_{2}(x_{1}, x_{3})}{p_{1}(x_{1})p_{1}(x_{2})p_{1}(x_{3})},  (29)

which is the Kirkwood superposition approximation. We see that the Kirkwood closure is a good approximation when at least two particles are distant, and independent. However, it fails when all three particles are close to each other.

Next, we find the n-level Kirkwood closure relation for the n + 1-particle pdf in terms of the n-particle pdf, using probabilistic considerations.

Proposition: The n-level Kirkwood closure relation is given by

p_{n}(x_{1}, x_{2}, …, x_{n}) = ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}}. (30)

Proof: We have already seen that the approximation holds for n = 2, 3. Assuming, by induction, that at least two particles are distant , and independent, we can assume without loss of generality, that particles 1 and n are far apart, , and independent, |x_{1} − x_{n}| > d. Using Bayes’ law, we find that

p_{n}(x_{1}, x_{2}, …, x_{n}) = p_{n−1}(x_{1}, x_{2}, …, x_{n−1})p_{n}(x_{n}|x_{1}, x_{2}, …, x_{n−1}).

Because particles 1 and n are far, [[Comment: introductory adverbial phrases and clauses require a comma. ‘Since implies a time relation in strict English; ‘because’ is preferred by all grammarians]] apart it follows that

p_{n}(x_{n}|x_{1}, x_{2}, …, x_{n−1}) = p_{n−1}(x_{n}|x_{2}, x_{3}, …, x_{n−1}) = \frac{p_{n−1}(x_{2}, x_{3}, …, x_{n})}{p_{n−2}(x_{2}, x_{3}, …, x_{n−1})}.

Hence,

p_{n}(x_{1}, x_{2}, …, x_{n}) = \frac{p_{n−1}(x_{1}, x_{2}, …, x_{n−1})p_{n−1}(x_{2}, x_{3}, …, x_{n})}{p_{n−2}(x_{2}, x_{3}, …, x_{n−1})}. (31)

It follows from the induction assumption that for every j = 2, 3, …, n − 1, with particles 1 and n far apart, we have

1 = p_{n−1}(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n})  (32)

× ∏_{k=1}^{n−2} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n, i_{l}≠j} p_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}}.

Multiplying eqs.(31) and (32) (for all j = 2, 3, …, n − 1) ends the proof. ⋄

Corollary: For n = 4 Kirkwood’s formula becomes

p_{4}(x_{1}, x_{2}, x_{3}, x_{4})  (33)

= \frac{p_{3}(x_{1}, x_{2}, x_{3})p_{3}(x_{1}, x_{2}, x_{4})p_{3}(x_{1}, x_{3}, x_{4})p_{3}(x_{2}, x_{3}, x_{4})p_{1}(x_{1})p_{1}(x_{2})p_{1}(x_{3})p_{1}(x_{4})}{p_{2}(x_{1}, x_{2})p_{2}(x_{1}, x_{3})p_{2}(x_{1}, x_{4})p_{2}(x_{2}, x_{3})p_{2}(x_{2}, x_{4})p_{2}(x_{3}, x_{4})}.

In the case of the entire space Ω = ℝ^{3}, we define

g_{n}(x_{1}, x_{2}, …, x_{n}) = lim_{Ω→ℝ^{3}} \frac{p_{n}(x_{1}, x_{2}, …, x_{n})}{∏_{j=1}^{n} p_{1}(x_{j})}.  (34)

Dividing equation (30) by ∏_{j=1}^{n} p_{1}(x_{j}) gives

\frac{p_{n}(x_{1}, x_{2}, …, x_{n})}{∏_{j=1}^{n} p_{1}(x_{j})} = ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} (\frac{p_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})}{∏_{j=1}^{k} p_{1}(x_{i_{j}})})^{(−1)^{n−1−k}} ,

(35)

where we used the combinatorial identity

∑_{k=1}^{n−1} \binom{n − 1}{k − 1}(−1)^{n−1−k} = 1.  (36)

Note that the k = 1 terms in the product of equation (35) cancel out, so the product may begin from k = 2

\frac{p_{n}(x_{1}, x_{2}, …, x_{n})}{∏_{j=1}^{n} p_{1}(x_{j})} = ∏_{k=2}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} (\frac{p_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})}{∏_{j=1}^{k} p_{1}(x_{i_{j}})})^{(−1)^{n−1−k}} .

(37)

Taking the limit Ω → ℝ^{3} we obtain the n-level Kirkwood closure relation

g_{n}(x_{1}, x_{2}, …, x_{n}) = ∏_{k=2}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} g_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}}. (38)

Examples:

• n = 3

g_{3}(x_{1}, x_{2}, x_{3}) = g_{2}(x_{1}, x_{2})g_{2}(x_{1}, x_{3})g_{2}(x_{2}, x_{3})

• n = 4

g_{4}(x_{1}, x_{2}, x_{3}, x_{4}) = \frac{g_{3}(x_{1}, x_{2}, x_{3})g_{3}(x_{1}, x_{2}, x_{4})g_{3}(x_{1}, x_{3}, x_{4})g_{3}(x_{2}, x_{3}, x_{4})}{g_{2}(x_{1}, x_{2})g_{2}(x_{1}, x_{3})g_{2}(x_{1}, x_{4})g_{2}(x_{2}, x_{3})g_{2}(x_{2}, x_{4})g_{2}(x_{3}, x_{4})}

6 High Level Entropy Closure

In this section we use the maximum entropy principle to derive the n-level closure relation, and compare the resulting closure relation with the n-level probabilistic Kirkwood closure of Section 5. The problem at level n (n ≥ 2) is to find an approximation for the n + 1-particle pdf in terms of the n-particle pdf. For example, the Kirkwood superposition approximation (20) closes the hierarchy at level n = 2. We use the principle of maximum entropy to obtain the closure relation. Similar to the derivation of Section 2 we assume that the n-particle pdf p_{n}(x_{1}, x_{2}, …, x_{n}) is known, and we search for the n + 1-particle pdf p_{n+1}(x_{1}, x_{2}, …, x_{n+1}) that maximizes the entropy

H = − ∫_{Ω^{n+1}} p_{n+1}(x_{1}, x_{2}, …, x_{n+1}) log p_{n+1}(x_{1}, x_{2}, …, x_{n+1}) dx_{1} dx_{2} ⋯ dx_{n+1},

(39)

with the n + 1 constraints that p_{n} is the marginal of p_{n+1}

p_{n}(x_{1}, x_{2}, …, x_{n}) = ∫_{Ω} p_{n+1}(x_{1}, …, x_{j}, …, x_{n+1}) dx_{j}, j = 1, 2, …, n + 1.

(40)

Since p_{2} is the marginal of p_{n} (n ≥ 2), it follows that p_{2} is also known. Therefore, for a pairwise additive potential, maximizing the Helmholtz free energy of the n + 1-particle system is equivalent to minimizing the entropy of the n + 1-particle system. Introducing the Lagrange multipliers λ_{j}(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n+1}), j = 1, 2, …, n+1, the Euler-Lagrange equation gives

− log p_{n+1} − 1 + ∑_{j=1}^{n+1} λ_{j} = 0.  (41)

Since the n particles are interchangeable,

p_{n+1}(x_{1}, …, x_{n+1}) = ∏_{j=1}^{n+1} γ(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n+1}).  (42)

Integration with respect to x_{n+1} yields

p_{n}(x_{1}, …, x_{n}) = γ(x_{1}, …, x_{n}) ∫_{Ω} ∏_{j=1}^{n} γ(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n+1}) dx_{n+1}.

(43)

Solving the non-linear integral equation (43) for γ and substituting in equation (42) gives the n-level closure relation of the n-level BBGKY hierarchy equation for p_{n}

0 = fex(x_{1})p_{n}(x_{1}, x_{2}, …, x_{n}) + ∑_{j=2}^{n} f(x_{j}, x_{1})p_{n}(x_{1}, x_{2}, …, x_{n})
−k_{B}T ∇_{x_{1}}p_{n}(x_{1}, x_{2}, …, x_{n})  (44)
+(N − n) ∫_{Ω} f(x_{3}, x_{1})p_{n+1}(x_{1}, x_{2}, …, x_{n+1}) dx_{n+1}.

6.1 Entire Space Systems

We have seen in Section 3 that for entire space systems Ω = ℝ^{3}, the maximum entropy principle yields the Kirkwood superposition approximation (n = 2). In this section we show that the maximum entropy principle results in the probabilistic Kirkwood closure (30) for all levels n ≥ 2.

First, let Ω ⊂ ℝ^{3} be a bounded domain. Let

h_{n}^{(Ω)}(x_{1}, x_{2}, …, x_{n}) = \frac{p_{n}^{(Ω)}(x_{1}, x_{2}, …, x_{n})}{∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}}} ,

(45)

and

δ^{(Ω)}(x_{1}, x_{2}, …, x_{n}) = \frac{γ^{(Ω)}(x_{1}, x_{2}, …, x_{n})}{∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k} \frac{n−k}{n−k+1}}} .

(46)

Dividing equation (43) by ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}} we obtain

h_{n}^{(Ω)}(x_{1}, x_{2}, …, x_{n}) = δ^{(Ω)}(x_{1}, x_{2}, …, x_{n}) ∫_{Ω} F (x_{1}, x_{2}, …, x_{n}, x_{n+1})

× ∏_{j=1}^{n} δ^{(Ω)}(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n+1}) dx_{n+1},

(47) where

F (x_{1}, x_{2}, …, x_{n}, x_{n+1}) = ∏_{k=1}^{n−2} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}(n−k−1)}

× ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k−1}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k−1}}, x_{n+1})^{(−1)^{n−1−k}(n−k)}.

For min_{1≤j≤n} |x_{n+1} − x_{j}| ≫ 1, the n + 1-th particle becomes independent of particles 1, 2, …, n, and we have

p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k−1}}, x_{n+1}) = p_{1}^{(Ω)}(x_{n+1})p_{k−1}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k−1}}) (1 + o(1)) ,

(48)

for all k = 1, 2, …, n and all sets of indices 1 ≤ i_{1} < i_{2} < … < i_{k−1} ≤ n. Therefore, to leading order

F (x_{1}, x_{2}, …, x_{n}, x_{n+1}) = ∏_{k=1}^{n−2} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n} p_{k}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−1−k}(n−k−1)}

× ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k−1}≤n} p_{k−1}^{(Ω)}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k−1}})^{(−1)^{n−1−k}(n−k)}

× ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k−1}≤n} p_{1}^{(Ω)}(x_{n+1})^{(−1)^{n−1−k}(n−k)}

= ∏_{k=1}^{n−1} ∏_{1≤i_{1}<i_{2}<...<i_{k−1}≤n} p_{1}^{(Ω)}(x_{n+1})^{(−1)^{n−1−k}(n−k)}

= p_{1}^{(Ω)}(x_{n+1})∑_{k=1}^{n−1}(−1)^{n−1−k}(n − k)\binom{n}{k − 1}

= p_{1}^{(Ω)}(x_{n+1}),  (49)

where we have used the combinatorial identity

∑_{k=1}^{n−1}(−1)^{n−1−k}(n − k)\binom{n}{k − 1} = 1.  (50)

We conclude that

h_{n}^{(Ω)}(x_{1}, x_{2}, …, x_{n}) = δ^{(Ω)}(x_{1}, x_{2}, …, x_{n}) ∫_{Ω} p_{1}^{(Ω)}(x_{n+1}) (1 + o(1))

× ∏_{j=1}^{n} δ^{(Ω)}(x_{1}, …, x_{j−1}, x_{j+1}, …, x_{n+1}) dx_{n+1}.

Therefore,

lim_{Ω→ℝ^{3}} δ^{(Ω)}(x_{1}, x_{2}, …, x_{n}) = lim_{Ω→ℝ^{3}} h_{n}^{(Ω)}(x_{1}, x_{2}, …, x_{n})  (51)

Substituting eq. (46) in eq. (42) while using the relations (45) and (51), and the definition (34) we find that

g_{n+1}(x_{1}, …, x_{n+1}) = ∏_{k=2}^{n} ∏_{1≤i_{1}<i_{2}<...<i_{k}≤n+1} g_{k}(x_{i_{1}}, x_{i_{2}}, …, x_{i_{k}})^{(−1)^{n−k}}. (52)

We observe that for entire space systems, the probabilistic Kirkwood closure

(38) agrees with the maximum entropy closure (52) for all orders n.

6.2 Closure at the highest level n = N − 1

Although the probabilistic Kirkwood closure and the maximum entropy closure agree in the case of entire space systems, they differ in the case of confined systems of finite number of particles N. It appears that the maximum entropy closure (42) is exact in a confined system when applied at the highest level n = N − 1, while the probabilistic Kirkwood closure

(30) is not exact. In other words, the maximum entropy closure relation yields the Boltzmann distribution (21). What appeared at first to be an approximation, turns out to be the exact result. Indeed, let

γ(x_{1}, x_{2}, …, x_{N−1})  (53)

= (\frac{1}{Z_{N}} )^{1/N} exp \{− [\frac{1}{N − 2} ∑_{1≤i<j≤N−1} U(x_{i}, x_{j}) + \frac{1}{N − 1} ∑_{j=1}^{N−1} Uex(x_{j})] /k_{B}T \} .

Clearly,

∏_{j=1}^{N} γ(x_{1}, …, x_{j−1}, x_{j}, …, x_{N}) = \frac{1}{Z_{N}} e^{−U(x_{1},...,x_{N})/k_{B}T} ,  (54)

which is the Boltzmann distribution. The Boltzmann distribution obviously satisfies the BBGKY equation (44). Therefore, we have found a solution to BBGKY equation which satisfies the closure relation (42). This solution coincides with the Boltzmann distribution, and so we conclude that it is the exact solution. For the N − 1-particle pdf p_{N−1} we have

p_{N−1}(x_{1}, …, x_{N−1}) = \frac{1}{Z_{N}} ∫_{Ω} e^{−U(x_{1},...,x_{N})/k_{B}T} dx_{N},  (55)

which is both the exact and “approximate” result.

The probabilistic Kirkwood approximation (30) gives however a different result. If it were to be exact, then its resulted N − 1-particle pdf must be [[COMMENT: You use the subjunctive voice correctly, but no one uses it any more at all, if they ever did!]]given by equation (55). Therefore, all lower level pdf’s p_{n} (n ≤ N − 1) are (multiple) integrals of the Boltzmann distribution. Therefore, by the closure (30), p_{N} should have been a multiplication of integrals of the Boltzmann distribution, which is a contradiction to the known form of the Boltzmann distribution (21).

Although the maximum entropy closure is exact at the highest order n = N − 1, it is not exact at lower orders. The observation [[Comment: I do not know what observation! I also need a description of the REASONS we believe .....”]] made here motivates us to believe that the maximum entropy closure will turn out to be more accurate (i.e., fit the experimental data better) then the probabilistic Kirkwood closure in confined systems, even when used at lower levels (n = 2, 3).

7 Confined Systems

Systems in bounded domains are particularly important, because only bounded domains can include the spatially nonuniform boundary conditions needed to describe devices, with spatially distinct inputs, outputs, and (sometimes) power supplies. A large fraction of electrochemistry involves such devices as batteries or concentration cells. A large fraction of molecular biology involves such devices as proteins that transport ions across otherwise impermeable membranes [7, 8].

In the general case, where Ω ⊂ ℝ^{3} is a bounded domain, there is no analytic solution to the system (6). We propose to solve this system by the following iterative scheme 1. Initial guess γ^{(0)}(x_{1}, x_{2}). Set i = 0.

2. Solve for p_{2}^{(i)}(x_{1}, x_{2}) the non-homogeneous linear equation

k_{B}T ∇_{x_{1}}p_{2}^{(i)}(x_{1}, x_{2}) − f(x_{2}, x_{1})p_{2}^{(i)}(x_{1}, x_{2}) − fex(x_{1})p_{2}^{(i)}(x_{1}, x_{2})

= (N − 2) ∫_{Ω} f(x_{3}, x_{1})γ^{(i)}(x_{1}, x_{2})γ^{(i)}(x_{2}, x_{3})γ^{(i)}(x_{1}, x_{3}) dx_{3}.

3. Solve the non-linear system for γ^{(i+1)}

p_{2}^{(i)}(x_{1}, x_{2}) = γ^{(i+1)}(x_{1}, x_{2}) ∫_{Ω} γ^{(i+1)}(x_{2}, x_{3})γ^{(i+1)}(x_{1}, x_{3}) dx_{3}. (56)

4. i ← i + 1. Return to step 2, until convergence.

The analysis of the previous section indicates that a good initial guess might be

γ^{(0)}(x_{1}, x_{2}) = \sqrt{p_{1}(x_{1})p_{1}(x_{2})} g_{2}(x_{1}, x_{2}),  (57)

where g_{2}(x_{1}, x_{2}) is the solution to the BBGKY equation with the Kirkwood superposition approximation in the entire space ℝ^{3}. This solution can be found rather easily using the inherited symmetries of the problem. For example, it is well known that if fex = 0, then g_{2}(x_{1}, x_{2}) = g_{2}(|x_{1} − x_{2}|), and the problem for g_{2} becomes one dimensional.

Step 2 requires the solution of a linear partial differential equation in a bounded region. This equation can be written in a gradient form

∇_{x_{1}} [e^{(U(x_{1},x_{2})+Uex(x_{1}))/k_{B}T} p_{2}(x_{1}, x_{2})] =

−\frac{N − 2}{k_{B}T}e^{(U(x_{1},x_{2})+Uex(x_{1}))/k_{B}T} ∫_{Ω} ∇_{x_{1}}U(x_{1}, x_{3})p_{3}(x_{1}, x_{2}, x_{3}) dx_{3}.

The identity ∇_{x_{1}} × ∇_{x_{1}}u(x_{1}) = 0, for all u, imposes a solvability condition for γ. Indeed, taking the curl of the last equation, together with the closure

(4) result in

∇_{x_{1}} × [e^{(U(x_{1},x_{2})+Uex(x_{1}))/k_{B}T} ∫_{Ω} ∇_{x_{1}}U(x_{1}, x_{3})γ(x_{1}, x_{2})γ(x_{1}, x_{3})γ(x_{2}, x_{3}) dx_{3}]

= 0.

In step 3 we solve a non-linear integral equation. We suggest solving the non-linear equation (5) by a Newton-Raphson iterative scheme. Let γ^{(n)}(x_{1}, x_{2}) be the n-th iteration. Define the operator Γ^{(n)} : Ω^{2} → Ω^{2} as follows

Γ^{(n)}u(x, z) = ∫_{Ω} γ^{(n)}(x, y)u(z, y) dy.  (58)

Let the operator S : Ω^{2} → Ω^{2} be the symmetrization operator

Su(x, y) = u(y, x).  (59)

The Newton-Raphson iteration scheme suggests

γ^{(n+1)} = γ^{(n)} + Δ(x_{1}, x_{2}),  (60)

where Δ(x_{1}, x_{2}) satisfies the linear integral equation

p_{2}(x_{1}, x_{2}) − p_{2}^{(n)}(x_{1}, x_{2}) = \frac{p_{2}^{(n)}(x_{1}, x_{2})}{γ^{(n)}(x_{1}, x_{2})}Δ(x_{1}, x_{2})
+γ^{(n)}(x_{1}, x_{2})Γ^{(n)}Δ(x_{1}, x_{2})
+γ^{(n)}(x_{1}, x_{2})SΓ^{(n)}Δ(x_{1}, x_{2}),

where

p_{2}^{(n)}(x_{1}, x_{2}) = γ^{(n)}(x_{1}, x_{2}) ∫_{Ω} γ^{(n)}(x_{1}, x_{3})γ^{(n)}(x_{2}, x_{3}) dx_{3}.  (61)

We may write the iteration equivalently as

γ^{(n+1)} = γ^{(n)} + (\frac{p_{2}^{(n)}}{γ^{(n)}} + γ^{(n)}Γ^{(n)} + γ^{(n)}SΓ^{(n)})^{−1} (p_{2} − p_{2}^{(n)}) .  (62)

The algorithm steps are performed until convergence is achieved.

We have yet to test our generalized Kirkwood closure in practice. The resulting pair correlation function should be compared with MD or MC simulations of particles in a confined region. We can only hope that it will outperform the regular Kirkwood superposition approximation in bounded domains. The observation of subsection 6.2 and the generality of the maximum entropy principle (minimum Helmholtz free energy) motivates this belief. We expect to see the difference between the two closures near the boundary walls of the domain. We expect to use our closure to determine the customary properties of nonequilbrium systems, e.g., current and conductance.

8 Mixtures

The maximum entropy principle can also be used to find closure relations of mixtures, both in confined domains and in the entire space. Suppose a mixture of S ≥ 2 species, with N_{α} (α = 1, 2, …, S) particles of each specie. Let N = ∑_{α=1}^{S} be the total number of particles of all specie. There are S^{2} 2-particle pdfs,

p_{2}^{αβ}(x_{1}, x_{2}) α, β ∈ {1, 2, …, S},

that exhibit the symmetry p_{2}^{αβ}(x_{1}, x_{2}) = p_{2}^{βα}(x_{2}, x_{1}). In this section we briefly discuss how to find the closure relation in the mixture problem.

In the maximum entropy approach, one is searching for S^{3} 3-particle pdfs p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}), α, β, γ ∈ {1, 2, …, S}, that bring the entropy

H = − ∑_{α,β,γ=0}^{S} \frac{N_{α}}{N} \frac{N_{β}}{N} \frac{N_{γ}}{N} ∫_{Ω} p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) log p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) dx_{1} dx_{2} dx_{3}

(63)

to maximum, with the 3S^{3} marginal constraints

p_{2}^{αβ}(x_{1}, x_{2}) = ∫_{Ω} p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) dx_{3},  (64)
p_{2}^{αγ}(x_{1}, x_{3}) = ∫_{Ω} p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) dx_{2},
p_{2}^{βγ}(x_{2}, x_{3}) = ∫_{Ω} p_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) dx_{1}.

This variational problem is solved using the Euler-Lagrange formulation similar to the derivation done in Section 2.

In the case of an entire space system, it can be shown in a similar way to the method of subsections 3.2 and 6.1, that the mixture entropy closure coincides with the probabilistic Kirkwood closure. For example, they both suggest XXXX that the triplets correlation functions are related to the pair correlation function via

g_{3}^{αβγ}(x_{1}, x_{2}, x_{3}) = g_{2}^{αβ}(x_{1}, x_{2})g_{2}^{αγ}(x_{1}, x_{3})g_{2}^{βγ}(x_{2}, x_{3}),  (65)

for α, β, γ ∈ {1, 2, …, S}. Similarly, one can obtain XXXX high order closures.

In confined systems, the Euler-Lagrange formulation leads to integral equation of the form (5). Note that since the entropy (63) depends on the particle fraction \frac{N_{α}}{N}, we expect the resulting confined system pair correlation will also depend on the particle fraction.

9 Discussion and Summary

We have used the maximum entropy principle to derive a closure relation for the BBGKY hierarchy. It is possible to consider functionals other than the entropy functional that will yield different (known) closures. In fact, a somewhat similar approach is exercised in the density functional theory (DFT). In the DFT problem setup, one assumes a functional of the pair correlation function (e.g., the Helmholtz free energy) to be maximized. The function that brings the given functional to its maximal value is the resulting pair correlation function. Using this method, one can recover some of the Ornstein-Zernike integral closures, that relate the direct and indirect correlation functions, like the PY closure for instance. Our approach is different from the DFT in that we find a relation between the probability correlation (density) functions of successive orders, rather than a relation between the direct and indirect correlation functions.

In this paper we have used the maximum entropy principle to derive a closure relation for the BBGKY hierarchy. This approach to the closure problem appears to be new. We proved that for entire space systems, the maximum entropy closure relation coincides with the probabilistic Kirkwood superposition approximation for all orders of the hierarchy. In confined systems, with finite number of particles, the maximum entropy closure differ from the Kirkwood superposition approximation. In particular, when applied to the highest level of the hierarchy, the maximum entropy closure is exact, while the probabilistic Kirkwood approximation is not. Together with the generality of the maximum entropy principle, we are lead to believe that the maximum entropy closure will outperform the Kirkwood superposition approximation, even for low order approximations. We expect to notice the differences between the resulting pair correlation function especially near the domain boundaries. We expect to use the closure to make estimates of the nonequilibrium properties that naturally arise in a system with spatially inomogeneous boundary conditions. The implementation and results of using the closure will be the subject of a separate paper.

10 Acknowledgments

Special thanks to Bob Eisenberg, Douglas Henderson and Roland Roth for reading the manuscript and for their useful comments.

References

[1] R. S. Berry, S. Rice, J. Ross, Physical Chemistry, 2nd edition. Oxford University Press, NY 2000.

[2] J.M.G. Barthel, A.H. Baumgartel (Editor), H. Krienke, Physical Chemistry of Electrolyte Solutions: Modern Aspects, Steinkopf, Dietrich Pub. 1998.

[3] H. T. Davis, Statistical Mechanics of Phases, Interfaces, and Thin Films, (Advances in Interfacial Engineering Series), John Wiley & Sons, NY, 1995.

[4] D. Henderson (Editor), Fundamentals of inhomogeneous fluids, Marcel Dekker, New York, 1992

[5] E. L. Lehmann, Testing Statistical Hypotheses, 2nd Ed., Springer Verlag, 1997.

[6] E. L. Lehmann, G. Casella, Theory of Point Estimation, 2nd Ed., Springer Verlag, 1998.

[7] R.S. Eisenberg, “Channels as enzymes.” J. Memb. Biol. 115, pp.112 (1990).

[8] R.S. Eisenberg, “Ion Channels as Devices”, Journal of Computational Electronics 2, pp.245-249 (2003).