\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage{latexsym}
\usepackage{amstext}
\usepackage{amssymb}
\usepackage{epsfig}
\usepackage{graphics}

\font \bb=msbm10 at 12pt
\newcommand{\mb}[1]{ \mbox{\boldmath$#1$} }
\newcommand{\ds}{\displaystyle}
\begin{document}

\title{Maximum Entropy Formulation of the Kirkwood Superposition Approximation}
\author{A. Singer\footnote{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},\ Z. Schuss\footnote{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}}
\maketitle

\begin{abstract}
Using a variational formulation, we derive the Kirkwood
superposition approximation for \textbf{[[ COMMENT: the reader
will not understand the phrase ``entire space"]]} systems at
equilibrium. We define the entropy of particle triplets,\textbf{[[
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. \textbf{[[ 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. \textbf{[[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 \textbf{many nonequilibrium systems and most}
electrochemical or biophysical device\textbf{s}, e.g., a
concentration cell, battery, or ionic channel. The generalized
closure is slightly different \textbf{[[ COMMENT: Grammatical
English requires ``different from" instead of ``different than"]]}
the Kirkwood closure, \textbf{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. }.
\end{abstract}

\section{Introduction}
The pair correlation function is one of the cornerstones in the
theory of simple liquids \cite{BRR,Davis,Barthel,Henderson}. Many
thermodynamic properties of the fluid can be derived from the pair
function. \textbf{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 $\ds \frac{dp}{d\rho}=0$, and hence \textbf{it
predicts solidification}. Neither the PY or HNC theories can do
this \textbf{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 \textsf{qualitative}
properties very different from the original function.
Scientifically, mathematical issues can sometimes be sidestepped,
if the   \textbf{social process} of scientific investigation
(`guess and check') converges to an agreed solution. If error
estimates are not available, the scientific process is
\textbf{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 \textbf{(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\textbf{ 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.
\textbf{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.}

\section{Maximum Entropy}\label{sec:entropy}
Let $\Omega \subset \mathbb{R}^3$. Suppose
$p_2(\mb{x}_1,\mb{x}_2)$ is a known symmetric probability
distribution function (pdf)
\begin{eqnarray*}
\int_{\Omega \times \Omega}
p_2(\mb{x}_1,\mb{x}_2)\,d\mb{x}_1\,d\mb{x}_2
&=& 1 \\
&&\nonumber\\
p_2(\mb{x}_1,\mb{x}_2) &=& p_2(\mb{x}_2,\mb{x}_1) \\
&&\nonumber\\
 p_2(\mb{x}_1,\mb{x}_2) &\geq& 0.
\end{eqnarray*}
\textbf{XXXX} $p_2(\mb{x}_1,\mb{x}_2)$ represents the joint pdf of
finding two particles at locations $\mb{x}_1$ and
$\mb{x}_2$,\textbf{ as usually defined in the statistical
mechanics of fluids.}

Consider the following optimization problem. Find a pdf
$p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)$ of the triplets, that satisfies
the constraints
\begin{eqnarray*}
\phi_1(p_3) &=& \int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_3 -
p_2(\mb{x}_1,\mb{x}_2) = 0\\
&&\\ \phi_2(p_3) &=& \int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_2 -
p_2(\mb{x}_1,\mb{x}_3) = 0\\
&&\\ \phi_3(p_3) &=& \int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_1 -
p_2(\mb{x}_2,\mb{x}_3) = 0
\end{eqnarray*}
which means that $p_2$ is the marginal of $p_3$, by the definition
of the marginal probability density, and brings the entropy
functional
\begin{equation}
H(p_3) = -\int_{\Omega \times \Omega \times \Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) \log
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) \,d\mb{x}_1\,d\mb{x}_2\,d\mb{x}_3
\end{equation}
to its maximum. This is a variational formulation of the closure
problem. To find its solution we use the Lagrange multipliers
method. Let $\lambda_1(\mb{x}_2,\mb{x}_3),
\lambda_2(\mb{x}_1,\mb{x}_3), \lambda_3(\mb{x}_1,\mb{x}_2)$ be the
multipliers of the constraints, and the functional
\begin{eqnarray*}
F(p_3, \lambda_1, \lambda_2, \lambda_3) &=&
H+\lambda_1\phi_1 + \lambda_2\phi_2 + \lambda_3\phi_3\\
&&\\ & = &  -\int_{\Omega \times \Omega \times \Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) \log
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) \,d\mb{x}_1\,d\mb{x}_2\,d\mb{x}_3
\nonumber \\
&&\\ & & +\lambda_1(\mb{x}_2,\mb{x}_3)\left(\int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_1 -
p_2(\mb{x}_2,\mb{x}_3) \right) \nonumber \\
&&\\
&&+\lambda_2(\mb{x}_1,\mb{x}_3)\left(\int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_2 -
p_2(\mb{x}_1,\mb{x}_3) \right) \nonumber \\
&&\\
&&+\lambda_3(\mb{x}_1,\mb{x}_2)\left(\int_{\Omega}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_3 -
p_2(\mb{x}_1,\mb{x}_2) \right). \nonumber
\end{eqnarray*}
Since $F$ is independent of the derivatives of $p_3$, the
Euler-Lagrange equations read
\begin{eqnarray*}
\nabla_{\mb{x}_1}\left[-\log p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) -1 +
\lambda_1(\mb{x}_2,\mb{x}_3) + \lambda_2(\mb{x}_1,\mb{x}_3) +
\lambda_3(\mb{x}_1,\mb{x}_2)\right] &=& 0,\nonumber \\
&&\\
 \nabla_{\mb{x}_2}\left[-\log p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)
-1 + \lambda_1(\mb{x}_2,\mb{x}_3) + \lambda_2(\mb{x}_1,\mb{x}_3) +
\lambda_3(\mb{x}_1,\mb{x}_2)\right] &=& 0,\nonumber \\
&&\\ \nabla_{\mb{x}_3}\left[-\log p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)
-1 + \lambda_1(\mb{x}_2,\mb{x}_3) + \lambda_2(\mb{x}_1,\mb{x}_3) +
\lambda_3(\mb{x}_1,\mb{x}_2)\right] &=& 0.\nonumber
\end{eqnarray*}
Therefore
\begin{equation}
-\log p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) -1 +
\lambda_1(\mb{x}_2,\mb{x}_3) + \lambda_2(\mb{x}_1,\mb{x}_3) +
\lambda_3(\mb{x}_1,\mb{x}_2) = C,
\end{equation}
where $C$ is a constant, or equivalently
\begin{equation}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) =
\gamma_1(\mb{x}_2,\mb{x}_3)\gamma_2(\mb{x}_1,\mb{x}_3)\gamma_3(\mb{x}_1,\mb{x}_2),
\end{equation}
where
\begin{eqnarray*}
\gamma_1(\mb{x}_2,\mb{x}_3) &=&
e^{\lambda_1(\mb{x}_2,\mb{x}_3)-(C+1)/3} \\
&&\\ \gamma_2(\mb{x}_1,\mb{x}_3) &=&
e^{\lambda_2(\mb{x}_1,\mb{x}_3)-(C+1)/3} \\
&&\\ \gamma_3(\mb{x}_1,\mb{x}_2) &=&
e^{\lambda_3(\mb{x}_1,\mb{x}_2)-(C+1)/3}
\end{eqnarray*}
Clearly, $\gamma_i \geq 0\, (i=1,2,3)$, and therefore $p_3 \geq
0$. Moreover, $\gamma_1=\gamma_2=\gamma_3$ because $p_2$ is
symmetric. Let $\gamma = \gamma_1 = \gamma_2 = \gamma_3$. Then
$p_3$ takes the form
\begin{equation}
\label{eq:closure} p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) =
\gamma(\mb{x}_1,\mb{x}_2)
\gamma(\mb{x}_2,\mb{x}_3)\gamma(\mb{x}_1,\mb{x}_3).
\end{equation}
We find $\gamma$ from the constraint that $p_2$ is the marginal of
$p_3$
\begin{equation}
\label{eq:int-equation} p_2(\mb{x}_1,\mb{x}_2) =
\gamma(\mb{x}_1,\mb{x}_2)
\int_{\Omega}\gamma(\mb{x}_1,\mb{x}_3)\gamma(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3.
\end{equation}
The symmetry of $p_2$, i.e.
$p_2(\mb{x}_1,\mb{x}_2)=p_2(\mb{x}_2,\mb{x}_1)$ follows that of
$\gamma$, $\gamma(\mb{x}_1,\mb{x}_2)=\gamma(\mb{x}_2,\mb{x}_1)$.

\textbf{XXXX The pdf $p_2$ is unknown. We know, however, that it
satisfies}  the BBGKY equation
\begin{eqnarray*}
0&=&\mb{f}_{\mbox{ex}}(\mb{x}_1)p_2(\mb{x}_1,\mb{x}_2)+\mb{f}(\mb{x}_2,\mb{x}_1)
p_2(\mb{x}_1,\mb{x}_2)- k_BT \nabla_{\mb{x}_1}
p_2(\mb{x}_1,\mb{x}_2)\nonumber
\\
&&\\
 &&+ (N-2) \displaystyle \int_{\Omega}
\mb{f}(\mb{x}_3,\mb{x}_1) p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)
\,d\mb{x}_3,\label{eq:BBGKY}
\end{eqnarray*}
where $\mb{f}(\mb{x}_2,\mb{x}_1)$ is the force \textbf{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
$\mb{x}_1$, by another particle located at $\mb{x}_2$, and
$\mb{f}_{\mbox{ex}}(\mb{x}_1)$ is an external force field acting
on a particle located at $\mb{x}_1$. Substituting the maximum
entropy closure (\ref{eq:closure}) into the BBGKY equation
(\ref{eq:BBGKY}), together with the integral equation
(\ref{eq:int-equation})\textbf{, produces [[COMMENT: `results in'
is ungrammatical]]} a system of integral equations for $p_2$ and
$\gamma$,
\begin{eqnarray}
0&=&\mb{f}_{\mbox{ex}}(\mb{x}_1)p_2(\mb{x}_1,\mb{x}_2)+\mb{f}(\mb{x}_2,\mb{x}_1) p_2(\mb{x}_1,\mb{x}_2)
- k_BT \nabla_{\mb{x}_1} p_2(\mb{x}_1,\mb{x}_2)\nonumber\\
&&\nonumber\\
 &&+ (N-2) \displaystyle \int_{\Omega}
\mb{f}(\mb{x}_3,\mb{x}_1) \gamma(\mb{x}_1,\mb{x}_2)
\gamma(\mb{x}_2,\mb{x}_3)\gamma(\mb{x}_1,\mb{x}_3) \,d\mb{x}_3
\nonumber \\
&&\nonumber\\
 p_2(\mb{x}_1,\mb{x}_2) &=& \gamma(\mb{x}_1,\mb{x}_2)
\int_{\Omega}\gamma(\mb{x}_1,\mb{x}_3)\gamma(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3\label{system}
\end{eqnarray}
Solving the system (\ref{system}) will produce the pair function
$p_2$.

\section{Two Examples}
\label{sec:examples} In this section we discuss two examples in
which the system (\ref{system}) can be solved or simplified.
\subsection{Non interacting particles in an external field}
\textbf{Non-interacting particles in an external field are
described by $\mb{f}\equiv 0$ and the system} (\ref{system}) is
reduced to
\begin{eqnarray}
0&=&\mb{f}_{\mbox{ex}}(\mb{x}_1)p_2(\mb{x}_1,\mb{x}_2)
- k_BT \nabla_{\mb{x}_1} p_2(\mb{x}_1,\mb{x}_2)\label{eq:bolz}\\
&&\nonumber\\
 p_2(\mb{x}_1,\mb{x}_2) &=& \gamma(\mb{x}_1,\mb{x}_2)
\int_{\Omega}\gamma(\mb{x}_1,\mb{x}_3)\gamma(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3\label{eq:int}
\end{eqnarray}
Solving equation (\ref{eq:bolz}) yields
\begin{equation}
p_2(\mb{x}_1,\mb{x}_2) = e^{-U_{\mbox{ex}}(\mb{x}_1)/k_BT}
h(\mb{x}_2),
\end{equation}
where $\mb{f}_{\mbox{ex}}(\mb{x}_1) =
-\nabla_{\mb{x}_1}U_{\mbox{ex}}(\mb{x}_1)$, and $h(\mb{x}_2)$ is
an arbitrary function (the integration constant). Since
$p_2(\mb{x}_1,\mb{x}_2)=p_2(\mb{x}_2,\mb{x}_1)$, we have
\begin{equation}
p_2(\mb{x}_1,\mb{x}_2) = C^{-1}
e^{-\left[U_{\mbox{ex}}(\mb{x}_1)+U_{\mbox{ex}}(\mb{x}_2)
\right]/k_BT},
\end{equation}
where $\ds C =
{\left[\int_{\Omega}e^{-U_{\mbox{ex}}(\mb{x})/k_BT}\,d\mb{x}\right]^2}$
is a normalization constant. As expected, we find that
$\mb{x}_1,\mb{x}_2$ are independent random variables,
$p_2(\mb{x}_1,\mb{x}_2) = p_1(\mb{x}_1)p_1(\mb{x}_2)$. In this
case, the solution to (\ref{eq:int}) is given by
\begin{equation}
\gamma(\mb{x}_1,\mb{x}_2) = \sqrt{p_2(\mb{x}_1,\mb{x}_2)} =
\sqrt{p_1(\mb{x}_1)}\sqrt{p_1(\mb{x}_2)}.
\end{equation}
Indeed,
\begin{eqnarray*}
&&\gamma(\mb{x}_1,\mb{x}_2)\int_{\Omega}
\gamma(\mb{x}_1,\mb{x}_3)\gamma(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3=
\\
&&\\
&&= \sqrt{p_1(\mb{x}_1)p_1(\mb{x}_2)} \int_{\Omega}
\sqrt{p_1(\mb{x}_1)p_1(\mb{x}_3)}
\sqrt{p_1(\mb{x}_2)p_1(\mb{x}_3)}\,d\mb{x}_3\nonumber \\
&&\\
 && =p_1(\mb{x}_1)p_1(\mb{x}_2) \int_{\Omega}
p_1(\mb{x}_3)\,d\mb{x}_3 = p_1(\mb{x}_1,\mb{x}_2).\nonumber
\end{eqnarray*}
\subsection{Entire Space Systems}
\label{sub:entire} In the limit process $\Omega \to \mathbb{R}^3$
all the pdfs tend to zero, so it is more convenient to work with
densities. \textbf{[[COMMENT: I know this will seem trivial to
you, but the scientist--\textsf{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 $\Omega \subset \mathbb{R}^3$ be a bounded domain. The
previous example of non-interacting particles \textbf{implies the
definition}
\begin{equation}
\label{eq:delta-def} \delta(\mb{x}_1,\mb{x}_1) =
\frac{\gamma(\mb{x}_1,\mb{x}_2)}{\sqrt{p_1(\mb{x}_1)p_1(\mb{x}_2)}},
\end{equation}
so equation (\ref{eq:int}) reads
\begin{equation}
\label{eq:delta}
\frac{p_2(\mb{x}_1,\mb{x}_2)}{p_1(\mb{x}_1)p_1(\mb{x}_2)}=\delta(\mb{x}_1,\mb{x}_2)
\int_{\Omega}
p_1(\mb{x}_3)\delta(\mb{x}_1,\mb{x}_3)\delta(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3.
\end{equation}
We rewrite equation (\ref{eq:delta}) as
\begin{eqnarray}
&&\frac{p_2^{(\Omega)}(\mb{x}_1,\mb{x}_2)}{p_1^{(\Omega)}(\mb{x}_1)p_1^{(\Omega)}(\mb{x}_2)}=
\label{eq:int-omega}\\
&&\nonumber\\
&& \delta^{(\Omega)}(\mb{x}_1,\mb{x}_2) \int_{\Omega}
p_1^{(\Omega)}(\mb{x}_3)\delta^{(\Omega)}(\mb{x}_1,\mb{x}_3)
\delta^{(\Omega)}(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3,\nonumber
\end{eqnarray}
where $p_2^{(\Omega)}(\mb{x}_1,\mb{x}_2)=p_2(\mb{x}_1,\mb{x}_2)$,
$\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2)=\delta(\mb{x}_1,\mb{x}_2)$,
$p_1^{(\Omega)}(\mb{x}_1)=p_1(\mb{x}_1)$, to emphasize their
dependency on the specific domain $\Omega$. Let
\begin{equation}
g_2(\mb{x}_1,\mb{x}_2) = \lim_{\Omega \to \mathbb{R}^3}
\frac{p_2^{(\Omega)}(\mb{x}_1,\mb{x}_2)}{p_1^{(\Omega)}(\mb{x}_1)p_1^{(\Omega)}(\mb{x}_2)}.
\end{equation}
For example, if the two particles become independent when they are
separated,
\begin{equation}
\label{eq:independent}
p_2^{(\Omega)}(\mb{x}_1,\mb{x}_2)=p_1^{(\Omega)}(\mb{x}_1)p_1^{(\Omega)}(\mb{x}_2)\left(1+
o(1)\right), \mbox{ for } |\mb{x}_1-\mb{x}_2| \gg 1,
\end{equation}
 then $\lim_{|\mb{x}_2|\to
\infty}g_2(\mb{x}_1,\mb{x}_2) = 1$.

Next, we show that $\lim_{\Omega \to
\mathbb{R}^3}\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2)=g_2(\mb{x}_1,\mb{x}_2)$,
given the assumption (\ref{eq:independent}). Indeed,
\begin{eqnarray}
\int_{\Omega}
p_1(\mb{x}_3)\frac{p_2(\mb{x}_1,\mb{x}_3)}{p_1(\mb{x}_1)p_1(\mb{x}_3)}
\frac{p_2(\mb{x}_2,\mb{x}_3)}{p_1(\mb{x}_2)p_1(\mb{x}_3)}\,d\mb{x}_3
& = & \int_{\Omega} p_1(\mb{x}_3) \left(1+o(1) \right)\,d\mb{x}_3
\nonumber \\
&&\nonumber\\ & = & 1 + o(1).
\end{eqnarray}
Taking the limit $\Omega \to \mathbb{R}^3$  \textbf{[[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 (\ref{eq:int-omega}) follows
\begin{equation}
\label{eq:delta-g2} \delta(\mb{x}_1,\mb{x}_2) =
g_2(\mb{x}_1,\mb{x}_2),
\end{equation}
as asserted.

We interpret equation (\ref{eq:delta-g2}) as the Kirkwood
superposition approximation. Equations (\ref{eq:closure}) and
(\ref{eq:delta-def}) implies that the triplet pdf satisfies
\begin{equation}
\frac{p_3^{(\Omega)}(\mb{x}_1,\mb{x}_2,\mb{x}_3)}
{p_1^{(\Omega)}(\mb{x}_1)p_1^{(\Omega)}(\mb{x}_2)p_1^{(\Omega)}(\mb{x}_3)}=
\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2)\delta^{(\Omega)}(\mb{x}_1,\mb{x}_3)\delta^{(\Omega)}(\mb{x}_2,\mb{x}_3).
\end{equation}
Taking the limit $\Omega \to \mathbb{R}^3$, using equation
(\ref{eq:delta-g2}), we obtain
\begin{eqnarray}
g_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)& =& \lim_{\Omega \to
\mathbb{R}^3}\frac{p_3^{(\Omega)}(\mb{x}_1,\mb{x}_2,\mb{x}_3)}
{p_1^{(\Omega)}(\mb{x}_1)p_1^{(\Omega)}(\mb{x}_2)p_1^{(\Omega)}(\mb{x}_3)}\nonumber\\
&&\nonumber\\
&=&
g_2(\mb{x}_1,\mb{x}_2)g_2(\mb{x}_1,\mb{x}_3)g_2(\mb{x}_2,\mb{x}_3),
\label{eq:Kirkwood}
\end{eqnarray}
which is the Kirkwood closure relation \textbf{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 \ref{sec:probabilistic}.\textbf{[[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
\cite{Lehmann1,Lehmann2} and beyond \textbf{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.

\section{Minimum Helmholtz Free Energy}
\textbf{XXX Elementary textbooks in statistical mechanics mention
that the Boltzmann distribution}
\begin{equation}
\label{eq:boltzmann} p_N(\mb{x}_1,\ldots,\mb{x}_N) = \frac{1}{Z_N}
e^{-U(\mb{x}_1,\ldots,\mb{x}_N)/k_BT},
\end{equation}
brings the Helmholtz free energy
\begin{eqnarray}
F(p) &=& U(p) - k_BTH(p) \label{eq:helmholtz} \\
&=& \int_{\Omega^N}
U(\mb{x}_1,\ldots,\mb{x}_N)p(\mb{x}_1,\ldots,\mb{x}_N)\,d\mb{x}_1\cdots\,d\mb{x}_N
\nonumber \\
&& + k_BT \int_{\Omega^N} p(\mb{x}_1,\ldots,\mb{x}_N)\log
p(\mb{x}_1,\ldots,\mb{x}_N)\,d\mb{x}_1\cdots\,d\mb{x}_N, \nonumber
\end{eqnarray}
to its minimum under the normalization constraint
\begin{equation}
\int_{\Omega^N}p(\mb{x}_1,\ldots,\mb{x}_N)\,d\mb{x}_1\cdots\,d\mb{x}_N
= 1.
\end{equation}
For a pairwise additive potential together with an external field
force
\begin{equation}
U(\mb{x}_1,\ldots,\mb{x}_N) = \sum_{1 \leq i < j \leq N}
U(\mb{x}_i,\mb{x}_j) + \sum_{j=1}^N U_{\mbox{ex}}(\mb{x}_j),
\end{equation}
the potential energy term $U(p)$ of the Helmholtz free energy
(\ref{eq:helmholtz}), takes the simple form
\begin{eqnarray}
U(p) & = & \int_{\Omega^N} \left(\sum_{1 \leq i < j \leq N}
U(\mb{x}_i,\mb{x}_j) + \sum_{j=1}^N U_{\mbox{ex}}(\mb{x}_j)\right)
p(\mb{x}_1,\ldots,\mb{x}_N)\,d\mb{x}_1\cdots\,d\mb{x}_N
\nonumber \\
& = & \frac{N(N-1)}{2} \int_{\Omega^2} U(\mb{x}_1,\mb{x}_2)
p_2(\mb{x}_1,\mb{x}_2)\,d\mb{x}_1\,d\mb{x}_2 + N\int_{\Omega}
U_{\mbox{ex}}(\mb{x}_1) p_1(\mb{x}_1)\,d\mb{x}_1,\nonumber
\end{eqnarray}
where
\begin{eqnarray}
p_2(\mb{x}_1,\mb{x}_2) &=&
\int_{\Omega^{N-2}}p_N(\mb{x}_1,\ldots,\mb{x}_N)\,d\mb{x}_1\cdots\,d\mb{x}_N,\nonumber\\
p_1(\mb{x}_1) &=& \int_{\Omega} p_2(\mb{x}_1,\mb{x}_2)\,d\mb{x}_2,
\nonumber
\end{eqnarray}
are the marginal distributions. If the pdf
$p_2(\mb{x}_1,\mb{x}_2)$ is assumed to be known, as in Section
\ref{sec:entropy}, 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(\mb{x}_1,\mb{x}_2)$ is known, is equivalent \textbf{to}
maximizing the entropy, since $U(p)$ is constant \textbf{during
the minimizing process. [[COMMENT: 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.]]}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Sunday Morning in Parkis%%%%%%%%%%%%%%%%%%

\section{Probabilistic Interpretation of the Kirkwood Closure}
\label{sec:probabilistic} The Kirkwood superposition approximation
(\ref{eq:Kirkwood}) 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 \textbf{of higher orders} of the BBGKY
hierarchy. The problem at level $n$ ($n\geq 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 (\ref{eq:Kirkwood}) 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 {\bf
REFERENCES}. However, the computational complexity increases
drastically with $n$.

First, consider the case $n=2$. We assume that particles become
independent as they are \textbf{separated} (\ref{eq:independent}),
\textbf{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 (\ref{eq:independent}) into an approximation, we assume
that there exists a distance $d>0$ such that
\begin{equation}
\label{eq:ind} p_2(\mb{x}_1,\mb{x}_2) =
p_1(\mb{x}_1)p_1(\mb{x}_2),
\end{equation}
for $|\mb{x}_1-\mb{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
\ref{f:configurations}).
\begin{figure}
\mbox{
\begin{minipage} {\textwidth}
\begin{center}
\begin{tabular}{c}
\psfig {figure=configurations.eps,width=7cm}\\
\end{tabular}
\end{center}
\end{minipage}}
\\
\caption{\small Four configurations of three particles. (a) no
intersections (b) one intersection (c) two intersections (d) three
intersections.} \label{f:configurations}
\end{figure}
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 $|\mb{x}_1-\mb{x}_3|>d$. Applying Bayes' law we have
\begin{equation}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) =
p_3(\mb{x}_3|\mb{x}_1,\mb{x}_2)p_2(\mb{x}_1,\mb{x}_2).
\end{equation}
By the \textbf{independence} assumption (\ref{eq:ind}) we have
\begin{equation}
p_3(\mb{x}_3|\mb{x}_1,\mb{x}_2) = p_2(\mb{x}_3|\mb{x}_2).
\end{equation}
\textbf{Therefore,}
\begin{equation}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) = p_2(\mb{x}_3|\mb{x}_2)
p_2(\mb{x}_1,\mb{x}_2) =
\frac{p_2(\mb{x}_2,\mb{x}_3)}{p_1(\mb{x}_2)}p_2(\mb{x}_1,\mb{x}_2).
\end{equation}
Multiplying by $\ds 1 =
\frac{p_2(\mb{x}_1,\mb{x}_3)}{p_1(\mb{x}_1)p_1(\mb{x}_3)}$ we
obtain
\begin{equation}
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) =
\frac{p_2(\mb{x}_1,\mb{x}_2)p_2(\mb{x}_2,\mb{x}_3)p_2(\mb{x}_1,\mb{x}_3)}{p_1(\mb{x}_1)p_1(\mb{x}_2)p_1(\mb{x}_3)},
\end{equation}
which is the Kirkwood superposition approximation. We see that the
Kirkwood closure is a good approximation when at least two
particles are distant\textbf{, 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.\\

\noindent {\bf Proposition:} {\em The $n$-level Kirkwood closure
relation is given by
\begin{equation}
\label{eq:prob-Kirkwood} p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n}
p_k(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}}.
\end{equation}
}\\
\noindent{\bf Proof:} We have already seen that the approximation
holds for $n=2,3$. Assuming, by induction, that at least two
particles are distant \textbf{, and independent, } we can assume
without loss of generality, that particles $1$ and $n$ are far
apart, \textbf{, and independent, } $|\mb{x}_1-\mb{x}_n|>d$. Using
Bayes' law, we find that
\[p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) = p_{n-1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n-1})
p_n(\mb{x}_n|\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n-1}).\]
\textbf{Because} particles $1$ and $n$ are far, \textbf{[[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(\mb{x}_n|\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n-1}) =
p_{n-1}(\mb{x}_n|\mb{x}_2,\mb{x}_3,\ldots,\mb{x}_{n-1}) =
\frac{p_{n-1}(\mb{x}_2,\mb{x}_3,\ldots,\mb{x}_n)}{p_{n-2}(\mb{x}_2,\mb{x}_3,
\ldots,\mb{x}_{n-1})}.\] Hence,
\begin{equation}
\label{eq:step0} p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\frac{p_{n-1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n-1})p_{n-1}(\mb{x}_2,\mb{x}_3,
\ldots,\mb{x}_n)}{p_{n-2}(\mb{x}_2,\mb{x}_3,\ldots,\mb{x}_{n-1})}.
\end{equation}
It follows from the induction assumption that for every $j=2,3,\ldots,n-1$, with
particles $1$ and $n$ far apart, we have
\begin{eqnarray}
1&=& p_{n-1}(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_n) \label{eq:prod1}\\
&& \times \prod_{k=1}^{n-2} \prod_{1 \leq
i_1 < i_2 < \ldots <i_k \leq n, \; i_l \neq j}
p_k(\mb{x}_{i_1},\mb{x}_{i_2},\ldots,\mb{x}_{i_k})^{(-1)^{n-1-k}}.\nonumber
\end{eqnarray}
Multiplying eqs.(\ref{eq:step0}) and (\ref{eq:prod1}) (for all $j=2,3,\ldots,n-1$) ends the
proof. $\diamond$\\

\noindent {\bf Corollary:} {\em For $n=4$ Kirkwood's formula
becomes}
\begin{eqnarray}
&&p_4(\mb{x}_1,\mb{x}_2,\mb{x}_3,\mb{x}_4) \\
&&=
\frac{p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)p_3(\mb{x}_1,\mb{x}_2,\mb{x}_4)
p_3(\mb{x}_1,\mb{x}_3,\mb{x}_4)p_3(\mb{x}_2,\mb{x}_3,\mb{x}_4)p_1(\mb{x}_1)p_1(\mb{x}_2)p_1(\mb{x}_3)p_1(\mb{x}_4)}
{p_2(\mb{x}_1,\mb{x}_2)p_2(\mb{x}_1,\mb{x}_3)p_2(\mb{x}_1,\mb{x}_4)p_2(\mb{x}_2,\mb{x}_3)
p_2(\mb{x}_2,\mb{x}_4)p_2(\mb{x}_3,\mb{x}_4)}. \nonumber
\end{eqnarray}
In the case of the entire space $\Omega=\mathbb{R}^3$, we define
\begin{equation}
\label{eq:g-def} g_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\lim_{\Omega \to \mathbb{R}^3}
\frac{p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)}{\prod_{j=1}^n
p_1(\mb{x}_j)}.
\end{equation}
Dividing equation (\ref{eq:prob-Kirkwood}) by $\prod_{j=1}^n
p_1(\mb{x}_j)$ gives
\begin{equation}
\label{eq:prob-1}
\frac{p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)}{\prod_{j=1}^n
p_1(\mb{x}_j)} = \prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots
< i_k \leq n} \left(\frac{p_k(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})}{\prod_{j=1}^k
p_1(\mb{x}_{i_j})}\right)^{(-1)^{n-1-k}},
\end{equation}
where we used the combinatorial identity
\begin{equation}
\sum_{k=1}^{n-1} {n-1 \choose k-1}(-1)^{n-1-k} = 1.
\end{equation}
Note that the $k=1$ terms in the product of equation
(\ref{eq:prob-1}) cancel out, so the product may begin from $k=2$
\begin{equation}
\frac{p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)}{\prod_{j=1}^n
p_1(\mb{x}_j)} = \prod_{k=2}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots
< i_k \leq n} \left(\frac{p_k(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})}{\prod_{j=1}^k
p_1(\mb{x}_{i_j})}\right)^{(-1)^{n-1-k}}.
\end{equation}
Taking the limit $\Omega \to \mathbb{R}^3$ we obtain the $n$-level
Kirkwood closure relation
\begin{equation}
\label{eq:prob-g} g_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\prod_{k=2}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n}
g_k(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}}.
\end{equation}
Examples:
\begin{itemize}
\item $n=3$
$$g_3(\mb{x}_1,\mb{x}_2,\mb{x}_3) = g_2(\mb{x}_1,\mb{x}_2)g_2(\mb{x}_1,\mb{x}_3)g_2(\mb{x}_2,\mb{x}_3)$$
\item $n=4$
$$g_4(\mb{x}_1,\mb{x}_2,\mb{x}_3,\mb{x}_4) = \frac{g_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)g_3(\mb{x}_1,\mb{x}_2,\mb{x}_4)
g_3(\mb{x}_1,\mb{x}_3,\mb{x}_4)g_3(\mb{x}_2,\mb{x}_3,\mb{x}_4)}{g_2(\mb{x}_1,\mb{x}_2)g_2(\mb{x}_1,\mb{x}_3)
g_2(\mb{x}_1,\mb{x}_4)g_2(\mb{x}_2,\mb{x}_3)g_2(\mb{x}_2,\mb{x}_4)g_2(\mb{x}_3,\mb{x}_4)}
$$
\end{itemize}

\section{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 \ref{sec:probabilistic}. The problem at level $n$ ($n\geq
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 (\ref{eq:Kirkwood}) 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 \ref{sec:entropy} we assume that the $n$-particle pdf
$p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)$ is known, and we search
for the $n+1$-particle pdf
$p_{n+1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n+1})$ that maximizes
the entropy
\begin{equation}
H = -\int_{\Omega^{n+1}}
p_{n+1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n+1})\log
p_{n+1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n+1}) \, d\mb{x}_1
\,d\mb{x}_2\cdots \,d\mb{x}_{n+1},
\end{equation}
with the $n+1$ constraints that $p_n$ is the marginal of $p_{n+1}$
\begin{equation}
p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\int_{\Omega}p_{n+1}(\mb{x}_1,\ldots,\mb{x}_j,\ldots,\mb{x}_{n+1})\,d\mb{x}_j,
\quad j=1,2,\ldots,n+1.
\end{equation}
Since $p_2$ is the marginal of $p_n$ ($n \geq 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
$\lambda_j(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_{n+1})$,
$j=1,2,\ldots,n+1$, the Euler-Lagrange equation gives
\begin{equation}
-\log p_{n+1} - 1 + \sum_{j=1}^{n+1} \lambda_j = 0.
\end{equation}
Since the $n$ particles are interchangeable,
\begin{equation}
\label{eq:p_n+1} p_{n+1}(\mb{x}_1,\ldots,\mb{x}_{n+1}) =
\prod_{j=1}^{n+1}
\gamma(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_{n+1}).
\end{equation}
Integration with respect to $\mb{x}_{n+1}$ yields
\begin{equation}
\label{eq:p_n-gamma}p_n(\mb{x}_1,\ldots,\mb{x}_n) =
\gamma(\mb{x}_1,\ldots,\mb{x}_n) \int_{\Omega}\prod_{j=1}^n
\gamma(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_{n+1})\,d\mb{x}_{n+1}.
\end{equation}
Solving the non-linear integral equation (\ref{eq:p_n-gamma}) for
$\gamma$ and substituting in equation (\ref{eq:p_n+1})
\textbf{gives }the $n$-level closure relation of the $n$-level
BBGKY hierarchy equation for $p_n$
\begin{eqnarray}
0&=&\mb{f}_{\mbox{ex}}(\mb{x}_1)p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)+\sum_{j=2}^n
\mb{f}(\mb{x}_j,\mb{x}_1)
p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)\nonumber \\
&&- k_BT \nabla_{\mb{x}_1}
p_n(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) \label{eq:BBGKY-n}\\
 &&+ (N-n) \displaystyle \int_{\Omega}
\mb{f}(\mb{x}_3,\mb{x}_1)
p_{n+1}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{n+1})
\,d\mb{x}_{n+1}.\nonumber
\end{eqnarray}

\subsection{Entire Space Systems}
\label{sub:entire2} We have seen in Section \ref{sec:examples}
that for entire space systems $\Omega=\mathbb{R}^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
(\ref{eq:prob-Kirkwood}) for all levels $n\geq 2$.

First, let $\Omega \subset \mathbb{R}^3$ be a bounded domain. Let
\begin{equation}
\label{eq:h_n} h_n^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\frac{p_n^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)}{\prod_{k=1}^{n-1}\prod_{1
\leq i_1 < i_2 < \ldots < i_k \leq n}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}}},
\end{equation}
and
\begin{equation}
\label{eq:delta-n}
\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\frac{\gamma^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)}{\prod_{k=1}^{n-1}\prod_{1
\leq i_1 < i_2 < \ldots < i_k \leq n}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}\frac{n-k}{n-k+1}}}.
\end{equation}
Dividing equation (\ref{eq:p_n-gamma}) by
$\prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}}$ we obtain
\begin{eqnarray}
h_n^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) &=&
\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) \int_{\Omega}
F(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n,\mb{x}_{n+1}) \nonumber \\
&& \times \prod_{j=1}^n
\delta^{(\Omega)}(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_{n+1})\,d\mb{x}_{n+1}
, \nonumber \\ \label{eq:g_n-delta}
\end{eqnarray}
where
\begin{eqnarray}
F(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n,\mb{x}_{n+1}) &=&
\prod_{k=1}^{n-2}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}(n-k-1)} \nonumber \\
&& \times \prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots <
i_{k-1} \leq n} p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_{k-1}},\mb{x}_{n+1})^{(-1)^{n-1-k}(n-k)}. \nonumber
\end{eqnarray}
For $\min_{1\leq j \leq n} |\mb{x}_{n+1}-\mb{x}_j| \gg 1$, the
$n+1$-th particle becomes independent of particles $1,2,\ldots,n$,
and we have
\begin{equation}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_{k-1}},\mb{x}_{n+1}) =
p_1^{(\Omega)}(\mb{x}_{n+1})p_{k-1}^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2},
\ldots ,\mb{x}_{i_{k-1}})\left(1+o(1)\right),
\end{equation}
for all $k=1,2,\ldots,n$ and all sets of indices $1\leq i_1 < i_2
< \ldots < i_{k-1} \leq n$. Therefore, to leading order
\begin{eqnarray}
&& F(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n,\mb{x}_{n+1}) =
\prod_{k=1}^{n-2}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n}
p_k^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots
,\mb{x}_{i_k})^{(-1)^{n-1-k}(n-k-1)} \nonumber \\
&& \quad \times \prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots
< i_{k-1} \leq n} p_{k-1}^{(\Omega)}(\mb{x}_{i_1},\mb{x}_{i_2},
\ldots
,\mb{x}_{i_{k-1}})^{(-1)^{n-1-k}(n-k)} \nonumber \\
&& \quad \times \prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots
< i_{k-1} \leq n} p_1^{(\Omega)}(\mb{x}_{n+1})^{(-1)^{n-1-k}(n-k)}
\nonumber \\
&& = \prod_{k=1}^{n-1}\prod_{1 \leq i_1 < i_2 < \ldots < i_{k-1}
\leq n} p_1^{(\Omega)}(\mb{x}_{n+1})^{(-1)^{n-1-k}(n-k)}
\nonumber \\
&& = p_1^{(\Omega)}(\mb{x}_{n+1})^{\ds \sum_{k=1}^{n-1}
(-1)^{n-1-k}(n-k){n \choose k-1}} \nonumber \\
&& = p_1^{(\Omega)}(\mb{x}_{n+1}),
\end{eqnarray}
where we have used the combinatorial identity
\begin{equation}
\sum_{k=1}^{n-1} (-1)^{n-1-k}(n-k){n \choose k-1} = 1.
\end{equation}
We conclude that
\begin{eqnarray}
h_n^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) &=&
\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) \int_{\Omega}
p_1^{(\Omega)}(\mb{x}_{n+1})\left(1+o(1) \right) \nonumber \\
&& \times \prod_{j=1}^n
\delta^{(\Omega)}(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_{j+1},\ldots,\mb{x}_{n+1})\,d\mb{x}_{n+1}
. \nonumber
\end{eqnarray}
Therefore,
\begin{equation}
\label{eq:entire-space} \lim_{\Omega \to \mathbb{R}^3}
\delta^{(\Omega)}(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n) =
\lim_{\Omega \to \mathbb{R}^3} h_n^{(\Omega)}
(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_n)
\end{equation}
Substituting eq. (\ref{eq:delta-n}) in eq. (\ref{eq:p_n+1}) while
using the relations (\ref{eq:h_n}) and (\ref{eq:entire-space}),
and the definition (\ref{eq:g-def}) we find that
\begin{equation}
\label{eq:entropy-g} g_{n+1}(\mb{x}_1,\ldots,\mb{x}_{n+1}) =
\prod_{k=2}^{n}\prod_{1 \leq i_1 < i_2 < \ldots < i_k \leq n+1}
g_k(\mb{x}_{i_1},\mb{x}_{i_2}, \ldots ,\mb{x}_{i_k})^{(-1)^{n-k}}.
\end{equation}
We observe that for entire space systems, the probabilistic
Kirkwood closure (\ref{eq:prob-g}) agrees with the maximum entropy
closure (\ref{eq:entropy-g}) for all orders $n$.

\subsection{Closure at the highest level $n=N-1$}
\label{sub: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 (\ref{eq:p_n+1}) is exact \textbf{in a confined system
when applied } at the highest level $n=N-1$, while the
probabilistic Kirkwood closure (\ref{eq:prob-Kirkwood}) is not
exact. In other words, the maximum entropy closure relation yields
the Boltzmann distribution (\ref{eq:boltzmann}). \textbf{What
appeared at first to be} an approximation, turns out to be the
exact result. Indeed, let
\begin{eqnarray}
&& \gamma(\mb{x}_1,\mb{x}_2,\ldots,\mb{x}_{N-1}) \label{eq:gamma-N}\\
&&\quad =\left(\frac{1}{Z_N}\right)^{1/N}
\exp\left\{-\left[\frac{1}{N-2}\sum_{1\leq i < j \leq N-1}
U(\mb{x}_i,\mb{x}_j)+\frac{1}{N-1}\sum_{j=1}^{N-1}U_{\mbox{ex}}(\mb{x}_j)\right]/k_BT\right\}.\nonumber
\end{eqnarray}
Clearly,
\begin{equation}
\prod_{j=1}^N
\gamma(\mb{x}_1,\ldots,\mb{x}_{j-1},\mb{x}_j,\ldots,\mb{x}_N) =
\frac{1}{Z_N}\, e^{-U(\mb{x}_1,\ldots,\mb{x}_N)/k_BT},
\end{equation}
which is the Boltzmann distribution. The Boltzmann distribution
obviously satisfies the BBGKY equation (\ref{eq:BBGKY-n}).
Therefore, we have found a solution to BBGKY equation which
satisfies the closure relation (\ref{eq:p_n+1}). 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
\begin{equation}
\label{eq:p_N-1} p_{N-1}(\mb{x}_1,\ldots,\mb{x}_{N-1}) =
\frac{1}{Z_N} \int_{\Omega}
e^{-U(\mb{x}_1,\ldots,\mb{x}_N)/k_BT}\,d\mb{x}_N,
\end{equation}
which is both the exact and \textbf{``approximate'' }result.

The probabilistic Kirkwood approximation (\ref{eq:prob-Kirkwood})
gives however a different result. If it were to be exact, then its
resulted $N-1$-particle pdf must \textbf{be [[COMMENT: You use the
subjunctive voice correctly, but no one uses it any more at all,
if they ever did!]]}given by equation (\ref{eq:p_N-1}). Therefore,
all lower level pdf's $p_n$ ($n\leq N-1$) are (multiple) integrals
of the Boltzmann distribution. Therefore, by the closure
(\ref{eq:prob-Kirkwood}), $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
(\ref{eq:boltzmann}).

Although the maximum entropy closure is exact at the highest order
$n=N-1$, it is not exact at lower orders. The observation
\textbf{[[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$).

\section{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 \cite{Eisenberg1, Eisenberg2}.

In the general case, where $\Omega \subset \mathbb{R}^3$ is a
bounded domain, there is no analytic solution to the system
(\ref{system}). We propose to solve this system by the following
iterative scheme
\begin{enumerate}
\item Initial guess $\gamma^{(0)}(\mb{x}_1,\mb{x}_2)$. Set $i=0$.
\item Solve for $p_2^{(i)}(\mb{x}_1,\mb{x}_2)$ the
non-homogeneous linear equation
\begin{eqnarray}
&&k_BT \nabla_{\mb{x}_1}
p_2^{(i)}(\mb{x}_1,\mb{x}_2)-\mb{f}(\mb{x}_2,\mb{x}_1)p_2^{(i)}(\mb{x}_1,\mb{x}_2)-
\mb{f}_{\mbox{ex}}(\mb{x}_1)
p_2^{(i)}(\mb{x}_1,\mb{x}_2) \nonumber \\
&&\nonumber\\
&&=(N-2) \displaystyle \int_{\Omega}
\mb{f}(\mb{x}_3,\mb{x}_1) \gamma^{(i)}(\mb{x}_1,\mb{x}_2)
\gamma^{(i)}(\mb{x}_2,\mb{x}_3)\gamma^{(i)}(\mb{x}_1,\mb{x}_3)
\,d\mb{x}_3.\nonumber
\end{eqnarray}
\item Solve the non-linear system for $\gamma^{(i+1)}$
\begin{equation}
p_2^{(i)}(\mb{x}_1,\mb{x}_2) =
\gamma^{(i+1)}(\mb{x}_1,\mb{x}_2)\int_{\Omega}
\gamma^{(i+1)}(\mb{x}_2,\mb{x}_3)\gamma^{(i+1)}(\mb{x}_1,\mb{x}_3)\,d\mb{x}_3.
\end{equation}
\item $i \leftarrow i+1$. Return to step 2, until convergence.
\end{enumerate}
The analysis of the previous section indicates that a good initial
guess might be
\begin{equation}
\gamma^{(0)}(\mb{x}_1,\mb{x}_2) =
\sqrt{p_1(\mb{x}_1)p_1(\mb{x}_2)}\, g_2(\mb{x}_1,\mb{x}_2),
\end{equation}
where $g_2(\mb{x}_1,\mb{x}_2)$ is the solution to the BBGKY
equation with the Kirkwood superposition approximation in the
entire space $\mathbb{R}^3$. This solution can be found rather
easily using the inherited symmetries of the problem. For example,
it is well known that if $\mb{f}_{\mbox{ex}}=0$, then
$g_2(\mb{x}_1,\mb{x}_2) = g_2(|\mb{x}_1-\mb{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
\begin{eqnarray}
&&\nabla_{\mb{x}_1}\left[e^{(U(\mb{x}_1,\mb{x}_2)+
U_{\mbox{ex}}(\mb{x}_1))/k_BT}p_2(\mb{x}_1,\mb{x}_2)
\right]=\nonumber\\
&&\nonumber\\
&&
-\frac{N-2}{k_BT}e^{(U(\mb{x}_1,\mb{x}_2)+U_{\mbox{ex}}(\mb{x}_1))/k_BT}
\int_{\Omega}\nabla_{\mb{x}_1}U(\mb{x}_1,\mb{x}_3)
p_3(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_3.\nonumber
\end{eqnarray}
The identity $\nabla_{\mb{x}_1} \times \nabla_{\mb{x}_1}
u(\mb{x}_1) = 0$, for all $u$, imposes a solvability condition for
$\gamma$. Indeed, taking the curl of the last equation, together
with the closure (\ref{eq:closure}) result in
\begin{eqnarray*}
&&\nabla_{\mb{x}_1} \times
\left[e^{(U(\mb{x}_1,\mb{x}_2)+U_{\mbox{ex}}(\mb{x}_1))/k_BT}
\int_{\Omega}\nabla_{\mb{x}_1}U(\mb{x}_1,\mb{x}_3)
\gamma(\mb{x}_1,\mb{x}_2)\gamma(\mb{x}_1,\mb{x}_3)\gamma(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3\right]\\
&&\\
&&=0.
\end{eqnarray*}

In step 3 we solve a non-linear integral equation. We suggest
solving the non-linear equation (\ref{eq:int-equation}) by a
Newton-Raphson iterative scheme. Let
$\gamma^{(n)}(\mb{x}_1,\mb{x}_2)$ be the $n$-th iteration. Define
the operator $\Gamma^{(n)}:\Omega^2\to\Omega^2$ as follows
\begin{equation}
\Gamma^{(n)} u(\mb{x},\mb{z}) =
\int_{\Omega}\gamma^{(n)}(\mb{x},\mb{y}) u(\mb{z},\mb{y})
\,d\mb{y}.
\end{equation}
Let the operator $S:\Omega^2\to\Omega^2$ be the symmetrization
operator
\begin{equation}
Su(\mb{x},\mb{y}) = u(\mb{y},\mb{x}).
\end{equation}
The Newton-Raphson iteration scheme suggests
\begin{equation}
\gamma^{(n+1)} = \gamma^{(n)}+\Delta(\mb{x}_1,\mb{x}_2),
\end{equation}
where $\Delta(\mb{x}_1,\mb{x}_2)$ satisfies the linear integral
equation
\begin{eqnarray*}
p_2(\mb{x}_1,\mb{x}_2)-p_2^{(n)}(\mb{x}_1,\mb{x}_2) &=&
\frac{p_2^{(n)}(\mb{x}_1,\mb{x}_2)}{\gamma^{(n)}(\mb{x}_1,\mb{x}_2)}\Delta(\mb{x}_1,\mb{x}_2)
\nonumber \\
&&\\
&& + \gamma^{(n)}(\mb{x}_1,\mb{x}_2)
\Gamma^{(n)}\Delta(\mb{x}_1,\mb{x}_2)\nonumber
\\
&&\\
 &&+\gamma^{(n)}(\mb{x}_1,\mb{x}_2)
S\Gamma^{(n)}\Delta(\mb{x}_1,\mb{x}_2),
\end{eqnarray*}
where
\begin{equation}
p_2^{(n)}(\mb{x}_1,\mb{x}_2) = \gamma^{(n)}(\mb{x}_1,\mb{x}_2)
\int_{\Omega}\gamma^{(n)}(\mb{x}_1,\mb{x}_3)\gamma^{(n)}(\mb{x}_2,\mb{x}_3)\,d\mb{x}_3.
\end{equation}
We may write the iteration equivalently as
\begin{equation}
\gamma^{(n+1)}=\gamma^{(n)}+\left(\frac{p_2^{(n)}}{\gamma^{(n)}}+\gamma^{(n)}\Gamma^{(n)}
+\gamma^{(n)}S\Gamma^{(n)}\right)^{-1}\left(p_2-p_2^{(n)} \right).
\end{equation}
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
\ref{sub:N-1} 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. \textbf{We expect to use our closure to
determine the customary properties of nonequilbrium systems, e.g.,
current and conductance.}

\section{Mixtures}
The maximum entropy principle \textbf{can also be used} to find
closure relations of mixtures, both in confined domains and in the
entire space. Suppose a mixture of $S\geq 2$ species, with
$N_\alpha$ ($\alpha=1,2,\ldots,S$) particles of each specie. Let
$N=\sum_{\alpha=1}^S$ be the total number of particles of all
specie. There are $S^2$ $2$-particle pdfs,
$$p_2^{\alpha \beta}(\mb{x}_1,\mb{x}_2)\quad
\alpha,\beta\in\{1,2,\ldots,S\},$$ that exhibit the symmetry
$p_2^{\alpha \beta}(\mb{x}_1,\mb{x}_2) = p_2^{\beta
\alpha}(\mb{x}_2,\mb{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^{\alpha \beta
\gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3)$, $\alpha,\beta,\gamma \in
\{1,2,\ldots,S\}$, that bring the entropy
\begin{equation}
\label{eq:entropy-mixture} H = -\sum_{\alpha,\beta,\gamma=0}^S
\frac{N_\alpha}{N} \frac{N_\beta}{N} \frac{N_\gamma}{N}
\int_{\Omega} p_3^{\alpha \beta
\gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3) \log p_3^{\alpha \beta
\gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_1\,d\mb{x}_2\,d\mb{x}_3
\end{equation}
to maximum, with the $3S^3$ marginal constraints
\begin{eqnarray}
p_2^{\alpha \beta} (\mb{x}_1,\mb{x}_2) &=& \int_{\Omega}
p_3^{\alpha \beta \gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_3,
\\
p_2^{\alpha \gamma} (\mb{x}_1,\mb{x}_3) &=& \int_{\Omega}
p_3^{\alpha \beta \gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_2,
\nonumber \\
p_2^{\beta \gamma} (\mb{x}_2,\mb{x}_3) &=& \int_{\Omega}
p_3^{\alpha \beta \gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3)\,d\mb{x}_1.
\nonumber
\end{eqnarray}
This variational problem is solved using the Euler-Lagrange
formulation similar to the derivation done in Section
\ref{sec:entropy}.

In the case of an entire space system, it can be shown in a
similar way to the method of subsections \ref{sub:entire} and
\ref{sub:entire2}, that the mixture entropy closure coincides with
the probabilistic Kirkwood closure. For example, they both suggest
\textbf{XXXX} that the triplets correlation functions are related
to the pair correlation function via
\begin{equation}
g_3^{\alpha \beta \gamma}(\mb{x}_1,\mb{x}_2,\mb{x}_3) =
g_2^{\alpha \beta}(\mb{x}_1,\mb{x}_2)g_2^{\alpha
\gamma}(\mb{x}_1,\mb{x}_3)g_2^{\beta \gamma}(\mb{x}_2,\mb{x}_3),
\end{equation}
for $\alpha,\beta,\gamma \in \{1,2,\ldots,S \}$. Similarly, one
can obtain \textbf{XXXX} high order closures.

In confined systems, the Euler-Lagrange formulation leads to
integral equation of the form (\ref{eq:int-equation}). Note that
since the entropy (\ref{eq:entropy-mixture}) depends on the
particle fraction $\ds \frac{N_\alpha}{N}$, we expect the
resulting confined system pair correlation \textbf{will also}
depend on the particle fraction.

\section{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. \textbf{We expect to use the closure to make estimates
of the nonequilibrium properties that naturally arise in a system
with spatially inomogeneous boundary conditions.} \textbf{The
implementation and results of using the closure will be the
subject of a separate paper.}

\section{Acknowledgments}
Special thanks to Bob Eisenberg, Douglas Henderson and Roland Roth
for reading the manuscript and for their useful comments.

\begin{thebibliography}{99}
\bibitem{BRR}
R. S. Berry, S. Rice, J. Ross, {\em Physical Chemistry}, 2nd
edition. Oxford University Press, NY 2000.

\bibitem{Barthel}J.M.G. Barthel, A.H. Baumgartel (Editor), H. Krienke, {\em
Physical Chemistry of Electrolyte Solutions: Modern Aspects},
Steinkopf, Dietrich Pub. 1998.

\bibitem{Davis}
H. T. Davis, {\em Statistical Mechanics of Phases, Interfaces, and
Thin Films}, (Advances in Interfacial Engineering Series), John
Wiley \& Sons, NY, 1995.

\bibitem{Henderson}
D. Henderson (Editor), {\em Fundamentals of inhomogeneous fluids},
Marcel Dekker, New York, 1992

\bibitem{Lehmann1}
E. L. Lehmann, {\em Testing Statistical Hypotheses}, 2nd Ed.,
Springer Verlag, 1997.

\bibitem{Lehmann2}
E. L. Lehmann, G. Casella, {\em Theory of Point Estimation}, 2nd
Ed., Springer Verlag, 1998.

\bibitem{Eisenberg1}
R.S. Eisenberg, ``Channels as enzymes." {\em J. Memb. Biol.}
{\bf115}, pp.1–12 (1990).

\bibitem{Eisenberg2}
R.S. Eisenberg, ``Ion Channels as Devices", {\em Journal of
Computational Electronics} {\bf2}, pp.245-249 (2003).

\end{thebibliography}
\end{document}
