# On the equivalence of self-consistent equations for nonuniform liquids: a
unified description of the various modifications

Hiroshi Frusawa111e-mail<EMAIL_ADDRESS>Laboratory of
Statistical Physics, Kochi University of Technology, Tosa-Yamada, Kochi
782-8502, Japan.

###### Abstract

A variety of self-consistent (SC) equations have been proposed for non-uniform
states of liquid particles under external fields, including adsorbed states at
solid substrates and confined states in pores. External fields represent not
only confining geometries but also fixed solutes. We consider SC equations
ranging from the modified Poisson-Boltzmann equations for the Coulomb
potential to the hydrostatic linear response equation for the equilibrium
density distribution of Lennard-Jones fluids. Here, we present a unified
equation that explains the apparent diversity of previous forms and proves the
equivalence of various SC equations. This unified description of SC equations
is obtained from a hybrid method combining the conventional density functional
theory and statistical field theory. The Gaussian approximation of density
fluctuations around a mean-field distribution is performed based on the
developed hybrid framework, allowing us to derive a novel form of the grand-
potential density functional that provides the unified SC equation for
equilibrium density.

##  1 Introduction

Many studies on dense uniform liquids have demonstrated the dominant role of
short-ranged harshly repulsive intermolecular forces in determining structural
correlations, which provides the insight that attractive intermolecular forces
essentially cancel due to the symmetric configurations of molecules in uniform
liquids [1-3]. In contrast, translational invariance of molecular arrangements
is broken in non-uniform liquids [4-12]. Accordingly, both attractive and
repulsive forces can, for instance, significantly influence the structure of a
Lennard-Jones (LJ) fluid next to a hard wall because the vector sums of the
long-range attractive forces do not cancel each other.

There are several situations in which such non-uniform states of liquid
molecules emerge. Inhomogeneities occur [4-12]: (i) at interfaces between
liquid-gas, liquid-liquid, and liquid-crystal phases; (ii) in the adsorption
of liquids at solid substrates or walls; (iii) for fluids in confining
geometries, such as slits and pores; and (iv) in the sedimentation equilibrium
of colloidal particles under gravity. Theoretically, a non-uniform state can
be created by adding an interaction potential term associated with an external
field. For example, an external field may represent the presence of walls,
slits, or pores [4-14]. We can also find adequate fields that are not only
responsible for shifts and rounding of bulk phase transitions, but also induce
wetting or drying transitions and related phenomena such as capillary
condensation. Furthermore, an external field is able to describe the
interactions of fixed solutes with solvent liquids that significantly perturb
the solvent density around solutes [4-14].

Based on this theoretical approach, we can evaluate the solvation free energy
per solute and effective interaction potential between solutes [1, 4-14].
Specifically, the solvation thermodynamics in charged and polar systems
underlies important processes in biology, including protein folding and ligand
binding [13, 14]. This indicates that the Coulomb potential as well as the LJ
potential is significant when considering the long-range part of external
field. In fact, among the theoretical issues related to external Coulomb
fields are the formulations of the following diverse phenomena: (i) hydration
thermodynamics of hydrophobic and ionic solutes [13]; (ii) solvation of highly
charged colloids [14]; (iii) unconventional phenomena such as Coulomb
attraction between like-charged surfaces and overcharging of solvated objects
due to the accumulation of counterions in solvents (see [15] for a review);
and (iv) novel properties observed around the interfaces between electrodes
and room temperature ionic liquids, including the crossover from overscreening
to crowding with an increase in applied voltage and the long-range
underscreening detected by surface force measurements [16-19].

All of these phenomena associated with Coulomb interactions stem from the
strong correlations between charged particles, which requires going beyond the
Poisson-Boltzmann (PB) equation, a typical mean-field equation for calculating
the Coulomb potential induced by charge-charge interactions [15]. A variety of
modified PB equations [14-34], proposed so far, represent a subset of the
self-consistent (SC) equations considering intermolecular correlations. Our
concern here is with the following two types of the modified PB equations. One
type is derived from the charge smearing model [14, 20-24, 35-40] and the
other adds higher-order gradient terms to the original PB equation [15-18,
25-30]. Both types of modifications have coexisted as complementary to each
other, and the direct connections of these different types of the modified PB
equations with the liquid-state theory, or the vast knowledge on correlation
functions [1, 41], remain to be clarified.

Meanwhile, for dense non-Coulombic liquids under the influence of external
fields, there is another set of SC equations for predicting the non-uniform
distribution of average density $\rho_{\mathrm{eq}}({\bf r})$ in equilibrium
[4, 5, 14, 41-47]. This type of SC theory can be traced back to the van der
Waals interface theory [5, 42]. The essential ingredient of the classical van
der Waals theory is that gas-liquid coexistence can be attributed to
intermolecular competition between short-range repulsions and long-range
attractions, thereby providing a qualitatively accurate description of slowly
varying interfaces. A recent SC equation modified the van der Waals equation
to quantitatively describe local interfacial structures on the molecular scale
[5].

The goal of this paper is to show that there is an underlying framework to
yield various SC equations, despite their apparently diverse forms. To this
end, we demonstrate the equivalence of the aforementioned SC equations for
both the equilibrium density and Coulomb potential, using a hybrid method [48,
49] that combines the density functional theory (DFT) [6-12] and statistical
field theory [15, 23, 31-34, 36, 50-52].

The DFT incorporates as input microscopic correlations due to the inter-
particle interactions [6-12], thus providing a fruitful approach to the
descriptions of either microscopic structures in inhomogeneous fluids or
states of various matters, including not only hard spheres, but also soft-core
particles, mixtures and anisotropic particles. The statistical field theory
[15, 23, 31-34, 36, 50-52], on the other hand, can incorporate fluctuations,
especially important in the vicinity of the critical point, into the
formulation in systematic ways using the Gaussian approximation, many-loop
expansion, variation method, and the renormalization group theory. Benefiting
from both of the DFT and statistical field theory, the hybrid theory [48, 49]
allows us to incorporate knowledge from the liquid-state theory [1, 41] into a
field-theoretic treatment of density fluctuations in inhomogeneous fluids.
Accordingly, a unified equation for $\rho_{\mathrm{eq}}({\bf r})$ can be
obtained from a novel form of the grand-potential density functional. The
unified equation corresponds to an extension from the equation based on the
Percus’ test particle method [47, 53] for uniform liquids to that for non-
uniform liquids. It will also be demonstrated that the single unified equation
for $\rho_{\mathrm{eq}}({\bf r})$ can yield a variety of existing SC equations
[4, 5, 14-34, 41-47].

In the remainder of this paper, we first present a review of previous SC
equations [4, 5, 14-34, 41-47] and our main results related to the unification
discussed above: sections 2 and 3 provide the formal background and essential
results for demonstrating the equivalence between previous SC equations. In
Sections 4 and 5, comparisons are made between previous SC equations and the
obtained forms derived from a single unified equation for
$\rho_{\mathrm{eq}}({\bf r})$: we investigate the SC equations for non-
Coulombic liquids [4, 5, 14, 41-47] in Section 4, whereas Section 5 presents a
unified equation for $\rho_{\mathrm{eq}}({\bf r})$ expressed by the Coulomb
potential for verifying the equivalence of modified PB equations [14-34].
After demonstrating the relevance of our hybrid framework [48, 49] for SC
equations, Section 6 outlines the derivation scheme for the grand-potential
density functional, a key functional to obtain the unified SC equation.
Summary and conclusions are given in Section 7. In the appendices, full
details of the formulations regarding not only the present results but also
the previous modified PB equations can be found.

##  2 Formal background

###  2.1 Target systems in terms of the potential separation [1-5, 14, 20, 21,
46, 47]

In non-uniform liquids, there are two types of interaction potentials: one is
the interaction potential between liquid particles, which produces
intermolecular forces in both uniform and non-uniform liquids, whereas the
other is a potential that creates a known external field felt by constituent
particles, irrespective of liquid states, due to the presence of fixed
solutes, confining geometries, and so on. The former interaction potential
$v({\bf r})$, as a function of the intermolecular distance $|{\bf r}|\equiv
r$, includes the hard-sphere (HS) potential, LJ potential, and Coulomb
potential. The latter external potential $J({\bf r})$, which is uniquely
determined by the location ${\bf r}$ of particles, is caused not only by the
existence of a hard wall, but also by long-range repulsive or attractive
particle-wall interactions, such as Coulomb and LJ interactions. It should be
noted that not only $v({\bf r})$ and $J({\bf r})$ but also the other energetic
quantities, appearing in this paper, are defined in the $k_{B}T$-unit.

Let us introduce the liquid systems through explaining the potential splitting
into two parts:

$\displaystyle v({\bf r})$ $\displaystyle=v_{0}({\bf r})+v_{1}({\bf r}),$ (1)
$\displaystyle J({\bf r})$ $\displaystyle=J_{0}({\bf r})+J_{1}({\bf r}).$ (2)

In general, $v_{0}$ and $J_{0}$ denote harshly repulsive parts and $v_{1}$ and
$J_{1}$ slowly varying parts. In a dense and uniform fluid, the subscript ”0”
indicates that the potentials are those in a reference system that has been
found sufficient for determining the essential structural properties [1-3]. As
the first approximation in treating $J_{0}({\bf r})$, we have estimated from
$J_{0}$ an effective radius $R$ of fixed object, such as a fixed solute, that
excludes constituent liquid particles (see eq. (16) for the definition of the
effective radius $R$) [5, 42].

We investigate the SC equations for HS fluids, attractive particle systems,
and point-charge systems as listed in Table 1. Correspondingly, we give
concrete forms of the interaction potential splitting ($v=v_{0}+v_{1}$) for
the HS fluids, LJ ones and point-charge systems. The so-called WCA separation,
proposed by Weeks, Chandler and Andersen, applies to attractive particle
systems including the LJ fluids [1-3].

First, in the HS fluids, $v_{0}$ is identified with the full potential
$v_{\mathrm{HS}}$ of HS interactions (i.e., $v_{0}=v_{\mathrm{HS}}$), and it
follows that $v_{1}\equiv 0$. Next, the WCA separation for the LJ fluids
creates the WCA potential $v_{\mathrm{WCA}}({\bf r})$ arising from the LJ
potential $v_{LJ}$ truncated at the minimum as follows:

$\displaystyle v_{0}({\bf r})=v_{\mathrm{WCA}}({\bf
r})\equiv\left\\{\begin{array}[]{l}v_{LJ}({\bf
r})+\epsilon_{LJ}\quad\mathrm{if}\>r\leq r_{m}\\\ \\\
0\quad\qquad\qquad\>\>\,\mathrm{if}\>r>r_{m},\end{array}\right.$ (6)

and

$\displaystyle v_{1}({\bf
r})=\left\\{\begin{array}[]{l}-\epsilon_{LJ}\quad\qquad\>\>\mathrm{if}\>r\leq
r_{m}\\\ \\\ v_{LJ}({\bf
r})\quad\qquad\mathrm{if}\>r>r_{m}.\end{array}\right.$ (10)

In eqs. (6) and (10), $r_{m}$ and $\epsilon_{LJ}$ denote the separation and
depth of potential minimum, respectively. Equation (6) is a typical form of
the WCA potential in attractive particle systems.

Lastly, the point-charge systems have the Coulomb interaction potential that
is always repulsive:

$v({\bf r})=v_{\mathrm{el}}({\bf r})\equiv\frac{z^{2}l_{B}}{r},$ (11)

where we consider particles carrying charges of $+ze$ and the Bjerrum length
$l_{B}$ is defined as the distance at which two unit charges have interaction
energy $k_{B}T$. Despite the long-range nature, there is no characteristic
length for the potential splitting; nevertheless, previous studies have
proposed the potential splitting [14, 20, 21] such that

$\displaystyle v_{\mathrm{el}}({\bf r})$ $\displaystyle=v_{0}({\bf
r})+v_{1}({\bf r}),$ (12) $\displaystyle v_{0}({\bf r})$
$\displaystyle=z^{2}l_{B}\frac{\mathrm{erfc}(r/\xi)}{r},$ (13) $\displaystyle
v_{1}({\bf r})$ $\displaystyle=z^{2}l_{B}\frac{\mathrm{erf}(r/\xi)}{r},$ (14)

which has been demonstrated to be useful in finding the non-uniform Coulomb
potential by adjusting the characteristic length $\xi$. General forms of
$v_{1}({\bf r})$ in point-charge systems are given in Appendix A.2 where it is
clarified that this type of potential splitting in the Coulomb potential is
equivalent to the introduction of charge smearing model (or Onsager charge
smearing) for strongly-coupled Coulomb systems [14, 20-24, 35-40].

Table 1: A list of previous self-consistent equations derived from a unified
equation (39) [EL$\,=\,$Euler-Lagrange; RY$\,=\,$Ramakrishnan-Yussouff;
HLR$\,=\,$hydrostatic linear response; PB$\,=\,$Poisson-Boltzmann].
Constituents |  | Self-consistently
---
determined variables
Original equation | Modified equation
Hard spheres | Equilibrium density |  | EL equation of
---
the RY functional
| HLR equation
---
(eq. (23))
Attractive particles |  | Equilibrium density
---
Reference density
| ditto
---
—
| ditto
---
Mean-field equation
|  |  |  (eq. (15))
Point charges | Coulomb potential | PB equation |  | Finite-spread PB equation
---
|  |  |  (eqs. (24), (121) and (122))
Point charges | Coulomb potential | PB equation |  | Higher-order PB equation
---
|  |  |  (eqs. (26) and (124))

###  2.2 A variety of previous self-consistent equations

We outline previous SC equations [4, 5, 14-34, 42-47] in the above three
systems, which this paper aims to derive in a unified manner. Table 1 lists
the previous SC equations for these systems. In Table 1, the SC equations are
classified into two groups: the former group can determine the equilibrium
density for the HS and LJ fluids [4, 5, 14, 42-47], whereas the latter group
modifies the PB equation to improve the mean-field solution of Coulomb
potential [14-34].

Conventional mean-field equation for the reference density
$\rho_{\mathrm{ref}}({\bf r})$ [4, 5, 14, 42-47].— The conventional mean-field
equation is used when considering a slowly varying interaction potential
$v_{1}({\bf r})$ in a non-uniform liquid with the external potential $J({\bf
r})$ applied. The solution of the mean-field equation has been referred to as
the reference density $\rho_{\mathrm{ref}}({\bf r})$. In the outside of the
excluded region with the effective radius of $R$ [5, 42], which is defined
below in eq. (16), $\rho_{\mathrm{ref}}({\bf r})$ satisfies

$\displaystyle\rho_{\mathrm{ref}}({\bf r})=e^{\mu-J_{1}({\bf r})-\int d{\bf
r}^{\prime}v_{1}({\bf r}-{\bf r}^{\prime})\rho_{\mathrm{ref}}({\bf
r}^{\prime})},$ (15)

where $\mu$ denotes the chemical potential in the $k_{B}T$-unit. It is found
from eq. (15) that the spatial dependence of $\rho_{\mathrm{ref}}({\bf r})$
for HS fluids arises only from $J_{1}({\bf r})$ because of $v_{1}({\bf
r})\equiv 0$.

The usual form relating $J_{0}({\bf r})$ to $R$ is

$R=\int_{0}^{\infty}dr\,\left\\{1-e^{-J_{0}({\bf r})}\right\\},$ (16)

which is a measure of the region $w$ that excludes the liquid particles due to
the harshly repulsive field $J_{0}({\bf r})$. Considering that the reference
density $\rho_{\mathrm{ref}}({\bf r})$ necessarily vanishes inside the region
$w$ due to the presence of $J_{0}({\bf r})$, the reference density is
expressed by

$\displaystyle\rho_{\mathrm{ref}}({\bf
r})=\left\\{\begin{array}[]{l}0;\qquad{\bf r}\in w,\\\ \\\
\rho_{\mathrm{ref}}^{0}\,e^{-\Delta\phi({\bf r})};\>{\bf r}\notin
w,\end{array}\right.$ (20)

where the reference density
$\rho_{\mathrm{ref}}^{0}\equiv\rho_{\mathrm{ref}}({\bf r}_{0})$ at a specified
position of ${\bf r}_{0}$ is treated separately, and the shifted chemical
potential $\mu_{\mathrm{ref}}$, defined by

$\displaystyle\rho_{\mathrm{ref}}^{0}$ $\displaystyle\equiv
e^{\mu_{\mathrm{ref}}}=e^{\mu-\phi({\bf r}_{0})},$ (21)
$\displaystyle\phi({\bf r}_{0})$ $\displaystyle=J_{1}({\bf r}_{0})+\int d{\bf
r}^{\prime}v_{1}({\bf r}_{0}-{\bf r}^{\prime})\rho_{\mathrm{ref}}({\bf
r}^{\prime}),$ (22)

is utilized for obtaining the equilibrium density, according to the
hydrostatic approximation [4, 5, 43-45]. Accordingly, the reference density at
other locations is given by eq. (20) using $\rho_{\mathrm{ref}}^{0}$ and the
external potential difference $\Delta\phi({\bf r})\equiv\phi({\bf
r})-\phi({\bf r}_{0})$ from the external field given by eq. (22).

Previous studies based on the hydrostatic approximation have demonstrated that
the mean-field like treatment works well for the purpose of obtaining the
equilibrium density distribution of a dense non-uniform liquid with the help
of the above reference density $\rho_{\mathrm{ref}}^{0}$.

The SC equation for the equilibrium density $\rho_{\mathrm{eq}}({\bf r})$ [4,
5, 42-47].— In attractive particle systems as well as in HS fluids, the
equilibrium density $\rho_{\mathrm{eq}}({\bf r})$ is obtained from evaluating
the deviation from the reference density $\rho_{\mathrm{ref}}^{0}$ at a
specified position of ${\bf r}_{0}$. The evaluation is based on the
approximation of the hydrostatic linear response (HLR) [4, 5, 14, 42-47],
which reads for the equilibrium density at ${\bf r}_{0}$

$\displaystyle\rho_{\mathrm{eq}}({\bf
r}_{0})=\rho_{\mathrm{ref}}^{0}\left[1+\int d{\bf
r}^{\prime}\,c_{\mathrm{0}}({\bf r}_{0}-{\bf
r}^{\prime})\left\\{\rho_{\mathrm{eq}}({\bf
r}^{\prime})-\rho_{\mathrm{ref}}^{0}\right\\}\right],$ (23)

using the direct correlation function (DCF) $c_{0}({\bf r})$ for a uniform
system of liquid particles interacting via the harshly repulsive potential
$v_{0}({\bf r})$. The DCF depends on the density; in eq. (23), however, we use
the DCF at a uniform density of $\rho^{0}_{\mathrm{ref}}$. The same equation
(23) applies to both of attractive particle systems and HS fluids because the
HLR equation [4, 5, 14, 42-47] is formulated for $v_{0}({\bf r})$,
irrespective of the systems considered. Equation (23) is reduced to the
closure in the hypernetted chain (HNC) approximation [1] when the density
difference $\rho_{\mathrm{eq}}({\bf r})-\rho_{\mathrm{ref}}^{0}$ is related to
the total correlation function (TCF) $h_{0}({\bf r})$ in the
$v_{0}$–interacting systems: $\rho_{\mathrm{eq}}({\bf
r})-\rho_{\mathrm{ref}}^{0}=\rho_{\mathrm{ref}}^{0}h_{0}({\bf r}-{\bf
r}_{0})$, according to the Percus’ test particle method [47, 53]. In eq. (23),
we take into account the non-uniformity by incorporating the contribution of
slowly varying potential $v_{1}({\bf r})$ into the mean-field solution
$\rho_{\mathrm{ref}}^{0}$ as given by eq. (15).

Modified PB equations for the Coulomb potential $\Psi({\bf r})$ [14-30].—
Lastly, we mention two types of modified PB equations (see Appendix A for the
details). Both modifications are relevant to describe the charge-charge
correlations between charged particles in the one component plasma (OCP) at
strong coupling, concentrated electrolytes, and room temperature ionic
liquids. It has also been found that these two modifications explain the
deviation from the PB solution for counterion distribution over a wide range
of from the intermediate to strong coupling [14-30].

One type is the finite-spread PB equation [14, 20-24] based on the charge
smearing model. In the uniform OCP, a smeared charge distribution is simply
characterized by the mean distance $a$ between adjacent charges and has been
found useful to represent the strong repulsive correlations between the
charges; for non-uniform systems, $a$ is supposed to represent the mean
separation in the ground state (see Appendix A.1 for the details of the
definition). In the local molecular field theory [14, 20, 21], the Gaussian
charge smearing model provides

$\displaystyle\nabla^{2}\Psi({\bf r})$ $\displaystyle=-4\pi
l_{B}\overline{\overline{\,e^{\mu-J_{1}({\bf r})-\Psi({\bf r})}}}$ (24)
$\displaystyle\overline{\overline{\,\mathcal{O}({\bf r})}}$
$\displaystyle=\frac{1}{\\{\pi(a/m)^{2}\\}^{3/2}}\int d{\bf
r}^{\prime}\exp\left\\{-\left(\frac{m|{\bf r}-{\bf
r}^{\prime}|}{a}\right)^{2}\right\\}\mathcal{O}({\bf r}^{\prime}),$ (25)

where a constant $m$ is an adjusting parameter (see also Appendix A.3); for
instance, $m=1/1.4$ for the uniform OCP and $m=1/0.6$ for the counterion
system in the presence of one charged wall, according to previous results at
strong coupling [14, 20, 21].

We would also like to investigate the other type of higher-order PB equations
[15-18, 25-30]. This type of modification adds higher-order gradient terms to
the PB equation such that [26]

$\displaystyle\nabla^{2}\left\\{1-\frac{a^{2}}{4}\nabla^{2}+\frac{a^{4}}{16}\nabla^{4}\right\\}\Psi({{\bf
r}})=-4\pi l_{B}z^{2}e^{\mu-J_{1}({\bf r})-\Psi({\bf r})};$ (26)

see Appendix A.4 for the details of higher-order PB equations [15-18, 25-30].
The other type of higher-order PB equation, which is truncated at the
biharmonic $\nabla^{2}\nabla^{2}$–term, has been referred to as the Bazant-
Storey-Kornyshev (BSK) equation [16-18] in binary systems such as concentrated
electrolytes and room temperature ionic liquids.

##  3 Central results for proving the equivalence of previous self-consistent
equations

###  3.1 Grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$

As detailed in Section 6, the density-functional integral representation of
the grand potential is a powerful tool to evaluate fluctuating density fields
around the reference density $\rho_{\mathrm{ref}}$ with the help of the
conventional DFT [6-12]. The hybrid method [48, 49] combining the DFT [6-12]
and statistical field theory [15, 23, 31-34, 36, 50-52] allows us to obtain
the grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ into which the additional
contribution due to density fluctuations is incorporated. The resulting form
is

$\displaystyle\widetilde{\Omega}[\rho_{\mathrm{ref}}]$
$\displaystyle=\widetilde{\Omega_{0}}[\rho_{\mathrm{ref}}]+U_{1}[\rho_{\mathrm{ref}}],$
(27) $\displaystyle\widetilde{\Omega_{0}}[\rho_{\mathrm{ref}}]$
$\displaystyle=\mathcal{F}_{0}^{\mathrm{ex}}[\rho^{0}_{\mathrm{ref}}]-\frac{1}{2}\iint
d{\bf r}d{\bf r}^{\prime}\Delta\rho_{\mathrm{ref}}({\bf r})h_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$
$\displaystyle\qquad\qquad+\int d{\bf r}\left[\rho_{\mathrm{ref}}({\bf
r})\ln\rho_{\mathrm{ref}}({\bf r})-\rho_{\mathrm{ref}}({\bf
r})-\rho_{\mathrm{ref}}({\bf r})\mu\right],$ (28) $\displaystyle
U_{1}[\rho_{\mathrm{ref}}]$ $\displaystyle=\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\rho_{\mathrm{ref}}({\bf r})v_{\mathrm{att}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})+\int d{\bf
r}\,\rho_{\mathrm{ref}}({\bf r})J_{1}({\bf r}).$ (29)

In eq. (28), $\mathcal{F}_{0}^{\mathrm{ex}}[\rho^{0}_{\mathrm{ref}}]$
represents the excess part of the intrinsic Helmholtz free energy,
$\Delta\rho_{\mathrm{ref}}({\bf r})=\rho_{\mathrm{ref}}({\bf
r})-\rho_{\mathrm{ref}}^{0}$ and the reference TCF $h_{0}({\bf r})$ is related
to the reference DCF $c_{0}({\bf r})$ through the Ornstein-Zernike equation
given below. The subscript ”0” in the TCF have double meanings as well as that
in the DCF: the TCF is not only the function at a uniform density, which is
equal to a reference density $\rho_{\mathrm{ref}}^{0}=\rho_{\mathrm{ref}}({\bf
r}_{0})$ at ${\bf r}_{0}$, but is also the function in repulsive
$v_{0}$–interaction systems. The subscript ”1” on the right hand side (rhs) of
eq. (29) is altered to the subscript ”att”, clarifying that the potential
splitting given by eq. (1) applies only to the attractive particle systems in
our results: $v=v_{0}+v_{\mathrm{att}}$.

Equation (28) indicates that the repulsive part of the grand potential
$\widetilde{\Omega_{0}}$ is composed of two parts: the interaction energy term
and the ideal-gas contribution. The second term on the rhs of eq. (28),
representing the interaction energy, is expressed using the TCF instead of the
DCF, other than the conventional DFT [6-12]. In eq. (29), on the other hand,
the interaction energy due to the slowly varying attractive part
$v_{\mathrm{att}}({\bf r})$ of interaction potential is evaluated in the mean-
field approximation. Such treatment is apparently a crude estimation but has
been found to be relevant for evaluating the long-range contribution to
$\widetilde{\Omega}$ in non-uniform liquids [14, 46, 47]; actually, it is
quite difficult to find density-density correlations by explicitly considering
the non-uniformity.

###  3.2 Reference and equilibrium densities

As seen from the expression (29), we regard only the contribution due to the
long-range attractive interactions as the additional part of the grand-
potential density functional $\widetilde{\Omega}[\rho_{\mathrm{ref}}]$.
Correspondingly, in this study, $v_{0}({\bf r})$ reads for the following three
systems:

$\displaystyle v_{0}({\bf r})=\left\\{\begin{array}[]{l}v_{\mathrm{WCA}}({\bf
r})\quad(\mathrm{Attractive\,particles})\\\ \\\ v_{\mathrm{HS}}({\bf
r})\quad\quad(\mathrm{Hard\,spheres})\\\ \\\ v_{\mathrm{el}}({\bf
r})\quad\quad\>(\mathrm{Point\,charges}).\end{array}\right.$ (35)

In what follows, the same rule of the subscript ”0” as above applies to the
subscript ”0” in the TCF ($h_{0}({\bf r})$) and DCF ($c_{0}({\bf r})$).
Equation (15) with the above expression (35) implies that the spatial
dependence of $\rho_{\mathrm{ref}}({\bf r})$ is determined solely by the
external potential $J_{1}({\bf r})$ for both of the HS fluids and point-charge
systems because of $v_{\mathrm{att}}({\bf r})\equiv 0$.

Once the external potential $J_{1}$ is treated as a variable field
$\mathcal{J}$, we can differentiate the grand potential with respect to
$\mathcal{J}$, thereby generating the equilibrium density
$\rho_{\mathrm{eq}}({\bf r})$ with the help of a variable reference density
$\rho_{\mathrm{ref}}^{\mathcal{J}}({\bf r})\equiv e^{\mu-\mathcal{J}({\bf
r})-\int d{\bf r}^{\prime}v_{\mathrm{att}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})}$ as follows:

$\displaystyle\rho_{\mathrm{eq}}({\bf r})$
$\displaystyle=\left.\frac{\delta\widetilde{\Omega}[\rho^{\mathcal{J}}_{\mathrm{ref}}]}{\delta\mathcal{J}({\bf
r})}\right|_{\mathcal{J}=J_{1}}$
$\displaystyle=\left.\frac{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}({\bf
r})}{\delta\mathcal{J}({\bf
r})}\right|_{\mathcal{J}=J_{1}}\left.\frac{\delta\widetilde{\Omega}[\rho^{\mathcal{J}}_{\mathrm{ref}}]}{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}({\bf
r})}\right|_{\rho^{\mathcal{J}}_{\mathrm{ref}}=\rho_{\mathrm{ref}}}$
$\displaystyle=-\rho_{\mathrm{ref}}({\bf
r})\left.\frac{\delta\widetilde{\Omega}[\rho^{\mathcal{J}}_{\mathrm{ref}}]}{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}({\bf
r})}\right|_{\rho^{\mathcal{J}}_{\mathrm{ref}}=\rho_{\mathrm{ref}}},$ (36)

where use has been made of the following relation in the last line of the
above equation:

$\displaystyle\left.\frac{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}({\bf
r})}{\delta\mathcal{J}({\bf r})}\right|_{\mathcal{J}=J_{1}}$
$\displaystyle=\left.\frac{\delta}{\delta\mathcal{J}}\left\\{e^{\mu-\mathcal{J}({\bf
r})-\int d{\bf r}^{\prime}v_{\mathrm{att}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf
r}^{\prime})}\right\\}\right|_{\mathcal{J}=J_{1}}$ $\displaystyle=-e^{\mu-
J_{1}({\bf r})-\int d{\bf r}^{\prime}v_{\mathrm{att}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})}.$ (37)

It also follows from eqs. (27) to (29) that the functional differentiation of
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ with respect to
$\rho_{\mathrm{ref}}({\bf r})$ yields

$\displaystyle\left.\frac{\delta\widetilde{\Omega}[\rho^{\mathcal{J}}_{\mathrm{ref}}]}{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}({\bf
r})}\right|_{\rho^{\mathcal{J}}_{\mathrm{ref}}=\rho_{\mathrm{ref}}}$
$\displaystyle=-\int d{\bf r}^{\prime}\,h_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$
$\displaystyle\qquad\underline{+\int d{\bf r}^{\prime}\,v_{\mathrm{att}}({\bf
r}-{\bf r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})+J_{1}({\bf
r})-\mu+\ln\rho_{\mathrm{ref}}({\bf r})}+\rho_{\mathrm{ref}}({\bf
r})\left.\frac{\delta\mathcal{J}}{\delta\rho^{\mathcal{J}}_{\mathrm{ref}}}\right|_{\rho^{\mathcal{J}}_{\mathrm{ref}}=\rho_{\mathrm{ref}}}$
$\displaystyle=-\int d{\bf r}^{\prime}\,h_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})-1,$ (38)

where the resulting form seen in the last line has been obtained from both
using eq. (37) and canceling the underlined terms due to the relation (15).
Combination of eqs. (36) and (38) provides the single unified equation for
$\rho_{\mathrm{eq}}({\bf r})$:

$\displaystyle\rho_{\mathrm{eq}}({\bf r})=\rho_{\mathrm{ref}}({\bf
r})\left\\{1+\int d{\bf r}^{\prime}\,h_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})\right\\},$ (39)

which is our main result. Plugging the Ornstein-Zernike equation into eq.
(39), the HLR equation (23) [4, 5, 14, 42-47] will be obtained in the next
section.

##  4 Comparison with the HLR equation (23) for the equilibrium density
$\rho_{\mathrm{eq}}({\bf r})$

We use the Ornstein-Zernike equation [1],

$\displaystyle h_{0}({\bf r}-{\bf r}^{\prime})=c_{0}({\bf r}-{\bf
r}^{\prime})+\int d{\bf r}"\,c_{0}({\bf r}-{\bf
r}")\rho_{\mathrm{ref}}^{0}\,h_{0}({\bf r}"-{\bf r}^{\prime}),$ (40)

for the reference density $\rho_{\mathrm{ref}}^{0}$ given by eq. (21) because
both of the DCF and TCF, $c_{0}({\bf r})$ and $h_{0}({\bf r})$, are defined
for uniform liquids with the density of $\rho_{\mathrm{ref}}^{0}$. Plugging
eq. (40) into eq. (39) at ${\bf r}_{0}$, we have

$\displaystyle\frac{\rho_{\mathrm{eq}}({\bf
r}_{0})-\rho_{\mathrm{ref}}^{0}}{\rho_{\mathrm{ref}}^{0}}=\int d{\bf
r}^{\prime}\,h_{0}({\bf r}_{0}-{\bf r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf
r}^{\prime})$ $\displaystyle=\int d{\bf r}^{\prime}\,c_{0}({\bf r}_{0}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})+\int d{\bf
r}"\,c_{0}({\bf r}_{0}-{\bf r}")\rho_{\mathrm{ref}}^{0}\int d{\bf
r}^{\prime}h_{0}({\bf r}"-{\bf r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf
r}^{\prime}).$ (41)

The main result (39) is again inserted into $\rho_{\mathrm{ref}}^{0}\int d{\bf
r}^{\prime}h_{0}({\bf r}"-{\bf r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf
r}^{\prime})$ of the above expression with the relation
$\rho_{\mathrm{ref}}({\bf
r}")=\rho_{\mathrm{ref}}^{0}+\Delta\rho_{\mathrm{ref}}({\bf r}")$, yielding

$\displaystyle\rho_{\mathrm{ref}}^{0}\int d{\bf r}^{\prime}h_{0}({\bf r}"-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$
$\displaystyle=\rho_{\mathrm{ref}}^{0}\left\\{\frac{\rho_{\mathrm{eq}}({\bf
r}")-\rho_{\mathrm{ref}}({\bf
r}")}{\rho_{\mathrm{ref}}^{0}+\Delta\rho_{\mathrm{ref}}({\bf r}")}\right\\}$
$\displaystyle\approx\rho_{\mathrm{eq}}({\bf r}")-\rho_{\mathrm{ref}}({\bf
r}")-\left\\{\rho_{\mathrm{eq}}({\bf r}")-\rho_{\mathrm{ref}}({\bf
r}")\right\\}\frac{\Delta\rho_{\mathrm{ref}}({\bf
r}")}{\rho_{\mathrm{ref}}^{0}}$ $\displaystyle\approx\rho_{\mathrm{eq}}({\bf
r}")-\rho_{\mathrm{ref}}({\bf r}").$ (42)

It will be seen in Section 6 and Appendix B.2 that the above approximation
used in eq. (42) is consistent with the approximation in deriving the grand-
potential density functional $\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ given
by eqs. (27) to (29). Thus, we find that
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ leads to the HLR equation (23) via
eq. (39) [4, 5, 14, 42-47]: it follows from eq. (42) that eq. (41) reduces to

$\displaystyle\frac{\rho_{\mathrm{eq}}({\bf
r}_{0})-\rho_{\mathrm{ref}}^{0}}{\rho_{\mathrm{ref}}^{0}}$
$\displaystyle\approx\int d{\bf r}^{\prime}\,c_{0}({\bf r}_{0}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})+\int d{\bf
r}"\,c_{0}({\bf r}_{0}-{\bf r}")\left\\{{\rho_{\mathrm{eq}}({\bf
r}")-\rho_{\mathrm{ref}}({\bf r}")}\right\\}$ $\displaystyle=\int d{\bf
r}^{\prime}\,c_{0}({\bf r}_{0}-{\bf
r}^{\prime})\left\\{\rho_{\mathrm{eq}}({\bf
r}^{\prime})-\rho_{\mathrm{ref}}^{0}\right\\},$ (43)

which is equal to the HLR equation (23).

##  5 Comparison with modified PB equations for the Coulomb potential
$\Psi({\bf r})$

Before comparing the unified equation (39) for $\rho_{\mathrm{eq}}({\bf r})$
with the modified PB equations such as eqs. (24) to (26), we would like to see
that the essential equation (39) implies the Boltzmann distribution of
$\rho_{\mathrm{eq}}({\bf r})$. Actually, in the approximation of $1+x\approx
e^{x}$, eq. (39) reads for point-charge systems

$\displaystyle\rho_{\mathrm{eq}}({\bf r})$
$\displaystyle\approx\rho_{\mathrm{ref}}({\bf r})e^{-\Psi({\bf r})},$ (44)
$\displaystyle\Psi({\bf r})$ $\displaystyle=-\int d{\bf
r}^{\prime}\,h_{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime}),$ (45)

where $h_{\mathrm{el}}({\bf r})$ denotes the electrostatic TCF of point
charges interacting via $v_{\mathrm{el}}({\bf r})$ without the potential
splitting. This section investigates the modified PB equations, or the SC
equations for the Coulomb potential $\Psi({\bf r})$ created by Coulomb
interactions between point charges. We consider two types of the modified PB
equations: the finite-spread PB equations [14, 20-24] based on the charge
smearing model and the higher-order PB equations [15-18, 25-30]. We show below
that both forms can be derived from the Poisson equation for $\Psi({\bf r})$
defined by eq. (45) using the electrostatic TCF $h_{\mathrm{el}}({\bf r})$.

###  5.1 Relation between the Coulomb potential $\Psi({\bf r})$ and the DCF
$c_{\mathrm{el}}({\bf r})$

Plugging the Ornstein-Zernike equation (40) into eq. (45) as before, we have

$\displaystyle\Psi({\bf r})$ $\displaystyle=-\int d{\bf
r}^{\prime}\,c_{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})-\iint d{\bf
r}^{\prime}d{\bf r}"\,c_{\mathrm{el}}({\bf r}-{\bf
r}")\rho_{\mathrm{ref}}({\bf r}")h_{\mathrm{el}}({\bf r}"-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$ $\displaystyle=-\int
d{\bf r}^{\prime}c_{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})-\int d{\bf
r}"\,c_{\mathrm{el}}({\bf r}-{\bf r}")\rho_{\mathrm{ref}}({\bf r}")\Psi({\bf
r}")-\Psi_{c}.$ (46)

In the last line of eq. (46), the constant potential $\Psi_{c}$ is related to
the isothermal compressibility $\kappa_{T}$ as
$\Psi_{c}=\left(k_{B}T\rho_{\mathrm{ref}}^{0}\kappa_{T}\right)^{-1}-1$ because
of the following relation:

$\displaystyle 1-\rho_{\mathrm{ref}}^{0}\int d{\bf r}\,c_{\mathrm{el}}({\bf
r})=\frac{1}{k_{B}T\rho_{\mathrm{ref}}^{0}\kappa_{T}}.$ (47)

The resulting expressions of the Poisson equation for the above potential
$\Psi$ vary, according to the treatments of the Laplace equation for the
electrostatic DCF $c_{\mathrm{el}}({\bf r})$.

To see the different expressions of the Laplace equation, it is convenient to
introduce the Fourier-transform of $c_{\mathrm{el}}({\bf r})$. In general,
$c_{\mathrm{el}}({\bf k})$ can be written as

$c_{\mathrm{el}}({\bf k})=-z^{2}l_{B}\frac{f({\bf k})}{k^{2}},$ (48)

where $f({\bf k})$ take various forms depending on the approximations adopted:

$\displaystyle f({\bf
k})=\left\\{\begin{array}[]{l}e^{-\frac{(ka)^{2}}{4m}}\quad\qquad(\mathrm{HNC}\>\cite[cite]{[\@@bibref{}{ng}{}{}]})\\\
\\\
\left\\{\frac{3}{ka}j_{1}(ka)\right\\}^{2}\quad(\mathrm{Soft\,MSA}\>[22,39,40])\\\
\\\
\frac{1}{1+(ka)^{2}/4}\quad\qquad(\mathrm{Bessel\,smoothed\,model}\>[35]),\end{array}\right.$
(54)

where $k=|{\bf k}|$ and $a$ is the mean separation between adjacent charges as
mentioned before eq. (24) (see also Appendix A.1 for the detailed definition
of $a$). In eq. (54), we use $m=1.08$ according to the HNC approximation [38],
the soft MSA [22, 39, 40] signifies the mean spherical approximation (MSA) [1,
41] for soft spheres, and the smeared interaction potential is identified with
the DCF in the Bessel smoothed model [35] that supposes a smeared charge
distribution expressed by the modified Bessel function $K_{1}$ (see Appendix
A.2 for the details).

###  5.2 Finite-spread PB equation [14, 20-24

In the real-space representation of the electrostatic DCF given by eq. (48),
we have

$\displaystyle-\nabla^{2}\int d{\bf r}^{\prime}c_{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})$ $\displaystyle=\iint d{\bf
s}d{\bf r}^{\prime}\,\nabla^{2}v_{\mathrm{el}}({\bf r}-{\bf s})f({\bf s}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime})$ $\displaystyle=-4\pi
z^{2}l_{B}\int d{\bf s}\,\delta({\bf r}-{\bf s})\left<\rho_{\mathrm{ref}}({\bf
s})\right>$ $\displaystyle=-4\pi z^{2}l_{B}\left<\rho_{\mathrm{ref}}({\bf
r})\right>,$ (55)

where we have introduced the smeared density, $\left<\rho_{\mathrm{ref}}({\bf
r})\right>$, defined by

$\displaystyle\left<\rho_{\mathrm{ref}}({\bf s})\right>$ $\displaystyle=\int
d{\bf r}^{\prime}\,f({\bf s}-{\bf r}^{\prime})\,\rho_{\mathrm{ref}}({\bf
r}^{\prime})$ $\displaystyle=\int d{\bf r}^{\prime}\,f({\bf s}-{\bf
r}^{\prime})\,e^{\mu-J_{1}({\bf r}^{\prime})}.$ (56)

From eqs. (46) and (55), we find

$\displaystyle\nabla^{2}\Psi({\bf r})$ $\displaystyle=-4\pi
z^{2}l_{B}\left<\rho_{\mathrm{ref}}({\bf r})\right>+4\pi
z^{2}l_{B}\left<\rho_{\mathrm{ref}}({\bf r})\right>\Psi({\bf r})$ (57)
$\displaystyle\approx-4\pi z^{2}l_{B}\left<\rho_{\mathrm{eq}}({\bf
r})\right>,$ (58)

where the approximation (44) has been applied to the above second line:

$\displaystyle\left<\rho_{\mathrm{eq}}({\bf r})\right>$
$\displaystyle=\left<\rho_{\mathrm{ref}}({\bf r})e^{-\Psi({\bf r})}\right>$
$\displaystyle=\int d{\bf s}\,f({\bf r}-{\bf s})\,\rho_{\mathrm{ref}}({\bf
s})e^{-\Psi({\bf s})};$ (59)

see also Appendix A.3 and A.5 for the details. It is found from eqs. (55) and
(56) that the introduction of smeared charge model is compatible with the use
of the DCF in the liquid-state theory, as has been verified for the strongly
coupled OCP (see Appendix A.5 for the details) [37, 39, 40]. Equations (56),
(58) and (59) are reduced to eqs. (24) and (25) in the HNC approximation that
provides $f({\bf r}-{\bf s})=\\{\pi(a/m)^{2}\\}^{-3/2}\,e^{-\left(m|{\bf
r}-{\bf s}|/a\right)^{2}}$; however, a discrepancy exists between the values
of $m$: $m=1.08$ for the HNC [38], whereas the local molecular field theory
[14, 20, 21] set that $m=1/0.6\approx 1.67$ for a counterion system as
mentioned after eq. (25). Considering that the effective mean distance between
point charges is diminished to $a/m$ by the factor $1/m$ in the Gaussian
charge smearing model, it is suggested from the fitting result ($m\approx
1.67$) of the local molecular theory that the effective mean distance $a/m$
should be shorter than that of the HNC approximation of the uniform OCP model
for a quantitative description of inhomogeneous counterion distribution in
perpendicular direction to oppositely charged surface.

###  5.3 Higher-order Poisson-Boltzmann equations [15-18, 25-30]

The higher-order PB equations add higher-order gradient terms to the Laplacian
appearing in the usual Poisson equation (see also Appendix A.4 for the
details). The Fourier-space representation of eqs. (58) and (59) is

$\displaystyle-k^{2}\Psi(k)=-4\pi l_{B}z^{2}f({\bf
k})\rho_{\mathrm{ref}}(-{\bf k})e^{-\Psi({\bf k})}.$ (60)

While $f({\bf k})\rho(-{\bf k})$ can be interpreted as a smeared density when
leaving $f({\bf k})$ on the right hand side, the higher-order PB equations are
generated by investigating the contribution of $f^{-1}({\bf k})$ on the left
hand side. In other words, this subsection focuses on how to handle
$k^{2}/f({\bf k})$.

In the long-wave approximation valid for $ka\ll 1$, we have

$\displaystyle\frac{1}{f({\bf
k})}\approx\left\\{\begin{array}[]{l}1+\frac{(ka)^{2}}{4m}\quad\qquad(\mathrm{HNC})\\\
\\\ 1+\frac{(ka)^{2}}{5}\quad\quad(\mathrm{Soft\,MSA})\\\ \\\
1+\frac{(ka)^{2}}{4}\quad\qquad(\mathrm{Bessel}),\end{array}\right.$ (66)

In the above results, $m=1.08$ as before, whereas the above approximate form
for the soft MSA [37, 39, 40] is derived as follows: we have used the
approximation that

$\left\\{\frac{3}{ka}j_{1}(ka)\right\\}^{-2}\approx 1+\frac{(ka)^{2}}{5},$
(67)

expanding the Bessel function $j_{1}(x)$ as

$\displaystyle j_{1}(x)$ $\displaystyle=\frac{\sin x-x\cos x}{x^{2}}$
$\displaystyle\approx\frac{x}{3}\left(1-\frac{x^{2}}{10}\right),$ (68)

where use has been made of the approximation,

$\displaystyle\sin x-x\cos x\approx\frac{x^{3}}{3}-\frac{x^{5}}{30},$ (69)

based on the expansions such that $\sin x\approx x-x^{3}/6+x^{5}/120$ and
$\cos x\approx 1-x^{2}/2+x^{4}/24$.

It follows from eq. (66) that

$\displaystyle-\frac{k^{2}}{f({\bf k})}\approx-\left\\{1+(\gamma
a)^{2}k^{2}\right\\}k^{2},$ (70)

with various values of $\gamma$:

$\displaystyle\gamma=\left\\{\begin{array}[]{l}\frac{1}{2\sqrt{m}}\quad(\mathrm{HNC})\\\
\\\ \frac{1}{\sqrt{5}}\quad(\mathrm{Soft\,MSA})\\\ \\\
\frac{1}{2}\quad(\mathrm{Bessel}).\end{array}\right.$ (76)

It is found from eqs. (60), (70) and (76) that eqs. (58) and (59) can be
reduced to

$\displaystyle\left\\{1-(\gamma a)^{2}\nabla^{2}\right\\}\nabla^{2}\Psi({{\bf
r}})=-4\pi l_{B}z^{2}\rho_{\mathrm{ref}}({\bf r})e^{-\Psi({\bf r})},$ (77)

using eq. (76) as the values of $\gamma$. Thus, the higher-order PB equation
of the BSK form [15-18] has been verified.

To validate another type of the higher-order PB equation given by eq. (26), we
need to perform the expansion of $f^{-1}({\bf k})$ up to the order $k^{4}$. In
the HNC approximation [38], for example, it is straightforward to show that

$\frac{1}{f({\bf k})}=e^{\frac{(ka)^{2}}{4m}}\approx
1+\frac{1}{4m}(ka)^{2}+\frac{1}{32m^{2}}(ka)^{4},$ (78)

yielding

$\displaystyle\nabla^{2}\left(1-\frac{a^{2}}{4m}\nabla^{2}+\frac{a^{4}}{32m^{2}}\nabla^{4}\right)\Psi({{\bf
r}})=-4\pi l_{B}z^{2}e^{\mu-J_{1}({\bf r})-\Psi({\bf r})},$ (79)

with $m=1.08$. In the parentheses on the left hand side of eq. (79), the
coefficient of $\nabla^{2}$–term is close to that of of eq. (26), and yet the
coefficient of the $\nabla^{4}$–term is somewhat smaller than the previous
factor in eq. (26).

##  6 Derivation scheme of the grand-potential density functional given by
eqs. (27) to (29)

In this section, we outline how to derive eqs. (27) to (29) starting with the
density functional integral representation of the grand potential. More
details of the following formulations are given in Appendix B.

###  6.1 Density-functional integral representation of the grand potential for
the attractive particle systems

Let $\Omega[v,J]$ be the grand potential of non-uniform liquids where
constituent particles in an external potential $J({\bf r})$ interact via an
interaction potential $v({\bf r})$. Once the density functional form of
$\Omega[v,J]$ is found for attractive particle systems, other grand-potential
density functionals for only repulsive systems (i.e., hard sphere fluids and
point-charge systems) is obtained by setting that $v_{\mathrm{att}}({\bf
r})\equiv 0$ and replacing $v_{0}$ with either hard sphere or Coulomb
interaction potential (i.e., $v_{\mathrm{HS}}$ or $v_{\mathrm{el}}$).

As detailed in Appendix B, the grand potential can be represented by the
functional integral over the density field as follows:

$\displaystyle e^{-\Omega[v,\,J]}$ $\displaystyle=\int
D\rho\,e^{-\Omega^{*}[\rho]}.$ (80)

Here the conditional grand potential $\Omega^{*}[\rho]$ of a given density
field $\rho({\bf r})$ is an extension of the grand-potential density
functional formulated in the conventional DFT [6-12]:

$\displaystyle\Omega^{*}[\rho]=\Omega_{\mathrm{V}}[\rho]+\Delta\Omega,$ (81)

where $\Omega_{\mathrm{V}}[\rho]$ denotes the variational grand potential,
which reduces to the equilibrium grand potential $\Omega$ when
$\rho=\rho_{\mathrm{eq}}$, and $\Delta\Omega$ corresponds to the deviation
free energy from the equilibrium grand potential due to the imposition of a
given density field $\rho({\bf r})$, instead of $\rho_{\mathrm{eq}}({\bf r})$.
In this subsection, resulting forms of both $\Omega_{\mathrm{V}}[\rho]$ and
$\Delta\Omega[\rho]$ are only presented (see Appendix B.1 for the detailed
derivation).

The variational grand potential $\Omega_{\mathrm{V}}[\rho]$ is given by

$\displaystyle\Omega_{\mathrm{V}}[\rho]$
$\displaystyle=\mathcal{F}_{0}[\rho]-\int d{\bf r}\,\rho({\bf
r})\mu+U_{1}[\rho],$ (82) $\displaystyle U_{1}[\rho]$
$\displaystyle=\frac{1}{2}\iint d{\bf r}d{\bf r^{\prime}}\rho({\bf
r})v_{\mathrm{att}}({\bf r}-{\bf r}^{\prime})\rho({\bf r}^{\prime})+\int d{\bf
r}\,\rho({\bf r})J_{1}({\bf r}).$ (83)

As expected from the subscript ”0”, the first term $\mathcal{F}_{0}[\rho]$ on
the rhs of eq. (82) represents the intrinsic Helmholtz free energy of the only
repulsive systems where particles in the absence of external field $J({\bf
r})\equiv 0$ interact via the interaction potential $v_{0}({\bf r})$. We can
see from eq. (83) that the attractive interaction energy is evaluated in the
mean-field approximation. The second term on the rhs of eq. (83), which
expresses the one-body potential energy, is coupled only with the slowly-
varying external field $J_{1}$ because the configurational integration range
is set outside the excluded region defined by eqs. (16) and (20).

Following the functional integral formulation combined with the conventional
DFT, the first approximation of the additional contribution $\Delta\Omega$
provides the logarithmic form in the Gaussian approximation of the fluctuating
one-body potential field, dual to a given density field $\rho({\bf r})$ (see
Appendix B.2 for the derivation):

$\displaystyle\Delta\Omega$
$\displaystyle=\frac{1}{2}\ln\left|\mathrm{det}\,G_{0}({\bf
r};\rho^{0}_{\mathrm{ref}})\right|$ (84) $\displaystyle G_{0}({\bf
r};\rho^{0}_{\mathrm{ref}})$
$\displaystyle=\rho^{0}_{\mathrm{ref}}\left\\{\delta({\bf
r})+\rho^{0}_{\mathrm{ref}}h_{0}({\bf r})\right\\},$ (85)

where we have neglected variation of the density-density correlation function
from $G_{0}({\bf r};\rho^{0}_{\mathrm{ref}})$ for a uniform density
$\rho^{0}_{\mathrm{ref}}$, which is consistent with the other approximations
of $\mathcal{F}_{0}[\rho]$ described below.

###  6.2 Reference density $\rho_{\mathrm{ref}}$ as the saddle-point density

Starting with the above density functional integral form (80), we would like
to derive the mean-field density $\rho_{\mathrm{ref}}({\bf r})$ where the
attractive part $v_{\mathrm{att}}({\bf r})$ of the interaction potential is
considered. To this end, we divide $\mathcal{F}_{0}[\rho]$ into the ideal gas
and excess contributions, $\mathcal{F}_{0}^{\mathrm{id}}[\rho]$ and
$\mathcal{F}_{0}^{\mathrm{ex}}[\rho]$:

$\mathcal{F}_{0}[\rho]=\mathcal{F}_{0}^{\mathrm{id}}[\rho]+\mathcal{F}_{0}^{\mathrm{ex}}[\rho],$
(86)

so that we may rewrite, for later convenience, the conditional grand potential
$\Omega^{*}[\rho]$ given by eqs. (81) to (85) as

$\displaystyle\Omega^{*}[\rho]$
$\displaystyle=\mathcal{F}_{0}^{\mathrm{ex}}[\rho]+\Omega_{\mathrm{att}}^{*}[\rho]+\Delta\Omega,$
(87) $\displaystyle\Omega_{\mathrm{att}}^{*}[\rho]$
$\displaystyle=\mathcal{F}_{0}^{\mathrm{id}}[\rho]+U_{1}[\rho]-\int d{\bf
r}\,\rho({\bf r})\mu.$ (88)

Accordingly, eq. (80) can read

$\displaystyle e^{-\Omega[v,\,J]}$ $\displaystyle=e^{-\Delta\Omega}\int
D\rho\,e^{-\Omega_{\mathrm{att}}^{*}[\rho]}\left(\frac{\int
D\rho\,e^{-\mathcal{F}_{0}^{\mathrm{ex}}[\rho]-\Omega_{\mathrm{att}}^{*}[\rho]}}{\int
D\rho\,e^{-\Omega_{\mathrm{att}}^{*}[\rho]}}\right),$ (89)

using eqs. (87) and (88). This form (89) reveals that the mean-field
approximation of the attractive interaction energy becomes equivalent to the
use of the saddle-point density determined by

$\left.\frac{\delta\Omega_{\mathrm{att}}^{*}[\rho]}{\delta\rho}\right|_{\rho=\rho_{\mathrm{ref}}}=0$
(90)

because combination of eqs. (88) and (90) yields

$\displaystyle\ln\rho_{\mathrm{ref}}({\bf r})+\int d{\bf
r}^{\prime}v_{\mathrm{att}}({\bf r}-{\bf r}^{\prime})\rho_{\mathrm{ref}}({\bf
r}^{\prime})+J_{1}({\bf r})-\mu=0,$ (91)

which is identical to the previous result (15).

###  6.3 Hybrid evaluation based on the mean-field and Gaussian approximations

Following the previous treatment, we approximate the functional integral (89)
by

$\displaystyle e^{-\Omega[v,\,J]}\approx
e^{-\widetilde{\Omega}[\rho_{\mathrm{ref}}]}\equiv
e^{-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]-\Delta\Omega}\left(\frac{\int
Dn\,e^{-\mathcal{F}_{0}[\rho_{\mathrm{ref}}+n]-U_{1}[\rho_{\mathrm{ref}}+n]}}{e^{-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]}}\right),$
(92)

extracting only the saddle-point path for the attractive part of the grand
potential. The approximate form (92) of $\Omega[v,\,J]$ has been referred to
as the grand-potential density functional which is denoted by
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ in eq. (27). We expand
$\mathcal{F}_{0}[\rho_{\mathrm{ref}}+n]+U_{1}[\rho_{\mathrm{ref}}+n]$ around
the reference energy up to the quadratic term with respect to fluctuating
density field $n({\bf r})$, yielding (see Appendix B.2 for the derivation)

$\displaystyle\mathcal{F}_{0}[\rho_{\mathrm{ref}}+n]+U_{1}[\rho_{\mathrm{ref}}+n]-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]$
$\displaystyle\quad\approx\mathcal{F}^{\mathrm{ex}}_{0}[\rho_{\mathrm{ref}}]-\iint
d{\bf r}d{\bf r}^{\prime}\,n({\bf r})\,c_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$
$\displaystyle\qquad\qquad+\frac{1}{2}\iint d{\bf r}d{\bf r}^{\prime}\,n({\bf
r})\,G_{0}^{-1}({\bf r}-{\bf r}^{\prime};\rho_{\mathrm{ref}}^{0})\,n({\bf
r}^{\prime}),$ (93)

supposing that the intrinsic Helmholtz free energy
$\mathcal{F}_{0}[\rho_{\mathrm{ref}}]$ for the $v_{0}$-interaction systems are
of the Ramakrishnan-Yussouff functional form [1, 6, 7]:

$\displaystyle\mathcal{F}_{0}[\rho_{\mathrm{ref}}]=\mathcal{F}_{0}[\rho^{0}_{\mathrm{ref}}]-\frac{1}{2}\iint
d{\bf r}d{\bf r^{\prime}}\Delta\rho_{\mathrm{ref}}({\bf r})c_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})+\int d{\bf
r}\rho_{\mathrm{ref}}({\bf r})\ln\left\\{\frac{\rho_{\mathrm{ref}}({\bf
r})}{\rho^{0}_{\mathrm{ref}}}\right\\},$ (94)

where the same approximation used in eqs. (42), (84), and (85) has been
applied to the last line of eq. (93) (see also eq. (138)).

Accordingly, eq. (92) is further reduced to

$\displaystyle e^{-\widetilde{\Omega}[\rho_{\mathrm{ref}}]}$
$\displaystyle=e^{-\mathcal{F}_{0}^{\mathrm{ex}}[\rho_{\mathrm{ref}}]-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]-\Delta\Omega}\left(\int
Dn\,e^{-\Delta\mathcal{F}_{0}[n]}\right)$ (95)
$\displaystyle\Delta\mathcal{F}_{0}[n]$ $\displaystyle=\frac{1}{2}\iint d{\bf
r}d{\bf r}^{\prime}\,n({\bf r})\,G_{0}^{-1}({\bf r}-{\bf
r}^{\prime};\rho_{\mathrm{ref}}^{0})n({\bf r}^{\prime})-\iint d{\bf r}d{\bf
r}^{\prime}\,n({\bf r})c_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime}),$ (96)

indicating that the remaining task is to perform the functional integral over
the fluctuating $n$-field. As proved in Appendix B.3, eq. (96) reads

$\displaystyle\Delta\mathcal{F}_{0}[n]$ $\displaystyle=\frac{1}{2}\iint d{\bf
r}d{\bf r}^{\prime}\,\widetilde{n}({\bf r})\,G_{0}^{-1}({\bf r}-{\bf
r}^{\prime};\rho_{\mathrm{ref}}^{0})\widetilde{n}({\bf r}^{\prime})$
$\displaystyle\quad-\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\Delta\rho_{\mathrm{ref}}({\bf r})\left\\{h_{0}({\bf r}-{\bf
r}^{\prime})-c_{0}({\bf r}-{\bf
r}^{\prime})\right\\}\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime}),$ (97)

due to the shift of fluctuating field from $n$ to

$\displaystyle\widetilde{n}({\bf r})$ $\displaystyle=n({\bf r})-\iint d{\bf
r}^{\prime}d{\bf r}"\>G_{0}({\bf r}-{\bf
r}^{\prime};\rho_{\mathrm{ref}}^{0})c_{0}({\bf r}^{\prime}-{\bf
r}")\Delta\rho_{\mathrm{ref}}({\bf r}")$ $\displaystyle=n({\bf r})-\iint d{\bf
r}^{\prime}d{\bf r}"\>\rho_{\mathrm{ref}}^{0}\left\\{\delta({\bf r}-{\bf
r}^{\prime})+\rho_{\mathrm{ref}}^{0}h_{0}({\bf r}-{\bf
r}^{\prime}))\right\\}c_{0}({\bf r}^{\prime}-{\bf
r}")\Delta\rho_{\mathrm{ref}}({\bf r}")$ $\displaystyle=n({\bf
r})-\rho_{\mathrm{ref}}^{0}\int d{\bf r}"h_{0}({\bf r}-{\bf
r}")\Delta\rho_{\mathrm{ref}}({\bf r}")$ $\displaystyle\approx n({\bf
r})-\left\\{\rho_{\mathrm{eq}}({\bf r})-\rho_{\mathrm{ref}}({\bf
r})\right\\},$ (98)

where the last approximation is the same as eq. (42). It follows from eqs.
(95) and (97) that the Gaussian integration over the $\widetilde{n}$–field
leads to

$\displaystyle e^{-\widetilde{\Omega}[\rho_{\mathrm{ref}}]}$
$\displaystyle=e^{-\mathcal{F}^{\mathrm{ex}}_{0}[\rho_{\mathrm{ref}}]-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]-\Delta\Omega+\frac{1}{2}\iint
d{\bf r}d{\bf r}^{\prime}\Delta\rho_{\mathrm{ref}}({\bf r})\left\\{h_{0}({\bf
r}-{\bf r}^{\prime})-c_{0}({\bf r}-{\bf
r}^{\prime})\right\\}\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\underline{\left\\{\int
D\widetilde{n}\,e^{-\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\widetilde{n}({\bf r})G_{0}^{-1}({\bf r}-{\bf
r}^{\prime})\widetilde{n}({\bf r}^{\prime})}\right\\}}$
$\displaystyle=e^{-\mathcal{F}^{\mathrm{ex}}_{0}[\rho_{\mathrm{ref}}]-\Omega_{\mathrm{att}}^{*}[\rho_{\mathrm{ref}}]+\frac{1}{2}\iint
d{\bf r}d{\bf r}^{\prime}\Delta\rho_{\mathrm{ref}}({\bf r})\left\\{h_{0}({\bf
r}-{\bf r}^{\prime})-c_{0}({\bf r}-{\bf
r}^{\prime})\right\\}\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})},$ (99)

where the logarithmic term $\Delta\Omega$ given by eq. (84) is canceled by the
result of Gaussian integration over the $\widetilde{n}$–field in the
underlined term. Combining eqs. (88), (94) and (99), we obtain eqs. (27) to
(29), the central results for $\widetilde{\Omega}[\rho_{\mathrm{ref}}]$.

##  7 Summary and conclusions

Figure 1:  A schematic flow chart of proving the equivalence of previous SC
equations [4, 5, 14-34, 42-47] based on the hybrid method [48, 49] that
incorporates the density functional theory (DFT) [6-12] into the statistical
field theory [15, 23, 31-34, 36, 50-52]. The saddle-point equation (15), or
eq. (90), determines the reference density $\rho_{\mathrm{ref}}({\bf r})$. The
Gaussian approximation of density fluctuations around
$\rho_{\mathrm{ref}}({\bf r})$ yields the grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ given by eqs. (27) to (29). The
equilibrium density $\rho_{\mathrm{eq}}({\bf r})$ is obtained from
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ using the relation (36). The single
unified equation (39) for $\rho_{\mathrm{eq}}({\bf r})$ produces a variety of
previous SC equations, including the hydrostatic linear response (HLR)
equation (23) [4, 5, 14, 42-47], the finite-spread Poisson-Boltzmann (PB)
equation (24) (or eqs. (A.15) to (A.18)) [14, 20-24], and the higher-order PB
equation (26) (or eq. (A.20)) [15-18, 25-30].

The flow chart depicted in Fig. 1 outlines the unified description of previous
SC equations covering the mean-field equation (15), the HLR equation (23) [4,
5, 14, 42-47], the finite-spread PB equation [14, 20-24] given by eq. (24),
and the higher-order PB equation [15-18, 25-30] expressed by eq. (26). The
first step is that the hybrid method [48, 49], as a kind of thermodynamic
perturbation theory, yields the grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ of the reference density
$\rho_{\mathrm{ref}}({\bf r})$ as a solution of the saddle-point equation
(90), or eq. (15). Next, a single unified equation (39) for
$\rho_{\mathrm{eq}}({\bf r})$ is obtained from
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ using the relation (36). The final
step is to demonstrate that the unified equation (39) is reduced to the three
types of SC equations [4, 5, 14-34, 42-47].

Following the process in Fig. 1, we can divide the formulations presented thus
far into the following four parts: (i) the mean-field equation (15) for the
reference density $\rho_{\mathrm{ref}}({\bf r})$, (ii) the grand-potential
density functional $\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ given by eqs.
(27) to (29), (iii) the HLR equation (23) for the equilibrium density
$\rho_{\mathrm{eq}}({\bf r})$ [4, 5, 14, 42-47], and (iv) the modified PB
equations (58), (59), (77), and (79) [14-30]. A detailed summary of each part
is provided below.

(i) The mean-field equation (15) for the reference density
$\rho_{\mathrm{ref}}({\bf r})$.— The proposed formulation clarifies that the
HLR theory [4, 5, 14, 42-47] proceeds in the opposite direction to the
conventional thermodynamic perturbation [1-3] theory that treats the slowly
varying attractive interaction potential as a perturbation to the repulsive
reference system. In contrast to the conventional thermodynamic perturbation
theory, our formulations based on the HLR treatment adopt, as a reference free
energy, the mean-field free energy of attractive particle system under an
external field in the absence of a repulsive interaction potential, i.e.
$v_{0}({\bf r})\equiv 0$. It is the aim of the saddle-point equation (90) for
$\rho_{\mathrm{ref}}({\bf r})$ to determine a reference state that reflects
the entire profile of the non-uniform density distribution using a simple and
precise method. Along with the HLR treatment, we investigated a perturbative
contribution to the reference system from energetic aspects when recovering
the repulsive part $v_{0}({\bf r})$ of the interaction potential.

(ii) The grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$.— The hybrid method leads to a novel
density-functional form of the grand potential given by eqs. (27) to (29),
thereby demonstrating that the same density-functional form validates the
previous SC equations [4, 5, 14-30, 42-47]. It has been shown that there are
two main features of the density-functional form. One is that
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$ is a functional of the reference
density evaluated in the mean-field approximation, instead of the equilibrium
density $\rho_{\mathrm{eq}}({\bf r})$ used in the conventional DFT. The other
is that the TCF, specifically $-h_{0}({\bf r})$, is treated as effective
interaction potential instead of the DCF. The TCF, which is generated from the
DCF due to the Ornstein-Zernike equation, emerges as a result of considering
fluctuations around $\rho_{\mathrm{ref}}({\bf r})$ when repulsive interactions
are inserted into the systems as described above. In more detail, the second
term on the rhs of eq. (93), which is the coupling term for the fluctuating
density $n({\bf r})$ with the reference density $\rho_{\mathrm{ref}}({\bf
r})$, plays a key role in generating the effective interaction energy
expressed by the TCF. The non-vanishing coupling term reveals that the
proposed formulation follows the core concept of the HLR theory [4, 5, 14,
42-47].

(iii) The HLR equation (23) for the equilibrium density
$\rho_{\mathrm{eq}}({\bf r})$.— The difference between the HLR method and our
theory is that we treat the slowly varying part of the repulsive interaction
potential, including the long-range part of the Coulomb interaction potential,
as a perturbative potential added in the reference system, even though the
other slowly varying parts are evaluated in the mean-field approximation
according to the previous SC equations reviewed in Section 2. The local
molecular field theory [14, 20, 21] starts with the a priori separation of
$v_{\mathrm{el}}({\bf r})$ into harshly repulsive and slowly varying parts,
thus providing the finite-spread PB equation [14, 20-24] or the aforementioned
mean-field equation based on a Gaussian charge smearing model with the harshly
repulsive part of $v_{\mathrm{el}}({\bf r})$ to be cut-off [14, 20-24, 35-40].
However, our goal is to clarify the underlying mechanisms that generate a
variety of SC equations from a unified perspective. Therefore, we investigated
the theoretical background of why the Gaussian charge smearing model is
relevant for describing charge-charge correlations.

As a result, we verified the HLR equation [4, 5, 14, 42-47] for any repulsive
interaction system using a unified framework based on the recently developed
density functional integral representation [48, 49]. The field-theoretic
formulation of the HLR equation also demonstrates the benefits of the
hydrostatic approximation [4, 5, 14, 42-47], which regards a uniform liquid
with a density of $\rho_{\mathrm{ref}}^{0}$, instead of the bulk density, as
another reference system for evaluating the grand-potential density functional
of non-uniform liquids.

(iv) A set of modified PB equations, eqs. (58), (59), (77) and (79).— In terms
of the Coulomb potential $\Psi({\bf r})$ defined as eq. (45) using the TCF, we
can see that the HLR equation (23), expressed by the TCF, forms the basis of
validating the Boltzmann distribution (44) for $\rho_{\mathrm{eq}}({\bf r})$
in both the finite-spread PB equation [14, 20-24] and higher-order PB equation
[15-18, 25-30]. These modifications of the PB equation are ascribed to the
Poisson equation (55) of the electrostatic DCF $c_{\mathrm{el}}({\bf r})$.
Furthermore, the Ornstein-Zernike equation yields the SC equation (46) for
$\Psi({\bf r})$ with the use of $c_{\mathrm{el}}({\bf r})$. Combining these
relations, eqs. (46) and (55), we obtain the linearized SC equation (57) for
$\Psi({\bf r})$, which transforms to alternative representations, the finite
and higher-order PB equations, depending on how the weight function $f({\bf
r})$ is approximated and/or interpreted. We stress again that these
modifications of the conventional PB equation emerge due to the electrostatic
DCF, and that the electrostatic DCF naturally appears when calculating
$\rho_{\mathrm{eq}}({\bf r})$ from the grand-potential density functional
$\widetilde{\Omega}[\rho_{\mathrm{ref}}]$.

In conclusion, we have proved the equivalence of previous SC equations [4, 5,
14-34, 42-47] by demonstrating that the apparent diversity of SC equations can
be explained by the identical form of the grand-potential density functional
given by eqs. (27) to (29). This theoretical finding implies that not only
improvements to the previous forms but also derivations of any new equations
can be performed in a systematic manner, by either exploring a new system to
which our formulations may apply or considering higher-order terms beyond the
Gaussian approximation.

##  Appendix A Summary of previous results related to modified PB equations
for point-charge systems

###  A.1 Coupling constant in non-uniform systems

Let $a$ and $\Gamma$ be the mean distance between point charges and the
coupling constant representing the strength of Coulomb interactions,
respectively. As is well-known, $\Gamma$ in the uniform OCP is given by

$\displaystyle\Gamma=z^{2}\frac{l_{B}}{a},$ (100)

when the valence of each charge is $z$. In non-uniform point-charge systems,
however, the distance between charges cannot be uniquely determined.
Correspondingly, it is not trivial to define $\Gamma$, unlike the conventional
uniform OCP. Two types of coupling constants have been defined so far for the
typical non-uniform systems of small counterions that are dissociated from
oppositely charged macroions, or macroscopic objects including charged plates:
while one uses the Gouy-Chapman length as a reference length $a$, the other
definition of $\Gamma$ is based on an analogy with the two-dimensional OCP.

We have adopted the latter definition as follows. When all of the conterions
are localized on the entire surface of macroions oppositely charged, the
distance $b$ between counterions on the surface is defined by

$\pi b^{2}\sigma=z,$ (101)

using the number density $\sigma$ per area of surface charges. Obviously, $b$
is the minimum of counterion-counterion separation: $b\leq a$. The second
definition of $\Gamma$ in counterion systems therefore implies that eq. (100)
leads to the unique definition of $\Gamma$ for non-uniform systems as long as
$a$ is identified with $b$. That is, we set that $a=b$ not only in eq. (100)
but also in the main text. Consequently, the relation $\Gamma\gg 1$ continues
to represent the strong charge-charge correlations in non-uniform systems.

###  A.2 Two types of weighted densities

There are two methods to formulate charge smearing models [14, 20-24, 35-40]
although not described explicitly. To see this, we first investigate the
relationship between three density functional forms of the slowly varying part
of the Coulomb interaction energy $U_{1}^{\mathrm{el}}$ that is associated
with the slowly varying part of the Coulomb interaction potential such as eq.
(14).

The first form expresses $U_{1}^{\mathrm{el}}$ as the sum of Coulomb
interactions between smeared charges with the interaction potential being
given by the original interaction potential $v_{\mathrm{el}}({\bf r})$:

$\displaystyle U_{1}^{\mathrm{el}}[\rho]$
$\displaystyle=\frac{z^{2}l_{B}}{2}\iint d{\bf r}d{\bf
r}^{\prime}\,\overline{\,\rho}({\bf r})v^{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\,\overline{\,\rho}({\bf r}^{\prime}),$ (102)
$\displaystyle\overline{\,\rho}({\bf r})$ $\displaystyle=\int d{\bf
s}\,\rho({\bf s})\,\omega({\bf s}-{\bf r}),$ (103)

where $\overline{\,\rho}({\bf r})$ is smeared due to the weighted function
$\omega({\bf s}-{\bf r})$; for instance, we have

$\displaystyle\omega({\bf r})=\left\\{\begin{array}[]{l}\frac{1}{\\{\pi
a^{2}/(2m^{2})\\}^{3/2}}e^{-2\left(\frac{mr}{a}\right)^{2}}\quad(\mathrm{Gaussian})\\\
\\\ \frac{3}{4\pi a^{3}}\Theta(a-r)\quad\>(\mathrm{Homogeneous})\\\ \\\
\frac{2}{\pi
a^{2}}\left(\frac{K_{1}(2r/a)}{r}\right)\quad\>\>\>\>\>(\mathrm{Bessel}),\end{array}\right.$
(109)

which are the models that consider the radial distribution given by the
modified Bessel function $K_{1}$ [26] as well as the Gaussian and homogeneous
charge distributions.

Equations (102) and (103) naturally lead to the second description of
$U_{1}^{\mathrm{el}}$ where the slowly varying interaction potential
$v_{1}({\bf r})$ is introduced without smearing the density:

$\displaystyle U_{1}^{\mathrm{el}}[\rho]$
$\displaystyle=\frac{z^{2}l_{B}}{2}\iint d{\bf s}d{\bf s}^{\prime}\rho({\bf
s})v_{1}^{\mathrm{el}}({\bf s}-{\bf s}^{\prime})\rho({\bf s}^{\prime}).$ (110)

Comparing eqs. (102), (103), and (110), it is found that

$\displaystyle v_{1}^{\mathrm{el}}({\bf s}-{\bf s}^{\prime})$
$\displaystyle=\iint d{\bf r}d{\bf r}^{\prime}\,\omega({\bf s}-{\bf
r})\,v^{\mathrm{el}}({\bf r}-{\bf r}^{\prime})\,\omega({\bf r}^{\prime}-{\bf
s}^{\prime}),$ (111)

which also reads

$\displaystyle v_{1}^{\mathrm{el}}({\bf s}-{\bf s}^{\prime})=\iint d{\bf
r}^{\prime}d{\bf r}"\,v^{\mathrm{el}}({\bf s}-{\bf r}")\,\omega({\bf r}"-{\bf
r}^{\prime})\,\omega({\bf r}^{\prime}-{\bf s}^{\prime})$ (112)

due to the coordinate change such that ${\bf r}"={\bf s}-{\bf r}+{\bf
r}^{\prime}$.

Equations (110) and (112) thus generate the third form:

$\displaystyle U_{1}^{\mathrm{el}}[\rho]$
$\displaystyle=\frac{z^{2}l_{B}}{2}\iint d{\bf r}d{\bf r}"\rho({\bf
r})v^{\mathrm{el}}({\bf r}-{\bf r}")\overline{\overline{\,\rho}}({\bf r}"),$
(113) $\displaystyle\overline{\overline{\,\rho}}({\bf r}")$
$\displaystyle=\iint d{\bf r}^{\prime}d{\bf s}^{\prime}\,\omega({\bf r}"-{\bf
r}^{\prime})\,\omega({\bf r}^{\prime}-{\bf s}^{\prime})\rho({\bf
s}^{\prime}).$ (114)

Equation (114) represents the second type of the weighted density used in the
main text; as will be seen below, the expression (114) yields the same form as
the distribution function given in eq. (25).

###  A.3 Finite-spread PB equations [14, 20-24]

There are two kinds of expressions for the Coulomb potentials created by
counterion-counterion interactions in correspondence with the above two
definitions of smeared density, i.e., $\overline{\,\rho_{\mathrm{ref}}}({\bf
r})$ and $\overline{\overline{\,\rho_{\mathrm{ref}}}}({\bf r})$. While a
finite-spread PB theory [22] uses the former definition,

$\displaystyle\psi({\bf r})=\int d{\bf r}^{\prime}v^{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\,\overline{\,\rho_{\mathrm{ref}}}({\bf r}^{\prime}),$ (115)

the local molecular field theory [14, 20, 21], another finite-spread PB
theory, adopts the latter:

$\displaystyle\Psi({\bf r})$ $\displaystyle=\int d{\bf
r}^{\prime}\,v^{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\overline{\overline{\,\rho_{\mathrm{ref}}}}({\bf r}^{\prime})$
(116) $\displaystyle=\int d{\bf r}^{\prime}\,v_{1}^{\mathrm{el}}({\bf r}-{\bf
r}^{\prime})\rho_{\mathrm{ref}}({\bf r}^{\prime}).$ (117)

Both potentials are related to each other as follows:

$\Psi({\bf r})=\int d{\bf r}^{\prime}\,\omega({\bf r}-{\bf
r}^{\prime})\psi({\bf r}^{\prime}),$ (118)

indicating that $\Psi({\bf r})$ corresponds to the weighted potential of
$\psi({\bf r})$.

We can see from eq. (117), the key representation, that the reference density
given by eq. (15) obeys the Boltzmann distribution such that

$\displaystyle\rho_{\mathrm{ref}}({\bf r})$ $\displaystyle=e^{\mu-J_{1}({\bf
r})-\int d{\bf r}^{\prime}\,\omega({\bf r}-{\bf r}^{\prime})\psi({\bf
r}^{\prime})}$ (119) $\displaystyle=e^{\mu-J_{1}({\bf r})-\Psi({\bf r})}.$
(120)

Furthermore, it is found from the above definitions given by eqs. (116) and
(117) that the Poisson equations of $\psi({\bf r})$ and $\Psi({\bf r})$ yield

$\displaystyle\nabla^{2}\psi({\bf r})$ $\displaystyle=-4\pi z^{2}l_{B}\int
d{\bf r}^{\prime}\,\omega({\bf r}-{\bf r}^{\prime})\,\rho_{\mathrm{ref}}({\bf
r}^{\prime}),$ (121) $\displaystyle\nabla^{2}\Psi({\bf r})$
$\displaystyle=-4\pi z^{2}l_{B}\iint d{\bf r}^{\prime}d{\bf r}"\,\omega({\bf
r}-{\bf r}^{\prime})\,\omega({\bf r}^{\prime}-{\bf
r}")\rho_{\mathrm{ref}}({\bf r}"),$ (122)

respectively. Combining eqs. (119) and (121), we obtain an expression of the
finite-spread PB equation. As mentioned at the end of Appendix A.3, on the
other hand, eq. (122) reads

$\displaystyle\nabla^{2}\Psi({\bf r})$ $\displaystyle=-4\pi
z^{2}l_{B}\left\\{\frac{8}{\\{\pi a^{2}/m^{2}\\}^{3}}\iint d{\bf
r}^{\prime}d{\bf r}"\,e^{-2|m({\bf r}-{\bf r}^{\prime})/a|^{2}}\,e^{-2|m({\bf
r}^{\prime}-{\bf r}")/a|^{2}}\rho_{\mathrm{ref}}({\bf r}")\right\\}$
$\displaystyle=-4\pi z^{2}l_{B}\left\\{\frac{1}{\\{\pi
a^{2}/m^{2}\\}^{3/2}}\int d{\bf r}"\,e^{-|m({\bf r}-{\bf
r}")/a|^{2}}\rho_{\mathrm{ref}}({\bf r}")\right\\}$ (123)

for Gaussian charge smearing model. Combining eqs. (120) and (123), we verify
another form of the finite-spread PB equation given by eqs. (24) and (25) used
in the local molecular field theory [14, 20, 21].

###  A.4 Higher-order PB equations [15-18, 25-30]

One type of higher-order PB equations is truncated at the biharmonic
$\nabla^{2}\nabla^{2}$-term, which can be written as

$\displaystyle\left\\{(\gamma a)^{2}\nabla^{2}-1\right\\}\nabla^{2}\Psi({{\bf
r}})=4\pi l_{B}z^{2}\rho_{\mathrm{ref}}({\bf r})e^{-\Psi({\bf r})},$ (124)

where $\gamma$ takes the value of either $\Gamma$ or $\sim 10^{0}$ [16-18, 25,
28]. This type of equation is the same as the BSK equation in binary systems
as mentioned in the main text. The BSK equation [16-18] has been found to be a
powerful tool in investigating the structural and dynamical properties of
concentrated electrolytes and room temperature ionic liquids where $a$ is
identified with the diameter of ions.

Meanwhile, the other type has been presented in eq (26) which has the
$\nabla^{2}\nabla^{4}$-term as the gradient term of the highest order.
Incidentally, in eq. (26), the coefficients of the $\nabla^{2}\nabla^{2}$–and
$\nabla^{2}\nabla^{4}$–terms use a previous result of $\sigma$, given that
$\Delta=0$ in eq. (26) of Ref [26].

###  A.5 Correspondence between the electrostatic DCF $c_{\mathrm{el}}({\bf
r})$ and the slowly varying part $v_{1}^{\mathrm{el}}({\bf r})$ of the Coulomb
interaction potential

The electrostatic DCF relevant in the strong coupling regime of the OCP can be
regarded as $v_{1}^{\mathrm{el}}({\bf r})$ when considering specific charge
smearing models. Comparison between eqs. (48) and (112) implies that the
equality,

$f({\bf r}-{\bf r}^{\prime})=\int d{\bf r}"\omega({\bf r}-{\bf
r}")\,\omega({\bf r}"-{\bf r}^{\prime}),$ (125)

holds between $f({\bf r})$ and $\omega({\bf r})$ for the strongly coupled OCP.
More specifically, $f({\bf r})$ in the HNC approximation [38] is expressed by
$\omega({\bf r})$ of the Gaussian charge smearing model, whereas $f({\bf r})$
in the soft MSA [37, 39, 40] by $\omega({\bf r})$ of the homogeneous charge
smearing model (see also eqs. (54) and (109)).

##  Appendix B Details of Section 6

###  B.1 Density functional form of the conditional grand potential
$\Omega^{*}[\rho]$: derivation of eqs. (80) to (83) [48, 49]

The configurational representation $\Omega[v,\,J]$ is

$\displaystyle e^{-\Omega[v,\,J]}=\mathrm{Tr}\,\exp\left\\{-\sum_{i<j}v({\bf
r}_{i}-{\bf r}_{j})-\sum_{i}J({\bf r}_{i})\right\\}$
$\displaystyle\mathrm{Tr}\equiv\sum_{N=0}^{\infty}e^{N\mu}\frac{1}{N!}\int
d{\bf r}_{1}\cdots\int d{\bf r}_{N},$ (126)

where the interaction energy $\sum_{i<j}v({\bf r}_{i}-{\bf r}_{j})$ appearing
in the above exponent can read, according to the potential splitting given by
eq. (1),

$\displaystyle\sum_{i<j}v({\bf r}_{i}-{\bf r}_{j})$
$\displaystyle=\sum_{i<j}v_{0}({\bf r}_{i}-{\bf r}_{j})+\frac{1}{2}\iint d{\bf
r}d{\bf r^{\prime}}\hat{\rho}({\bf r})v_{\mathrm{att}}({\bf r}-{\bf
r}^{\prime})\hat{\rho}({\bf r}^{\prime})$ (127)

where $\hat{\rho}({\bf r})=\sum_{i=1}^{N}\delta({\bf r}-{\bf r}_{i})$
represents an instantaneous density of $N$ spheres located at $\\{{\bf
r}_{i}\\}\,(i=1,\cdots,N)$. Bearing in mind the expression (127), the density-
functional integral representation can be incorporated into eq. (126) via the
following relation:

$\displaystyle e^{-\Omega[v,\,J]}$ $\displaystyle=\int D\rho\,\prod_{{\bf
r}}\delta\left[\hat{\rho}({\bf r})-\rho({\bf
r})\right]e^{-\Omega[v_{0},0]-U_{1}[\rho]}$ $\displaystyle=\int
D\rho\left(\int D\psi\,e^{\int d{\bf r}\,\,i\psi({\bf r})\\{\hat{\rho}({\bf
r})-\rho({\bf r})\\}}e^{-\Omega[v_{0},0]-U_{1}[\rho]}\right)$
$\displaystyle=\int D\rho\,e^{-\Omega^{*}[\rho]},$ (128) $\displaystyle
U_{1}[\rho]$ $\displaystyle=\frac{1}{2}\iint d{\bf r}d{\bf
r^{\prime}}\rho({\bf r})v_{\mathrm{att}}({\bf r}-{\bf r}^{\prime})\rho({\bf
r}^{\prime})+\int d{\bf r}\,\rho({\bf r})J({\bf r}),$ (129)

where $\prod_{{\bf r}}\delta\left[\hat{\rho}({\bf r})-\rho({\bf r})\right]$
represents the constraint that a density field $\rho({\bf r})$ must be the
same as an instantaneous density distribution $\hat{\rho}({\bf r})$ at every
position in the system and an identity, $1=\int D\rho\,\prod_{{\bf
r}}\delta\left[\hat{\rho}({\bf r})-\rho({\bf r})\right]$, has been inserted
into the first line of the above equations.

The purpose of this subsection is to find the density-functional form of the
conditional grand potential $\Omega^{*}[\rho]$ with the help of the
conventional DFT [6-12]. To this end, it is necessary to manipulate the grand
potential $\Omega[v_{0},-i\psi]$ under the complex external field $-i\psi({\bf
r})$ because eq. (128) is converted to

$\displaystyle e^{-\Omega^{*}[\rho]}$ $\displaystyle=\int
D\psi\,e^{-U_{1}[\rho]-i\int d{\bf r}\,\rho({\bf r})\psi({\bf
r})}\,\mathrm{Tr}\,\,e^{-\sum_{i<j}v_{0}({\bf r}_{i}-{\bf
r}_{j})+\sum_{i}i\psi({\bf r}_{i})}$
$\displaystyle=e^{-U_{1}[\rho]}\,\underline{\left(\int
D\psi\,e^{-\Omega[v_{0},-i\psi]-i\int d{\bf r}\,\rho({\bf r})\psi({\bf
r})}\right)},$ (130)

implying that we can obtain the density functional from evaluating the
underlined $\psi$–integral.

Let the complex potential field $\psi$ be decomposed into two parts of the
fluctuating field $\phi$ and the saddle-point field $i\psi^{*}$:
$\psi=\phi+i\psi^{*}$. On the one hand, the saddle-point path provides the
first Legendre transform of $\Omega[v_{0},-i\psi]$:

$\displaystyle\left.\frac{\delta\left(\Omega[v_{0},-i\psi]\right)}{\delta\psi({\bf
r})}\right|_{\psi=i\psi^{*}}=-i\rho({\bf r}),$ (131)
$\displaystyle\mathcal{F}_{0}[\rho]=\Omega[v_{0},\psi^{*}]-\int d{\bf
r}\,\rho({\bf r})\left\\{\psi^{*}({\bf r})-\mu\right\\},$ (132)

where $\mathcal{F}_{0}[\rho]$ corresponds to the intrinsic Helmholtz free
energy of the $v_{0}$-interaction system in terms of the conventional DFT
[6-12]. On the other hand, the functional integration over fluctuating
potential field $\phi$ around $i\psi^{*}$ yields

$\displaystyle e^{-\Delta\Omega[\rho]}=\int
D\phi\,e^{-\Delta\mathcal{A}[\rho,\phi]},$ (133)

where

$\displaystyle\Delta\mathcal{A}[\rho,\phi]$
$\displaystyle=\Omega[v_{0},-i\phi+\psi^{*}]-\Omega_{0}[v_{0},\psi^{*}]+\int
d{\bf r}\,i\rho({\bf r})\phi({\bf r}).$ (134)

In the Gaussian approximation, the additional functional $\Delta\mathcal{A}$
around the saddle-point path $i\psi^{*}$ is written as

$\displaystyle\Delta\mathcal{A}[\rho,\phi]$ $\displaystyle=\frac{1}{2}\iint
d{\bf r}d{\bf r}^{\prime}\,\left.\phi({\bf
r})\frac{\delta^{2}\Omega_{0}}{\delta\psi({\bf r})\delta\psi({\bf
r}^{\prime})}\right|_{\psi=i\psi^{*}}\phi({\bf r}^{\prime})$
$\displaystyle=\frac{1}{2}\iint d{\bf r}d{\bf r}^{\prime}\,\phi({\bf
r})\,G_{0}({\bf r}-{\bf r}^{\prime};\rho({\bf r}))\,\phi({\bf r}^{\prime}),$
$\displaystyle G_{0}({\bf r}-{\bf r}^{\prime};\rho({\bf r}))$
$\displaystyle=\rho({\bf r})\rho({\bf r}^{\prime})h_{0}({\bf r}-{\bf
r}^{\prime})+\rho({\bf r})\delta({\bf r}-{\bf r}^{\prime}).$ (135)

It follows from eqs. (133) and (135) that

$\Delta\Omega=\frac{1}{2}\ln\left|\mathrm{det}\,G_{0}({\bf r}-{\bf
r}^{\prime};\rho({\bf r}))\right|.$ (136)

We thus obtain eqs. (80) to (83) from combining eqs. (130), (132) and (136).

###  B.2 Quadratic expansion of $\Omega^{*}[\rho]$ around the reference
density $\rho_{\mathrm{ref}}({\bf r})$: derivation of eqs. (84), (85), and
(93)

The functional Taylor expansion leads to

$\displaystyle\mathcal{F}_{0}[\rho_{\mathrm{ref}}+n]+U_{1}[\rho_{\mathrm{ref}}+n]$
$\displaystyle=\mathcal{F}_{0}[\rho_{\mathrm{ref}}]+U_{1}[\rho_{\mathrm{ref}}]+\int
d{\bf
r}\left.\left\\{\frac{\delta\mathcal{F}_{0}^{\mathrm{ex}}[\rho]}{\delta\rho({\bf
r})}\right|_{\rho=\rho_{\mathrm{ref}}}+\left.\frac{\delta\Omega_{\mathrm{att}}^{*}[\rho]}{\delta\rho({\bf
r})}\right|_{\rho=\rho_{\mathrm{ref}}}\right\\}n({\bf r})$
$\displaystyle\qquad\quad+\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\left\\{\left.\frac{\delta^{2}\mathcal{F}_{0}[\rho]}{\delta\rho({\bf
r})\delta\rho({\bf
r}^{\prime})}\right|_{\rho=\rho_{\mathrm{ref}}}+\left.\frac{\delta^{2}U_{1}[\rho]}{\delta\rho({\bf
r})\delta\rho({\bf
r}^{\prime})}\right|_{\rho=\rho_{\mathrm{ref}}}\right\\}n({\bf r})n({\bf
r}^{\prime})$
$\displaystyle\approx\mathcal{F}_{0}[\rho_{\mathrm{ref}}]+U_{1}[\rho_{\mathrm{ref}}]+\int
d{\bf r}\left.\frac{\delta\mathcal{F}_{0}^{\mathrm{ex}}[\rho]}{\delta\rho({\bf
r})}\right|_{\rho=\rho_{\mathrm{ref}}}n({\bf r})+\frac{1}{2}\iint d{\bf
r}d{\bf
r}^{\prime}\left.\frac{\delta^{2}\mathcal{F}_{0}[\rho]}{\delta\rho({\bf
r})\delta\rho({\bf r}^{\prime})}\right|_{\rho=\rho_{\mathrm{ref}}}n({\bf
r})n({\bf r}^{\prime}).$ (137)

When reduced to the last line of the above equation (137), we have used the
saddle-point equation (90), in addition to eqs. (81) to (88), or the equality,
$\mathcal{F}_{0}[\rho]+U_{1}[\rho]-\int d{\bf r}\rho({\bf
r})\mu=\mathcal{F}_{0}^{\mathrm{ex}}[\rho]+\Omega_{\mathrm{att}}[\rho]$, and
have also neglected the quadratic term in the $n$-expansion of
$U_{1}[\rho_{\mathrm{ref}}+n]$ due to the mean-field approximation of
attractive interaction energy. Plugging the RY functional (94) into eq. (137),
eq. (93) is verified with the approximation that

$\displaystyle\left.\frac{\delta^{2}\mathcal{F}_{0}[\rho]}{\delta\rho({\bf
r})\delta\rho({\bf
r}^{\prime})}\right|_{\rho=\rho_{\mathrm{ref}}}=G_{0}^{-1}({\bf r}-{\bf
r}^{\prime};\rho_{\mathrm{ref}}({\bf r}))$ $\displaystyle=-c_{0}({\bf r}-{\bf
r}^{\prime})+\frac{\delta({\bf r}-{\bf r}^{\prime})}{\rho_{\mathrm{ref}}({\bf
r})}$ $\displaystyle\approx-c_{0}({\bf r}-{\bf r}^{\prime})+\frac{\delta({\bf
r}-{\bf r}^{\prime})}{\rho_{\mathrm{ref}}^{0}}$ $\displaystyle\approx
G_{0}^{-1}({\bf r}-{\bf r}^{\prime};\rho_{\mathrm{ref}}^{0}),$ (138)

where the approximation used in the last line is the same as the approximation
in eq. (42). By applying the same approximation as eq. (138) to eq. (136), the
approximate logarithmic functional given by eqs. (84) and (85) can also be
validated.

###  B.3 Transforming eq. (96) to eq. (97) due to the shift of fluctuating
density field

We express $\Delta\mathcal{F}_{0}[n]$ given by eq. (96) as

$\displaystyle\Delta\mathcal{F}_{0}[n]=\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\,n({\bf r})\,G_{0}^{-1}({\bf r}-{\bf
r}^{\prime};\rho_{\mathrm{ref}}^{0})n({\bf r}^{\prime})-\iint d{\bf r}d{\bf
r}^{\prime}\,n({\bf r})c_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime}),$
$\displaystyle=\frac{1}{2}\int d{\bf 1}\,n({\bf 1})\left[\int d{\bf
2}\,G_{0}^{-1}({\bf 1}-{\bf 2})\underline{\left\\{n({\bf 2})-2\iint d{\bf
3}\,d{\bf 4}\>G_{0}({\bf 2}-{\bf 3})c_{0}({\bf 3}-{\bf
4})\Delta\rho_{\mathrm{ref}}({\bf 4})\right\\}}\right],$ (139)

where $G_{0}({\bf r}-{\bf r}^{\prime};\rho_{\mathrm{ref}}^{0})$ has been
abbreviated as $G_{0}({\bf r}-{\bf r}^{\prime})$ in the last equality. The
above underlined terms in the last line of eq. (139) implies the necessity to
shift the fluctuating density field from $n({\bf r})$to $\widetilde{n}({\bf
r})$ defined by eq. (98) for performing the Gaussian integration over the
fluctuating density field.

To see this, we rewrite the first term on the rhs of eq. (97) using the
$n$–field as follows:

$\displaystyle\frac{1}{2}\iint d{\bf 1}d{\bf 2}\,\widetilde{n}({\bf
1})G_{0}^{-1}({\bf 1}-{\bf 2})\widetilde{n}({\bf 2})$
$\displaystyle=\frac{1}{2}\iint d{\bf 1}d{\bf 2}\,n({\bf 1})G_{0}^{-1}({\bf
1}-{\bf 2})n({\bf 2})$ $\displaystyle-\frac{1}{2}\iiiint d{\bf 1}\,d{\bf
2}\,d{\bf 3}\,d{\bf 4}\>G_{0}^{-1}({\bf 1}-{\bf 2})$
$\displaystyle\qquad\times\left\\{n({\bf 1})G_{0}({\bf 2}-{\bf 3})c_{0}({\bf
3}-{\bf 4})\Delta\rho_{\mathrm{ref}}({\bf 4})+n({\bf 2})G_{0}({\bf 1}-{\bf
3})c_{0}({\bf 3}-{\bf 4})\Delta\rho_{\mathrm{ref}}({\bf 4})\right\\}$
$\displaystyle+\frac{1}{2}\iiint d{\bf 1}\,d{\bf 2}\,d{\bf 3}\>G_{0}^{-1}({\bf
1}-{\bf 2})$ $\displaystyle\qquad\times\left[\int d{\bf 4}\,G_{0}({\bf 2}-{\bf
3})c_{0}({\bf 3}-{\bf 4})\Delta\rho_{\mathrm{ref}}({\bf 4})\right]\left[\int
d{\bf 5}\,G_{0}({\bf 1}-{\bf 3})c_{0}({\bf 3}-{\bf
5})\Delta\rho_{\mathrm{ref}}({\bf 5})\right]$ $\displaystyle=\frac{1}{2}\iint
d{\bf r}d{\bf r}^{\prime}n({\bf r})G_{0}^{-1}({\bf r}-{\bf r}^{\prime})n({\bf
r}^{\prime})-\iint d{\bf r}d{\bf r}^{\prime}n({\bf r})c_{0}({\bf r}-{\bf
r}^{\prime})\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime})$
$\displaystyle\qquad+\frac{1}{2}\iiiint d{\bf 1}\,d{\bf 3}\,d{\bf 4}\,d{\bf
5}\Delta\rho_{\mathrm{ref}}({\bf 4})c_{0}({\bf 4}-{\bf 1})G_{0}({\bf 1}-{\bf
3})c_{0}({\bf 3}-{\bf 5})\Delta\rho_{\mathrm{ref}}({\bf 5})$
$\displaystyle=\Delta\mathcal{F}_{0}[n]+\frac{1}{2}\iint d{\bf r}d{\bf
r}^{\prime}\Delta\rho_{\mathrm{ref}}({\bf r})\left\\{h_{0}({\bf r}-{\bf
r}^{\prime})-c_{0}({\bf r}-{\bf
r}^{\prime})\right\\}\Delta\rho_{\mathrm{ref}}({\bf r}^{\prime}).$ (140)

In obtaining the above last line of eq. (140), use has been made of the
following derivation due to the Ornstein-Zernike equation (40):

$\displaystyle\iint d{\bf 1}\,d{\bf 3}\,c_{0}({\bf 4}-{\bf 1})G_{0}({\bf
1}-{\bf 3})c({\bf 3}-{\bf 5})$ $\displaystyle=\iint d{\bf 1}\,d{\bf
3}\,c_{0}({\bf 4}-{\bf 1})\left\\{\delta({\bf 1}-{\bf
3})+\rho_{\mathrm{ref}}^{0}h_{0}({\bf 1}-{\bf
3})\right\\}\rho_{\mathrm{ref}}^{0}c_{0}({\bf 3}-{\bf 5})$ $\displaystyle=\int
d{\bf 3}\,h_{0}({\bf 4}-{\bf 3})\rho_{\mathrm{ref}}^{0}c_{0}({\bf 3}-{\bf 5})$
$\displaystyle=h_{0}({\bf 4}-{\bf 5})-c_{0}({\bf 4}-{\bf 5}).$ (141)

Combination of eqs. (139) and (140) proves eq. (97).

## References

## References

  * [1] Hansen J P and Mcdonald I R 2013 Theory of Simple Liquids (Academic Press, London), and references therein.
  * [2] Widom B 1967 Science 157 375-382
  * [3] Chandler D, Weeks J D and Andersen H C 1983 Science 220 787-794
  * [4] Weeks J D, Vollmayr K and Katsov K 1997 Physica A 244 461-475
  * [5] Weeks J D 2002 Ann. Rev. Phys. Chem. 53 533-562
  * [6] Evans R 1979 Adv. Phys. 28 143-200.
  * [7] Evans R 1992 Fundamentals of inhomogeneous fluids (Marcel Dekker, New York) 85-175
  * [8] Likos C N 2001 Phys. Rep. 348 267-439
  * [9] Löwen H 2002 J. Phys.: Condens. Matt. 14 11897
  * [10] Wu J and Li Z 2007 Annu. Rev. Phys. Chem. 58 85-112
  * [11] Budkov Y A, Kolesnikov A L, Ivlev D V, Kalikin N N and Kiselev M G 2019 J. Mol. Liq. 276 801-805
  * [12] Budkov Y A and Kolesnikov A L 2020 J. Phys.: Condens. Mat. 32 365001
  * [13] Chandler D 2005 Nature 437 640-647
  * [14] Remsing R C, Liu S, Weeks J D 2016 P. Natl. Acad. Sci. USA 113 2819
  * [15] Naji A, Kanduc M, Forsman J and Podgornik R 2013 J. Chem. Phys 139 150901
  * [16] Bazant M Z, Storey B D and Kornyshev A A 2011 Phys. Rev. Lett. 106 046102
  * [17] Misra R P, de Souza J P, Blankschtein D and Bazant M Z 2019 Langmuir 35 11550-11565
  * [18] de Souza J P, Goodwin Z A H, McEldrew M, Kornyshev A A and Bazant M Z 2020 Phys. Rev. Lett. 125 116001
  * [19] Coles S W, Park C, Nikam R, Kanduc M, Dzubiella J and Rotenberg B 2020 J. Phys. Chem. B 124 1778-1786
  * [20] Chen Y G, Kaur C and Weeks J D 2004 J. Phys. Chem. B 108 19874
  * [21] Chen Y G and Weeks J D 2006 P. Natl. Acad. Sci. USA 103 7560
  * [22] Frydel D 2016 Adv. Chem. Phys. 160 209-260
  * [23] Frusawa H 2018 Phys. Rev. E 98 052130, and references therein
  * [24] Brilliantov N V, Rubí J M and Budkov Y A 2020 Phys. Rev. E 101 042135
  * [25] Santangelo C D 2006 Phys. Rev. E 73 041512
  * [26] Hatlo M M and Lue L 2009 Soft Matter 5 125-133.
  * [27] Hatlo M M and Lue L 2010 Europhys. Lett. 89 25002
  * [28] Blossey R, Maggs A C and Podgornik R 2017 Phys. Rev. E 95 060602
  * [29] de Souza J P and Bazant M Z 2020 J. Phys. Chem. C 124 11414-11421
  * [30] Avni Y, Adar R M and Andelman, D 2020 Phys. Rev. E 101 010601
  * [31] Netz R R and Orland H 2003 Eur. Phys. J. E 11 301-311
  * [32] Buyukdagli S, Achim C V and Ala-Nissila T 2011 J. Stat. Mech.: Theory E. 5 P05033
  * [33] Buyukdagli S and Blossey R 2014 J. Chem. Phys. 140 234903
  * [34] Buyukdagli S and Blossey R 2016 J. Phys.: Condens. Mat. 28 343001
  * [35] Warren P B and Vlasov A 2014 J. Chem. Phys 140 084904
  * [36] Patsahan O 2015 Condens. Matt. Phys. 18 24601-24601
  * [37] Frydel D 2016 J. Chem. Phys. 145 184703
  * [38] Ng K C 1974 J. Chem. Phys. 61 2680
  * [39] Rosenfeld Y and Gelbart WM 1984 J. Chem. Phys. 81 4574
  * [40] de Carvalho R L, Evans R and Rosenfeld Y 1999 Phys. Rev. E 59 1435
  * [41] Caccamo C, Pellicane G, Costa D, Pini D and Stell G 1999 Phys. Rev. E 60 5533
  * [42] Lum K, Chandler D and Weeks J D 1999 J. Phys. Chem. B 103 4570-4577
  * [43] Weeks J D, Selinger R L and Broughton J Q 1995 Phys. Rev. Lett. 74 2694
  * [44] Katsov K and Weeks J D 2001 Phys. Rev. Lett. 864 440
  * [45] Vollmayr-Lee K, Katsov K and Weeks J D 2001 J. Chem. Phys. 114 416-425
  * [46] Archer A J and Evans R 2013 J. Chem. Phys 138 014502
  * [47] Archer A J, Chacko B and Evans R 2017 J. Chem. Phys 147 034501
  * [48] Frusawa H 2018 J. Phys. A: Math. Theor. 52 065003
  * [49] Frusawa H 2020 Entropy 22 34
  * [50] Brilliantov N V 1998 Phys. Rev. E 58 2628
  * [51] Caillol J M, Patsahan O and Mryglod I 2006 Physica A 368 326-344
  * [52] Budkov Y A 2020 Phys. Chem. Chem. Phys. 22 14756-14772
  * [53] Percus J K 1962 Phys. Rev. Lett. 8 462